Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve functions to retrieve fields of certain type from introspection #5643

Merged
merged 2 commits into from
May 30, 2024
Merged

Improve functions to retrieve fields of certain type from introspection #5643

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ee12034
Add functions to retrieve fields of certain type
gassmoeller May 24, 2024
fe04d81
Check that composition type is valid
gassmoeller May 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 12 additions & 11 deletions benchmarks/entropy_adiabat/plugins/entropy_advection.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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

Expand Down
62 changes: 33 additions & 29 deletions include/aspect/introspection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Comment on lines 58 to 69
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you actually make use of these numeric values? My preference is to treat enums as symbolic names independent of any numerical representation.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I use these values to index into an std::vector below. That is why I explicitly specified them to make sure I can uniquely do the indexing. If I would just treat the enums as symbolic names I would create a utility function with a switch statement to return an unsigned int depending on type. However, I think according to the standard I am actually allowed to implicitly cast the enums to the unsigned ints.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just dislike this from a conceptual perspective. It requires keeping track of things in different places.

What would you think of keeping the enum as just a list of symbolic names, and storing data instead of a fixed-size vector into a std::map<enum_type, value_type>?

Copy link
Member Author

@gassmoeller gassmoeller May 29, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not happy about that solution because I do not know the performance implications. The cost of a lookup into a std::map already suprised us in the particle code before and the reason for this PR was to allow this lookup to be fast enough to be done multiple times per cell. What do you think about the alternative of keeping the enum as just a list of names, and having a helper function in an anonymous namespace in introspection.cc that is essentially:

unsigned int
index_of_type (enum type)
{
switch (type):
  case type_a:
    return 0;
  case type_b:
    return 1;
...
}

This way no code outside of Introspection can see that we translate the enum values into indices of a std::vector. However, it requires to always change the function if we change the enum. The beaty of the current solution is that no other code is affected if we add another enum option, everything will just work like before.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

But you still have to have two places that you keep in sync: The translation (whether in a function or by providing initializers) and the place where you store the number of enum elements.

I can be ok with the current state. Let's not make it more difficult!


/**
* 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.
Expand Down Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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> &
Comment on lines -541 to +560
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That was a bug to begin with. The fact that we marked these return types as const suggests that we always intended for them to be references.

get_names_for_fields_of_type (const CompositionalFieldDescription::Type &type) const;

/**
Expand Down Expand Up @@ -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;
};
}

Expand Down
3 changes: 2 additions & 1 deletion source/material_model/diffusion_dislocation.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down
3 changes: 2 additions & 1 deletion source/material_model/multicomponent.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion source/material_model/rheology/friction_models.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
12 changes: 7 additions & 5 deletions source/material_model/steinberger.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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. "
Expand All @@ -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
Expand Down
9 changes: 4 additions & 5 deletions source/material_model/visco_plastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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());
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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);
Expand Down
5 changes: 3 additions & 2 deletions source/material_model/viscoelastic.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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);

Expand Down
Loading
Loading