prepare property lookup

source/simulator/initial_conditions.cc

@@ -284,36 +284,51 @@ namespace aspect
  // need to write into it and we can not
  // write into vectors with ghost elements
 
- const Particle::Interpolator::Interface<dim> &particle_interpolator = particle_worlds[0]->get_interpolator();
- const Particle::Property::Manager<dim> &particle_property_manager = particle_worlds[0]->get_property_manager();
+ // For each particle world store the pairs of corresponding advection field index (first) and particle property index (second)
+ std::vector<std::vector<std::pair<unsigned int, unsigned int>>> particle_property_indices;
+ // A property component mask indicating for each particle world which particle properties need to be interpolated
+ std::vector<ComponentMask> property_mask;
 
- std::vector<unsigned int> particle_property_indices;
- ComponentMask property_mask (particle_property_manager.get_data_info().n_components(),false);
-
- for (const auto &advection_field: advection_fields)
+ for (unsigned int world_index = 0; world_index < particle_worlds.size(); ++world_index)
  {
- if (parameters.mapped_particle_properties.size() != 0)
- {
- const std::pair<std::string,unsigned int> particle_property_and_component = parameters.mapped_particle_properties.find(advection_field.compositional_variable)->second;
-
- const unsigned int particle_property_index = particle_property_manager.get_data_info().get_position_by_field_name(particle_property_and_component.first)
- + particle_property_and_component.second;
+ const Particle::Property::Manager<dim> &particle_property_manager = particle_worlds[world_index]->get_property_manager();
 
- particle_property_indices.push_back(particle_property_index);
+ particle_property_indices.push_back(std::vector<std::pair<unsigned int, unsigned int>>());
+ property_mask.push_back(ComponentMask(particle_property_manager.get_data_info().n_components(),false));
 
- property_mask.set(particle_property_index,true);
- }
- else
+ for (unsigned int advection_field=0; advection_field<advection_fields.size(); ++advection_field)
  {
- const unsigned int particle_property_index = std::count(introspection.compositional_field_methods.begin(),
- introspection.compositional_field_methods.begin() + advection_field.compositional_variable,
- Parameters<dim>::AdvectionFieldMethod::particles);
- AssertThrow(particle_property_index <= particle_property_manager.get_data_info().n_components(),
- ExcMessage("Can not automatically match particle properties to fields, because there are"
- "more fields that are marked as particle advected than particle properties"));
-
- particle_property_indices.push_back(particle_property_index);
- property_mask.set(particle_property_index,true);
+ if (parameters.mapped_particle_properties.size() != 0)
+ {
+ const std::pair<std::string,unsigned int> particle_property_and_component = parameters.mapped_particle_properties.find(advection_fields[advection_field].compositional_variable)->second;
+
+ // Check if the required particle property exists in the current particle world.
+ // If not: assume we find it in another world.
+ if (particle_property_manager.get_data_info().fieldname_exists(particle_property_and_component.first))
+ {
+ const unsigned int particle_property_index = particle_property_manager.get_data_info().get_position_by_field_name(particle_property_and_component.first)
+ + particle_property_and_component.second;
+
+ particle_property_indices[world_index].push_back({advection_field, particle_property_index});
+ property_mask[world_index].set(particle_property_index,true);
+ }
+ }
+ else
+ {
+ Assert(particle_worlds.size() == 1,
+ ExcMessage("Automatically mapping particle properties to compositional fields is only supported if there is exactly one set of particles. "
+ "Please specify the particle properties manually in the parameter file using the parameter 'Compositional Fields/Mapped particle properties'."));
+
+ const unsigned int particle_property_index = std::count(introspection.compositional_field_methods.begin(),
+ introspection.compositional_field_methods.begin() + advection_fields[advection_field].compositional_variable,
+ Parameters<dim>::AdvectionFieldMethod::particles);
+ AssertThrow(particle_property_index <= particle_property_manager.get_data_info().n_components(),
+ ExcMessage("Can not automatically match particle properties to fields, because there are"
+ "more fields that are marked as particle advected than particle properties"));
+
+ particle_property_indices[world_index].push_back({advection_field,particle_property_index});
+ property_mask[world_index].set(particle_property_index,true);
+ }
  }
  }
 
@@ -357,10 +372,10 @@ namespace aspect
  try
  {
  particle_properties =
- particle_interpolator.properties_at_points(particle_worlds[world_index]->get_particle_handler(),
- quadrature_points,
- property_mask,
- cell);
+ particle_worlds[world_index]->get_interpolator().properties_at_points(particle_worlds[world_index]->get_particle_handler(),
+  quadrature_points,
+  property_mask[world_index],
+  cell);
  }
  // interpolators that throw exceptions usually do not result in
  // anything good, because they result in an unwinding of the stack
@@ -390,14 +405,14 @@ namespace aspect
  // to the particle field interpolated at these points
  cell->get_dof_indices (local_dof_indices);
  const unsigned int n_dofs_per_cell = finite_element.base_element(base_element_index).dofs_per_cell;
- for (unsigned int j=0; j<advection_fields.size(); ++j)
+ for (unsigned int j=0; j<particle_property_indices[world_index].size(); ++j)
  for (unsigned int i=0; i<n_dofs_per_cell; ++i)
  {
  const unsigned int system_local_dof
- = finite_element.component_to_system_index(advection_fields[j].component_index(introspection),
+ = finite_element.component_to_system_index(advection_fields[particle_property_indices[world_index][j].first].component_index(introspection),
  /*dof index within component=*/i);
 
- particle_solution(local_dof_indices[system_local_dof]) = particle_properties[i][particle_property_indices[j]];
+ particle_solution(local_dof_indices[system_local_dof]) = particle_properties[i][particle_property_indices[world_index][j].second];
  }
  }
  }