Problem with update strain rheology

ASPECT
1

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
image990×383 63.3 KB

image
image985×383 64.8 KB

image
image980×368 59.8 KB

image
image971×378 64.1 KB

2

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 (590.8 KB)
log.txt (16.9 KB)

3

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.

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 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 </home/projects/shk00051/candi_dealii9.4.2/tmp/unpack/deal.II-v9.4.2/include/deal.II/dofs/dof_accessor.templates.h> in function
void dealii::DoFCellAccessor<dim, spacedim, lda>::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

4

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
image830×546 79.9 KB

image
image826×502 77.3 KB

5

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

6

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

7

Hi John,

Sure. Here is the test prm file I used.

nohealing.prm (4.3 KB)

Cheers,
Sibiao

8

Great! Then I also need to update my deal.ii version.

9

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
image755×173 49.7 KB

10

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. However, my pull request Optimize strain dependent rheology by gassmoeller · Pull Request #5335 · geodynamics/aspect · GitHub 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 a try and let us know if that fixes the issue?

Best,
Rene

11

Hi all,

FYI, the PR 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

12

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

13

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 FEFieldFunctionandFEPointEvaluation` 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

14

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 (65.8 KB)
initial_plastic_strain_1.txt (35.7 KB)

Cheers,
Sibiao

15

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