Add viscosity additional outputs

[lhy11009]

Pull Request Checklist. Please read and check each box with an X. Delete any part not applicable. Ask on the forum if you need help with any step.

Add outputs of diffusion and dislocation viscosity outputs from the "visco plastic" rheology module

For all pull requests:

• I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• I have tested my new feature locally to ensure it is correct.
• I have created a testcase for the new feature/benchmark in the tests/ directory.
• I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[lhy11009]

I have talked with @alarshi before, so I assume she would be willing to take a look. Just to brief on it, I largely model it from the "PlasticAdditionalOutputs" class that exists in the "visco plastic" rheology module and is created by the "visco plastic" material model.

[tjhei]

can you document what i is?

[lhy11009]

Here I changed it to "const unsigned int point_index", following the example of the fill_plastic_outputs

[lhy11009]

Thanks @tjhei for looking into this. I have:

• addressed the previous conversations.
• updated tests
• added changelog

And below are a few comments and questions from me.

Here, I checked all the tests that should be updated for the "visco plastic" material model with "named additional outputs" as outputs. Therefore, they should be the one fails after I run the updated code.

visco_plastic_additional_plastic_outputs
visco_plastic_adiabatic_pressure_in_plasticity
visco_plastic_derivatives_2d
visco_plastic_derivatives_3d
visco_plastic_phases_plasticity
visco_plastic_total_strain_healing_temperature_dependent
visco_plastic_track_noninitial_plastic_strain
visco_plastic_viscous_strain_weakening
visco_plastic_yield_noninitial_plastic_strain_particles
visco_plastic_yield_plastic_strain_weakening
visco_plastic_yield_plastic_viscous_strain_healing
visco_plastic_yield_plastic_viscous_strain_weakening
visco_plastic_yield_plastic_viscous_strain_weakening_particles (including visco_plastic_yield_plastic_viscous_strain_weakening in prm)
visco_plastic_yield_strain_weakening_compositions
visco_plastic_yield_strain_weakening_full_strain_tensor_3d
visco_plastic_yield_strain_weakening_full_strain_tensor
visco_plastic_yield_strain_weakening
visco_plastic_yield_strain_weakening_defect_corr_stokes (including visco_plastic_yield_strain_weakening in prm)
visco_plastic_yield_strain_weakening_no_advect_no_stokes (including visco_plastic_yield_strain_weakening in prm)
visco_plastic_yield_strain_weakening_particles (including visco_plastic_yield_strain_weakening in prm)
visco_plastic_yield_total_strain_healing
visco_plastic_yield_viscous_strain_weakening
visco_plastic_yield_viscous_strain_weakening_particles (including visco_plastic_yield_viscous_strain_weakening in prm)
visco_plastic_yield_temperature_activated_viscous_strain_weakening (including visco_plastic_yield_viscous_strain_weakening in prm)

Tests that failed in docker run matched all the expected ones.

    1151 - visco_plastic_additional_plastic_outputs (Failed)
    1154 - visco_plastic_adiabatic_pressure_in_plasticity (Failed)
    1183 - visco_plastic_phases_plasticity (Failed)
    1192 - visco_plastic_total_strain_healing_temperature_dependent (Failed)
    1193 - visco_plastic_track_noninitial_plastic_strain (Failed)
    1203 - visco_plastic_viscous_strain_weakening (Failed)
    1207 - visco_plastic_yield_noninitial_plastic_strain_particles (Failed)
    1209 - visco_plastic_yield_plastic_strain_weakening (Failed)
    1210 - visco_plastic_yield_plastic_strain_weakening_particles (Failed)
    1211 - visco_plastic_yield_plastic_viscous_strain_healing (Failed)
    1212 - visco_plastic_yield_plastic_viscous_strain_weakening (Failed)
    1213 - visco_plastic_yield_plastic_viscous_strain_weakening_particles (Failed)
    1214 - visco_plastic_yield_strain_weakening (Failed)
    1215 - visco_plastic_yield_strain_weakening_compositions (Failed)
    1216 - visco_plastic_yield_strain_weakening_defect_corr_stokes (Failed)
    1217 - visco_plastic_yield_strain_weakening_full_strain_tensor (Failed)
    1218 - visco_plastic_yield_strain_weakening_full_strain_tensor_3d (Failed)
    1220 - visco_plastic_yield_strain_weakening_no_advect_no_stokes (Failed)
    1221 - visco_plastic_yield_strain_weakening_particles (Failed)
    1222 - visco_plastic_yield_temperature_activated_viscous_strain_weakening (Failed)
    1223 - visco_plastic_yield_total_strain_healing (Failed)
    1224 - visco_plastic_yield_viscous_strain_weakening (Failed)
    1225 - visco_plastic_yield_viscous_strain_weakening_particles (Failed)

Following this step, I fixed the test outputs from the newly generated files.

[lhy11009]

I also removed the new test I created in the draft, as the additional viscosity outputs are well tested in all the tests mentioned before.
@tjhei , this is ready for a second round of review.
See also my question below:

[lhy11009]

Currently, I assign the value to max limit, if the viscosity is not computed by the viscous flow law (e.g. the dislocation creep would not be calculated if the flow law is diffusion). I'm not entirely sure this is the best approach. I have tested that assigning values to nan doesn't work and will trigger error in the post-processor.

The values will be written as "inf" in outputs.

[alarshi]

Hmm, I see. Does the issue persist after #6973 was merged?
But in either case, I am not sure if having both dislocation and diffusion viscosity outputs for each viscous flow law is the right approach. See my comment above.

[lhy11009]

Oh, this and 6973 are separated.

[lhy11009]

However, in the calculation step, it seems the nan value won't trigger issue, as indicated by the examples of other variables.

[tiannh7]

Hi @lhy11009, thanks for working on this!

I had a similar need about a year ago and implemented a prototype that also added diffusion_viscosity and dislocation_viscosity as named additional outputs. In that version, I also put the selection of named outputs into a subsection of the NamedAdditionalOutputs visualization postprocessor, so users can explicitly choose which fields to write.

For reference, the commit is here:
tiannh7@176e2fe

I think this kind of user-side selection could be useful here as well, especially for controlling creep/plastic/elastic-related outputs and avoiding unnecessary output/test changes when new named outputs are added.

[lhy11009]

Hi, @tiannh7 , thanks for looking into my PR! It's great to know that this will be useful to all of us.

I agreed to your comments. After looking into your commits, I think the key is to have this block as your commits:

            prm.enter_subsection("Named additional outputs");
            {
              const std::string pattern_of_names =
                "current cohesions|"
                "current friction angles|"
                "current yield stresses|"
                "plastic yielding|"
                "diffusion viscosity|"
                "dislocation viscosity|"
                "elastic shear modulus";

              prm.declare_entry("List of named outputs", "",
                                Patterns::MultipleSelection(pattern_of_names),
                                "A comma-separated list of named additional outputs to be visualized. "
                                "If left empty, all available named outputs will be output.\n\n"
                                "The following named outputs are available:\n\n" +
                                pattern_of_names);
            }

This is indeed a better behavior than the current PR. The current PR will just try to output everything under the "named additional outputs" without letting the user select what we want.
On the other hand, implementing something like this does require a bit of extra work. The main thing is to check if the selected pattern_of_names matches what is provided by the material model (e.g. diffusion and dislocation viscosity would be present in the "visco plastic" model, but may not be present, say, the "simple" model).
I think we can indeed introduce utilities like this in future PRs, and maybe keep the contents of the current PR.
Let me know what you think.

Hi @lhy11009, thanks for working on this!

I had a similar need about a year ago and implemented a prototype that also added diffusion_viscosity and dislocation_viscosity as named additional outputs. In that version, I also put the selection of named outputs into a subsection of the NamedAdditionalOutputs visualization postprocessor, so users can explicitly choose which fields to write.

For reference, the commit is here: tiannh7@176e2fe

I think this kind of user-side selection could be useful here as well, especially for controlling creep/plastic/elastic-related outputs and avoiding unnecessary output/test changes when new named outputs are added.

[alarshi]

@lhy11009 : Thank you for adding this functionality and investigating and fixing the failing tests!

Sorry for the delay in getting back to you. I added some comments as I was reviewing your PR, but before you address them, it would be good to get someone else's feedback as well.

I misunderstood our initial discussion. I thought you wanted to do something similar to what is done in grain_size material model (see disl_viscosities_out in https://github.com/geodynamics/aspect/blob/main/source/material_model/grain_size.cc), that has both diffusion and dislocation mechanisms so the respective outputs help visualize where the two regimes are dominant.

But, looking at your PR, you want to output pure dislocation and diffusion creep viscosities for all compositions, right? Since visco_plastic material model will output both dislocation and diffusion viscosities even when they are not computed, e.g., diffusion viscosity when the chosen flow law is dislocation creep, you assigned infinite values to diffusion and dislocation viscosity when they are not computed. Hmm, I think your reasoning is correct: if we want the material to deform only following dislocation creep, we can set a very high diffusion creep viscosity. I am not sure how to think about the infinite values for the FK flow law, as the viscosity is computed differently and does not use the Arrhenius equation.

With my limited understanding of the visco_plastic material model, I only suggested some implementation changes and not what the outputs are telling us.

[alarshi]

Here, it would be helpful to add the details. As I understand, these are the diffusion viscosities that are computed in Rheology::DiffusionCreep::compute_viscosity() function before yielding. The values are only relevant when diffusion creep is present, i.e., viscous flow law is either diffusion or composite.
Or something along these lines.

[alarshi]

Same as above.

[alarshi]

Hmm, I see. Does the issue persist after #6973 was merged?
But in either case, I am not sure if having both dislocation and diffusion viscosity outputs for each viscous flow law is the right approach. See my comment above.

[alarshi]

Since you are doing almost the same operation in the each case, you could replace ln 1057-1114 by initializing the viscosities using a default value (std::max in your case) and a flag to decide whether to output diffusion viscosities or/and dislocation viscosities.
Something like
bool compute_dislocation_viscosity = (viscous_flow_law != dislocation && viscous_flow_law != frank_kamenetskii); and

if (compute_dislocation_viscosity)
     viscosity_out->dislocation_viscosities[i]
                    = MaterialUtilities::average_value(
                        volume_fractions,
                        isostrain_viscosities.dislocation_viscosities,
                        viscosity_averaging);

and similarly for the diffusion_viscosities output.

[lhy11009]

@lhy11009 : Thank you for adding this functionality and investigating and fixing the failing tests!

Sorry for the delay in getting back to you. I added some comments as I was reviewing your PR, but before you address them, it would be good to get someone else's feedback as well.

I misunderstood our initial discussion. I thought you wanted to do something similar to what is done in grain_size material model (see disl_viscosities_out in https://github.com/geodynamics/aspect/blob/main/source/material_model/grain_size.cc), that has both diffusion and dislocation mechanisms so the respective outputs help visualize where the two regimes are dominant.

But, looking at your PR, you want to output pure dislocation and diffusion creep viscosities for all compositions, right? Since visco_plastic material model will output both dislocation and diffusion viscosities even when they are not computed, e.g., diffusion viscosity when the chosen flow law is dislocation creep, you assigned infinite values to diffusion and dislocation viscosity when they are not computed. Hmm, I think your reasoning is correct: if we want the material to deform only following dislocation creep, we can set a very high diffusion creep viscosity. I am not sure how to think about the infinite values for the FK flow law, as the viscosity is computed differently and does not use the Arrhenius equation.

With my limited understanding of the visco_plastic material model, I only suggested some implementation changes and not what the outputs are telling us.

Yes, this is a very helpful summary. Ideally, we should have just the one viscosity (diffusion or dislocation) output if that is what get calculated, but both of them are created in the same class, so I am not sure how to assign value to just one and make the post-processer not output the other one (post-processer would throw error if I don't initialize both of them)

[lhy11009]

Thanks, @alarshi for the comments. I think this would be very helpful. I would look for a more generic way to just include the ones computed. I would bring up the question in the next user meeting.