Add functions to retrieve fields of certain type

benchmarks/entropy_adiabat/plugins/entropy_advection.cc

@@ -41,7 +41,7 @@ namespace aspect
  const FiniteElement<dim> &fe = this->get_fe();
 
  const typename Simulator<dim>::AdvectionField advection_field = *scratch.advection_field;
- const std::vector<CompositionalFieldDescription> composition_descriptions = this->introspection().get_composition_descriptions();
+ const std::vector<CompositionalFieldDescription> &composition_descriptions = this->introspection().get_composition_descriptions();
  if (!advection_field.is_temperature()
  && composition_descriptions[advection_field.compositional_variable].type != CompositionalFieldDescription::entropy)
  return;
@@ -184,7 +184,7 @@ namespace aspect
  const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
  std::vector<double> residuals(n_q_points,0.0);
 
- const std::vector<CompositionalFieldDescription> composition_descriptions = this->introspection().get_composition_descriptions();
+ const std::vector<CompositionalFieldDescription> &composition_descriptions = this->introspection().get_composition_descriptions();
  if (!advection_field.is_temperature()
  && composition_descriptions[advection_field.compositional_variable].type != CompositionalFieldDescription::entropy)
  return residuals;
@@ -277,15 +277,16 @@ namespace aspect
  assemblers.advection_system[temperature_index].emplace_back (std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>());
  assemblers.advection_system_assembler_properties[temperature_index].needed_update_flags = update_hessians;
 
- const std::vector<CompositionalFieldDescription> composition_descriptions = simulator_access.introspection().get_composition_descriptions();
- for (unsigned int c=0; c<composition_descriptions.size(); ++c)
- if (composition_descriptions[c].type == CompositionalFieldDescription::entropy)
- {
- const unsigned int entropy_index = c + 1;
- assemblers.advection_system[entropy_index].clear();
- assemblers.advection_system[entropy_index].emplace_back (std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>());
- assemblers.advection_system_assembler_properties[entropy_index].needed_update_flags = update_hessians;
- }
+ if (simulator_access.introspection().composition_type_exists(CompositionalFieldDescription::entropy))
+ {
+ // Find the index of the entropy field and replace the assembler for it.
+ // The index of the entropy field is its index in the compositional fields plus one (for the temperature field).
+ const unsigned int entropy_index = simulator_access.introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy)[0]
+ + 1;
+ assemblers.advection_system[entropy_index].clear();
+ assemblers.advection_system[entropy_index].emplace_back (std::make_unique<Assemblers::EntropyAdvectionSystem<dim>>());
+ assemblers.advection_system_assembler_properties[entropy_index].needed_update_flags = update_hessians;
+ }
  }
 } // namespace aspect
 

include/aspect/introspection.h

@@ -57,17 +57,22 @@ namespace aspect
  */
  enum Type
  {
- chemical_composition,
- stress,
- strain,
- grain_size,
- porosity,
- density,
- entropy,
- generic,
- unspecified
+ chemical_composition = 0,
+ stress = 1,
+ strain = 2,
+ grain_size = 3,
+ porosity = 4,
+ density = 5,
+ entropy = 6,
+ generic = 7,
+ unspecified = 8
  } type;
 
+ /**
+ * The number of different types defined in Type.
+ */
+ constexpr static unsigned int n_types = 9;
+
  /**
  * This function translates an input string into the
  * available enum options for the type of compositional field.
@@ -480,20 +485,29 @@ namespace aspect
  /**
  * A function that returns the names of
  * compositional fields that correspond to chemical compositions.
+ *
+ * This function is shorthand for
+ * get_names_for_fields_of_type(CompositionalFieldDescription::chemical_composition).
  */
  const std::vector<std::string> &
  chemical_composition_field_names () const;
 
  /**
  * A function that returns the indices of
  * compositional fields that correspond to chemical compositions.
+ *
+ * This function is shorthand for
+ * get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition).
  */
  const std::vector<unsigned int> &
  chemical_composition_field_indices () const;
 
  /**
  * A function that returns the number of
  * compositional fields that correspond to chemical compositions.
+ *
+ * This function is shorthand for
+ * get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition).
  */
  unsigned int
  n_chemical_composition_fields () const;
@@ -509,7 +523,6 @@ namespace aspect
  bool
  composition_type_exists (const CompositionalFieldDescription::Type &type) const;
 
-
  /**
  * A function that gets the type of a compositional field as an input
  * parameter and returns the index of the first compositional field of
@@ -522,7 +535,6 @@ namespace aspect
  unsigned int
  find_composition_type (const CompositionalFieldDescription::Type &type) const;
 
-
  /**
  * A function that gets the name of a compositional field as an input
  * parameter and returns if the compositional field is used in this
@@ -538,14 +550,14 @@ namespace aspect
  * Get the indices of the compositional fields which are of a
  * particular type (chemical composition, porosity, etc.).
  */
- const std::vector<unsigned int>
+ const std::vector<unsigned int> &
  get_indices_for_fields_of_type (const CompositionalFieldDescription::Type &type) const;
 
  /**
  * Get the names of the compositional fields which are of a
  * particular type (chemical composition, porosity, etc.).
  */
- const std::vector<std::string>
+ const std::vector<std::string> &
  get_names_for_fields_of_type (const CompositionalFieldDescription::Type &type) const;
 
  /**
@@ -580,26 +592,18 @@ namespace aspect
  std::vector<CompositionalFieldDescription> composition_descriptions;
 
  /**
- * A vector that stores the names of the compositional fields
- * corresponding to chemical compositions that will
- * be used in the simulation.
- */
- std::vector<std::string> chemical_composition_names;
-
- /**
- * A vector that stores the indices of the compositional fields
- * corresponding to chemical compositions that will
- * be used in the simulation.
+ * A vector of vectors of composition names that stores the
+ * names of the compositional fields corresponding to each field type
+ * given in CompositionalFieldDescription.
  */
- std::vector<unsigned int> chemical_composition_indices;
+ std::vector<std::vector<std::string>> composition_names_for_type;
 
  /**
- * The number of compositional fields
- * that correspond to chemical compositions that will
- * be used in the simulation.
+ * A vector of vectors of composition indices that stores the
+ * indices of the compositional fields corresponding to each field type
+ * given in CompositionalFieldDescription.
  */
- const unsigned int n_chemical_compositions;
-
+ std::vector<std::vector<unsigned int>> composition_indices_for_type;
  };
 }
 

source/material_model/diffusion_dislocation.cc

@@ -454,7 +454,8 @@ namespace aspect
 
  // Rheological parameters for chemical compositions
  // increment by one for background:
- const unsigned int n_chemical_composition_fields = this->introspection().n_chemical_composition_fields() + 1;
+ const unsigned int n_chemical_composition_fields =
+ this->introspection().n_chemical_composition_fields() + 1;
 
  // Diffusion creep parameters
  diffusion_creep.initialize_simulator (this->get_simulator());

source/material_model/multicomponent.cc

@@ -42,7 +42,8 @@ namespace aspect
  // The (incompressible) Boussinesq approximation treats the
  // buoyancy term as Delta rho[i] * C[i], which implies that
  // compositional fields are given as volume fractions.
- const std::vector<double> volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i], this->introspection().chemical_composition_field_indices());
+ const std::vector<double> volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i],
+ this->introspection().chemical_composition_field_indices());
 
  equation_of_state.evaluate(in, i, eos_outputs);
  out.viscosities[i] = MaterialUtilities::average_value(volume_fractions, viscosities, viscosity_averaging);

source/material_model/rheology/friction_models.cc

@@ -97,7 +97,7 @@ namespace aspect
  // If chemical fields are present, volume_fractions will be of size 1+n_chemical_composition_fields.
  // The size of chemical_composition_field_indices will be one less.
  unsigned int index = 0;
- if (this->introspection().n_chemical_composition_fields() > 0)
+ if (this->introspection().composition_type_exists(CompositionalFieldDescription::chemical_composition))
  index = this->introspection().chemical_composition_field_indices()[volume_fraction_index-1];
  double friction_from_function =
  friction_function->value(Utilities::convert_array_to_point<dim>(point.get_coordinates()),index);

source/material_model/steinberger.cc

@@ -291,7 +291,8 @@ namespace aspect
  else
  {
  // We only want to compute mass/volume fractions for fields that are chemical compositions.
- mass_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i], this->introspection().chemical_composition_field_indices());
+ mass_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i],
+ this->introspection().chemical_composition_field_indices());
 
  // The function compute_volumes_from_masses expects as many mass_fractions as densities.
  // But the function compute_composition_fractions always adds another element at the start
@@ -534,11 +535,12 @@ namespace aspect
  equation_of_state.parse_parameters(prm);
 
  // Assign background field and do some error checking
+ const unsigned int n_chemical_composition_fields = this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition);
  has_background_field = ((equation_of_state.number_of_lookups() == 1) ||
- (equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields() + 1));
+ (equation_of_state.number_of_lookups() == n_chemical_composition_fields + 1));
  AssertThrow ((equation_of_state.number_of_lookups() == 1) ||
- (equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields()) ||
- (equation_of_state.number_of_lookups() == this->introspection().n_chemical_composition_fields() + 1),
+ (equation_of_state.number_of_lookups() == n_chemical_composition_fields) ||
+ (equation_of_state.number_of_lookups() == n_chemical_composition_fields + 1),
  ExcMessage("The Steinberger material model assumes that either there is a single material "
  "in the simulation, or that all compositional fields of the type "
  "chemical composition correspond to mass fractions of different materials. "
@@ -547,7 +549,7 @@ namespace aspect
  "or one additional file (if a background field is used). You have "
  + Utilities::int_to_string(equation_of_state.number_of_lookups())
  + " material data files, but there are "
- + Utilities::int_to_string(this->introspection().n_chemical_composition_fields())
+ + Utilities::int_to_string(n_chemical_composition_fields)
  + " fields of type chemical composition."));
 
  // Parse the Composition viscosity prefactors parameter

source/material_model/visco_plastic.cc

@@ -133,7 +133,7 @@ namespace aspect
  evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
  MaterialModel::MaterialModelOutputs<dim> &out) const
  {
- EquationOfStateOutputs<dim> eos_outputs (this->introspection().n_chemical_composition_fields()+1);
+ EquationOfStateOutputs<dim> eos_outputs (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)+1);
  EquationOfStateOutputs<dim> eos_outputs_all_phases (n_phases);
 
  std::vector<double> average_elastic_shear_moduli (in.n_evaluation_points());
@@ -176,7 +176,8 @@ namespace aspect
  n_phase_transitions_for_each_chemical_composition,
  eos_outputs);
 
- const std::vector<double> volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i], this->introspection().chemical_composition_field_indices());
+ const std::vector<double> volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i],
+ this->introspection().chemical_composition_field_indices());
 
  // not strictly correct if thermal expansivities are different, since we are interpreting
  // these compositions as volume fractions, but the error introduced should not be too bad.
@@ -385,16 +386,14 @@ namespace aspect
 
  std::vector<unsigned int> n_phases_for_each_composition = phase_function.n_phases_for_each_composition();
 
- const std::vector<unsigned int> indices = this->introspection().chemical_composition_field_indices();
-
  // Currently, phase_function.n_phases_for_each_composition() returns a list of length
  // equal to the total number of compositions, whether or not they are chemical compositions.
  // The equation_of_state (multicomponent incompressible) requires a list only for
  // chemical compositions.
  std::vector<unsigned int> n_phases_for_each_chemical_composition = {n_phases_for_each_composition[0]};
  n_phase_transitions_for_each_chemical_composition = {n_phases_for_each_composition[0] - 1};
  n_phases = n_phases_for_each_composition[0];
- for (auto i : indices)
+ for (auto i : this->introspection().chemical_composition_field_indices())
  {
  n_phases_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1]);
  n_phase_transitions_for_each_chemical_composition.push_back(n_phases_for_each_composition[i+1] - 1);

source/material_model/viscoelastic.cc

@@ -34,7 +34,7 @@ namespace aspect
  evaluate(const MaterialModel::MaterialModelInputs<dim> &in,
  MaterialModel::MaterialModelOutputs<dim> &out) const
  {
- EquationOfStateOutputs<dim> eos_outputs (this->introspection().n_chemical_composition_fields()+1);
+ EquationOfStateOutputs<dim> eos_outputs (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)+1);
 
  // Store which components to exclude during volume fraction computation.
  ComponentMask composition_mask(this->n_compositional_fields(), true);
@@ -50,7 +50,8 @@ namespace aspect
  {
  const std::vector<double> composition = in.composition[i];
 
- const std::vector<double> volume_fractions = MaterialUtilities::compute_only_composition_fractions(composition, this->introspection().chemical_composition_field_indices());
+ const std::vector<double> volume_fractions = MaterialUtilities::compute_only_composition_fractions(composition,
+ this->introspection().chemical_composition_field_indices());
 
  equation_of_state.evaluate(in, i, eos_outputs);