Merge pull request #345 from FluidityProject/distinct_particle_attrib…

…_on_spawn ENH: Allow for independent particle attribute values on initialisation and advection
FluidityProject · Jun 17, 2023 · aedd3df · aedd3df
2 parents a8d2e6f + d70044d
commit aedd3df
Show file tree

Hide file tree

Showing 33 changed files with 1,359 additions and 1,062 deletions.
diff --git a/assemble/Adapt_State.F90 b/assemble/Adapt_State.F90
@@ -1006,7 +1006,7 @@ subroutine adapt_state_first_timestep(states)
 
       !Set particle attributes and dependent fields
       call get_option("/timestepping/current_time", current_time)
-      call update_particle_attributes_and_fields(states, current_time, dt)
+      call update_particle_attributes_and_fields(states, current_time, dt, initial=.true.)
       call calculate_diagnostic_fields_from_particles(states)
 
       ! Form the new metric

diff --git a/assemble/Particle_Diagnostics.F90 b/assemble/Particle_Diagnostics.F90
@@ -48,16 +48,15 @@ module particle_diagnostics
   use detector_tools, only: temp_list_insert, insert, allocate, deallocate, temp_list_deallocate, &
        & remove, temp_list_remove
   use particles, only : get_particle_arrays, initialise_constant_particle_attributes, &
-       & update_particle_attributes_and_fields, get_particle_arrays, particle_lists
+       & particle_lists
   use multimaterial_module, only: calculate_sum_material_volume_fractions
 
   implicit none
 
   private
 
-  public :: update_particle_attributes_and_fields, calculate_particle_material_fields, &
-       & calculate_diagnostic_fields_from_particles, initialise_constant_particle_diagnostics, &
-       & particle_cv_check
+  public :: calculate_particle_material_fields, calculate_diagnostic_fields_from_particles, &
+       & initialise_constant_particle_diagnostics, particle_cv_check
 
   contains
 
@@ -95,9 +94,12 @@ subroutine initialise_constant_particle_diagnostics(state)
        do j = 1,particle_arrays(i)
           particle => particle_lists(list_counter)%first
           subgroup_path = trim(group_path) // "/particle_subgroup["//int2str(j-1)//"]"
-          if (option_count(trim(subgroup_path) // "/attributes/scalar_attribute/constant") + &
-               & option_count(trim(subgroup_path) // "/attributes/vector_attribute/constant") + &
-               & option_count(trim(subgroup_path) // "/attributes/tensor_attribute/constant")>0) then
+          if (option_count(trim(subgroup_path) // "/attributes/scalar_attribute/initial_attribute_value/constant") + &
+               & option_count(trim(subgroup_path) // "/attributes/scalar_attribute/attribute_value/constant") + &
+               & option_count(trim(subgroup_path) // "/attributes/vector_attribute/initial_attribute_value/constant") + &
+               & option_count(trim(subgroup_path) // "/attributes/vector_attribute/attribute_value/constant") + &
+               & option_count(trim(subgroup_path) // "/attributes/tensor_attribute/initial_attribute_value/constant") + &
+               & option_count(trim(subgroup_path) // "/attributes/tensor_attribute/attribute_value/constant")>0) then
              call initialise_constant_particle_attributes(state, subgroup_path, particle_lists(list_counter))
           end if
           list_counter = list_counter + 1
@@ -126,24 +128,24 @@ subroutine calculate_diagnostic_fields_from_particles(state)
 
     !Initialise diagnostic fields generated from particles
     do i = 1,size(state)
-       do k = 1,scalar_field_count(state(i))
-          s_field => extract_scalar_field(state(i),k)
-          if (have_option(trim(s_field%option_path)//"/diagnostic/algorithm::from_particles")) then
-             call get_option(trim(s_field%option_path)//"/name", name)
-             if (name=="MaterialVolumeFraction") cycle
-             call calculate_field_from_particles(state, i, s_field)
-          end if
-       end do
+      do k = 1,scalar_field_count(state(i))
+        s_field => extract_scalar_field(state(i),k)
+        if (have_option(trim(s_field%option_path)//"/diagnostic/algorithm::from_particles")) then
+          call get_option(trim(s_field%option_path)//"/name", name)
+          if (name=="MaterialVolumeFraction") cycle
+          call calculate_field_from_particles(state, i, s_field)
+        end if
+      end do
     end do
 
     !Initialise fields based on number of particles if present
     do i = 1,size(state)
-       do k = 1,scalar_field_count(state(i))
-          s_field => extract_scalar_field(state(i),k)
-          if (have_option(trim(s_field%option_path)//"/diagnostic/algorithm::number_of_particles")) then
-             call calculate_field_from_particles(state, i, s_field)
-          end if
-       end do
+      do k = 1,scalar_field_count(state(i))
+        s_field => extract_scalar_field(state(i),k)
+        if (have_option(trim(s_field%option_path)//"/diagnostic/algorithm::number_of_particles")) then
+          call calculate_field_from_particles(state, i, s_field)
+        end if
+      end do
     end do
 
   end subroutine calculate_diagnostic_fields_from_particles

diff --git a/docker/actions/Dockerfile.actions.bionic b/docker/actions/Dockerfile.actions.bionic
@@ -8,6 +8,8 @@ RUN apt-get -y update && \
       rm -rf /var/cache/apt/archives && \
       rm -rf /var/lib/apt/lists
 
+RUN python3 -m pip install --upgrade numpy scipy
+
 RUN adduser fluidity sudo
 RUN echo '%sudo ALL=(ALL) NOPASSWD:ALL' >> /etc/sudoers
 

diff --git a/examples/particle_rayleigh_taylor_mu10/particle-rayleigh-taylor-mu10-24.flml b/examples/particle_rayleigh_taylor_mu10/particle-rayleigh-taylor-mu10-24.flml
@@ -125,9 +125,11 @@
         </initial_position>
         <attributes>
           <scalar_attribute name="Chem1">
-            <constant>
-              <real_value rank="0">0</real_value>
-            </constant>
+            <attribute_value>
+              <constant>
+                <real_value rank="0">0</real_value>
+              </constant>
+            </attribute_value>
           </scalar_attribute>
         </attributes>
       </particle_subgroup>
@@ -156,9 +158,11 @@
         </initial_position>
         <attributes>
           <scalar_attribute name="Chem1">
-            <constant>
-              <real_value rank="0">1</real_value>
-            </constant>
+            <attribute_value>
+              <constant>
+                <real_value rank="0">1</real_value>
+              </constant>
+            </attribute_value>
           </scalar_attribute>
         </attributes>
       </particle_subgroup>

diff --git a/examples/particle_rayleigh_taylor_mu10/particle-rayleigh-taylor-mu10-48.flml b/examples/particle_rayleigh_taylor_mu10/particle-rayleigh-taylor-mu10-48.flml
@@ -125,9 +125,11 @@
         </initial_position>
         <attributes>
           <scalar_attribute name="Chem1">
-            <constant>
-              <real_value rank="0">0</real_value>
-            </constant>
+            <attribute_value>
+              <constant>
+                <real_value rank="0">0</real_value>
+              </constant>
+            </attribute_value>
           </scalar_attribute>
         </attributes>
       </particle_subgroup>
@@ -156,9 +158,11 @@
         </initial_position>
         <attributes>
           <scalar_attribute name="Chem1">
-            <constant>
-              <real_value rank="0">1</real_value>
-            </constant>
+            <attribute_value>
+              <constant>
+                <real_value rank="0">1</real_value>
+              </constant>
+            </attribute_value>
           </scalar_attribute>
         </attributes>
       </particle_subgroup>

diff --git a/femtools/CGAL_Tools.F90 b/femtools/CGAL_Tools.F90
@@ -2,8 +2,8 @@
 module cgal_tools
 
   use fldebug
-  use fields
   use element_numbering
+  use fields
   implicit none
 
   interface

diff --git a/femtools/Detector_Tools.F90 b/femtools/Detector_Tools.F90
@@ -784,7 +784,7 @@ end subroutine set_particle_scalar_attribute_from_python
   !! specified in the string func at that location.
   subroutine set_particle_scalar_attribute_from_python_fields(particle_list, state, positions, lcoords, ele, natt, &
        attributes, old_attr_names, old_attr_counts, old_attr_Dims, old_attributes, field_names, field_counts, old_field_names, &
-       old_field_counts, func, time, dt, is_array)
+       old_field_counts, func, time, dt, is_array, first_newly_init_part)
     !! Particle list for which to evaluate the function
     type(detector_linked_list), intent(in) :: particle_list
     !! Model state structure
@@ -828,6 +828,8 @@ subroutine set_particle_scalar_attribute_from_python_fields(particle_list, state
     real, intent(in) :: dt
     !! Whether this is an array-valued attribute
     logical, intent(in) :: is_array
+    !> Pointer to the first newly initialised particle
+    type(detector_type), pointer, optional :: first_newly_init_part
 
     ! locals
     integer :: i
@@ -864,7 +866,12 @@ subroutine set_particle_scalar_attribute_from_python_fields(particle_list, state
          field_vals, dim)
 
     ! copy old fields off particles
-    particle => particle_list%first
+    if (present(first_newly_init_part)) then
+      particle => first_newly_init_part
+    else
+      particle => particle_list%first
+    end if
+
     do i = 1, nparts
       old_field_vals(:,i) = particle%old_fields(:)
       particle => particle%next
@@ -874,9 +881,9 @@ subroutine set_particle_scalar_attribute_from_python_fields(particle_list, state
          lvx, lvy, lvz, time, dt, field_counts, field_names, field_vals, old_field_counts, old_field_names, &
          old_field_vals, old_attr_counts, old_attr_names, old_attr_dims, old_attributes, is_array, attributes, stat)
     if (stat/=0) then
-       ewrite(-1, *) "Python error, Python string was:"
-       ewrite(-1 , *) trim(func)
-       FLExit("Dying")
+      ewrite(-1, *) "Python error, Python string was:"
+      ewrite(-1 , *) trim(func)
+      FLExit("Dying")
     end if
 
     deallocate(field_vals)
@@ -943,7 +950,7 @@ end subroutine set_particle_vector_attribute_from_python
   !! specified in the string func at that location.
   subroutine set_particle_vector_attribute_from_python_fields(particle_list, state, positions, lcoords, ele, natt, &
        attributes, old_attr_names, old_attr_counts, old_attr_dims, old_attributes, field_names, field_counts, old_field_names, &
-       old_field_counts, func, time, dt, is_array)
+       old_field_counts, func, time, dt, is_array, first_newly_init_part)
     !! Particle list for which to evaluate the function
     type(detector_linked_list), intent(in) :: particle_list
     !! Model state structure
@@ -987,6 +994,8 @@ subroutine set_particle_vector_attribute_from_python_fields(particle_list, state
     real, intent(in) :: dt
     !! Whether this is an array-valued attribute
     logical, intent(in) :: is_array
+    !> Pointer to the first newly initialised particle
+    type(detector_type), pointer, optional :: first_newly_init_part
 
     ! locals
     integer :: i
@@ -1023,7 +1032,12 @@ subroutine set_particle_vector_attribute_from_python_fields(particle_list, state
          field_vals, dim)
 
     ! copy old fields off particles
-    particle => particle_list%first
+    if (present(first_newly_init_part)) then
+      particle => first_newly_init_part
+    else
+      particle => particle_list%first
+    end if
+
     do i = 1, nparts
       old_field_vals(:,i) = particle%old_fields(:)
       particle => particle%next
@@ -1109,7 +1123,7 @@ end subroutine set_particle_tensor_attribute_from_python
   !! specified in the string func at that location.
   subroutine set_particle_tensor_attribute_from_python_fields(particle_list, state, positions, lcoords, ele, natt, &
        attributes, old_attr_names, old_attr_counts, old_attr_dims, old_attributes, field_names, field_counts, old_field_names, &
-       old_field_counts, func, time, dt, is_array)
+       old_field_counts, func, time, dt, is_array, first_newly_init_part)
     !! Particle list for which to evaluate the function
     type(detector_linked_list), intent(in) :: particle_list
     !! Model state structure
@@ -1153,6 +1167,8 @@ subroutine set_particle_tensor_attribute_from_python_fields(particle_list, state
     real, intent(in) :: dt
     !! Whether this is an array-valued attribute
     logical, intent(in) :: is_array
+    !> Pointer to the first newly initialised particle
+    type(detector_type), pointer, optional :: first_newly_init_part
 
     ! locals
     integer :: i
@@ -1190,7 +1206,12 @@ subroutine set_particle_tensor_attribute_from_python_fields(particle_list, state
          field_vals, dim)
 
     ! copy old fields off particles
-    particle => particle_list%first
+    if (present(first_newly_init_part)) then
+      particle => first_newly_init_part
+    else
+      particle => particle_list%first
+    end if
+
     do i = 1, nparts
       old_field_vals(:,i) = particle%old_fields(:)
       particle => particle%next