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

Add compositional dependence on viscosity in the Steinberger material model #5260

Conversation

poulamiiroy
Copy link
Contributor

This pull request supersedes the PR #5248.

Added compositional dependence on viscosity in the Steinberger material model. The original Steinberger material model uses one viscosity profile for all of the compositions. So, I have added prefectors which are multiplied by original viscosity based on volume fraction.

For all pull requests:

Copy link
Contributor

@naliboff naliboff left a comment

Choose a reason for hiding this comment

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

@poulamiiroy - Thank you for the contribution and breaking the original PR into smaller contributions! I think with a few moderate adjustments the PR will be in good shape. Please let me know if you have any questions about my comments, and thank you again for the contribution!

@@ -274,6 +274,7 @@ namespace aspect
std::vector<double> conductivity_reference_temperatures;
std::vector<double> conductivity_exponents;
std::vector<double> saturation_scaling;
std::vector<double> prefactors;
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
std::vector<double> prefactors;
std::vector<double> viscosity_prefactors;

Can you also move this variable down a few lines and add documentation that described what it represents. The variable is currently lumped with other variables related to thermal conductivity.

const SymmetricTensor<2,dim> &,
const Point<dim> &position) const
{
const double depth = this->get_geometry_model().depth(position);
const double adiabatic_temperature = this->get_adiabatic_conditions().temperature(position);
// Calculate volume fractions from mass fractions
Copy link
Contributor

Choose a reason for hiding this comment

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

Can the following four lines be removed?

@@ -189,11 +194,14 @@ namespace aspect
// For an explanation on this formula see the Steinberger & Calderwood 2006 paper
const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
// Limit the lateral viscosity variation to a reasonable interval

Copy link
Contributor

Choose a reason for hiding this comment

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

I would remove this extra line, and instead add it between lines 195 and 196.

const double vis_lateral = std::max(std::min(std::exp(vis_lateral_exp),max_lateral_eta_variation),1/max_lateral_eta_variation);

const double vis_radial = radial_viscosity_lookup->radial_viscosity(depth);
const double vis_compositional = MaterialUtilities::average_value (volume_fractions, prefactors, viscosity_averaging);
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
const double vis_compositional = MaterialUtilities::average_value (volume_fractions, prefactors, viscosity_averaging);
// Volume fraction averaged value of the composition-dependent viscosity prefactors
const double vis_prefactor = MaterialUtilities::average_value (volume_fractions, prefactors, viscosity_averaging);

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I wanted to keep it as vis_compositional since I am already renaming "prefactors" as "viscosity_prefactors" as you have suggested. What do you think?

@@ -268,11 +276,16 @@ namespace aspect

// Evaluate the equation of state properties over all evaluation points
equation_of_state.evaluate(eos_in, eos_outputs);
// std::vector<Vector<double>> viscosities(in.n_evaluation_points(), this->n_compositional_fields());
Copy link
Contributor

Choose a reason for hiding this comment

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

Can you remove this and additional commented lines below if they are no longer needed?


return std::max(std::min(vis_lateral * vis_radial,max_eta),min_eta);
return std::max(std::min(vis_lateral * vis_radial * vis_compositional,max_eta),min_eta);
Copy link
Contributor

Choose a reason for hiding this comment

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

I confess I don't understand this equation, as it looks like viscosity values are being multiplied together. Would it be possible to add the equation here (or further above) in a comment, along with a brief description?

Copy link
Member

Choose a reason for hiding this comment

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

Agreed, as written it looks like multiplication of different viscosities. Would you mind renaming

  • vis_lateral to temperature_prefactor
  • vis_compositional to compositional_prefactor
  • vis_radial to eta_ref?

Copy link
Contributor

Choose a reason for hiding this comment

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

I second @bobmyhill renaming scheme for the values

"geometric, or maximum composition.");
prm.declare_entry ("Prefactors", "1.e21",
Patterns::Anything(),
"List of viscosities for background mantle and compositional fields,"
Copy link
Contributor

Choose a reason for hiding this comment

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

Just to be sure, are the "Prefactors" in fact viscosity values, or are they dimensionless composition-dependent factors that the viscosity is multiplied by?

Copy link
Contributor

Choose a reason for hiding this comment

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

Follow up - depending on the answer, I would adjust the parameter name and description to be a bit more specific.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Prefactors are the dimensionless composition-dependent factors, for example, if one composition needs to have 10 times higher viscosity than ambient, then prefactor for this composition would be 10.


viscosity_averaging = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme",
prm);
// viscosities = Utilities::parse_map_to_double_array (prm.get("Viscosities"),
Copy link
Contributor

Choose a reason for hiding this comment

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

Remove commented lines?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Thanks a lot John for reviewing my pull request. I will follow what you have suggested. Please let me know for further discussion.

Copy link
Contributor

@jdannberg jdannberg left a comment

Choose a reason for hiding this comment

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

This looks goof to me, thank you for contributing this! Can you address the few small comments I left, and then I think this is good to be merged!

Comment on lines 284 to 285
// if (in.requests_property(MaterialProperties::viscosity))
// out.viscosities[i] = viscosity(in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
Copy link
Contributor

Choose a reason for hiding this comment

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

If these lines are not needed anymore, simply remove them.

Comment on lines 320 to 321
// out.viscosities[i] = MaterialUtilities::average_value (volume_fractions[i], viscosities, viscosity_averaging);
// out.viscosities[i] = viscosities[i];
Copy link
Contributor

Choose a reason for hiding this comment

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

Same here, these lines can be removed.

Comment on lines 456 to 471
"model. Larger values will be cut off.");
"model. Larger values will be cut off."),
Copy link
Contributor

Choose a reason for hiding this comment

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

I do not think you want to change this. This actually needs to be a semicolon.

Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"),
"When more than one compositional field is present at a point "
"with different viscosities, we need to come up with an average "
"viscosity at that point. Select a weighted harmonic, arithmetic, "
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
"viscosity at that point. Select a weighted harmonic, arithmetic, "
"viscosity at that point. Select a weighted harmonic, arithmetic, "

Comment on lines 514 to 517
// viscosities = Utilities::parse_map_to_double_array (prm.get("Viscosities"),
// list_of_composition_names,
// has_background_field,
// "Viscosities");
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
// viscosities = Utilities::parse_map_to_double_array (prm.get("Viscosities"),
// list_of_composition_names,
// has_background_field,
// "Viscosities");

I think John is just asking to remove these lines.

subsection Steinberger model
set Use lateral average temperature for viscosity = false
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
# set Derivatives file names = /home/poulami/software/aspect_30_09_2022/aspect/data/material-model/steinberger/test-steinberger-compressible/constant_material_small_hefesto_derivatives.txt
Copy link
Contributor

Choose a reason for hiding this comment

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

If you're not using that line, simply remove it.

set Use lateral average temperature for viscosity = false
set Data directory = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
# set Derivatives file names = /home/poulami/software/aspect_30_09_2022/aspect/data/material-model/steinberger/test-steinberger-compressible/constant_material_small_hefesto_derivatives.txt
set Radial viscosity file name = radial-visc.txt #, radial-visc.txt, radial-visc.txt ## lithosphere at top
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
set Radial viscosity file name = radial-visc.txt #, radial-visc.txt, radial-visc.txt ## lithosphere at top
set Radial viscosity file name = radial-visc.txt # lithosphere at top

Same here, remove everything that is not used.

set Data directory = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
# set Derivatives file names = /home/poulami/software/aspect_30_09_2022/aspect/data/material-model/steinberger/test-steinberger-compressible/constant_material_small_hefesto_derivatives.txt
set Radial viscosity file name = radial-visc.txt #, radial-visc.txt, radial-visc.txt ## lithosphere at top
set Lateral viscosity file name = temp-viscosity-prefactor.txt ##, temp-viscosity-prefactor.txt, temp-viscosity-prefactor.txt
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
set Lateral viscosity file name = temp-viscosity-prefactor.txt ##, temp-viscosity-prefactor.txt, temp-viscosity-prefactor.txt
set Lateral viscosity file name = temp-viscosity-prefactor.txt

Copy link
Contributor

Choose a reason for hiding this comment

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

It does not look like you are actually using this file. Is that intentional? If you do not want to use it in your test case below, simply remove it from your pull request.

Comment on lines 192 to 193
//The lateral variation of viscosity due to lateral temperature (Visc_lT)
//Visc_lT = exp [(-1)*(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT)], Eq. 6 of the paper
Copy link
Contributor

Choose a reason for hiding this comment

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

I think what @naliboff was asking for was a more detailed explanation of why these things are multiplied together. I have made a suggestion below:

Suggested change
//The lateral variation of viscosity due to lateral temperature (Visc_lT)
//Visc_lT = exp [(-1)*(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT)], Eq. 6 of the paper
// We here compute the lateral variation of viscosity due to temperature (vis_lateral) as
// V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))] as in Eq. 6 of the paper.
// We get H/nR from the lateral_viscosity_lookup->lateral_viscosity function.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

@jdannberg Thanks a lot for your review! I made the changes.

Copy link
Contributor

@jdannberg jdannberg left a comment

Choose a reason for hiding this comment

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

It looks like many of the tests that use the steinberger material model are failing now. They all have the error message

An error occurred in line <943> of file </__w/aspect/aspect/source/material_model/utilities.cc> in function
    double aspect::MaterialModel::MaterialUtilities::average_value(const std::vector<double, std::allocator<double> >&, const std::vector<double, std::allocator<double> >&, const aspect::MaterialModel::MaterialUtilities::CompositionalAveragingOperation&)
The violated condition was: 
    volume_fractions.size() == parameter_values.size()
Additional information: 
    The volume fractions and parameter values vectors used for averaging
    have to have the same length!

To figure out what might be the problem, I would suggest to run these tests locally on your machine. For example, you could run them in the debugger to see what the two values are (volume fractions and viscosity prefactors) and which of these actually has the wrong length.

Copy link
Contributor

Choose a reason for hiding this comment

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

You will need to update this test output with the one from the official tester (and same for the statistics file). The file that you added here still, for example, includes the output path on your own computer (Writing graphical output: /home/poulami/software/Aspect_globalscale/aspect/tests/multicomponent_steinberger/solution/solution-00001)

@jdannberg
Copy link
Contributor

/rebuild

Copy link
Member

@bobmyhill bobmyhill left a comment

Choose a reason for hiding this comment

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

A couple of suggested changes to variable names that will make code comprehension easier.


return std::max(std::min(vis_lateral * vis_radial,max_eta),min_eta);
return std::max(std::min(vis_lateral * vis_radial * vis_compositional,max_eta),min_eta);
Copy link
Member

Choose a reason for hiding this comment

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

Agreed, as written it looks like multiplication of different viscosities. Would you mind renaming

  • vis_lateral to temperature_prefactor
  • vis_compositional to compositional_prefactor
  • vis_radial to eta_ref?

@@ -187,13 +189,22 @@ namespace aspect
delta_temperature = temperature-adiabatic_temperature;

// For an explanation on this formula see the Steinberger & Calderwood 2006 paper
// We here compute the lateral variation of viscosity due to temperature (vis_lateral) as
// V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))] as in Eq. 6 of the paper.
// We get H/nR from the lateral_viscosity_lookup->lateral_viscosity function.
const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
Copy link
Member

Choose a reason for hiding this comment

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

Please rename vis_lateral_exp to log_temperature_prefactor

Copy link
Contributor

Choose a reason for hiding this comment

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

A follow-up question on these lines.

The equation in the comments is written as V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))], while the actual equation on line 195 has delta_temperature/(temperature*adiabatic_temperature).

Is the equation in the comments referring to the same equation as the one on line 195? If yes, can you confirm which one is correct or expand the documentation a bit further to explain the differences between the two equations?

Copy link
Contributor

Choose a reason for hiding this comment

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

@naliboff I suggested to write the equation like this, this is how it is written in the Steinberger and Calderwood paper.
temperature is the same as T_adiabatic + dT

Copy link
Contributor

@naliboff naliboff left a comment

Choose a reason for hiding this comment

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

@poulamiiroy - A few additional comments on top of those from @bobmyhill.


return std::max(std::min(vis_lateral * vis_radial,max_eta),min_eta);
return std::max(std::min(vis_lateral * vis_radial * vis_compositional,max_eta),min_eta);
Copy link
Contributor

Choose a reason for hiding this comment

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

I second @bobmyhill renaming scheme for the values

@@ -187,13 +189,22 @@ namespace aspect
delta_temperature = temperature-adiabatic_temperature;

// For an explanation on this formula see the Steinberger & Calderwood 2006 paper
// We here compute the lateral variation of viscosity due to temperature (vis_lateral) as
// V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))] as in Eq. 6 of the paper.
// We get H/nR from the lateral_viscosity_lookup->lateral_viscosity function.
const double vis_lateral_exp = -1.0*lateral_viscosity_lookup->lateral_viscosity(depth)*delta_temperature/(temperature*adiabatic_temperature);
Copy link
Contributor

Choose a reason for hiding this comment

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

A follow-up question on these lines.

The equation in the comments is written as V_lT = exp [-(H/nR)*dT/(T_adiabatic*(T_adiabatic + dT))], while the actual equation on line 195 has delta_temperature/(temperature*adiabatic_temperature).

Is the equation in the comments referring to the same equation as the one on line 195? If yes, can you confirm which one is correct or expand the documentation a bit further to explain the differences between the two equations?

Copy link
Contributor

@jdannberg jdannberg left a comment

Choose a reason for hiding this comment

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

Okay, the code looks much better now and most of the tests are succeeding as well! 👍

But there is still one test that is failing (steinberger_phase_fractions_background). The test finishes now, but the results are different, so you'll need to figure out why the output is different (i.e. which column in the graphical output is different? Is it the viscosity? What has changed compared to before?).

In addition, you have accidentally added some lines in your last commit that are just printing statements you commented out.

Comment on lines +203 to +204
// std::cout<<"size vf"<<volume_fractions.size()<<std::endl;
// std::cout<<"size vis_pre"<<viscosity_prefactors.size()<<std::endl;
Copy link
Contributor

Choose a reason for hiding this comment

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

You added some lines again here that are code you commented out and that is no longer needed. Please remove them.

// std::cout<<"size vf"<<volume_fractions.size()<<std::endl;
// std::cout<<"size vis_pre"<<viscosity_prefactors.size()<<std::endl;
const double compositional_prefactor = MaterialUtilities::average_value (volume_fractions, viscosity_prefactors, viscosity_averaging);
// std::cout<<"size vf"<<volume_fractions.size()<<std::endl;
Copy link
Contributor

Choose a reason for hiding this comment

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

Same here, please remove this line.

@@ -303,7 +305,7 @@ namespace aspect
chemical_compositions.push_back(in.composition[i][c]);

mass_fractions = MaterialUtilities::compute_composition_fractions(chemical_compositions, *composition_mask);

// std::cout<<"mf"<<mass_fractions<<std::endl;
Copy link
Contributor

Choose a reason for hiding this comment

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

Please remove this line as well.

@bobmyhill
Copy link
Member

@poulamiiroy Ignore my last; thank you for making the variable name changes we suggested.
Do you need any help with figuring out why the tests are failing?

@@ -551,7 +580,15 @@ namespace aspect
+ " fields of type chemical composition."));

has_background_field = (equation_of_state.number_of_lookups() == n_chemical_fields + 1);

std::vector<std::string> list_of_composition_names = this->introspection().get_composition_names();
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 want fields that correspond to chemical compositions, rather than all compositional fields:

Suggested change
std::vector<std::string> list_of_composition_names = this->introspection().get_composition_names();
std::vector<std::string> list_of_composition_names = this->introspection().chemical_composition_field_names();

This might fix a couple of the tests.

if (equation_of_state.number_of_lookups() == 1)
list_of_composition_names.resize(1, "background");
if ((equation_of_state.number_of_lookups() == 1) && (has_background_field))
list_of_composition_names.resize(0, "background");
Copy link
Member

Choose a reason for hiding this comment

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

Are you sure that these lines cover all the potential use cases of the material model?

The way you have implemented this, if the number of lookups is equal to 1, the vector of composition names is replaced by {"background"}. If additionally a background field is implied, the vector is empty {}.

What happens when the number of lookups is not equal to 1? You should test that parse_map_to_double_array does the correct thing in these cases.

You might be also be interested in mirroring the changes I made here, using the new map parser:
https://github.com/geodynamics/aspect/pull/5039/files#diff-cfe301f322bebee87432e53137677e6372743e5e142733d8fb63b5599c9322cb

@poulamiiroy
Copy link
Contributor Author

Thank you for your comments! I am working on it. I'll ping you when the PR needs a new review. :)

// if ((equation_of_state.number_of_lookups() > 1) && (has_background_field))
// list_of_composition_names.resize(equation_of_state.number_of_lookups(), "background");
if (equation_of_state.number_of_lookups() > 1)
list_of_composition_names.resize(n_chemical_fields, "background");
Copy link
Member

Choose a reason for hiding this comment

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

The line that you've added here appends to the list_of_composition_names vector with entries called "background" (e.g. {"a", "b", "background", "background", "background", "background", ...}.

I doubt that is what you want to do. The background field is the first compositional field, not the last. And there is only ever one of them.

Have another look at the link above. I suspect "insert" will be more useful to you than "resize".

Copy link
Member

Choose a reason for hiding this comment

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

If you're ever unsure of what a C++ function does, you can always look it up: https://cplusplus.com/reference/vector/vector/resize/
https://cplusplus.com/reference/vector/vector/insert/

@bobmyhill
Copy link
Member

Closed in favour of #5362

@bobmyhill bobmyhill closed this Aug 14, 2023
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.

None yet

4 participants