aspect
aspect copied to clipboard
geodynamics
Enable use of Iterative Advection schemes in combination with strain-dependent rheologies
This PR enables the use of iterative Advection schemes in the strain rheology module, after testing and discussions with @anne-glerum, @jdannberg, and others confirmed that nonlinear iterative Advection schemes should not produce unrealistic accumulation of strain within compositional fields updated with reactions terms.
- [X] I have followed the instructions for indenting my code.
- [X] I have tested my new feature locally to ensure it is correct.
- [X] I have created a testcase for the new feature/benchmark in the tests/ directory.
- [X] I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.
So everything was already working correctly? Good!
Yep, it appears so. @jdannberg write out a nice explanation of why it should work elsewhere, which may be worth adding to this PR for documentation. Short answer - during each nonlinear iteration the reaction term is being added to the old solution and not the current linearization point.
@naliboff It looks like Anne already did a review, so here just for completeness is my explanation for why it should already work as intended:
If you use reaction terms, then they’re assembled in the equation. So if we simplify everything to a first-order scheme, it would look something like
(C_new – C_old)/Delta t + u grad C = R/Delta t
where R is the reaction term. If we make this super simple and ignore advection for now, we can get the new C as
C_new = C_old + R
C_old is the old solution, which is not updated within nonlinear iterations (unless you use advection methods like prescribed fields with diffusion, for them it’s different). So even if R depends on the solution (here, I believe, it is something like strain_rate * Delta t), you get a reasonable C_new, because R is added to C_old, which does not change. The only problem would be if you had something like R = -C (using the current linearization point for R). In this case, your first iteration would set C_new to C_old – C_old = 0, but then in the next one, R would be zero because C_new is zero, so C_new = C_old, so then in the next iteration R = - C_old, and so on. You basically jump between the values. That is the problem I had in the melt models, but I do not see this as a problem here.
And I think the idea of using reaction_terms = max(update, - C) is that you can never make the composition smaller than zero with your reaction terms (you can never subtract anything more than the current value of the composition). I think that makes sense too.
@anne-glerum - thanks for the review, and I've addressed all of the comments you made above (one more note below). @jdannberg - Thanks for investigating and posting the notes above.
@anne-glerum - I tried to run the same test with the Newton Stokes solver, but always converged within one iteration after the first time step. At least during the first step the strain values stayed at 0 during the nonlinear iterations. Next, I ran the continental_extension cookbook test with the iterative Advection scheme and there the change in the max/min strain values during each nonlinear iteration was much smaller than the magnitude of the strain values. Would you like me to add this test to the PR?
During t0 the reaction terms are not computed, so no changes in the strain would be accumulated anyway. If the Picard Stokes solver does several nonlinear iterations, is there a way to make the Newton Stokes solver do several as well? By lowering the switch tolerance for example so that it does more Picard solves? What about the defect correction solver, does that do more iterations?