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]

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]

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?

[bangerth]

Did you have the time to look into this? In fact, I'm surprised we ever had the deviator here. We multiply the right hand side by the test function, not the deviator of anything related to the test function. I see no real good reason to take the deviator here.

[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.

[bangerth]

I read through the entire history of this PR again -- what a heroic effort :-) Is it correct? I'm not sure any of are quite sure about this, but I think it goes in the right direction.

At its core, this patch does two things:

• It fixes the issue with the definition of the deviator in 2d. I think we're all convinced that that is independent right, and if you want a suggestion for future patches: try to break things into smaller bits -- the introduction and use of the consistent deviator by themselves are probably not controversial, and could have been merged much quicker.
• The changes to the Newton solver. That one is perhaps not as clear. I think what it comes down to is that it is an issue of documentation: What exactly is it that a function takes as argument. For example, if we have f(s) and document that s is a stress tensor, then if the function internally uses the deviatoric stress tensor D s (with the deviator tensor), say f(s) = (D s) : (D s), then if we ask for the derivative df/ds, the function needs to return (D s) : D + D : (D s). If we say that the function takes the deviatoric stress as argument, then the function would compute its value as f(s) = s:s and its derivative is s : I + I : s. But since in 3d, we have that DD=D, the function could also return (D s) : D + D : (D s) and that would also be correct: wherever the derivative will be used, it will be used on deviatoric tensors, and the multiplication with D will not make a difference. I think we may simply have been sloppy in defining what functions take as arguments, and consequently what we expect the derivative to be. I imagine that we just got lucky in the past because DD=D, and so whether or not the function took one or the other as input, and whether it had D in the return value did not matter. It's just that that is now no longer the case, and we need to be more concerned about correctness. I don't think anyone has looked as carefully at the code as @YiminJin , but I also am pretty sure that there are more places in the code that are either outright wrong, or at least ambiguous because incompletely documented.

I'm in favor of merging this patch if you go through the remaining open comments by scrolling top to bottom through https://github.com/geodynamics/aspect/pull/6471/files As mentioned, I'm not completely sure this is all correct, but I think it is clearly better than it was before.

[bangerth]

Can you undo this change?

[bangerth]

Did you have the time to look into this? In fact, I'm surprised we ever had the deviator here. We multiply the right hand side by the test function, not the deviator of anything related to the test function. I see no real good reason to take the deviator here.

[bangerth]

Did you omit the case (velocity_block_stabilization & Newton::Parameters::Stabilization::PD) != Newton::Parameters::Stabilization::none on purpose? The condition is now no longer tested.

[YiminJin]

I think we should always use the full strain rate here. In the test results of the Spiegelman benchmark (as shown above), using full strain rate is faster than using deviatoric strain rate in most cases. Besides, I cannot see the reason from the formulas to use deviatoric stress when doing the PD stabilization (maybe I missed something?)

[bangerth]

You probably should apply the indentation script to undo this change.

[bangerth]

If I see this right, you only use these variables when you do averaging, right? If so, perhaps declare these variables here as empty arrays, and then resize them to the correct size inside the if condition immediately below -- this way, you don't pay the price of memory allocation if you don't actually need the variables.

[YiminJin]

Yes, I think it is better to allocate the space when necessary. There are 3 more places with the same problem (NewtonStokesPreconditioner, NewtonStokesCompressiblePreconditioner, NewtonStokesIncompressibleTerms) and I have modified them all.

By doing so I found a small problem: there are two booleans rebuild_stokes_matrix and rebuild_newton_stokes_matrix in struct internal::assembly::Scratch::StokesSystem. In NewtonStokesIncompressibleTerms we use rebuild_newton_stokes_matrix to determine whether to assemble the matrix or not, but in NewtonStokesCompressibleStrainRateViscosityTerm we use rebuild_stokes_matrix (in the main branch). In fact, I think we do not need two booleans for this. Shall I open an issue for this?

[gassmoeller]

By doing so I found a small problem: there are two booleans rebuild_stokes_matrix and rebuild_newton_stokes_matrix in struct internal::assembly::Scratch::StokesSystem. In NewtonStokesIncompressibleTerms we use rebuild_newton_stokes_matrix to determine whether to assemble the matrix or not, but in NewtonStokesCompressibleStrainRateViscosityTerm we use rebuild_stokes_matrix (in the main branch). In fact, I think we do not need two booleans for this. Shall I open an issue for this?

This may indeed be an issue (I didnt check), but yes, please open a separate issue for it. Maybe @MFraters can look into that, but it shouldnt become part of this PR.

[bangerth]

Separately, (i) you'll need to run the indentation script, (ii) could you reduce the number of commits to either one, or to a set of commits that are individually self-contained (say, one for the new deviator functions, one that uses the new functions, one that updates the Newton scheme, one for the updated tests, or some such)?

[gassmoeller]

@YiminJin will you have time to come back to this PR in the next few weeks? We would like to start preparing the next release towards the end of this month and I would prefer to have this patch in there (at the very least the part about the 2D consistent deviator).

[YiminJin]

@bangerth @gassmoeller I shall apologize for being absent for a long time. I thought I could come back to this after finishing my own tasks during Christmas, but I felt so exhausted that I took a few days' break.

@bangerth Thank you for going through this PR. I have fixed the space-allocating problem in a new commit (there is a small problem. Please take a look at my reply to your comment). About the format problems: I have fixed the problems you pointed out, but I believe there are more problems (like extra white spaces) that cannot pass the format check. However, if I use astyle to reformat the file, there would be some annoying changes. For example, line 89-90 in source/simulator/assemblers/newton_stokes.cc is:

      const std::shared_ptr<const MaterialModel::MaterialModelDerivatives<dim>> derivatives
        = scratch.material_model_outputs.template get_additional_output_object<MaterialModel::MaterialModelDerivatives<dim>>();

After reformatting, it becomes

      const std::shared_ptr<const MaterialModel::MaterialModelDerivatives<dim>> derivatives
                                                                             = scratch.material_model_outputs.template get_additional_output_object<MaterialModel::MaterialModelDerivatives<dim>>();

I think it is because >> is identified as a in/out signal. Are there some better ways to reformat the modified files?

[tjhei]

I think it is because >> is identified as a in/out signal. Are there some better ways to reformat the modified files?

You could use > > and see if that helps.

[YiminJin]

@tjhei Thanks for the suggestion. It circumvents the astyle problem, but I searched all the source files in ASPECT and cannot find a single case of "> >" (other places with ">>" passed the format check, I do not know why). So I downloaded the "changes-astyle.diff" file and made modifications manually.

[tjhei]

That's not how this works. 😀 The diff file is provided for you to make your formatting match the official formatting given by astyle.

I would suggest you accept the fact that one line of code is not formatted exactly the way you would want it to be. After all, the functionality is what counts, not the looks. 😉

[gassmoeller]

I shall apologize for being absent for a long time.

No need to apologize, we all have a lot of other things to do, and we understand if someone needs to focus on other things (work or vacation) for a while. I just wanted to check if we can plan with this PR for the release or not.

It looks like we are almost ready. Like Timo said, dont worry about the indentation for now. Our indentation rule at the moment is: we use whatever astyle 2.04 produces. This looks awkward in a few places, but we just leave it as is for now.

Could you:

• Update the test results as explained at the end of this page: https://aspect-documentation.readthedocs.io/en/latest/user/extending/testing/writing-tests.html
• Squash the number of commits to just one or a few (logically organized) commits

Then we just need to check with @bangerth if his comments were addressed.

[YiminJin]

I have applied the .diff files and now I can see in the summary that .diff files are still there but empty. But there are still failing checks. Why is that?

[YiminJin]

I located the problem under the help of @alarshi : the nonlinear channel flow test failed. Its top-left block has a negative eigenvalue even if the SPD stabilization is applied.

So it comes back to the problem in the channel flow test: I cannot understand its material model. The problem 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}$
Do you have time to take a look at it @MFraters ? I am not sure if I have missed something.