Rebasing of merged_particles branch

assemble/Adapt_State.F90

@@ -44,6 +44,7 @@ module adapt_state_module
   use parallel_fields
   use intersection_finder_module
   use fields
+  use profiler
   use state_module
   use vtk_interfaces
   use halos
@@ -62,6 +63,7 @@ module adapt_state_module
   use adaptive_timestepping
   use detector_parallel
   use diagnostic_variables
+  use particle_diagnostics, only: initialise_constant_particle_diagnostics, initialise_particle_diagnostics
   use checkpoint
   use edge_length_module
   use boundary_conditions_from_options
@@ -987,9 +989,15 @@ subroutine adapt_state_first_timestep(states)
       call copy_to_stored_values(states,"Iterated")
       call relax_to_nonlinear(states)
 
+      !Set constant particle attributes and MVF fields based on particles
+      call initialise_constant_particle_diagnostics(states)
+
       call calculate_diagnostic_variables(states)
       call calculate_diagnostic_variables_new(states)
 
+      !Set particle attributes and dependent fields
+      call initialise_particle_diagnostics(states)
+
       call enforce_discrete_properties(states)
       if(have_option("/timestepping/adaptive_timestep/at_first_timestep")) then
         ! doing this here helps metric advection get the right amount of advection
@@ -1165,6 +1173,16 @@ subroutine adapt_state_internal(states, metric, initialise_fields)
       ! been deallocated
       call tag_references()
 
+      !Check if particle lists are initialised
+      if (get_num_detector_lists()>0) then
+        call get_registered_detector_lists(detector_list_array)
+        ! Check if particle elements exist on this processor, pack to other processor if not
+        do j = 1, size(detector_list_array)
+           call l2_halo_detectors(detector_list_array(j)%ptr, old_positions, states(1))
+        end do
+
+      end if
+
       ! Generate a new mesh field based on the current mesh field and the input
       ! metric
       if (.not. vertical_only) then
@@ -1211,11 +1229,11 @@ subroutine adapt_state_internal(states, metric, initialise_fields)
 
       if (get_num_detector_lists()>0) then
         ! Update detector element and local_coords for every detector in all lists
-        call get_registered_detector_lists(detector_list_array)
+        call profiler_tic("find_particles_mesh_adapt")
         do j = 1, size(detector_list_array)
            call search_for_detectors(detector_list_array(j)%ptr, new_positions)
         end do
-
+        call profiler_toc("find_particles_mesh_adapt")
 #ifdef DDEBUG
         ! Sanity check that all local detectors are owned
         call get_registered_detector_lists(detector_list_array)

assemble/Diagnostic_fields_wrapper.F90

@@ -76,6 +76,7 @@ subroutine calculate_diagnostic_variables(state, exclude_nonrecalculated)
     type(scalar_field), pointer :: s_field
     type(vector_field), pointer :: v_field
     logical :: diagnostic
+    integer :: diagnostic_particles
 
     ! An array of submaterials of the current phase in state(istate).
     type(state_type), dimension(:), pointer :: submaterials
@@ -361,8 +362,10 @@ subroutine calculate_diagnostic_variables(state, exclude_nonrecalculated)
 
        s_field => extract_scalar_field(state(i), "MaterialVolumeFraction", stat)
        if(stat == 0) then
-         diagnostic = have_option(trim(s_field%option_path)//"/diagnostic")
-         if(diagnostic .and. .not.(aliased(s_field))) then
+         diagnostic = have_option(trim(s_field%option_path)//"/diagnostic/algorithm::Internal")
+         !Check if MaterialVolumeFraction field is set from particles, if so don't calculate internal MVF here
+         diagnostic_particles = option_count("material_phase/scalar_field::MaterialVolumeFraction/diagnostic/algorithm::from_particles")
+         if(diagnostic .and. .not.(aliased(s_field) .and. diagnostic_particles.lt.1)) then
            if(recalculate(trim(s_field%option_path))) then
              call calculate_sum_material_volume_fractions(state, s_field)
              call scale(s_field, -1.0)

assemble/Makefile.dependencies

@@ -37,7 +37,8 @@ Adapt_State.o ../include/adapt_state_module.mod: Adapt_State.F90 \
    ../include/metric_assemble.mod ../include/metric_tools.mod \
    ../include/mpi_interfaces.mod ../include/node_boundary.mod \
    ../include/parallel_fields.mod ../include/parallel_tools.mod \
-   ../include/pickers.mod ../include/populate_state_module.mod \
+   ../include/particle_diagnostics.mod ../include/pickers.mod \
+   ../include/populate_state_module.mod ../include/profiler.mod \
    ../include/project_metric_to_surface_module.mod ../include/quadrature.mod \
    ../include/reference_counting.mod ../include/reserve_state_module.mod \
    ../include/sam_integration.mod ../include/sparse_tools.mod \
@@ -628,6 +629,20 @@ OceanSurfaceForcing.o ../include/oceansurfaceforcing.mod: \
    ../include/parallel_fields.mod ../include/parallel_tools.mod \
    ../include/state_module.mod ../include/transform_elements.mod
 
+../include/particle_diagnostics.mod: Particle_Diagnostics.o
+	@true
+
+Particle_Diagnostics.o ../include/particle_diagnostics.mod: \
+   Particle_Diagnostics.F90 ../include/detector_data_types.mod \
+   ../include/detector_tools.mod ../include/elements.mod ../include/fdebug.h \
+   ../include/field_options.mod ../include/fields.mod \
+   ../include/fields_base.mod ../include/fields_calculations.mod \
+   ../include/fldebug.mod ../include/futils.mod \
+   ../include/global_parameters.mod ../include/halos.mod \
+   ../include/multimaterial_module.mod ../include/parallel_fields.mod \
+   ../include/parallel_tools.mod ../include/particles.mod \
+   ../include/pickers.mod ../include/profiler.mod ../include/state_module.mod
+
 ../include/petsc_solve_state_module.mod: Petsc_Solve_State.o
 	@true
 

assemble/Multimaterials.F90

@@ -226,13 +226,13 @@ subroutine calculate_diagnostic_material_volume_fraction(state)
     type(scalar_field), pointer :: sfield
     logical :: diagnostic
 
-    ! How many diagnostic MaterialVolumeFraction fields do we have in state?
+    ! How many diagnostic internal MaterialVolumeFraction fields do we have in state?
     ! Note that state contains all the submaterials of the current phase, including the phase itself.
     ! Therefore, if the only material is the phase itself, diagnostic_count should be 0. Otherwise,
     ! it should be 1.
     diagnostic_count = 0
     do i = 1, size(state)
-       if(have_option(trim(state(i)%option_path)//"/scalar_field::MaterialVolumeFraction/diagnostic")) then
+       if(have_option(trim(state(i)%option_path)//"/scalar_field::MaterialVolumeFraction/diagnostic/algorithm::Internal")) then
           diagnostic_count = diagnostic_count + 1
           ! Record the index of the state containing the diagnostic MaterialVolumeFraction field
           diagnostic_state_index = i 
@@ -241,7 +241,7 @@ subroutine calculate_diagnostic_material_volume_fraction(state)
 
     if(diagnostic_count>1) then
       ewrite(0,*) diagnostic_count, ' diagnostic MaterialVolumeFractions'
-      FLExit("Only one diagnostic MaterialVolumeFraction permitted.")
+      FLExit("Only one internal diagnostic MaterialVolumeFraction permitted.")
     end if
 
     if(diagnostic_count==1) then
@@ -253,7 +253,7 @@ subroutine calculate_diagnostic_material_volume_fraction(state)
 
       do i = 1,size(state)
         sfield=>extract_scalar_field(state(i),'MaterialVolumeFraction',stat)
-        diagnostic=(have_option(trim(sfield%option_path)//'/diagnostic'))
+        diagnostic=(have_option(trim(sfield%option_path)//'/diagnostic/algorithm::Internal'))
         if ( (stat==0).and.(.not. aliased(sfield)).and.(.not.diagnostic)) then
           call addto(sumvolumefractions, sfield)
         end if
@@ -807,9 +807,9 @@ subroutine calculate_sum_material_volume_fractions(state,sumvolumefractions)
     !locals
     integer :: i, stat
     type(scalar_field), pointer :: sfield
-    logical :: prognostic, diagnostic, prescribed
+    logical :: prognostic, diagnostic, prescribed, diagnostic_particles
 
-    diagnostic = have_option(trim(sumvolumefractions%option_path)//"/diagnostic")
+    diagnostic = have_option(trim(sumvolumefractions%option_path)//"/diagnostic/algorithm::Internal")
     if(.not.diagnostic) return
 
     call zero(sumvolumefractions)
@@ -819,7 +819,8 @@ subroutine calculate_sum_material_volume_fractions(state,sumvolumefractions)
       if(stat==0) then
         prognostic = have_option(trim(sfield%option_path)//"/prognostic")
         prescribed = have_option(trim(sfield%option_path)//"/prescribed")
-        if ((.not.aliased(sfield)).and.(prognostic.or.prescribed)) then
+        diagnostic_particles = have_option(trim(sfield%option_path)//"/diagnostic/algorithm::from_particles")
+        if ((.not.aliased(sfield)).and.(prognostic.or.prescribed.or.diagnostic_particles)) then
           call addto(sumvolumefractions, sfield)
         end if
       end if

assemble/Sam_integration.F90

@@ -1500,8 +1500,7 @@ subroutine transfer_detectors(old_mesh, new_positions, node_ownership)
           ! Pack this node for sending
 
           ! Integer data
-          isend_data(owner)%ptr(data_index(owner) * idata_size + 1) = node%type
-          isend_data(owner)%ptr(data_index(owner) * idata_size + 2) = node%id_number
+          isend_data(owner)%ptr(data_index(owner) * idata_size + 1) = node%id_number
           ! Real data
           rsend_data(owner)%ptr(data_index(owner) * rdata_size + 1:data_index(owner) * rdata_size + new_positions%dim) = node%position
 
@@ -1574,8 +1573,7 @@ subroutine transfer_detectors(old_mesh, new_positions, node_ownership)
         allocate(node)
 
         ! Integer data
-        node%type = ireceive_data(i)%ptr((j - 1) * idata_size + 1)
-        node%id_number = ireceive_data(i)%ptr((j - 1) * idata_size + 2)
+        node%id_number = ireceive_data(i)%ptr((j - 1) * idata_size + 1)
 
         ! Real data
         allocate(node%position(new_positions%dim))

assemble/Zoltan_integration.F90

@@ -2059,7 +2059,7 @@ subroutine update_particle_list_element(detector_list_array)
     integer :: old_local_element_number, new_local_element_number, old_universal_element_number
     integer, allocatable :: ndets_being_sent(:)
     real, allocatable :: send_buff(:,:), recv_buff(:,:)
-    logical do_broadcast
+    logical :: do_broadcast, sent
     integer, dimension(3) :: attribute_size
     integer :: total_attributes
     type(element_type), pointer :: shape
@@ -2081,7 +2081,6 @@ subroutine update_particle_list_element(detector_list_array)
              attribute_size(1)=size(detector%attributes)
              attribute_size(2)=size(detector%old_attributes)
              attribute_size(3)=size(detector%old_fields)
-             total_attributes=sum(attribute_size)
           end if
        end if
        do while (associated(detector))
@@ -2130,6 +2129,17 @@ subroutine update_particle_list_element(detector_list_array)
           deallocate(ndets_being_sent)
           cycle
        end if
+       i=1
+       sent=.false.
+       do while (sent.eqv..false.)
+          if (ndets_being_sent(i)>0) then
+             call mpi_bcast(attribute_size, 3, getPINTEGER(), i-1, MPI_COMM_FEMTOOLS, ierr)
+             assert(ierr == MPI_SUCCESS)
+             total_attributes=sum(attribute_size)
+             sent=.true.
+          end if
+          i=i+1
+       end do
        ewrite(2,*) "Broadcast required, initialising..."
 
        ! Allocate memory for all the particles you're going to send
@@ -2323,12 +2333,7 @@ subroutine transfer_fields(zz)
        add_detector => detector
        detector => detector%next
 
-       ! update detector name if names are present on the list, otherwise det%name=id_number
-       if (allocated(detector_list_array(add_detector%list_id)%ptr%detector_names)) then
-          add_detector%name=detector_list_array(add_detector%list_id)%ptr%detector_names(add_detector%id_number)
-       else
-          add_detector%name=int2str(add_detector%id_number)
-       end if
+       add_detector%name=int2str(add_detector%id_number)
 
        ! move detector to the correct list
        call move(add_detector, zoltan_global_unpacked_detectors_list, detector_list_array(add_detector%list_id)%ptr)

diagnostics/Diagnostic_Fields_New.F90.in

@@ -308,7 +308,10 @@ contains
         if(.not. recalculate_diagnostic_variable(s_field, exclude_nonrecalculated = exclude_nonrecalculated)) cycle
         ! and whether this field has already been calculated (as a dependency)
         if(has_scalar_field(calculated_states(i), trim(s_field%name))) cycle
-
+        ! Check if the field should be calculated from particles instead
+        if(have_option(trim(s_field%option_path) // "/diagnostic/algorithm::from_particles")) cycle
+        if(have_option(trim(s_field%option_path) // "/diagnostic/algorithm::number_of_particles")) cycle
+
         ! Calculate dependencies
         call calculate_dependencies(states, i, s_field, &
           & dep_states_mask = calculated_states, exclude_nonrecalculated = exclude_nonrecalculated)
@@ -842,6 +845,12 @@ SINGLE_STATE_SCALAR_DIAGNOSTICS
       ! Diagnostic algorithms that require the full states array
 MULTIPLE_STATE_SCALAR_DIAGNOSTICS
 
+      case("from_particles")
+         ! Not handled here
+
+      case("number_of_particles")
+         ! Not handled here
+
       case default
         if(present(stat)) then
           stat = 1

femtools/Detector_Data_Types.F90

@@ -38,9 +38,7 @@ module detector_data_types
 
   public :: detector_type, rk_gs_parameters, detector_linked_list, &
             detector_list_ptr, stringlist, attr_names_type, attr_write_type, field_phase_type, &
-            STATIC_DETECTOR, LAGRANGIAN_DETECTOR, allocate, deallocate
-
-  integer, parameter :: STATIC_DETECTOR=1, LAGRANGIAN_DETECTOR=2
+            allocate, deallocate
 
   type stringlist
     !!< Container type for a list of strings.
@@ -88,8 +86,6 @@ module detector_data_types
      integer :: element
      !! Local coordinates of the detector in that element.
      real, dimension(:), allocatable :: local_coords
-     !! Whether the detector is static or Lagrangian.
-     integer :: type = STATIC_DETECTOR
      !! Identification number indicating the order in which the detectors are read
      integer :: id_number
      !! Identification number indicating parent processor when detector was created
@@ -111,6 +107,8 @@ module detector_data_types
      !! Pointers for detector linked lists
      TYPE (detector_type), POINTER :: next=> null()
      TYPE (detector_type), POINTER :: previous=> null()
+     TYPE (detector_type), POINTER :: temp_next => null()
+     TYPE (detector_type), POINTER :: temp_previous => null()
   end type detector_type
 
   ! Parameters for lagrangian detector movement