[WIP] Use correct stresses in visco-elasto-plastic rheology

[anne-glerum]

This PR is at the moment for discussing with everybody interested the changes to the elasticy and visco-plastic plugins that @bobmyhill, @naliboff and I discussed yesterday. @dansand, your input would be great to have!

FYI @EstherHeck

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.

[naliboff]

@anne-glerum - Thanks for starting this discussion! I've added comments under most of your discussion points, but will wait to add further comments until @bobmyhill and @dansand have a chance to chime in.

[naliboff]

More thoughts about where we might need to use different stress values.

[naliboff]

@anne-glerum - Partial review at the end of the today. Bravo for getting this all implemented and particularly the comments, which very succinctly and clearly outlines the "old versus new" issue in the context of ASPECTs solvers (+ how things are being updated)! 👏 Independent of elasticity, the comments are an excellent resource for anyone trying to understand the solver schemes and solution updates!

[naliboff]

I initially thought we would want to do this at t=0, but now I'm not sure it makes sense to do so if we have not actually computed a numerical time step yet? However, I think we would still definitely want to add the force term to the RHS of the momentum equation at t=0.

[anne-glerum]

At the moment, we use the fixed_elastic_timestep when no numerical timestep is available yet, so that's probably fine. If the initial fields are zero, the reaction terms will be zero. Thinking about it again, I think the reaction terms should be computed at t0 just in case the initial stresses on the fields are not zero.

The force terms are already added at t0.

[naliboff]

I think these lines should be removed from the code, but the key points added into the ASPECT manual when we put in a section for the VEP rheology.

[gassmoeller]

Can you explain again why these lines are not needed anymore? If they are unnecessary or wrong, they should be deleted (no need to keep them). If they are a feature that is unsupported with the new way we do elastic stressses, we could keep a note about this (which feature, which PR it was removed), but still remove the code to avoid outdated code and confusion (we can always get it back from the git history).

[naliboff]

I don't understand what the Assert here would be checking. Can you elaborate?

[naliboff]

I'm still confused about what this assert would be checking.

Anne: @bobmyhill do you remember?

[anne-glerum]

Thanks for the comments @naliboff , I've addressed them where possible. Maybe @bobmyhill , @jdannberg or @gassmoeller can also give this a look?

[naliboff]

@anne-glerum - I added a few more suggestions and questions to the PR, nearly entirely just regarding documentation. Overall, the PR is in great shape and really well documented (thank you!). I did not get through the full PR, but did review the key updates to elasticity.cc, etc.

A few additional general questions:

1. Can the stress_component_statistics still be removed, per a comment in the .cc file? If it does need to be included, I will review that as well.
2. Are the updates to stress.cc implementing the same changes as proposed in fixed stress visualization for elastic materials #4375?

[naliboff]

Maybe uncomment these lines for now, but add in a comment above saying that this operation will never be done due to an assert (or something to this effect)?

[naliboff]

I'm still confused about what this assert would be checking.

Anne: @bobmyhill do you remember?

[anne-glerum]

1. Can the stress_component_statistics still be removed, per a comment in the .cc file? If it does need to be included, I will review that as well.

I think the statistics postprocessor could still be useful, as it gives the statistics of the full stress at time t, instead of the the rotated and advected, but not updated, stress at time t. I've changed the comments to reflect this.

2. Are the updates to stress.cc implementing the same changes as proposed in #4375?

No. The changes in #4375 are valid for the VEP implementation in main, the changes in stress.cc are for the VEP implementation in this PR.

[gassmoeller]

I will need to finish this later, but I have a number of small items and questions at this stage.

[gassmoeller]

This is potentially a problem. By first interpolating the particle property to the field and then assigning back to the particle we introduce a lot of numerical diffusion. Why does the approach above not work? Because of the nonlinear solver?

[anne-glerum]

I think my reasoning was as follows:

• In the operator splitting step, we update the solution (the stress compositional fields) of the previous timestep.
• After that we advance the particles and evaluate the updated solution at their location
• We then update the particle properties, in this case the stress. If we use the stress values on the particles, then these do not include the operator splitting update. If we instead use the stress values that were updated and interpolated to the particle locations, then we include the operator splitting update.

[anne-glerum]

Update: now we use two particle properties, composition and elastic_stress. composition updates the particle properties with the new field solutions that include the operator splitting step. elastic_stress then applies the reaction term. So now we can stick with
particle->get_properties()[data_position + i] += material_outputs.reaction_terms[0][i];

[gassmoeller]

Can you explain again why these lines are not needed anymore? If they are unnecessary or wrong, they should be deleted (no need to keep them). If they are a feature that is unsupported with the new way we do elastic stressses, we could keep a note about this (which feature, which PR it was removed), but still remove the code to avoid outdated code and confusion (we can always get it back from the git history).

[gassmoeller]

A few more comments. Let me know if I should take a look at the postprocessors as well.

[gassmoeller]

@anne-glerum I tried to have a look at this PR, but you have a lot of unrelated changes in your branch (see https://github.com/geodynamics/aspect/pull/4370/files). It looks like you may have merged the 2.4.0 release branch, or some other changes. Can you rebase and try to remove the unrelated changes from your branch?

[anne-glerum]

Hi @gassmoeller , yes I can see that makes reviewing very complicated. I rebased to the release because that is the version we wanted to use for the paper.

I can make a copy of this branch (based on the release) and then rebase this PR to main. Would it be okay if I squash my commits before I rebase? I remember rebasing to the release being very laborious because of a multitude of commits giving merge conflicts.

Alternatively, I update the copy of this branch I still have from before the rebase to the release with the changes I have made in the meantime and then rebase that to main. I'd then still prefer to squash the commits first. I'm not sure what is preferable or whether it makes a difference.

[gassmoeller]

I can make a copy of this branch (based on the release) and then rebase this PR to main. Would it be okay if I squash my commits before I rebase?

Yes, I think that is the best approach. Keep a copy based on the release (for reference in the paper), but rebase to main so that we can merge it. Squashing the commits is fine.

[anne-glerum]

I can make a copy of this branch (based on the release) and then rebase this PR to main. Would it be okay if I squash my commits before I rebase?

Yes, I think that is the best approach. Keep a copy based on the release (for reference in the paper), but rebase to main so that we can merge it. Squashing the commits is fine.

Okay, thanks. I will do that and let you know when I'm done :)

[anne-glerum]

Hi @gassmoeller, I have rebased to main. I also ran the third benchmark (2D square under simple shear for half of the model time) again; I summarize the results below. Maybe you could keep the issues listed below in mind when going through the code.

1. Results are independent of resolution for both fields and particles. 👍
2. Results are independent of the number of nonlinear iterations allowed for fields. 👍
3. The error during shearing decreases for both fields and particles with smaller timesteps. 👍
4. The step in error when shearing is stopped decreasing for smaller timesteps for both fields and particles 👍 , but only up to a timestep size of 0.0025 s. 👎
   
5. After shearing has stopped, the particle runs use all the allowed nonlinear iterations (NI). The more iterations are allowed, the larger the error. 👎
   
6. The error is not exactly the same between fields and particles 👎 , but the difference decreases with timestep size 👍 . After shearing has stopped, the difference is largely due to the step in error for particle runs.
   

[anne-glerum]

ad. 6. Could one difference between the particle and the field runs lie in the reaction_term? It computes the spin tensor W from the solution of the previous timestep, but for particles this solution is evaluated at the new location of the particles I think? But this shouldn't matter, as the spin tensor is constant over the domain.

ad. 5. The reaction_term used to update the stresses on the particles changes linearly over the NI within one timestep instead of becoming constant. The figure shows the reaction_term value for a specific particle for a run with 100 NI allowed from the moment shearing is stopped .

EDIT:
Actually, the composition property plugin also applies its update every nonlinear iteration, increasing the stress value in every nonlinear iteration. This happens before the elastic_stresses property plugin requests the reaction term. The composition plugin is only needed to apply the update from the operator splitting to the particle properties, so it should only act once per timestep.

EDIT 2:
Only applying the composition property plugin during the first nonlinear iteration leads to a slower decrease in the stresses after shearing has stopped that leads to a large error.

[gassmoeller]

I started a review, but didnt get through the important part yet, I will continue tomorrow or Wednesday. Just some comments for now.

[gassmoeller]

Hi @anne-glerum, I think I understand the problem (after about 2 hours of reconstructing the order of operations 😄) and have an idea for a fix, but it is a bit complex:
Problem:
By using the combination of particle properties composition and elastic_stress, and interpolating the resulting property back into a compositional field, we are circumventing the particle reset that happens at the start of every nonlinear iteration. Here is why: At the start of a NI, particle properties are reset. Then the particle property stress is overwritten from the solution vector in the composition property plugin. Then the elastic_stress plugin modifies that value. Then that value is interpolated back to the field into the solution vector. Now the cycle starts again, we reset particles, but then overwrite the properties with the solution vector (which was not reset). Therefore stress is infinitely accumulated in every iteration.
Solution:
I think this will be tricky, I dont see an easy solution. We cannot limit updating the particle properties from the field to the first nonlinear iteration, because in subsequent NI's the particles will be reset to their original values from the previous timestep again, so we would be missing out on the update that happened in the operator splitting step. I think we need to discuss three avenues and decide which one is most promising:

1. break the cycle of constant update somehow, e.g. by letting the elastic_stress particle property have its own set of properties (instead of reusing composition). I dont quite see if this will work, but it feels like it should be possible somehow.
2. Make operator splitting work within nonlinear iterations. Then we could move the operator splitting step into the loop, and presumably would reset the fields as well at the start of each NI (then we could use our current workflow, and updates would not accumulate over NIs).
3. Implement operator splitting on particles. We dont really need a full operator splitting, just a possibility to update particle properties before they are advected (e.g. via a signal). Then we could do the work we currently do in the operator splitting step directly in the particle property, and that would mean we dont need to update particle properties from the fields (and therefore the reset at the beginning of the NI would work correctly). This would mean elastic stress would be completely stored and computed on particles, and the compositional field would be a read-only convenience copy for the Stokes equation only.

Each solution is a bit tricky, but i think we need one of them. I prefer 2 or 3, because they also add additional functionality to ASPECT, and are not limited to our current problem. Maybe something we need to work on during the hackathon?

Could you try the problematic benchmark with single Advection, iterated Stokes? It may not be very accurate, but the solution for particles and fields should be much closer together?

[anne-glerum]

Ohhh brilliant @gassmoeller! Thank you for taking the time to figure this out!

Concerning option 3, we would need access to the old_solution at the location of particles to compute the operator splitting update, instead of the solution the plugins get now, but that's easily done.

It seems a bit like duplicating functionality (having an operator splitting for fields and for particles), but completely storing and computing stresses on the particles is a clean option. And, like you say, it's an update to a property value, it doesn't need the additional functionality of solving the ODEs over multiple timesteps.

Maybe we can discuss it tomorrow. And yes, this does sound like a hack project.

Could you try the problematic benchmark with single Advection, iterated Stokes? It may not be very accurate, but the solution for particles and fields should be much closer together?

Using iterated Advection and Stokes with max one NI or single Advection, iterated Stokes gives the same result. It is indeed much closer to the fields, but still a bit off:

(compare dark blue solid, dashed and dotted lines above)

(compare pink and blue line with vertical dash above)

[gassmoeller]

A first batch of comments, I have to split this into several reviews because of the length of the PR.

[gassmoeller]

Please include the first line in the if-else as in:

if (this->get_parameters().enable_elasticity == false)
          heating_model_outputs.heating_source_terms[q] = stress * deviatoric_strain_rate;
else
{
  ...
}

[gassmoeller]

Is it worth to have this complicated estimate instead of just returning 1.0 (= elastic_timestep() / elastic_timestep())? It will likely only matter for the zeroeth timestep in models that set a maximum time step. But we should not accumulate any elastic stresses in time step 0 anyway.

[anne-glerum]

It's a bit complicated, but in all our benchmarks we set the maximum timestep to either equal the elastic timestep or to be a specific factor smaller than the elastic timestep. If the initial stress is nonzero, the timestep ratio matters for the force term and viscous dissipaton, and the timestep ratio also feeds into the elastic viscosity.

[gassmoeller]

Next batch of comments mostly concerning the particle property and the postprocessors.

[gassmoeller]

I do not understand this part. Is it equivalent to:

for (n=0; n<2*n_independent_components)
  composition[n] = particle_property[data_position + n];
for (n=2*n_independent_components; n<n_compositional_fields)
  composition[n] = old_solution[compositional_index[n]];

If so, maybe it is easier to write it that way (e.g. flip the order of the loops and do not overwrite entries that have just been set).

Also just because I see it here: The code assumes that all stress fields are the first fields listed, is this correct? I thought I saw somewhere that we try to move away from that restriction and only require that stress fields are contiguous (but not necessarily the first compositional fields). We can discuss if that is necessary for this PR or something for later, but I feel uncomfortable introducing more features that require a certain position of compositional fields now that we have the compositional_field_type infrastructure.

[anne-glerum]

Yes, I've started to move away from the stress fields having to be the first fields. I'll make sure that is consistently done everywhere.

I now only fill the non-stress fields with the old_solution, so there is no overwriting with the particle values anymore.

[gassmoeller]

I vaguely remember that this was discussed somewhere, but please remind me. The reason the sign is flipped between the elastic and non-elastic case is that the stress stored in the ve_stress fields already has the correct sign convention? Maybe add this to the comment in line 139.

[anne-glerum]

The implementation in main takes the stress directly from the stored fields and multiplies it with -1 to get the right sign convention. This branch's implementation builds the full shear stress from the fields stress_0 and stress_old (not negating them) and then multiplies the whole shear stress with -1. So the signs are the same?

[gassmoeller]

Next batch of comments, this time mostly for the material models.

[gassmoeller]

With #5643 you can simplify this code to:

Suggested change

	const std::vector<unsigned int> stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress);
	for (auto it = stress_field_indices.begin(); it != stress_field_indices.end(); ++it)
	composition_mask.set(*it, false);
	}
	const std::vector<unsigned int> &stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress);
	for (unsigned int field_index: stress_field_indices)
	composition_mask.set(field_index, false);
	}

But one additional question: The comment says to set a mask representing the viscoelastic stress tensor components, but to me it looks like the code is doing the opposite, it sets a mask that excludes all the stress tensor components. Do I misunderstand the code or does the comment need an update?

[anne-glerum]

The strain rheology mask function above starts with trues and for each strain field sets the true to false. The corresponding stress tensor components are also set to false here. I'll make the comment more clear.

This whole function is actually no longer used since
#5542 and could be removed, also from main.

[gassmoeller]

This condition could use an update to the compositional type machinery: this->introspection().composition_type_exists(CompositionalFieldDescription::stress)

[anne-glerum]

I think the check can be removed. If elasticity is used with the VE or VP material model, we already do this check in parse_parameters in elasticity.cc. I did add the check in.requests_property(MaterialProperties::reaction_rates)) in fill_reaction_rates.

[gassmoeller]

you can use the introspection to ask for material type here (see my comment in visco_plastic).

[anne-glerum]

same here, I've removed this check.

[gassmoeller]

I didnt make it through all files, but here are some smaller suggestions. I think this is getting quite close to merge. Let's talk a bit more about the field duplication.

[bobmyhill]

Awesome work. A few very minor suggestions, but this looks almost good to go.

[bobmyhill]

Can all of the code in this if statement be moved into elasticity.cc?

[anne-glerum]

It could be, but the stress tensors are used in a later if statement to compute the effective strain rate. That function call could get the MaterialInputs in and unwrap the tensors from the fields itself. It would do it for each volume fraction though.

[bobmyhill]

Might we want to have a separate shear_modulus_averaging parameter declared in the elastic rheology module?

[gassmoeller]

@anne-glerum let us know when you want us to take another look at this.

[anne-glerum]

@gassmoeller I will, I want to test it again first on our cluster.