Hi,
I am running 2D visco-elastic-plastic subduction models. Currently, my goal is to be able to run a dynamic subduction model forward and let the elastic stresses build-up, but then after a certain point in time, reset the elastic stresses to 0 (either during the simulation or at a checkpoint restart). Since the elastic stresses are tracked as the first three compositional fields, ve_stress_xx, ve_stress_yy, ve_stress_xy, I am not sure if changing something where the elastic stresses are calculated or where the compositional fields in general are is a better approach? I had first tried making a change to /source/material_model/rheology/elasticity.cc where the stresses are calculated, but this only reduces the elastic stresses, does not make them zero, as I realized now this only affects the reaction_term:
Original snippet of elasticity.cc:
My edit trying to set the stresses to 0 at 1000 years:
If resetting compositional field values is a better approach, I am wondering about methods people have tried in the past. I did find this old thread on the forum with some ideas, https://community.geodynamics.org/t/resetting-composition-from-checkpoint-restart/1398 , but I am curious as to any issues with these in particular I may need to consider for the elasticity.
Thank you!
Rebecca:
the issue is that what you are trying to do (setting accumulated stresses to zero) is not something that is physical. ASPECT was built to simulate physical processes, and so there is not currently a mechanism to do what you want.
The options I see are this:
-
You can just “hack” ASPECT. In your case, one might put a statement at an appropriate place – say at the beginning or end of the code that executes a time step – that checks whether it is the right time to reset the stresses, and then does so. Since you are storing stresses in compositional fields, this could be as easy as just setting certain vectors to zero. This is obviously an ugly solution, and a change that cannot be packaged into a pull request and merged into the development sources of ASPECT – i.e., you have to port forward these changes every time you upgrade to a newer version of ASPECT – but it is a workable solution.
-
ASPECT has a system of “signals”. These are, in essence, places where we call a list of functions that have registered an interest in being notified when ASPECT is in specific places. There are, for example, signals that are executed just before or just after the mesh is refined, or just after the initial conditions have been set. The idea is that user plugins can register for these signals and a plugin function is then called every time the signal is executed. One way forward for you would be to introduce in ASPECT a signal that is executed at the beginning or end of a time step (whatever is most convenient for you) and that you can then use to register a function in your own plugin. This function would check whether it is the right time, and if so zero out the vectors that store the stresses. Adding the signal is something we’d be happy to do since by default it doesn’t do anything at all and so has no detrimental effect on ASPECT. The unphysical part of the action is something that you would then implement in your own plugin.
Does this help?
Best
Wolfgang
Hi Wolfgang,
We know that this is physical, but we need it in order to demonstrate that starting our models with no stresses does not impact the results. This is an issue we face in modeling the real Earth often - we want to investigate the state of deformation in a slab at the present day but we have no way of including the history of deformation leading up to the present day - so our models start out with stress equal to zero. To test whether this is an issue in interpreting the results of our models we want to create a synthetic test in which we run a model forward in time so it an acquire the elastic stresses physically. We then want to compare that output, to restarting the model from the final state and building up the elastic stresses from that point. This will only be a robust test if everything else (Temperature, geometry of other compositional fields) is unchanged. So, doing the synthetic test is what requires us to do something “unphysical” with code, like setting the stresses to zero.
I think we understand our options, but are also limited somewhat by not being expert C++ programmers. Becky and I have no idea how to use signals. Also, most of our experience with modifying Aspect is limited to the material models, so we don’t have a clear roadmap of the rest of the code. We have tried to trace back in the code where the stress components of in.composition are assigned from the FE solution vector but we have not yet been successful.
Maybe the biggest help you can give us, is to point us to the correct part of the code where in.composition is assigned from the solution vector - and then we can hack that??). We did figure out that we can’t change in.composition elsewhere (like in the Material model) because it is a constant (rightly so).
I also talk with Haoyuan about this and he has had similar test cases in the past where being able to change part of the solution from the checkpoint would have been useful, he is interested in following up on signal option during the hackathon, as it seems that there may be other cases where having this option for cases like these would be useful.
Thanks for all your help!
Magali
