Velocity and Pressure L2 norm doesn't converge with mesh refinement

Hi all,

I am trying to implement a 3D Cartesian exponentially varying viscosity benchmark on ASPECT. I did get the code running. However, my velocity and pressure L2 norm stay almost the same with increasing refinements (only uniform applied). I am wondering if anyone has any suggestion. Here I also attach some of my data.

refinement Velocity L2 norm Pressure L2 norm
3 2.832740e-02 1.420663e+04
4 2.271762e-02 1.420674e+04
5 2.261151e-02 1.420675e+04
6 2.260930e-02 1.420675e+04

I have a feeling that if my implementations of benchmark and boundary condition are correct, then the velocity L2 norm could be due to the solver while the pressure L2 norm could be due to the pressure normalization.

Thanks in advance,

Hi Ray,

Thanks for posting this to the forum!

I have a feeling that if my implementations of benchmark and boundary condition are correct, then the velocity L2 norm could be due to the solver while the pressure L2 norm could be due to the pressure normalization.

Can you describe the benchmark a bit further, in particular the variations in viscosity (magnitude, spatial changes), boundary conditions, solver tolerances, pressure normalization, etc?

Are you seeing any odd behavior in the linear solvers and does the problem require a non-linear solve?

There have been quite a few studies looking at how the velocity and pressure L2 norm vary as function of grid resolution, so in principal I agree you should be seeing better improvement with increasing refinement level.


It wouldn’t be the first time that I’ve seen this sort of thing in my own work
:slight_smile: A good first step would be to visualize the solution you get and compare
it (visually) with what you expect. That might help you decide which of these
might the the case:

  • The exact solution you implemented is wrong
  • The computed solution is wrong
    Both of course could be the case!

If you find that the computed solution is wrong, pay attention to where it
is wrong:

  • If time dependent, are the initial conditions right at least?
  • Are the boundary conditions correct?
    If both of these seem correct, let ASPECT also output the coefficients you are
    using for you, and again compare that against what you think they should be.


Hi John and Wolfgang,

Thanks for the messages! The benchmark I have been implementing is based on a paper by Popov etal (2014). I have been using 1e-12 as my linear solver tolerance, but the result doesn’t change much when I increase to 1e-7. Also the pressure normalization is done throughout the volume (I also tried surface pressure normalization, but it didn’t do much either).

The problem does require non-linear solver (no Advection, iterated Stokes). In addition, I have a bit silly question, that the output often shows the message like this:

Rebuilding Stokes preconditioner…
Solving Stokes system… 69+0 iterations.

I am wondering what the two numbers represent (iteration cheap and expensive solver?).

I will try with Wolfgang’s approach and hopefully find some clues.


Yes. The solver converged with the cheap solver already and didn’t need any
iterations with the expensive solver.

I happen to have reviewed this paper. Popov et al, Solid Earth 5, 2014. Given the complexity of the calculations in the first 5 pages or so, I did not re-derive every single equation of the manuscript myself, but instead assumed that the convergence rates they showed could not be a happy coincidence. However, I remember vividly sending lots of comments back about the quality of the plots and errors on plot labels, captions, and much more. However it was nearly published as is. How can a serious journal allow figures such as this ?

Having said this, there is already in ASPECT a 3D exponential viscosity benchmark (the ‘Burstedde’ one), so what is it in this benchmark that is of interest to you? Also, have you verified that the matlab codes in the supplementary material correspond 100% to the equations in the paper?

Hi Cedric,

Thanks for your message. I happened to come across this paper and thought it would be interesting to see how ASPECT recovers it. I never thought this would trigger such hidden informations.

Indeed, the supplementary matlab code is quite inconsistent with the paper in terms of choices of parameters. I may need to look more into this to see if this is actually due to the inaccuracy of analytical solutions.



Just a brief update. For the analytical solution, it seems like the solution doesn’t subtract the constant that is big enough so the mean pressure can go to zero (difficult to compare with the numerical solution after pressure normalization, but the numerical result is still a bit different without normalization). I haven’t figured out the velocity yet, but I would suspect something is off in the derivations (i.e. integration constant).

I also have another question. I noticed that for benchmark templates provided by ASPECT, the solutions tend to have an extra component such that:

  ComponentSelectFunction<dim> comp_u(std::pair<unsigned int, unsigned int>(0,dim),
  ComponentSelectFunction<dim> comp_p(dim, dim+2);

I am always curious why the number of components is ‘dim + 2’ instead of ‘dim + 1’, as I tend to think that velocity has 3 components ( in 3D case) and pressure has 1. Thus where is this extra component coming from?

Thanks in advance,

There is a parameter that determines how the pressure is to be normalized. By
default, the average pressure on the surface of the model is set to zero. If
you want the average pressure in the domain of the model to be zero, you need
to change that parameter.