diff --git a/benchmarks/entropy_adiabat/plugins/entropy_advection.cc b/benchmarks/entropy_adiabat/plugins/entropy_advection.cc index c7256af0c1b..368d91378d8 100644 --- a/benchmarks/entropy_adiabat/plugins/entropy_advection.cc +++ b/benchmarks/entropy_adiabat/plugins/entropy_advection.cc @@ -41,7 +41,7 @@ namespace aspect const FiniteElement &fe = this->get_fe(); const typename Simulator::AdvectionField advection_field = *scratch.advection_field; - const std::vector composition_descriptions = this->introspection().get_composition_descriptions(); + const std::vector &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 residuals(n_q_points,0.0); - const std::vector composition_descriptions = this->introspection().get_composition_descriptions(); + const std::vector &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.advection_system_assembler_properties[temperature_index].needed_update_flags = update_hessians; - const std::vector composition_descriptions = simulator_access.introspection().get_composition_descriptions(); - for (unsigned int c=0; c>()); - assemblers.advection_system_assembler_properties[entropy_index].needed_update_flags = update_hessians; - } + if (this->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 = this->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.advection_system_assembler_properties[entropy_index].needed_update_flags = update_hessians; + } } } // namespace aspect diff --git a/include/aspect/introspection.h b/include/aspect/introspection.h index 2224fd8a0c0..531dfc9a10f 100644 --- a/include/aspect/introspection.h +++ b/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,6 +485,9 @@ 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 & chemical_composition_field_names () const; @@ -487,6 +495,9 @@ namespace aspect /** * 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 & chemical_composition_field_indices () const; @@ -494,6 +505,9 @@ namespace aspect /** * 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 + const std::vector & 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 + const std::vector & get_names_for_fields_of_type (const CompositionalFieldDescription::Type &type) const; /** @@ -580,26 +592,18 @@ namespace aspect std::vector 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 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 chemical_composition_indices; + std::vector> 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> composition_indices_for_type; }; } diff --git a/source/material_model/diffusion_dislocation.cc b/source/material_model/diffusion_dislocation.cc index 34d589725aa..093fea7880a 100644 --- a/source/material_model/diffusion_dislocation.cc +++ b/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()); diff --git a/source/material_model/multicomponent.cc b/source/material_model/multicomponent.cc index 70a88b9fa0d..66040fe62ac 100644 --- a/source/material_model/multicomponent.cc +++ b/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 volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i], this->introspection().chemical_composition_field_indices()); + const std::vector 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); diff --git a/source/material_model/rheology/friction_models.cc b/source/material_model/rheology/friction_models.cc index d5dc45fd9f0..bb2d49361bb 100644 --- a/source/material_model/rheology/friction_models.cc +++ b/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(point.get_coordinates()),index); diff --git a/source/material_model/steinberger.cc b/source/material_model/steinberger.cc index bc355464696..061f7c649e3 100644 --- a/source/material_model/steinberger.cc +++ b/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 diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index 92f3a1ce198..cef81e360e1 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -133,7 +133,7 @@ namespace aspect evaluate(const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out) const { - EquationOfStateOutputs eos_outputs (this->introspection().n_chemical_composition_fields()+1); + EquationOfStateOutputs eos_outputs (this->introspection().get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)+1); EquationOfStateOutputs eos_outputs_all_phases (n_phases); std::vector 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 volume_fractions = MaterialUtilities::compute_only_composition_fractions(in.composition[i], this->introspection().chemical_composition_field_indices()); + const std::vector 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,8 +386,6 @@ namespace aspect std::vector n_phases_for_each_composition = phase_function.n_phases_for_each_composition(); - const std::vector 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 @@ -394,7 +393,7 @@ namespace aspect std::vector 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); diff --git a/source/material_model/viscoelastic.cc b/source/material_model/viscoelastic.cc index af359ba94b2..f25c6375a4a 100644 --- a/source/material_model/viscoelastic.cc +++ b/source/material_model/viscoelastic.cc @@ -34,7 +34,7 @@ namespace aspect evaluate(const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out) const { - EquationOfStateOutputs eos_outputs (this->introspection().n_chemical_composition_fields()+1); + EquationOfStateOutputs 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 composition = in.composition[i]; - const std::vector volume_fractions = MaterialUtilities::compute_only_composition_fractions(composition, this->introspection().chemical_composition_field_indices()); + const std::vector volume_fractions = MaterialUtilities::compute_only_composition_fractions(composition, + this->introspection().chemical_composition_field_indices()); equation_of_state.evaluate(in, i, eos_outputs); diff --git a/source/simulator/introspection.cc b/source/simulator/introspection.cc index b7ee02694e6..e7c0e0d220d 100644 --- a/source/simulator/introspection.cc +++ b/source/simulator/introspection.cc @@ -256,10 +256,16 @@ namespace aspect compositional_field_methods(parameters.compositional_field_methods), composition_names(parameters.names_of_compositional_fields), composition_descriptions(parameters.composition_descriptions), - chemical_composition_names (get_names_for_fields_of_type(CompositionalFieldDescription::chemical_composition)), - chemical_composition_indices (get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition)), - n_chemical_compositions (get_number_of_fields_of_type(CompositionalFieldDescription::chemical_composition)) - {} + composition_names_for_type(CompositionalFieldDescription::n_types), + composition_indices_for_type(CompositionalFieldDescription::n_types) + { + // Set up the indices for the different types of compositional fields + for (unsigned int c=0; c & Introspection::chemical_composition_field_names () const { - return chemical_composition_names; + return get_names_for_fields_of_type(CompositionalFieldDescription::chemical_composition); } @@ -381,7 +387,7 @@ namespace aspect const std::vector & Introspection::chemical_composition_field_indices () const { - return chemical_composition_indices; + return get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition); } @@ -390,7 +396,7 @@ namespace aspect unsigned int Introspection::n_chemical_composition_fields () const { - return n_chemical_compositions; + return get_number_of_fields_of_type(CompositionalFieldDescription::Type::chemical_composition); } @@ -399,10 +405,7 @@ namespace aspect bool Introspection::composition_type_exists (const CompositionalFieldDescription::Type &type) const { - for (unsigned int c=0; c 0; } @@ -411,9 +414,9 @@ namespace aspect unsigned int Introspection::find_composition_type (const typename CompositionalFieldDescription::Type &type) const { - for (unsigned int c=0; c 0) + return composition_indices_for_type[type][0]; + return composition_descriptions.size(); } @@ -433,31 +436,19 @@ namespace aspect template - const std::vector + const std::vector & Introspection::get_indices_for_fields_of_type (const CompositionalFieldDescription::Type &type) const { - std::vector indices; - - for (unsigned int i=0; i - const std::vector + const std::vector & Introspection::get_names_for_fields_of_type (const CompositionalFieldDescription::Type &type) const { - std::vector names; - - for (unsigned int i=0; i::get_number_of_fields_of_type (const CompositionalFieldDescription::Type &type) const { - unsigned int n = 0; - - for (unsigned int i=0; i