Add compositional dependence on viscosity in the Steinberger material model

[bobmyhill]

This PR duplicates #5260, modifying the logic used to read the viscosity prefactors from the prm, and tidying the code.

This PR slightly changes the test output for steinberger_phase_fractions_background (velocities and Vp anomalies change by 1e-8). I have checked that "compositional_prefactor" for this simulation always evaluates to 1 (or sometimes to 1.0000000000000002).

[poulamiiroy]

@bobmyhill, Perfect! Thank you for your work on this! We can use this PR.

[jdannberg]

Thank you for finishing this and cleaning up the code!
I have only one question (about the background field), but otherwise this looks good to me!

[jdannberg]

Would it make more sense to call this something like

Suggested change

	prm.declare_entry ("Viscosity prefactors", "1",
	prm.declare_entry ("Compositional viscosity prefactors", "1",

to be more precise about what the parameter is for? Who knows what others will want to add in the future (just thinking about phase prefactors, etc.)
Or is this name to be consistent with another material model?

[bobmyhill]

This is the name given also in the "latent heat" and "ascii reference profile" material models. I think we should probably keep it the same.

[jdannberg]

I just checked latent heat, and it has a Composition viscosity prefactor.

aspect/source/material_model/latent_heat.cc

Lines 256 to 258 in 8946fb9

	prm.declare_entry ("Composition viscosity prefactor", "1.0",
	Patterns::Double (0.),
	"A linear dependency of viscosity on composition. Dimensionless prefactor.");

The "composition reaction" model are has Composition viscosity prefactor 1 and Composition viscosity prefactor 2. So I think it would make sense to make it consistent with that.

The Viscosity prefactors in the "ascii reference profile" are actually for phase transitions, not for compositions. So keeping that name would be very confusing.

[jdannberg]

I don't understand this logic. What if a model has one compositional field and one look-up? Then it wouldn't have a background field, right?

If I understand the code below correctly, in this case you would erase the name of this field from the list of allowed keys, so if one tried to use the name of this field to specify the viscosity, it would not work. (And I admit that it's unlikely anyone would do that, since it doesn't make a lot of sense to specify a Compositional prefactor if there's just one field, I just thought it was confusing.)

Also what happens if there are 5 chemical fields and one look-up? Do we want to allow the user to specify 5 different compositional prefactors in that case? You could of course always get the same behaviour by loading the same look-up table 5 times, but that does not seem very efficient. But I am also fine with not allowing that right now (this is something that can be added in the future without any incompatible changes).

[bobmyhill]

Yes, this is also confusing to me, but is required to make two of the existing tests work.

The problematic tests are steinberger_phases.prm and steinberger_several_fields.prm. Here's the blurb:
This is a copy of the steinberger compressible test. It is intended to test the ability of the steinberger material model to use the first data file as a background composition, and ignore other compositional fields if only a single material file is provided.

These are tests that you wrote in 2021. As an alternative to the complicated logic in the parse_parameters function, we could instead define the types of compositional fields in those prms.

[bobmyhill]

@jdannberg I've decided to simplify the logic in the source code and require that the type is specified for fields that aren't of type "chemical composition".

[jdannberg]

I am not really happy with this solution. It is not just because of not being backwards compatible, but it would also be inconsistent with all of our other parameters (which @gassmoeller pointed out when we discussed this). For all other parameters, you can just specify one value, and this value will be used for all fields (of whatever field type). So you should be able to just specify one table, which will then be used for all fields (independent of field type). For tables this is actually even more important than for other parameters, because it may require a significant amount of memory to read in the same table several times.

Is there no way to do this right with the original logic of having a background field?

[bobmyhill]

Ok, first let me describe what I think steinberger_several_fields.prm should do:

• There are three compositional fields and one lookup.
• The compositional fields don't correspond to chemical compositions (values from 0 to 3e7). They represent some other tracked variable.
• The field type for these fields is "chemical composition" (the default value), but on main branch this is ignored in steinberger.cc because there is just one lookup.

Now I'll describe what I originally implemented (https://github.com/geodynamics/aspect/pull/5362/files/06ce80544b7e64ea1025552d25178cf202f77b03) and how that maintained the status-quo:

• On main branch, there is a background field if there is one look-up and no compositional fields, but not if there is one look-up and some compositional fields.
• In this PR, I changed that so that there is always a background field if there is one look-up. That is in better agreement with the descriptions in the two tests: It is intended to test the ability of the steinberger material model *to use the first data file as a background composition*, and ignore other compositional fields if only a single material file is provided.
• This doesn't change any existing functionality, but it does mean that users cannot use those fields to modify viscosities.

To answer your questions:

If I understand the code below correctly, in this case you would erase the name of this field from the list of allowed keys, so if one tried to use the name of this field to specify the viscosity, it would not work.

Yes, that is (I think) the correct behaviour according to the test descriptions, which suggest that if there is one lookup, compositional fields don't correspond to different material proportions.

Also what happens if there are 5 chemical fields and one look-up? Do we want to allow the user to specify 5 different compositional prefactors in that case?

No, that would be inconsistent with the current test descriptions, and it would also lead to undefined behaviour. How would we know which compositional fields correspond to chemical compositions (given that the behaviour on main branch explicitly ignores all fields, including those of type chemical composition)?

So you should be able to just specify one table, which will then be used for all fields (independent of field type). For tables this is actually even more important than for other parameters, because it may require a significant amount of memory to read in the same table several times.

On main branch, if there's only one lookup table, all compositional fields are explicitly ignored by the steinberger material model.

[bobmyhill]

I agree that the solution I implemented after your initial comments is further away from what you envisage. Sorry, I misunderstood your message. I'm happy to remove that commit - we could probably merge the PR then and discuss further changes separately.

One way to implement what you have in mind (one lookup, multiple viscosity prefactors) would be to add a further parameter cross-referencing the lookups to the fields of type "chemical composition" (i.e. a vector of unsigned ints or strings corresponding to the field names). I think the number of viscosity prefactors (after parsing) should always be equal to the number of fields of type chemical composition, so steinberger_phases.prm and steinberger_several_fields.prm would still need to be modified to specify the field types, and that would still break backwards compatibility.

I can't think of another way to do what you want without serious opportunity for confusion, so I'd appreciate any ideas you might have!

[jdannberg]

I think allowing different viscosity prefactors if there is only one look-up is a separate issue we don’t have to discuss now, since adding functionality later in another PR is always fine, so I totally agree with you on that. (And I can think about ideas for how to do that.)

If everything works as before when you just remove the “remove one lookup case from steinberger“ commit, I think that would be a good solution.

That still leaves the change with the name of the background field, i.e., if it is called background or whatever name the user gave it. I think it is okay to change the name of the key from how it is now on main (name of the compositional field) to background, and I agree that is more consistent. The case where it would break input files is so special that I just hope it won't be an issue.

The part that was important to me is keeping the option to have several fields (compositional fields, chemical fields, whatever) and to just read in one table for all of them. And from your second message, it seems like that would work when removing the “remove one lookup case from steinberger“ commit?

Just to explain, this is intended for the material model, since the documentation of the “Material file names” parameter says

The file names of the material data (material data is assumed to be in order with the ordering of the compositional fields). Note that there are three options on how many files need to be listed here: 1. If only one file is provided, it is used for the whole model domain, and compositional fields are ignored. 2. If there is one more file name than the number of compositional fields, then the first file is assumed to define a `background composition' that is modified by the compositional fields. If there are exactly as many files as compositional fields, the fields are assumed to represent the fractions of different materials and the average property is computed as a sum of the value of the compositional field times the material property of that field.

[bobmyhill]

Thank you for the review :)
I'll try to get round to making the changes tomorrow!

[bobmyhill]

@cmills1095, @gassmoeller: any ideas why readthedocs is failing? This PR passed until today, and I don't think I changed anything significant with respect to the docs.

[jdannberg]

I just checked latent heat, and it has a Composition viscosity prefactor.

aspect/source/material_model/latent_heat.cc

Lines 256 to 258 in 8946fb9

	prm.declare_entry ("Composition viscosity prefactor", "1.0",
	Patterns::Double (0.),
	"A linear dependency of viscosity on composition. Dimensionless prefactor.");

The "composition reaction" model are has Composition viscosity prefactor 1 and Composition viscosity prefactor 2. So I think it would make sense to make it consistent with that.

The Viscosity prefactors in the "ascii reference profile" are actually for phase transitions, not for compositions. So keeping that name would be very confusing.

[jdannberg]

Thank you for working to fix this!

I think there should still be a way to allow for the original logic of having a background field. It will just make reading in the prefactors a little bit more complicated. I think you just have to have a special case when there is one chemical field and one look-up table.

[bobmyhill]

@jdannberg thanks! I removed the offending commit, rebased and squashed. The testers pass, so this should be good to go if you're happy.

[jdannberg]

Yes, thank you for all your work on this!

[gassmoeller]

@gassmoeller: any ideas why readthedocs is failing? This PR passed until today, and I don't think I changed anything significant with respect to the docs.

Sorry, just saw this question. Sometimes the build of readthedocs fails randomly, I think it is caused by whatever cloud they are using to build the documentation. Usually a rerun fixes it (triggered manually by me, or by just pushing a new commit or amending an existing commit).