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.

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

Multicomponent entropy method - average the looked-up material properties #5622

Merged

Conversation

RanpengLi
Copy link
Contributor

Following #5615 which allows the entropy model to read in multiple files, I make the entropy model loop over look-up tables for all components, then average the material properties according to the components' volume or mass fractions.

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

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

Nice, this moves in the right direction. I think the functionality you added is fine, but I have some suggestions for the code structure. Take a look if my suggestions work for you.

/**
* The thermodynamic lookup equation of state.
*/
EquationOfState::ThermodynamicTableLookup<dim> equation_of_state;
Copy link
Member

Choose a reason for hiding this comment

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

So here you introduce a new member object into the EntropyModel class. But I dont see where you actually use this member. Is it necessary? The ThermodynamicTableLookup duplicates a lot of the functionality of the EntropyModel class, and eventually we may want to replace some parts of EntropyModel by using ThermodynamicTableLookup instead, but that doesnt seem like it happened in this PR.

@@ -163,6 +164,9 @@ namespace aspect
const unsigned int projected_density_index = this->introspection().compositional_index_for_name("density_field");
const unsigned int entropy_index = this->introspection().compositional_index_for_name("entropy");

std::vector<EquationOfStateOutputs<dim>> eos_outputs (in.n_evaluation_points(), material_file_names.size());
Copy link
Member

Choose a reason for hiding this comment

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

I think you do not need a std::vector of EquationOfStateOutputs. If you check the for-loop below, you only ever access eos_outputs[i] at a time, you never need all of eos_outputs. So you may as well make this simply a EquationOfStateOutputs<dim> eos_outputs(material_file_names.size()) and then access that single object in the loop below. This is different for the Steinberger material model (it accesses the whole vector of eos_outputs in line 278).

The same is true for the volume_fractions object in the next line.

Comment on lines 294 to 301
for (unsigned int j=0; j<material_file_names.size(); ++j)
{
eos_outputs[i].vp[j] = entropy_reader[j]->seismic_vp(entropy, pressure);
eos_outputs[i].vs[j] = entropy_reader[j]->seismic_vs(entropy,pressure);
}

seismic_out->vp[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].vp, MaterialUtilities::arithmetic);
seismic_out->vs[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].vs, MaterialUtilities::arithmetic);
Copy link
Member

Choose a reason for hiding this comment

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

So this is the place why you introduced vp and vs into EquationOfStateOutputs. But you should be able to write this just as well as:

Suggested change
for (unsigned int j=0; j<material_file_names.size(); ++j)
{
eos_outputs[i].vp[j] = entropy_reader[j]->seismic_vp(entropy, pressure);
eos_outputs[i].vs[j] = entropy_reader[j]->seismic_vs(entropy,pressure);
}
seismic_out->vp[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].vp, MaterialUtilities::arithmetic);
seismic_out->vs[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].vs, MaterialUtilities::arithmetic);
std::vector<double> vp (material_file_names.size());
std::vector<double> vs (material_file_names.size());
for (unsigned int j=0; j<material_file_names.size(); ++j)
{
vp[j] = entropy_reader[j]->seismic_vp(entropy,pressure);
vs[j] = entropy_reader[j]->seismic_vs(entropy,pressure);
}
seismic_out->vp[i] = MaterialUtilities::average_value (volume_fractions[i], vp, MaterialUtilities::arithmetic);
seismic_out->vs[i] = MaterialUtilities::average_value (volume_fractions[i], vs, MaterialUtilities::arithmetic);

Can you give this a try?

source/material_model/equation_of_state/interface.cc Outdated Show resolved Hide resolved
eos_outputs[i].densities,
true);

out.densities[i] = MaterialUtilities::average_value (volume_fractions[i], eos_outputs[i].densities, MaterialUtilities::arithmetic);
Copy link
Member

Choose a reason for hiding this comment

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

Yes, this is exactly what you want and should work just fine. Can you already add a test case that shows the functionality is working, or are there any pieces missing? I would guess that in the current state the entropy to temperature conversion would only happen correctly for the first composition table, but all other material properties are already averaged correctly?

@RanpengLi RanpengLi changed the title Multicomponent entropy method - average the looked-up material properties [WIP]Multicomponent entropy method - average the looked-up material properties May 30, 2024
@RanpengLi RanpengLi changed the title [WIP]Multicomponent entropy method - average the looked-up material properties Multicomponent entropy method - average the looked-up material properties May 30, 2024
@RanpengLi
Copy link
Contributor Author

This pull request is ready for another review. Thank you!

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

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

Looks closer. Some comments left. Let me know when you have addressed them.

source/adiabatic_conditions/compute_entropy_profile.cc Outdated Show resolved Hide resolved
source/adiabatic_conditions/compute_entropy_profile.cc Outdated Show resolved Hide resolved
source/material_model/entropy_model.cc Outdated Show resolved Hide resolved
source/material_model/equation_of_state/interface.cc Outdated Show resolved Hide resolved
source/material_model/entropy_model.cc Outdated Show resolved Hide resolved
source/material_model/entropy_model.cc Outdated Show resolved Hide resolved
@RanpengLi RanpengLi force-pushed the read-multi-table-eos-output branch 2 times, most recently from 038b214 to 731d9c3 Compare May 31, 2024 23:24
@RanpengLi
Copy link
Contributor Author

RanpengLi commented May 31, 2024

This pull request is ready for another review. Thank you again!

@gassmoeller
Copy link
Member

You need to run make indent

@RanpengLi RanpengLi force-pushed the read-multi-table-eos-output branch from 4333d13 to 52e41df Compare May 31, 2024 23:59
@RanpengLi
Copy link
Contributor Author

Done!

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

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

Ok, one more small change that I marked up (sorry I didnt see that in the earlier reviews). Afterwards this is ready to merge.

@@ -161,7 +164,18 @@ namespace aspect
MaterialModel::MaterialModelOutputs<dim> &out) const
{
const unsigned int projected_density_index = this->introspection().compositional_index_for_name("density_field");
const unsigned int entropy_index = this->introspection().compositional_index_for_name("entropy");
//TODO : need to make it work for more than one field
const std::vector<unsigned int> entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
const std::vector<unsigned int> entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);
const std::vector<unsigned int> &entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);

//TODO : need to make it work for more than one field
const std::vector<unsigned int> entropy_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::entropy);
const unsigned int entropy_index = entropy_indices[0];
const std::vector<unsigned int> composition_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition);
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
const std::vector<unsigned int> composition_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition);
const std::vector<unsigned int> &composition_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::chemical_composition);

@RanpengLi RanpengLi force-pushed the read-multi-table-eos-output branch from 52e41df to e0e6d8b Compare June 1, 2024 02:16
@gassmoeller gassmoeller merged commit ecc67a9 into geodynamics:main Jun 1, 2024
8 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants