Make deviator() and second_invariant() consistent with plane strain assumption in 2D

[YiminJin]

This PR modifies the functions deviator() and second_invariant() to make them consistent with the plane strain assumption in 2D. The modified functions are placed in namespace Utilities::Tensors. For details of the mathematics, please refer to #6434 and #6459. A direct impact of this modification is that it sharpen the shear bands in plastic models with nonzero dilation angle.

This PR is a follow-up of #6373, but it will affect other material models. Also, we need several test problems to show the plane strain functions work well (or to find out if we can do better), which have not been implemented. Please make comments if you have suggestions or questions about this PR.

[bangerth]

Others know this area much better than me, but in any case, this deserves a changelog entry.

[gassmoeller]

Thank you for tackling these questions. I left a few comments. I think this is the right approach, but there are still some places that compute the deviatoric strain rate as strain rate - 1/dim * trace, please fix those as well.

In light of the discussion in #6434 I think this is the right path forward, but like Wolfgang commented, please add a changelog entry. And we will likely have to update a lot of test results.

[gassmoeller]

here we still compute the shear strain rate as strain_rate -1./dim *trace. Does this need to be changed as well? If so, there are several places in the code base where we compute the shear strain rate in this way. Please search for /dim in ASPECT's directories and check which places need fixes.

[YiminJin]

Thank you for pointing this out. I found two more 1/dim * trace expressions in material models, and they are all replaced by Utilities::Tensors::deviator(). Also, a changelog is added for these modifications.

[gassmoeller]

why do we have a different way to compute the second invariant here? This way of computing the invariant also appears in the spiegelman benchmark cases. We should figure out in which way this is identical or different to the previous form.

[bangerth]

It's worth mentioning that @cedrict and I looked into this at some point in the past and found that people have all sorts of incompatible definitions of the second invariant. Sometimes it included a factor of 2, sometime it didn't, and similar shenanigans. I would suggest sticking closely to the definition we had previously used, and only change 1/dim to 1/3.

[YiminJin]

Sorry, this is my miss.

[gassmoeller]

you probably didnt mean to rename this function, right?

[YiminJin]

No, it is my mistake. I just used a script to replace the string second_invariant by Utilities::Tensors::deviatoric_tensor_inv2 in all the files.

Do you think the second invariant of strain rate tensor should be changed here? I do not know if Spiegelman use a special function here (the function calculates the norm of the strain rate tensor).

[gassmoeller]

Just a few more thoughts I had after the review.

[gassmoeller]

Add one sentence that in 3D this function simply returns the usual second invariant.

[gassmoeller]

Add one sentence that in 3D this returns the usual definition of a deviatoric tensor.

[YiminJin]

The comments for the two functions are revised accordingly.

[gassmoeller]

I am unsure about this name, maybe second_invariant_plane_strain would explain better what the purpose of this function is?

[YiminJin]

I do not know. I agree that plane_strain should be added to the name, but deviatoric should also be added, because the function returns the correct result only when the input tensor is deviatoric. Yet, we are not able to check if the input tensor is deviatoric, for the deviator of 2D plane strain tensor is not "deviatoric" in the common sense...

[YiminJin]

Completion of this PR is harder than expected. Currently there are two major problems:

1. The names of the plane-strain-consistent functions. After discussing with @gassmoeller and @bangerth , I decide to use consistent_deviator and consistent_second_invariant_of_deviatoric_tensor temporarily. Apparently, the latter is not a good name. I want to include both "consistent" and "deviatoric tensor" in the name to remind the users that this function only works for deviatoric tensors consistent with plane strain assumption in 2D. Could anyone help me find a better name for this function?
2. The viscosity derivative w.r.t. strain rate in material models should be modified. Currently, we calculate the derivative w.r.t. deviatoric strain rate $\varepsilon'$. In order to do so, we need to:
   (a) calculate the deviatoric strain tensor, $\varepsilon'=dev(\varepsilon)$;
   (b) add a small increment, $\tilde\varepsilon' = \varepsilon' + d\varepsilon'$;
   (c) feed it to the function that calculate the viscosity, $\eta = \eta(\tilde\varepsilon')$.
   However, in most cases, the viscosity is dependent on the deviatoric strain tensor, i.e. the function $\eta(\cdot)$ will call the function $dev(\cdot)$ again. This will cause problem, because the deviator of a deviatoric tensor in the plane-strain sense is no longer itself again:
   $dev(dev(\varepsilon)) \neq dev(\varepsilon)$.
   Now I have only modified the visco plastic model, and the Newton solver works. There are other material models, such as the drucker prager model, the grain size model, and some models in tests and benchmarks need to be modified.
   @MFraters Could you please take a look at this comment and check if it is correct? If it is correct, could you help me modify the Spiegelman benchmark? I think Spiegelman intended to use plane strain assumption in his original paper (Eq. (4)).

[gassmoeller]

Just a few small comments I saw while reading through the PR.

[gassmoeller]

was this the place where you removed the deviator on purpose?

[YiminJin]

Yes. I think we need to be more careful when using consistent_deviator(), since we cannot apply it more than once to a symmetric tensor. I move the deviator here to the place where it is called --- MaterialModel::Rheology::ViscoPlastic::compute_isostrain_viscosities (`source/material_model/rheology/visco_plastic.cc, line 310). But I need to think this over, because under this change the Newton assembler gets the full strain rate instead of deviatoric strain rate.

[bangerth]

By the way, is the name of the variable now wrong?

[MFraters]

So instead of strain_rate it should be something like consistent_deviator_of_strain_rate or just deviatoric_strain_rate?

[gassmoeller]

I think this comment is not addressed yet. The name of the variable is edot_deviator, but it now uses the full strain rate. Does that name of the variable have to be changed?

[gassmoeller]

add another empty line here

[gassmoeller]

and another empty line here

[YiminJin]

Sure. I also add another empty line between to_voigt_stiffness_vector and levi_civita to make the format consistent.

[gassmoeller]

why is this no longer the deviatoric force?

[YiminJin]

Actually I do not know. But the Stokes also uses the full elastic force, and the deviatoric elastic force does not produce sharp shear bands in the strip footing test. I need some time to figure out why.

[MFraters]

What did you find?

[gassmoeller]

A plan forward for this PR that I just discussed with @YiminJin:

1. We agree that this change fundamentally is the right path forward.
2. @YiminJin will upload the changed test results to this PR so we have an overview of the changes caused by this PR at the moment.
3. Material models should not hand over deviatoric strain rates to subfunctions, because applying the deviator again changes the results as @YiminJin described above. If a strain rate is expected as input of a function, it needs to be the full strain rate. This will affect the viscosity derivative computation in some material models like visco_plastic.
4. We need to test what the change in 3. will mean for the convergence rate and stability of the Newton solver. Benchmarks to re-run would be the Spiegelman benchmark and the nonlinear channel flow benchmark. We should also check if we have a compressible benchmark with the Newton solver and if that still gives the correct results.

We want to make sure this PR is tested and merged before the next release, since it fixes an important bug in the current implementation.

[YiminJin]

@gassmoeller The reasons for removing the deviator when calculating the viscoelastic strain rate and assembling the system rhs are:

1. When the material model is incompressible, the top-left block of the Stokes system should be
   $\int_\Omega2\eta\nabla_S\varphi : \varepsilon d\Omega$, (1)
   not
   $\int_\Omega\nabla_S\boldsymbol\varphi : dev(\varepsilon) d\Omega$, (2)
   because if we use (2), then the matrix becomes
   $\int_\Omega2\eta\nabla_S\varphi : \left[\nabla_S\varphi - \frac{1}{3}(\nabla\cdot\varphi)I\right] d\Omega$,
   which is not symmetric. Therefore, when we use Newton method, the system Jacobian should be
   $\int_\Omega(2\eta\nabla_S\varphi : \varepsilon + 2\nabla_S\varphi : \varepsilon\otimes\frac{\partial\eta}{\partial\varepsilon} : \nabla_S\varphi)d\Omega$,
   where $\varepsilon$ is the full strain rate.

I know it contradicts my modifications on the Newton assembler last year, but I cannot remember why I use deviatoric strain rate at that time...

2. In source/material_model/rheology/elasticity.cc, deviator is never applied to the stress tensor, for it is deviatoric itself. Although there would be numerical errors if we use field method to advect the deviatoric stress fields, I think it is better to keep the consistency between the codes material model and assembler.

The comments above are only simple analyses. I think our final choice should depend on experiments. I am working on it.

[YiminJin]

Sorry, I was wrong again:
$\int_\Omega\nabla_S\varphi : (\nabla_S\varphi - \frac{1}{3}\nabla\cdot\varphi I)d\Omega = \int_\Omega[\nabla_S\varphi : \nabla_S\varphi - \frac{1}{3}(\nabla\cdot\varphi)(\nabla\cdot\varphi)]d\Omega.$
The matrix is symmetric.

[bangerth]

You already found out that
$\int_\Omega\nabla_S\varphi : (\nabla_S\varphi - \frac{1}{3}\nabla\cdot\varphi I)d\Omega = \int_\Omega[\nabla_S\varphi : \nabla_S\varphi - \frac{1}{3}(\nabla\cdot\varphi)(\nabla\cdot\varphi)]d\Omega$

Moreover, we have
$\int_\Omega(\nabla\cdot\varphi I) : (\nabla_S\varphi - \frac{1}{3}\nabla\cdot\varphi I)d\Omega = \int_\Omega[ (\nabla\cdot\varphi I) (I : \nabla_S\varphi - \frac{1}{3}\nabla\cdot\varphi I : I) ]d\Omega = \int_\Omega[ (\nabla\cdot\varphi I) (\nabla\cdot\varphi - \nabla\cdot\varphi) ]d\Omega = 0$
and as a consequence,
$\int_\Omega\nabla_S\boldsymbol\varphi : dev(\varepsilon) d\Omega = \int_\Omega dev(\varepsilon) : dev(\varepsilon) d\Omega$.

Both lines of the argument indicate that the system matrix is symmetric.

[YiminJin]

I am trying out the benchmark in Spiegelman et al. (2016) with the modifications in this PR and #6546 . As a first step, I use the visco plastic model and the same material parameters as in benchmarks/newton_newton_solver_benchmark_set/spiegelman_et_al_2016/input.prm (except for the "regularization" viscosity). The prm file is:
spiegelman_et_al_2016.prm.txt
The results:

where $\eta^{vp}$ is the viscosity of the viscous damper. In all these tests, I cannot make the Newton solver converge to $10^{-4}$ within 100 iterations, and the Newton solver does not converge when Use Newton residual scaling method is off. This is apparently an unacceptable result. However, when I run the original benchmark with the main branch, the result looks like

which is quite different from the figures in Spiegelman et al. (2016). Do we have model settings that can reproduce the results in the paper with a good convergence rate? @gassmoeller @MFraters

[MFraters]

hmm, that means they have been broken somewhere along the way. Just checking, have you tried running the original files and material model here: https://github.com/geodynamics/aspect/tree/main/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016?

[YiminJin]

Just checking, have you tried running the original files and material model here: https://github.com/geodynamics/aspect/tree/main/benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016

Yes. I only changed the refinement level from 4 to 6.

[YiminJin]

I derived the analytical expression of the system Jacobian under plane-strain assumption for DP rheology, and I think I have found the reason for the different convergence behaviors with and without stabilization shown in #6160 .

The derivation is provided in the following document:
plane_strain_DP.pdf

In one word, the reason for the better performance of using deviatoric strain rate when SPD stabilization is turned on is that: when using the full strain rate, the Jacobian matrix is not semi-positive-definite, and we need a smaller scaling factor to restore the positive-definiteness. However, I think we should use the full strain rate even if it slows down the convergence rate, because it is correct (if we use the deviatoric strain rate, it would be equivalent to modifying the rheological model).

I made some modifications in source/simulator/assemblers/newton_stokes.cc and benchmarks/newton_solver_benchmark_set/spiegelman_et_al_2016/drucker_prager_compositions.cc accordingly. When using 8 levels of refinements (1024 * 256 cells), the result now looks like:

The shear bands are still not as sharp as those in the original paper of Fraters et al. (2019), but are similar to those in Spiegelman et al. (2016). Surprisingly, the nonlinear solver converges to $10^{-14}$ with only 13 iterations.

I also plotted the convergence curves with low-resolution models:

The convergence curves of Newton solver are very close to those plotted in #6160. Furthermore, these curves are produced by models with boundary velocity of 5 mm/y (I forgot to restore the boundary velocity after the high-resolution experiment), so the convergence behavior is actually better than before (in the original paper, the nonlinear solver fails to converge to $10^{-14}$ with Vel = 5 mm/y).

What is your opinion about my modifications @gassmoeller @MFraters ? If they are acceptable, I think it is time to apply the differences in test results and prepare for the final merge.

[YiminJin]

Corrections:

1. the parameters I used to reproduce the convergence curves were not the same as in the original paper. Here are the results using the same parameters ($\eta_{ref}$ = 1e23, Vel = 2.5 mm/y, mLT = 1e-8):

The curves are almost the same as those plotted in #6160 (this time we use full strain rate when assembling the system Jacobian).

2. The reference viscosity I used to reproduce the high-resolution result was one magnitude smaller than that in the original paper. Here is the result using the same parameters ($\eta_{ref}$ = 1e24, Vel = 5 mm/y, mLT=1e-8):

It converges to about 3e-6. The strain rate field now is very similar to Fig. 3 in the original paper.

[MFraters]

Thanks for doing all this work! That is great.

Do I understand it correctly that the convergence plot you show should be the same as the bottom left one in figure 5?

The stabilization definitely looks better, even though, looking at the Picard iteration the problems is more difficult, since that one is not converging.

Do you have the same plots for the more difficult cases? It would be interesting to see how it behaves there.

[MFraters]

why not simply t[0][0]*t[0][0]?

Same in other places.

[YiminJin]

I thought fixed_power() saves an addition operation for $(t_{11}+t_{22})^2$, and I apply it to other places just for a neat format...

But I think you are right. I read the source codes of fixed_power() and now I think it may cost more (there are if-else statements, and I am not sure if it will be optimized by the compiler). So I changed to the simpler form, as you suggested.

[YiminJin]

@MFraters The convergence behaviors for more difficult cases are shown below:

It can be seen that the convergence behavior of the current code is better indeed: the unstabilized Newton solver converges for the difficult cases with second-order rate. However, without stabilization the linear solver takes more iterations to converge, and in the most difficult case (Vel = 12.5 mm/y) the cheap iterations are even exhausted. So it is still worth using the stabilization in larger models.

[YiminJin]

Supplement to my comments on my derivation of the analytical formula of system Jacobian:

The derivation is under plane-strain assumption, but the conclusion also applies to the original 2D formulation. In the original 2D formulation, we have
$(2\Vert\boldsymbol\varepsilon\Vert\Vert\boldsymbol\varepsilon'\Vert)^2 - (\Vert\boldsymbol\varepsilon\Vert^2 + \boldsymbol\varepsilon:\boldsymbol\varepsilon')^2 = -\dfrac{1}{4}(\varepsilon_{xx} + \varepsilon_{yy})^4$.
Hence, $\mathbb{E}$ still has a negative eigenvalue if we use full strain rate in the assembler, and the convergence rate is slowed down by PD stabilization.

[YiminJin]

Eq. (9) in my derivation is calculated with the help of Sympy. The python codes are shown below:

from sympy import symbols, pprint, simplify

# Define the components of 2D strain rate as independent variables
e_xx, e_yy, e_xy = symbols('epsilon_11, epsilon_22, epsilon_12')

# Dalculate the deviatoric strain rate under plane-strain assumption
d_xx = 2 * e_xx / 3 - e_yy / 3
d_yy = 2 * e_yy / 3 - e_xx / 3
d_zz = -(e_xx + e_yy) / 3
d_xy = e_xy

# The minimum eigenvalue of the tangent operator is proportional to
# ||e|| ||d|| - ||e||^2 - e : d
# To determine whether it is positive or negative, we can compare
# the values of ||e|| ||d|| and (||e||^2 + e : d)
e_norm2 = e_xx**2 + e_yy**2 + 2 * e_xy**2
d_norm2 = d_xx**2 + d_yy**2 + d_zz**2 + 2 * d_xy**2
ed = e_xx * d_xx + e_yy * d_yy + 2 * e_xy * d_xy

print('With plane-strain:')
pprint(simplify(4 * e_norm2 * d_norm2 - (e_norm2 + ed)**2))

# Now do the same for deviatoric strain rate without plane-strain assumption
d_xx = e_xx / 2 - e_yy / 2
d_yy = e_yy / 2 - e_xx / 2
d_xy = e_xy

e_norm2 = e_xx**2 + e_yy**2 + 2*e_xy**2
d_norm2 = d_xx**2 + d_yy**2 + 2*d_xy**2
ed = e_xx * d_xx + e_yy * d_yy + 2 * e_xy * d_xy

print('Without plane-strain:')
pprint(simplify(4 * e_norm2 * d_norm2 - (e_norm2 + ed)**2))

The outputs are:

With plane-strain:
     4      3           2     2           3      4
  ε₁₁    4⋅ε₁₁ ⋅ε₂₂   2⋅ε₁₁ ⋅ε₂₂    4⋅ε₁₁⋅ε₂₂    ε₂₂ 
- ─── - ────────── - ────────── - ───────── - ───
   9        9             3           9        9  
Without plane-strain:
     4                 2    2                4
  ε₁₁      3        3⋅ε₁₁ ⋅ε₂₂         3   ε₂₂ 
- ──── - ε₁₁ ⋅ε₂₂ - ────────── - ε₁₁⋅ε₂₂ - ───
   4                    2                  4  

It is easy to see that the expressions with and without plane-strain are equal to
$-\dfrac{1}{9}(\varepsilon_{11} + \varepsilon_{22})^4$
and
$-\dfrac{1}{4}(\varepsilon_{11} + \varepsilon_{22})^4$
respectively. These values characterize the minimum eigenvalues of the tangent operator $\mathbb D = 2\eta\mathbb I + \mathbb E$ with and without plane-strain assumption. They are both negative, but the absolute value of $\mathbb D$ with plane-strain asumption is smaller than that without plane-strain assumption, so the convergence rate should be less affected by the PD stabilization.

[MFraters]

Generally looks good to me. A few small comments to clear things up.

I don't see a significant difference in the unstabilized convergence when it converges (maybe I am missing something), but the stabilized one is definitely better.

Can you fix the tests? Than we can also see what difference it makes there.

@bangerth, can you have a look at the derivations?

[MFraters]

I am personally not really a fan of using functions like this for something so simple. You are basically inlining the following code here: https://github.com/dealii/dealii/blob/93160909dbf3bbfe986ad5320b675737f89d6e00/include/deal.II/base/utilities.h#L944-L961

Since it is an inline constexpr, and with compiler optimizations, it is probably almost as fast, or as fast as writing it out yourself. You can argue both ways about which one is more clear to read, so it is fine to keep it, but I did want to mention it.

[YiminJin]

I did not read the source code of Utilities::fixed_power() when I made the change. Now I think you are right: using the function may not be faster. I have changed it to explicit form in a new commit.

[MFraters]

I don't think I fully understand this. Currently the Newton solver is only stabilized for the incompressible case. Did you mean incompressible or do you say it is now also stabilized for the compressible case?

[YiminJin]

I forget why I wrote this. I deleted the second sentence in a new commit.

[MFraters]

So instead of strain_rate it should be something like consistent_deviator_of_strain_rate or just deviatoric_strain_rate?

[MFraters]

What did you find?

[YiminJin]

I tried the default model runs (using run.sh) for the nonlinear channel flow benchmark and the Spiegelman et al. (2016) benchmark. Here are the results:

1. Nonlinear channel flow

a. Always use full strain rate (FSR):

b. Use full strain rate when not stabilized, and use deviatoric strain rate (DSR) when stabilized:

The figures correspond to Fig. 2 in the original paper:

The results can be summarized as follows:

	without stabilization	with stabilization
without line search	Almost the same in both cases	Harder to converge in both cases
with line search	Slightly better in both cases	Faster with FSR, slightly slower with DSR

2. Spiegelman et al. (2016)

a. Always use FSR:

b. Use FSR when not stabilized, and use DSR when stabilized:

The figures correspond to Fig. 4 in the original paper:

	without stabilization	with stabilization
Vel: 2.5 cm/y, $\eta_{ref}$: 1e23	Slightly faster in both cases	Faster in both cases
Vel: 5 cm/y, $\eta_{ref}$: 1e24	No better, no worse	Much faster in both cases, FSR slightly better
Vel: 12.5 cm/y, $\eta_{ref}$: 5e24	Easier to converge in both cases	No better, no worse

The results show that:

1. The convergence behavior is not affected by changing from 2D tensor to plane-strain tensor;
2. The new scaling factor (merged last year) does lead to faster convergence rate when using PD stabilization; meanwhile, it might also reduce the stability a little bit;
3. Using deviatoric strain rate in the system Jacobian does not improve convergence rate or stability when using PD stabilization; instead, it leads to slower convergence rate in the channel flow benchmark.

I have not figured out why using deviatoric strain rate results in slower convergence in the channel flow benchmark. But, according to the results and the mathematical work, using full strain rate seems to be the right choice. That is also the reason I changed deviatoric strain rate to full strain the in elastic rheology.

There is an issue to be noticed: when using 2D tensor, one can use deviator() as a filter, i.e. apply it repeatedly and obtain the same result; but when using plane-strain tensor, one can only apply consistent_deviator() once, otherwise the resulting tensor would be incorrect. Currently there is no guarantee for this.

What do you think about the results @MFraters @bangerth @gassmoeller ? Do I need to carry out experiments with elastic rheology? (I have tested the full strain rate version with associated-plastic-flow tests (strip-footing, pure shear), but not the deviatoric strain rate version).

[YiminJin]

Sorry, I made a big mistake in the last comment: I forgot that the 2D tensors in the channel flow benchmark have not been changed to plane-strain tensors.

When I tried to do the modification, I found something that looks weird: line 183 and line 189 of benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/simple_nonlinear.cc calculate the viscosity and the viscosity derivatives according to the first two equations in Appendix B1 of the original paper:
$\eta = \eta_0^{-1/n}\varepsilon_{II}^{1/n-1}$,
$\dfrac{\partial\eta}{\partial\boldsymbol\varepsilon} = \eta\bigg(\dfrac{1}{n}-1\bigg)\dfrac{\boldsymbol\varepsilon}{\Vert\boldsymbol\varepsilon\Vert^2}$
However, line 180 multiplies edot_ii by 2, which leads to
$\eta = \eta_0^{-1/n}(2\varepsilon_{II})^{1/n-1}$
Is there something I am missing? @MFraters

[YiminJin]

I also made a mistake when computing the deviatoric strain rates in system Jacobian: I forgot to change deviator() to Utilities::Tensors::consistent_deviator(). Here are the correct results for DSR in Spiegelman et al. (2016) benchmark:

It is almost the same as the full strain rate case. I also compared the distribution of SPD factor, and the results are again the same. I think this is because the material model is incompressible, hence $\dfrac{\eta}{\boldsymbol a : \boldsymbol b}$ is very close to -1, which leads to
$\alpha_{SPD} = -c_{safety}\dfrac{\eta}{\boldsymbol a : \boldsymbol b}\approx c_{safety}$
in the yielding regions.

@gassmoeller Do you remember how we came to the conclusion that using deviatoric strain rate when stabilized is faster than using full strain rate?

[YiminJin]

I have tried to test the VEP rheology with plastic dilation, and I found some problems in source/simulator/assemblers/newton_stokes.cc:

1. Function NewtonStokesCompressibleStrainRateTerms::execute() does not multiply phi_p with pressure_scaling;
2. The function does not symmetrize the top left block when Stabilization::symmetric is set;
3. The function does not take the spd factor into account.

I also found that the Newton Stokes assembler uses StokesCompressiblePreconditioner, which is inconsistent with the system Jacobian.

After fixing all the problems, and adding a NewtonStokesCompressiblePreconditioner, the solving speed of the strip footing problem improved about 40%.

And, from the compressible case, I confirmed the argument that we should use full strain rate in NewtonStokesIncompressibleTerms. But I cannot be sure, since I keep making mistakes on this issue.

[gassmoeller]

@YiminJin Concerning the deviatoric strain rate: I think the relevant PR where I tested this was #6159 (comment), and in particular this code: https://github.com/geodynamics/aspect/pull/6159/files#diff-67a939b674da49b34479b331895052a08d8712da7befddadaba9390cd782ee05R183. As far as I can tell we did not decide to use the deviatoric strain rate for the stabilized version because it is faster, it seems like there was a theoretical argument that the stabilization factor was only safe if applied to a deviatoric strain rate. In other words the guarantee that the stabilization makes the matrix positive definite is only valid if the used strain rate is deviatoric. But I have to say I dont remember the details or if I ever made that derivation myself (I would guess I accepted it from the paper or some earlier PR discussion). From looking at the earlier PR my main point was that we need to use the full strain rate for the unstabilized version in order to guarantee the expected convergence rate of the Newton method. If your version here works with the full strain rate for both versions so much the better.

[YiminJin]

I apologize for making so many mistakes during the discussion.

During the morning meeting, the documentation and the conclusion were correct, but I failed to explain it. I added some comments to the derivation to make it clearer:
fsr_vs_dsr.pdf
The point is, we have removed the volumetric components in class NewtonStokesCompressibleStrainRateTerms, so we should use full strain rate in class NewtonStokesIncompressibleTerms.

Yet, I understand that we should use full strain rate does not mean that using full strain rate is faster, because using full strain rate breaks the symmetry and semi-positive-definiteness of the top-left block. I have tested the Spiegelman benchmark and the result show that using full strain rate (without stabilization) leads to a little more linear iterations but much fewer nonlinear iterations (similar to the results shown in #6160 ). But one benchmark is not enough.

@MFraters Could you help me with the nonlinear channel flow benchmark? I do not know if I misunderstand your implementation.

[YiminJin]

I plot the number of cheap linear iterations in each Newton iteration for the Spiegelman benchmark:

1. Use full strain rate in NewtonStokesIncompressibleTerms:

a. residuals:

b. number of cheap iterations:

2. Use deviatoric strain rate in NewtonStokesIncompressibleTerms:

a. residuals:

b. number of cheap iterations:

Remark:

1. The cheap iterations are limited by 1000;
2. Five cases report linear solver errors (The iterative (top left) solver did not converge):
   Vel = 12.5 mm/yr, eta_ref = 5e24, mLT = 1e-8, 0 Picard, unstabilized, with full strain rate;
   Vel = 12.5 mm/yr, eta_ref = 5e24, mLT = 1e-8, 5 Picard, unstabilized, with full strain rate;
   Vel = 12.5 mm/yr, eta_ref = 5e24, mLT = 1e-8, 25 Picard, unstabilized, with full strain rate;
   Vel = 12.5 mm/yr, eta_ref = 5e24, mLT = 1e-8, 0 Picard, unstabilized, with deviatoric strain rate;
   Vel = 12.5 mm/yr, eta_ref = 5e24, mLT = 1e-8, 5 Picard, unstabilized, with deviatoric strain rate;

It can be seen that:

1. When not stabilized, using full strain rate generally takes fewer Newton iterations than using deviatoric strain rate, but the convergence is more difficult (takes more cheap linear iterations, and easier to fail);
2. When stabilized, the convergence behaviors of full strain rate and deviatoric strain rate are similar.

Fact 1 confirms the point that using full strain rate breaks the symmetry and semi-positive-definiteness of the top-left block. In the easiest case (Vel = 2.5 mm/yr, eta_ref = 1e23), the convergence curve of stabilized and unstabilized cases are almost the same when using deviatoric strain rate, which suggests that using deviatoric strain rate is equivalent to symmetrizing the matrix (changing the descending direction, slowing down the convergence, but restoring the symmetry and makes the linear solver stabler).

The reason that using deviatoric strain rate or not affects the convergence behavior when using incompressible model is unclear. I guess it comes from numerical error. The distribution of strain rate and velocity divergence in case {Vel = 5 mm/yr, eta_ref = 1e24, mLT = 1e-8} are:

It can be seen that the numerical error (velocity divergence) is in the same magnitude as the strain rate. So the absolute value of the negative eigenvalue of the tangent operator when using full strain rate is large at some nodes.

[MFraters]

I apologize for making so many mistakes during the discussion.

That is normal. I really appreciate all the work you put into this!

During the morning meeting, the documentation and the conclusion were correct, but I failed to explain it. I added some comments to the derivation to make it clearer: fsr_vs_dsr.pdf The point is, we have removed the volumetric components in class NewtonStokesCompressibleStrainRateTerms, so we should use full strain rate in class NewtonStokesIncompressibleTerms.

Yet, I understand that we should use full strain rate does not mean that using full strain rate is faster, because using full strain rate breaks the symmetry and semi-positive-definiteness of the top-left block. I have tested the Spiegelman benchmark and the result show that using full strain rate (without stabilization) leads to a little more linear iterations but much fewer nonlinear iterations (similar to the results shown in #6160 ). But one benchmark is not enough.

This is a difficult one. In the end the Newton derivatives do not decide whether you converge to the correct solution, only how fast. So in the end for us, it doesn't matter too much what the "correct" method of computing the derivative is, but what derivative results in the the fastest and most reliable convergence.

The problem you state is that method which results in the fastest convergence (full strain-rate) is not the same as the method which results in the most reliable convergence (deviatoric strain rate) in the unstabilized case.

When not stabilized, using full strain rate generally takes fewer Newton iterations than using deviatoric strain rate, but the convergence is more difficult (takes more cheap linear iterations, and easier to fail);

The difference in the amount of failure cases (2 vs 3) seems minimal to me. I think there is a third thing we should look at, which is not converging to the tolerance. Based on your figures, it seems to me that using the full strain-rate has a higher chance of actually converging to the required tolerance.

So, based on these results, I would say my preference is on the full strain-rate.

The reason that using deviatoric strain rate or not affects the convergence behavior when using incompressible model is unclear. I guess it comes from numerical error.

Could be, although it seems systematic. Maybe it could have something to do with whether every point is the cell is guaranteed to be divergence free, or just the cell as a whole. @bangerth, do you have an idea?

@MFraters Could you help me with the nonlinear channel flow benchmark? I do not know if I misunderstand your implementation.

Sure, could you be a bit more specific what the issue is? In general, there are two version, one with a prescribed velocity (input_v.prm) and one with a prescribed traction (input_t.prm). You need to run cmake to compile the material model and then you should be able to use the run script.

[YiminJin]

@MFraters Thank you for the detailed comment! I prefer the full strain rate too.

The issue with the nonlinear channel flow benchmark (with prescribed velocity) is that:

line 183 and line 189 of benchmarks/newton_solver_benchmark_set/nonlinear_channel_flow/simple_nonlinear.cc calculate the viscosity and the viscosity derivatives according to the first two equations in Appendix B1 of the original paper:
$\eta = \eta_0^{-1/n}\varepsilon_{II}^{1/n-1}$,
$\dfrac{\partial\eta}{\partial\boldsymbol\varepsilon} = \eta\bigg(\dfrac{1}{n}-1\bigg)\dfrac{\boldsymbol\varepsilon}{\Vert\boldsymbol\varepsilon\Vert^2}$
However, line 180 multiplies edot_ii by 2, which leads to
$\eta = \eta_0^{-1/n}(2\varepsilon_{II})^{1/n-1}$
Is there something I am missing?

[YiminJin]

I tested the strip footing model:

The Newton solver parameters are set as follows:

Parameter	Value
Max nonlinear iterations	100
Nonlinear solver tolerance	1e-5
Max pre-Newton nonlinear iterations	3
Max Newton line search iterations	3
Maximum linear Stokes solver tolerance	1e-2
Use Newton residual scaling method	false
Stabilization preconditioner	SPD
Stabilization velocity block	SPD

The convergence behaviors of the full strain rate case and the deviatoric strain rate case are as follows:

Clearly, the full strain rate prevails. It can be explained by the fact that in compressible/prescribed-dilation models, the assembler NewtonStokesCompressibleStrainRateTerms has already removed the volumetric part, so if we use deviatoric strain rate in NewtonStokesIncompressibleTerms, then we actually minus the volumetric part twice, which does not matter in 2D or 3D cases, but is incorrect in plane-strain cases.

[gassmoeller]

Thanks for putting all the effort into this! I will take another look at the main changes again, but here are a few small comments for now.

[gassmoeller]

Can you add a sentence about which models will be affected? E.g.:

Suggested change

	Changed: The deviator and second invariant of symmetric tensors are modified to
	be consistent with the plane strain assumption in 2D.
	Fixed: The deviator and second invariant of symmetric tensors were modified to
	be consistent with the plane strain assumption in 2D. This fixes a bug
	in the strain rate and strain rate invariant computation of many material
	models for compressible 2D models. It also results in better
	nonlinear solver convergence due to some fixes to the Newton solver.

[YiminJin]

Thank you for providing a nice example! However, there are some statements I am not quite sure of:

1. I do not know if the original implementation should be called a ``bug'', since ASPECT has never declared that it applies plane-strain assumption in 2D.
2. The modifications do not only changes compressible models, but also models including prescribed dilation. Basically, all the material models that calls deviator() are affected.
   Based on these points, I modified the statements a little in a new commit. Please help me revise them if the words are inappropriate.

[gassmoeller]

I think this comment is not addressed yet. The name of the variable is edot_deviator, but it now uses the full strain rate. Does that name of the variable have to be changed?

[gassmoeller]

Since you opened that pull request we also added the file stokes_matrix_free_global_coarsening. Could you fix it there as well?

[YiminJin]

Thanks for finding the flaws. The variable name has been fixed.

The GMG solver requires more attention. Besides changing deviatoric strain rate to full strain rate, I found two other problems:

1. The Newton derivatives related to the volumetric part is missing;
2. When material averaging is set to none, the Newton scheme does not work because the uninitialized viscosity_derivative_averaging_weights are assembled into the system.

I fixed the problems in a new commit. However, I cannot make the strip footing model run with GMG solver. The Kaus benchmark (with $\phi=\psi=30^\circ$) works, but it also requires Newton residual scaling method to achieve convergence. Below are the convergence behaviors with and without the Newton derivatives related to the volumetric part (only nonlinear iterations 10 to 19 are presented):

The prm file is:
kaus_2010_extension.prm.txt

It can be seen that adding the volumetric part does not lead to a big improvement in convergence rate (except for the 18th nonlinear iteration). Nevertheless, I feel it is safer to follow the ``correct'' formula.

Please let me know if there are problems in my codes, or more benchmarks to run.

[gassmoeller]

This looks generally in very good shape and is close to merging. Thanks for working through this complicated change and fixing other bugs on the way 👍.

I have one remaining question about something in the code, and a few minor questions on the test results, but before you work through the tests, check my last comment. It looks to me like you updated the test results with a wrong deal.II version, and this could explain some of the changes.

Let me know when you want me to take another look.

[gassmoeller]

There is something here I dont understand. As far as I understand newton_factor_wrt_strain_rate_table(cell,q) contains the full value if the newton factors are not averaged, but if the newton factors are averaged, it contains newton_factor[q] * averaging_weight[q] (see here). This is why with averaging here the contributions of all quadrature points are added up to the total contribution, while without averaging here only the full value of the newton_factor_wrt_strain_rate_table is used to compute deta_deps_times_sym_grad_u. If I understand this correctly then I am confused why here, where we compute velocity_terms we simply use newton_factor_wrt_strain_rate_table in either case. Are we not using a too small number here if the newton factors are averaged?

If I am correct and this is an oversight I think there are three possible solutions:

1. Either we need to compute a correctly cell-averaged value of newton_factor_wrt_strain_rate_table(cell,q) earlier in this function, and then use this value to compute both deta_deps_times_sym_grad_u and velocity_terms.
2. Or we need to always store the full newton_factor_wrt_strain_rate_table in the variable of this name, and store the averaging weights in a separate variable. Then we use newton_factor_wrt_strain_rate_table(cell,q) * averaging_weights(cell,q) wherever necessary.
3. Or we can perform the averaging already in stokes_matrix_free_local_smoothing.cc::evaluate_material_model, i.e. if we average we store the same (averaged) number in all entries of newton_factor_wrt_strain_rate_table(cell,q).

If we only need the averaged value of newton_factor_wrt_strain_rate_table I would prefer to choose option 3, i.e. perform the averaging immediately where these values are computed. However, if we also need the not averaged value, we will have to pick either 1 or 2.

Did I understand this correctly? I am not an expert in this code and might have misunderstood something. Maybe @MFraters or @tjhei also know more about this.

[YiminJin]

Thanks for finding this! It is a bug indeed, and it comes from a trick I used to reduce the computation. Let me break it down:

The top-left block of the system Jacobian w.r.t. incompressible terms can be expressed as
$$2\sum_qw_qJ_q\nabla_s\boldsymbol\varphi_{Iq} : \boldsymbol\varepsilon_q\otimes\dfrac{\partial\eta_q}{\partial\boldsymbol\varepsilon_q} : \nabla_s\boldsymbol\varphi_{Jq}$$
When using cellwise averaging, the expression becomes
$$2\sum_qw_qJ_q\nabla_s\boldsymbol\varphi_{Iq} : \boldsymbol\varepsilon_q\otimes\sum_r\tilde w_r\dfrac{\partial\eta_r}{\partial\boldsymbol\varepsilon_r} : \nabla_s\boldsymbol\varphi_{Jr}$$
where $\tilde w_r$ is the averaging weight. Now do the symmetrizing:
$$\sum_qw_qJ_q\nabla_s\boldsymbol\varphi_{Iq} : \boldsymbol\varepsilon_q\otimes\sum_rw_r\dfrac{\partial\eta_r}{\partial\boldsymbol\varepsilon_r} : \nabla_s\boldsymbol\varphi_{Jr} + \sum_r\tilde w_r\nabla_s\boldsymbol\varphi_{Ir} : \dfrac{\partial\eta_r}{\partial\boldsymbol\varepsilon_r}\otimes\sum_qw_qJ_q\boldsymbol\varepsilon_q : \nabla_s\boldsymbol\varphi_{Jq}$$
Let us focus on the second term. It can be rewritten as
$$\sum_q\sum_rw_qJ_q\tilde w_r\nabla_s\boldsymbol\varphi_{Ir} : \dfrac{\partial\eta_r}{\partial\boldsymbol\varepsilon_r}\otimes\boldsymbol\varepsilon_q : \nabla_s\boldsymbol\varphi_{Jq}$$
I thought that since we sum over $q$ and $r$, we can swap the two subscripts:
$$\sum_q\sum_r\color{red}w_qJ_q\color{black}\tilde w_q\nabla_s\boldsymbol\varphi_{Iq} : \dfrac{\partial\eta_q}{\partial\boldsymbol\varepsilon_q}\otimes\boldsymbol\varepsilon_r : \nabla_s\boldsymbol\varphi_{Jr}$$
This is the current implementation. However, I forgot that the subscripts of $w_qJ_q$ should also change. In the matrix-free framework, there is no way to do this. So the trick is wrong.
Moreover, I forgot this issue for the compressible terms. I think these are the reasons for the GMG solver to be unstable right now.
I will try your suggestions to fix this problem.

[YiminJin]

After further study I found that the current implementation is not just a "trick", but the only way I found to do the symmetrization.

First, we cannot pre-compute the averaged values in StokesMatrixFreeHandlerLocalSmoothingImplementation::evaluate_material_model(), because the quantity to be averaged is the product of $\partial\eta/\partial\boldsymbol\varepsilon$ and $\nabla_s\boldsymbol\varphi_J$, and in the matrix-free framework $\boldsymbol\varphi_J$ is bound to the solution vector (which changes after each linear iteration). Thus, we must do the averaging in local_apply() and try to separate $w_qJ_q\nabla_s\boldsymbol\varphi_{Iq}$ from the rest terms.

Currently, the only way I find is to rewrite the second term by
$\sum_q\sum_r w_qJ_q \dfrac{w_rJ_r}{w_qJ_q}\tilde w_q\nabla_s\boldsymbol\varphi_{Iq} : \dfrac{\partial\eta_q}{\partial\boldsymbol\varepsilon_q}\otimes\boldsymbol\varepsilon_r : \nabla_s\boldsymbol\varphi_{Jr}$
and the symmetric tensor to submit (to be multiplied by $w_qJ_q\nabla_s\boldsymbol\varphi_{Iq}$) for each quadrature point q is
$\tilde w_q\dfrac{\partial\eta_q}{\partial\boldsymbol\varepsilon_q}\dfrac{\sum_rw_rJ_r\boldsymbol\varepsilon_r : \nabla_s\boldsymbol\varphi_{Jr}}{w_qJ_q}$
Similarly, the term corresponding to compressibility is
$-\dfrac{\tilde w_q}{3}\dfrac{\partial\eta_q}{\partial\boldsymbol\varepsilon_q}\dfrac{\sum_rw_rJ_r(\nabla\cdot\boldsymbol u_r)(\nabla\cdot\boldsymbol\varphi_{Jr})}{w_qJ_q}$

I tried this, but the convergence behavior is not improved. I still cannot run the strip footing model even with Newton residual scaling method.

Are there mistakes in my derivation @bangerth @tjhei ?

The pdf below is the original documentation for the viscosity averaging issues in Newton solver:
cell_averaging.pdf

[YiminJin]

My apologies: I found the mistake myself. It is not in the math, but in the code: I mistakenly used trace(sym_grad_u) to represent the velocity divergence, but it should be trace(cell_data->strain_rate_table(cell,q)). After fixing this, the GMG Newton solver converges for the strip footing model (although not as fast as the AMG solver), even without cellwise averaging.

I will deal with the other problems later.

[gassmoeller]

This is the part I am concerned about. If we do not average, newton_factor_wrt_strain_rate_table contains the full value. If we average, it only contains the contribution of quadrature point q (viscosity_derivative_wrt_strain_rate[q] * viscosity_derivative_averaging_weights[q]). I think this is dangerous because it changes the meaning of the variable without changing the name. If we can change it, we should perform the averaging here, and store the same value in every quadrature point.

[gassmoeller]

See my comment below about the derivative with respect to velocity, the same actually applies to pressure. I dont think this has the same problem in the matrix free operator, because this variable seems to be only ever used to compute deta_dp_times_p and not by itself, but maybe check again.

[gassmoeller]

I dont understand why this test output changed. As far as I can see this test should use particles (which is also indicated by the output for the earlier timestep). Why is it now solving composition equations for the fields?

[gassmoeller]

Something in this test went wrong. Before this PR the analytical and finite difference derivative matched pretty accurately, but now there is a significant difference (factor ~1.5x). Is maybe the analytical formula that is in the test .cc file wrong?

[gassmoeller]

Where does this nan come from? It seems suspicious.

[gassmoeller]

nan here as well.

Oh, this is the ODE solver iterations. This was fixed in deal.II 9.6. Did you by any chance update the test results with a deal.II version 9.5? The correct way to do this would be to download the suggested test results from our reference tester with deal.II 9.6. It is not shown at the bottom of your PR at the moment, because of the conflict. But if you resolve the conflict, look for the tester geodynamics/aspect-tester:jammy-dealii-9.6-v1 and use its output to update the test results.

Disregard my other comments on the test results until you have updated the test results from that tester.