Choosing a Material Model for Boussinesq


I’m an undergrad student trying to set up a relatively simple 3D spherical shell computation using the Boussinesq approximation and non-dimensional input parameters to reproduce a published paper’s results. The paper implements a variable, temperature-dependent viscosity (Frank-Kamenetskii rheology, if that helps describe my case) that I need to work out in ASPECT. However, I cannot find an appropriate material model that I can both: 1) fit the viscosity scaling to, and 2) apply the paper’s parameters as inputs.

In the ASPECT manual, the “nondimensional” model’s viscosity equation can be made to fit the viscosity I want, and I seem to be able to provide the input parameters. However, although in the “model name” subsection of the manual it mentions the model being built for Boussinesq, ALA or TALA approximations, there is no input parameter to specify to use a Boussinesq equation of state (only ALA or TALA). Thus, there is an unwanted depth-dependence in the density formula used by the model.

The density formula I want is the same as in the “simpler” model.

How can I set up the material model I want (Boussinesq state, temperature-dependent viscosity)? Any help would be greatly appreciated.


Hi Ian,

thank you for using ASPECT and welcome to the community!
There is a separate part of the input file that lets you specify what approximation of the equations you want to use, so this is not part of the material model. In order to use the Boussinesq approximation, you would add

subsection Formulation
  set Formulation = Boussinesq approximation

to your input file. This is independent of the material model, so in principle, you can use any material model you want with that. It sounds like the nondimensional model would be the easiest way forward (but you could also use the simpler model if that works better). For a Boussinesq model, Di should be zero, and this should remove the pressure-dependence in the density.

Hope that helps!

Hi Juliane,

Thank-you, and I appreciate the welcome. Your advice was indeed helpful to move forward with my project.

I ran a test calculation and encountered this error saying the Stokes iterative solver didn’t converge. Could I get some help in identifying what I did that ASPECT doesn’t like? I looked at the convergence history and the change after each iteration at the end was sub-decimal place.

The following lines are from the standard output error:

An error occurred in line <1264> of file </scinet/niagara/software/2018a/opt/intel-2018.2-intelmpi-2018.2/aspect/2.0.1/
deal.II-v9.0.0/include/deal.II/lac/solver_gmres.h> in function
void dealii::SolverFGMRES<VectorType>::solve(const MatrixType &, VectorType &, const VectorType &, const Preconditi
onerType &) [with MatrixType = aspect::internal::StokesBlock, PreconditionerType = aspect::internal::BlockSchurPrecondi
tioner<dealii::TrilinosWrappers::PreconditionAMG, dealii::TrilinosWrappers::PreconditionBase>, VectorType = dealii::Tri
The violated condition was:
Additional information:
Iterative method reported convergence failure in step 1000. The residual in the last step was 488.312.

This error message can indicate that you have simply not allowed a sufficiently large number of iterations for your ite
rative solver to converge. This often happens when you increase the size of your problem. In such cases, the last resid
ual will likely still be very small, and you can make the error go away by increasing the allowed number of iterations
when setting up the SolverControl object that determines the maximal number of iterations you allow.

The other situation where this error may occur is when your matrix is not invertible (e.g., your matrix has a null-spac
e), or if you try to apply the wrong solver to a matrix (e.g., using CG for a matrix that is not symmetric or not posit
ive definite). In these cases, the residual in the last iteration is likely going to be large.

The parameters I ran were altered from a successful test-run using the material model ‘simpler’. The changes I made were:

  • Material model to ‘nondimensional’
  • inner/outer radii to 0.6 and 1.6, to describe a system with depth =1
  • boundary temperatures to 1 and 0
  • gravity to equal the Rayleigh number, compensating for the reduction to 1 of other values

Do I perhaps need to lower the CFL number? I’m honestly stuck, since I don’t see how my changes relate to the error information. I hope this isn’t too much to ask about.


Hi Ian,

you are right, if you get an error message that the solver didn’t converge, it may be hard to tell what caused the problem. As you have a setting that worked for a test run, I would go back to the setting that worked, and then make the changes you listed below one by one so that you can find out which one caused the problem. You can run these models in a low resolution.

If that doesn’t help, the other strategy I use to spot what’s wrong is to look at the graphical output of the temperature/density/viscosity etc. to see if something looks weird. Does your model crash in the very first time step before it produces any output? In that case, you can use the Nonlinear solver scheme = single Advection, no Stokes. This means you skip the part that produces the error messages, so you don’t get the right solution, but you can at least generate graphical output to check everything is correct with your material properties and initial conditions.

Good luck!