# Problem with update strain rheology **URL:** \ **Category:** ASPECT\ **Created:** [August 1, 2023, 7:41pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057 "2023-08-01T19:41:25Z")\ **Posts on this page:** 15\ **Page:** 1
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 1, 2023, 7:41pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/1 "2023-08-01T19:41:25Z")
Hello, I am currently testing the ‘update strain rheology logic’ provided by @jbnaliboff. I’d like to verify the accuracy of the old\_plastic\_strain calculation when computing the reaction term in the “plastic\_strain” compositional field. To give you some context, I have set the initial plastic strain to 1, and the delta\_e\_ii\_plastic (= edot\_ii \* t, with no strain healing) remains constant at 0.01. Ideally, the plastic strain should increase by 0.01 at each time step, following this pattern: 1 → 1.01 (first timestep) → 1.02 (second timestep) → 1.03 (third timestep)… Therefore, the old\_plastic strain values should be 1 (timestep 1) → 1.01 (timestep 2) → 1.02 (timestep 3)… However, the results of my test show that the old\_plastic\_strain values do not match this pattern. In the attached snapshots, you can see that the old\_plastic\_strain values are 1 (timestep 1, _correct_), 1.00778 (timestep 2, _wrong_), 1.01641 (timestep 3, _wrong_)… My question is how the old\_plastic\_strain is calculated? Could the inconsistency in the test results be attributed to the iterative solution of the plastic\_strain system? Thanks and looking forward to your suggestions! Best, Sibiao ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/f/f36282b5b0efe8719962fe78fbe4d674a57c515d.jpeg) ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/5/5bca6976d3b55074016eda5ac4b006c64805e637.jpeg) ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/5/535fbb7d4aa186836dfaac5a56f52797f79fa858.jpeg) ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/b/bc7f4341f40ff6cd1509a6da6db6bb2293065a61.jpeg) ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 2, 2023, 3:37pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/2 "2023-08-02T15:37:00Z")
Hi John @jbnaliboff and all, I also tested the continental\_extension case in the cookbooks and it failed. Here I attach the error message and the log file. [error.txt](https://community.geodynamics.org/uploads/short-url/lZ1OHAmW4zuPtsX8RjZMG0tYrt9.txt) (590.8 KB) [log.txt](https://community.geodynamics.org/uploads/short-url/cuM5Rw9FDjTbKT088tdIAhLpqtt.txt) (16.9 KB) ---
**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\ **Post date:** [August 5, 2023, 6:36pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/3 "2023-08-05T18:36:08Z")
Hi @sibiao - Thanks for posting these questions and tests to the forum. First, I think the good news is that your test shows that the strain magnitude updates with each time step (0.01) are correct. Do you agree? > My question is how the old\_plastic\_strain is calculated? Could the inconsistency in the test results be attributed to the iterative solution of the plastic\_strain system? The explanation below largely follows from the discussion from [pull request 4336](https://github.com/geodynamics/aspect/pull/4336). To start, `in.composition` is an extrapolation of the old solution at the current time step. Thus, when we were originally subtracting the values from in.composition to update the reaction terms, there was always a small error that only became noticeable in models with very small time steps. The` old_plastic_strain` (or any other of the other “old” strain values) are the old solution (previous time step) extracted at the current cell. I imagine this difference you see between the expected and observed old values are created during this extraction (interpolation?), but this is then canceled out somewhere else during the updates. I’m not sure exactly where it is canceled out, but I imagine it may be during the reaction term updates where `in.composition` is used. For example, can you print out the value of `in.composition[i][this->introspection().compositional_index_for_name("plastic_strain")` in strain\_dependent.cc to see if it also has similar deviations from the expected values? > I also tested the continental\_extension case in the cookbooks and it failed. Here I attach the error message and the log file. That’s weird. We have a [test case](https://github.com/geodynamics/aspect/blob/main/tests/continental_extension.prm) for that cookbook, that also runs for two time steps. Your model crashed on the second timestep, with an error I haven’t seen before: > * * * > > An error occurred in line \<2721\> of file \ in function > void dealii::DoFCellAccessor\::get\_dof\_values(const InputVector&, ForwardIterator, ForwardIterator) const [with InputVector = dealii::TrilinosWrappers::MPI::Vector; ForwardIterator = double\*; int dimension\_ = 2; int space\_dimension\_ = 2; bool level\_dof\_access = false] > The violated condition was: > this-\>is\_artificial() == false > Additional information: > Can’t ask for DoF indices on artificial cells. What version of deal.II and ASPECT did you use to run this model? I’m unfortunately not on a computer where I can easily test the most recent ASPECT master, but can do so in a day or two. Cheers, John ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 7, 2023, 7:58am UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/4 "2023-08-07T07:58:12Z")
Hi @jbnaliboff, Thanks for the reply! The extraction or interpolation of old\_plastic\_strain is what I also suspect is causing such difference. I have printed out the value of in.composition[plastic\_strain]. See the attached snapshots. The extracted value of old\_plastic\_strain is the value of current plastic strain (= in.composition[I][…“plastic\_strain”] from the previous timestep. You are right. The reaction term should be somewhere canceled out during the updates. Maybe by the advection term? Regarding the continental\_extension test, I used deal.ii 9.4.2 and ASPECT version 2.6.0-pre (master, 69faf2563). Update: If I use free-slip BC at the top instead of free surface, it works. Cheers, Sibiao ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/1/13a7bb27dc1720a152f5e2b982f286b3aa4d44b7.jpeg) ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/0/0f06aaadc902191755a49def84d32f78e1e3b2c3.jpeg) ---
**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\ **Post date:** [August 7, 2023, 7:05pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/5 "2023-08-07T19:05:58Z")
Hi @sibiao, 1. Just to make sure the correct values are being computed for the plastic strain accumulation, can you confirm that the output values of the plastic strain (i.e., what is output by postprocessors) do indeed increase by a factor of 0.1? If not, can you post the prm file you used? 2. FYI, I was able to reproduce the error with the continental extension cookbook on with the same aspect commit and deal.II 9.4.0. I’m going to test the prm file with deal.II 9.5.0, and then from there try to sort out what commit(s) produced the error. Cheers, John ---
**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\ **Post date:** [August 7, 2023, 7:14pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/6 "2023-08-07T19:14:45Z")
Hi Sibiao, Update - the continental extension cookbook runs successfully with deal.II 9.5.0. For now, I suggest updating to 9.5.0 and we will try to sort out what the issue is with earlier versions of deal.II soon. Cheers, John ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 7, 2023, 7:35pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/7 "2023-08-07T19:35:15Z")
Hi John, Sure. Here is the test prm file I used. [nohealing.prm](https://community.geodynamics.org/uploads/short-url/zYgnsfVuE6LoOJikXZbGnBZHqV7.prm) (4.3 KB) Cheers, Sibiao ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 7, 2023, 7:37pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/8 "2023-08-07T19:37:01Z")
> [@jbnaliboff](#): > > Hi Sibiao, > > Update - the continental extension cookbook runs successfully with deal.II 9.5.0. For now, I suggest updating to 9.5.0 and we will try to sort out what the issue is with earlier versions of deal.II soon. > > Cheers, > John Great! Then I also need to update my deal.ii version. ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 9, 2023, 7:51pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/9 "2023-08-09T19:51:47Z")
Hi John @jbnaliboff, Update: I have installed the latest deai.II version: 9.5.1, and run the continental\_extension case. Unfortunately, I encountered the same error as before. Since I run it on HLRN, I am not sure if such issue is related to the different systems. Hi Anne @anne-glerum, could I kindly request you to spare a moment at your convenience to run the same case from your HLRN account? Thank you so much! For your reference, the last commit of ASPECT is 69faf2563. Cheers, Sibiao ![image](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/b/ba082d1a01c06767397eb1a8647d06e28bb687e3.png) ---
**Author:** ![gassmoeller](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/gassmoeller/32/11_2.png) [@gassmoeller](https://community.geodynamics.org/u/gassmoeller)\ **Post date:** [August 11, 2023, 5:19pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/10 "2023-08-11T17:19:46Z")
Hi @sibiao, I dont have much to contribute to this discussion at this point, except I noted the following point: Your exception is thrown inside the `FEFieldFunction` class. As you noted this class was introduced in [Update strain rheology by naliboff · Pull Request #4336 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/pull/4336). However, my pull request [Optimize strain dependent rheology by gassmoeller · Pull Request #5335 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/pull/5335) will remove this class again in favor of a faster implementation in FEPointEvaluation. My PR should not change any functionality, but it should affect the error message you get (or maybe even fix it, if the problem is how we use the FEFieldFunction class). So maybe give the branch in [Optimize strain dependent rheology by gassmoeller · Pull Request #5335 · geodynamics/aspect · GitHub](https://github.com/geodynamics/aspect/pull/5335) a try and let us know if that fixes the issue? Best, Rene ---
**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\ **Post date:** [August 11, 2023, 6:19pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/11 "2023-08-11T18:19:57Z")
Hi all, FYI, the [PR](https://github.com/geodynamics/aspect/pull/5335) from @gassmoeller has been merged into the main branch, and indeed the continental extension cookbook now runs successfully when ASPECT was compiled with deal.II 9.4.0. @sibiao - can you try the ASPECT main branch on your system to confirm? @gassmoeller - this error deriving `FEFieldFunction` class must not have appeared in the melt materials as those are not typically used in combination with a free surface. I’ll make a note in the github issue that we should probably update the code there as well. Cheers, John ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 11, 2023, 8:25pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/12 "2023-08-11T20:25:29Z")
Hi @jbnaliboff and @gassmoeller, Yes, I tried the new version with merging of the PR from Rene. The continental\_extension cookbook also works on my system, for both deal.ii 9.4.2 and 9.5.1! Furthermore, I revisited the plastic strain test as mentioned earlier. Notably, when I initialized the plastic strain value to 0.01, which is the same as the value of plastic strain increment or reaction term. I noticed that the expected value and observed old solution are generally consistent now. I think this update should be working without problems. I would like to thank both of you for your support and time spent on this issue! Cheers, Sibiao ---
**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\ **Post date:** [August 15, 2023, 3:00am UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/13 "2023-08-15T03:00:48Z")
Hi @sibiao, > Yes, I tried the new version with merging of the PR from Rene. The continental\_extension cookbook also works on my system, for both deal.ii 9.4.2 and 9.5.1! Thanks for confirming! I will go back and check to see whether the continental extension cookbook test would crash using the `FEFieldFunction` class if we ran it for a few more time steps. This may be useful for future error-proofing. > Furthermore, I revisited the plastic strain test as mentioned earlier. Notably, when I initialized the plastic strain value to 0.01, which is the same as the value of plastic strain increment or reaction term. I noticed that the expected value and observed old solution are generally consistent now. Interesting, can you post the output from these tests? I’m not entirely sure what is causing this differences between the FEFieldFunction`and`FEPointEvaluation` approaches, but it would be good to pinpoint the exact change. I also think it may be worth updating (simplifying) the strain rheology tests to ensure they can be easily compared with analytical solution. > I think this update should be working without problems. I would like to thank both of you for your support and time spent on this issue! Likewise, thanks for highlighting the issue and helping with the testing. Cheers, John ---
**Author:** ![sibiao](https://avatars.discourse-cdn.com/v4/letter/s/49beb7/32.png) [@sibiao](https://community.geodynamics.org/u/sibiao)\ **Post date:** [August 15, 2023, 5:17pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/14 "2023-08-15T17:17:09Z")
Hi John, Sure. Here I attach two outputs of setting initial plastic strain to 1 and 0.01. Correction: There is no difference between the FEFieldFunction and FEPointEvaluation approaches, the results of models using them are the same. But the model with FEPointEvaluation runs significantly faster than another one. [initial\_plastic\_strain\_0.01.txt](https://community.geodynamics.org/uploads/short-url/ekKlYBI0qywNCBMm66K01tnKoHR.txt) (65.8 KB) [initial\_plastic\_strain\_1.txt](https://community.geodynamics.org/uploads/short-url/rtpZMBqvG0ISTx1cjVI7PCxGraG.txt) (35.7 KB) Cheers, Sibiao ---
**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\ **Post date:** [August 16, 2023, 9:31pm UTC](https://community.geodynamics.org/t/problem-with-update-strain-rheology/3057/15 "2023-08-16T21:31:47Z")
Hi Sibiao, Great to hear both methods are giving the same results for your test, but the FePointEvaluation is much faster (thanks @gassmoeller!). Cheers, John