geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
224 stars 235 forks source link

The definition of Stokes residual output? #2269

Open Shangxin-Liu opened 6 years ago

Shangxin-Liu commented 6 years ago

When I encounter convergence failure, there is an error message shown in my file saying (I took one of my convergence failure tests for example):

"Additional information: Iterative method reported convergence failure in step 10000. The residual in the last step was 6.61287e+17."

I notice that the residual 6.61287e+17 shown here is also the residual shown at the last iteration step in the file solver_history.txt:

9997 6.61287e+17 9998 6.61287e+17 9999 6.61287e+17 10000 6.61287e+17

I wonder what's the definition of the Stokes residual here? I cannot find it in the manual.

Is it the absolute Stokes residual like || M x - F ||(x the velocity and pressure solution) or normalized relative Stokes residual? How is it computed from velocity and pressure solution at each Stokes iteration step?

I feel that I cannot get a sense of how hard my problem is to converge without knowing the definition of the Stokes residual here.

gassmoeller commented 6 years ago

Hi Shangxing, this residual is the norm of the residual of the full Stokes equations, it is not normalized or relative. In addition, it is not "only" the Stokes residual, it is the residual of the reformulated equation, for which we scale the pressure equation to make it of a comparable magnitude to the velocity equation (see Section 3.2.4 of Kronbichler(2012) for a description of why and how). Thus, the number you see is not helpful at all to tell you "how hard" your problem is. The only thing you can see from it is how far the solver is converging (i.e. how much did the residual decrease since the first iteration, and is it close to the linear solver tolerance you selected in the input file). If your residual started at e+23 and ends at e+17 with a tolerance of 1e-7 you are nearly converged and the solver just does not manage the last order of magnitude, if you started at e+18 then you barely reduced the residual at all, and you should look at which component of your model makes it so hard to solve (i.e. viscosity contrast, resolution).

Shangxin-Liu commented 6 years ago

Thanks for the explanation, Rene. I will also take a look at the paper.

bangerth commented 6 years ago

@Shangxin-Liu -- is the error message you show in your original post the only thing you get to see? I thought that the exception also shows the first residual, so that you can compute by how much the residual has already been reduced (as @gassmoeller already mentioned).

Shangxin-Liu commented 6 years ago

@bangerth Yes. I only showed part of the solver history in the first email. The first row residual at 0 iteration step before the 200 cheap iterations is 1.24211e+20 and I set the linear solver tolerance to 1e-3.

0 1.24211e+20 1 1.2421e+20 2 1.2421e+20 3 1.24209e+20 4 1.24205e+20 5 1.24202e+20 .......

The final residual is 6.61287e+17 so the residual has been reduced by over 2 order of magnitudes, as Rene mention. But one thing I noticed is that in the expensive iterations, from the iteration steps 5091, the residual maintains constantly at the value of 6.61287e+17 until the final 10000 iterations. So it seems that this is the minimum residual the solver can reach in this problem.

So is the first row residual at "0 cheap iteration step" computed from the initial guess of the velocity and pressure, i.e., zero velocity and the pressure profile set up from adiabatic conditions model?

bangerth commented 6 years ago

What I had meant is that I had remembered that deal.II doesn't only print the last, but also the first residual of an iterative method so that one can easily see by how much it has been reduced. But I looked it up in deal.II/include/deal.II/lac/solver_control.h, and it does not in fact do that. But you found a way to output the solver history which also contains this.

The information in the solver history file is produces exclusively in GMRES. So yes, the one for iteration zero is really computed from the initial guess of velocity and pressure. If you are in time step zero and nonlinear iteration zero, this is indeed using a zero velocity and the (scaled) adiabatic pressure. In later time steps of nonlinear iterations, we use the solution from the previous time step(s) to get a good starting point if I recall correctly.