Dear users of ASPECT

I am a new user who tried to learn how to use the program both by reading the manual and by doing simulations of my interest. I am interested in performing simulations regarding the convective motion of the mantle by adding a layer of hyperdense material to see its behavior and evolution during the simulation.

This is well documented in the various cookbooks (for example, active composition), but when I want to do more realistic simulations, for example when I want to simulate the origin and evolution of LLSVP (Large Low Shear Velocity Provinces) alias Superplumes, it becomes impossible to run my simulations that always end in error, for one reason and another (especially my lack of general familiarity with the program).

In particular, I was interested in reproducing the simulations present in the following studies:

“Vastly Different Heights of LLVPs Caused by Different Strengths of Historical Slab Push”. Yuan et al. (2022)

Internal structure of ultralow-velocity zones consistent with origin from a basal magma ocean. Pachhai et al. 2022

In the studies there are the videos of what I more or less want to achieve.

To do this, I tried several ways. First, I tried to modify the cookbook “onset of convection” by inserting more realistic parameters. I couldn’t add a density gradient, but I was able to add the hyperdense layer to the base.

The simulation goes on for about 16 steps, then crashes for the following error, this is what the terminal tells me.

```
--------------------------------------------------------------------------
francyrad@MacBook-Pro-di-Francesco convection_2d % mpirun -np 8 aspect onset_of_convection.prm
-----------------------------------------------------------------------------
-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
-- . version 2.5.0-pre
-- . using deal.II 9.4.2
-- . with 32 bit indices and vectorization level 0 (64 bits)
-- . using Trilinos 12.18.1
-- . using p4est 2.3.2
-- . running in DEBUG mode
-- . running with 8 MPI processes
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
The output directory <convection/> provided in the input file appears not to exist.
ASPECT will create it for you.
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-- For information on how to cite ASPECT, see:
-- https://aspect.geodynamics.org/citing.html?ver=2.5.0-pre&sha=&src=code
-----------------------------------------------------------------------------
Number of active cells: 768 (on 5 levels)
Number of degrees of freedom: 16,838 (6,402+833+3,201+3,201+3,201)
*** Timestep 0: t=0 years, dt=0 years
Solving temperature system... 0 iterations.
Solving C_1 system ... 0 iterations.
Solving C_2 system ... 0 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 27+0 iterations.
Postprocessing:
RMS, max velocity: 0.692 m/year, 2.92 m/year
Writing graphical output: convection/solution/solution-00000
Temperature min/avg/max: 158.7 K, 1309 K, 4000 K
Compositions min/max/mass: 0/2.9e+06/9.082e+18 // 0/2.9e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 5.288e+11 W, -6.07e+04 W
*** Timestep 1: t=21669.5 years, dt=21669.5 years
Solving temperature system... 12 iterations.
Solving C_1 system ... 14 iterations.
Solving C_2 system ... 9 iterations.
Solving Stokes system... 29+0 iterations.
Postprocessing:
RMS, max velocity: 0.16 m/year, 0.614 m/year
Writing graphical output: convection/solution/solution-00001
Temperature min/avg/max: 158.7 K, 1316 K, 4000 K
Compositions min/max/mass: -1.286e+04/3.018e+06/9.989e+18 // -5.7/2.9e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 4.585e+10 W, -7.837e+04 W
*** Timestep 2: t=108380 years, dt=86710.9 years
Solving temperature system... 104 iterations.
Solving C_1 system ... 113 iterations.
Solving C_2 system ... 13 iterations.
Solving Stokes system... 28+0 iterations.
Postprocessing:
RMS, max velocity: 1.67 m/year, 5.81 m/year
Writing graphical output: convection/solution/solution-00002
Temperature min/avg/max: 158.8 K, 1292 K, 4000 K
Compositions min/max/mass: -2.133e+05/3.87e+06/1.175e+19 // -295.8/2.9e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 1.855e+12 W, -7.831e+04 W
*** Timestep 3: t=119769 years, dt=11388.4 years
Solving temperature system... 14 iterations.
Solving C_1 system ... 16 iterations.
Solving C_2 system ... 12 iterations.
Solving Stokes system... 28+0 iterations.
Postprocessing:
RMS, max velocity: 1.07 m/year, 3.65 m/year
Writing graphical output: convection/solution/solution-00003
Temperature min/avg/max: 158.8 K, 1297 K, 4000 K
Compositions min/max/mass: -1.893e+05/3.672e+06/1.187e+19 // -891.5/2.901e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 6.881e+11 W, -7.817e+04 W
*** Timestep 4: t=138750 years, dt=18981.3 years
Solving temperature system... 13 iterations.
Solving C_1 system ... 14 iterations.
Solving C_2 system ... 13 iterations.
Solving Stokes system... 28+0 iterations.
Postprocessing:
RMS, max velocity: 1.5 m/year, 4 m/year
Writing graphical output: convection/solution/solution-00004
Temperature min/avg/max: 158.8 K, 1302 K, 4000 K
Compositions min/max/mass: -1.362e+05/3.317e+06/1.197e+19 // -1.932e+04/2.912e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 4.979e+11 W, -6.859e+04 W
*** Timestep 5: t=159254 years, dt=20503.6 years
Solving temperature system... 27 iterations.
Solving C_1 system ... 24 iterations.
Solving C_2 system ... 22 iterations.
Solving Stokes system... 33+0 iterations.
Postprocessing:
RMS, max velocity: 19.1 m/year, 51.9 m/year
Writing graphical output: convection/solution/solution-00005
Temperature min/avg/max: 158.7 K, 1305 K, 4000 K
Compositions min/max/mass: -8.555e+04/3.102e+06/1.199e+19 // -1.873e+05/3.018e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 4.422e+11 W, 1.407e+06 W
*** Timestep 6: t=161000 years, dt=1745.85 years
Solving temperature system... 24 iterations.
Solving C_1 system ... 23 iterations.
Solving C_2 system ... 19 iterations.
Solving Stokes system... 33+0 iterations.
Postprocessing:
RMS, max velocity: 3.89 m/year, 18.2 m/year
Writing graphical output: convection/solution/solution-00006
Temperature min/avg/max: 158.7 K, 1304 K, 4000 K
Compositions min/max/mass: -7.687e+04/3.092e+06/1.199e+19 // -5.924e+04/3.016e+06/9.082e+18
Heat fluxes through boundary parts: 0 W, 0 W, 2.626e+11 W, 1.499e+06 W
*** Timestep 7: t=165991 years, dt=4991.41 years
Solving temperature system... retrying linear solve with different preconditioner...
---------------------------------------------------------
TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped
to avoid a possible deadlock.
---------------------------------------------------------
----------------------------------------------------
Exception 'ExcMessage (exception_message.str())' on rank 0 on processing:
--------------------------------------------------------
An error occurred in line <2779> of file </Users/francyrad/Documents/aspect_new/source/utilities.cc> in function
void aspect::Utilities::linear_solver_failed(const std::string &, const std::string &, const std::vector<SolverControl> &, const std::exception &, const MPI_Comm &, const std::string &)
The violated condition was:
false
Additional information:
The iterative advection solver in Simulator::solve_advection did not
converge.
The initial residual was: 2.972608e+22
The final residual is: 2.972608e+22
The required residual for convergence is: 3.580309e+12
See convection/solver_history.txt for the full convergence history.
The solver reported the following error:
--------------------------------------------------------
An error occurred in line <1076> of file
</Users/francyrad/dealii-candi/deal.II-v9.4.2/include/deal.II/lac/solver_gmres.h>
in function
void
dealii::SolverGMRES<dealii::TrilinosWrappers::MPI::Vector>::solve(const
MatrixType &, VectorType &, const VectorType &, const
PreconditionerType &) [VectorType =
dealii::TrilinosWrappers::MPI::Vector, MatrixType =
dealii::TrilinosWrappers::SparseMatrix, PreconditionerType =
dealii::TrilinosWrappers::PreconditionILU]
The violated condition was:
iteration_state == SolverControl::success
Additional information:
Iterative method reported convergence failure in step 1000. The
residual in the last step was 2.97261e+22.
This error message can indicate that you have simply not allowed a
sufficiently large number of iterations for your iterative solver to
converge. This often happens when you increase the size of your
problem. In such cases, the last residual 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-space), 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 positive definite). In these cases, the
residual in the last iteration is likely going to be large.
Stacktrace:
-----------
#0 2 aspect 0x00000001046dfbe0
_ZN6dealii11SolverGMRESINS_16TrilinosWrappers3MPI6VectorEE5solveINS1_12SparseMatrixENS1_15PreconditionILUEEEvRKT_RS3_RKS3_RKT0_
+ 2896: 2 aspect 0x00000001046dfbe0
_ZN6dealii11SolverGMRESINS_16TrilinosWrappers3MPI6VectorEE5solveINS1_12SparseMatrixENS1_15PreconditionILUEEEvRKT_RS3_RKS3_RKT0_
#1 3 aspect 0x00000001046dee28
_ZN6aspect9SimulatorILi2EE15solve_advectionERKNS1_14AdvectionFieldE +
1916: 3 aspect 0x00000001046dee28
_ZN6aspect9SimulatorILi2EE15solve_advectionERKNS1_14AdvectionFieldE
#2 4 aspect 0x00000001046e3e10
_ZN6aspect9SimulatorILi2EE30assemble_and_solve_temperatureEbPd + 292:
4 aspect 0x00000001046e3e10
_ZN6aspect9SimulatorILi2EE30assemble_and_solve_temperatureEbPd
#3 5 aspect 0x00000001046e4afc
_ZN6aspect9SimulatorILi2EE36solve_single_advection_single_stokesEv +
32: 5 aspect 0x00000001046e4afc
_ZN6aspect9SimulatorILi2EE36solve_single_advection_single_stokesEv
#4 6 aspect 0x00000001045d8348
_ZN6aspect9SimulatorILi2EE14solve_timestepEv + 184: 6 aspect
0x00000001045d8348 _ZN6aspect9SimulatorILi2EE14solve_timestepEv
#5 7 aspect 0x00000001045d7440
_ZN6aspect9SimulatorILi2EE3runEv + 900: 7 aspect
0x00000001045d7440 _ZN6aspect9SimulatorILi2EE3runEv
#6 8 aspect 0x000000010520126c
_Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
+ 480: 8 aspect 0x000000010520126c
_Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
#7 9 aspect 0x0000000105200ac0 main +
968: 9 aspect 0x0000000105200ac0 main
#8 10 dyld 0x000000018fd4fe50 start +
2544: 10 dyld 0x000000018fd4fe50 start
--------------------------------------------------------
Stacktrace:
-----------
#0 2 aspect 0x000000010467b0f8 _ZN6aspect9Utilities20linear_solver_failedERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_RKNS1_6vectorIN6dealii13SolverControlENS5_ISC_EEEERKSt9exceptionRKP19ompi_communicator_tS9_ + 916: 2 aspect 0x000000010467b0f8 _ZN6aspect9Utilities20linear_solver_failedERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_RKNS1_6vectorIN6dealii13SolverControlENS5_ISC_EEEERKSt9exceptionRKP19ompi_communicator_tS9_
#1 3 aspect 0x00000001046def00 _ZN6aspect9SimulatorILi2EE15solve_advectionERKNS1_14AdvectionFieldE + 2132: 3 aspect 0x00000001046def00 _ZN6aspect9SimulatorILi2EE15solve_advectionERKNS1_14AdvectionFieldE
#2 4 aspect 0x00000001046e3e10 _ZN6aspect9SimulatorILi2EE30assemble_and_solve_temperatureEbPd + 292: 4 aspect 0x00000001046e3e10 _ZN6aspect9SimulatorILi2EE30assemble_and_solve_temperatureEbPd
#3 5 aspect 0x00000001046e4afc _ZN6aspect9SimulatorILi2EE36solve_single_advection_single_stokesEv + 32: 5 aspect 0x00000001046e4afc _ZN6aspect9SimulatorILi2EE36solve_single_advection_single_stokesEv
#4 6 aspect 0x00000001045d8348 _ZN6aspect9SimulatorILi2EE14solve_timestepEv + 184: 6 aspect 0x00000001045d8348 _ZN6aspect9SimulatorILi2EE14solve_timestepEv
#5 7 aspect 0x00000001045d7440 _ZN6aspect9SimulatorILi2EE3runEv + 900: 7 aspect 0x00000001045d7440 _ZN6aspect9SimulatorILi2EE3runEv
#6 8 aspect 0x000000010520126c _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb + 480: 8 aspect 0x000000010520126c _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
#7 9 aspect 0x0000000105200ac0 main + 968: 9 aspect 0x0000000105200ac0 main
#8 10 dyld 0x000000018fd4fe50 start + 2544: 10 dyld 0x000000018fd4fe50 start
--------------------------------------------------------
Aborting!
----------------------------------------------------
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 1.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
```

The problem would then be to adjust the max number of iterations. I looked in the manual and tried to insert fixes, but without any success. Do you know how I can solve this problem?

The problem would then be to adjust the max number of iterations. I looked in the manual and tried to insert fixes, but without any success. Do you know how I can solve this problem? This is the code:

```
# This setup is a copy of the
# benchmarks/onset-of-convection/convection-box-base.prm
# contributed by Max Rudolph, with the difference that
# parameter values are specified explicitly in the input file
# (rather than through an ipython notebook).
set Dimension = 2
set Output directory = convection
subsection Formulation
set Formulation = Boussinesq approximation
end
set Pressure normalization = surface
set Surface pressure = 0
set Use conduction timestep = true
subsection Geometry model
set Model name = box
subsection Box
set X extent = 9.424777e6 # pi * 3e6
set Y extent = 2.9e6
set X repetitions = 3
set Y repetitions = 1
end
end
subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,z
set Function constants = p=1, L=9.424777e6, H=3.0e6, pi=3.1415926536, k=2
set Function expression = 2500 * (1.0-z/H) - p*cos(k*pi*x/L)*sin(pi*z/H)
end
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 4000
set Left temperature = 0
set Right temperature = 0
set Top temperature = 273
end
end
subsection Boundary velocity model
set Tangential velocity boundary indicators = left, right, bottom, top
end
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10.0
end
end
subsection Compositional fields
set Number of fields = 2
end
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = if(y<341176, 2.9e6, 0) ; if(y>2558824, 2.9e6, 0)
end
end
subsection Material model
set Model name = simple
subsection Simple model
set Reference density = 3300
set Reference specific heat = 1250
set Reference temperature = 273
set Thermal conductivity = 4.0
set Thermal expansion coefficient = 3e-5
set Viscosity = 1e23
set Density differential for compositional field 1 = 1.5
end
end
subsection Mesh refinement
set Initial global refinement = 4
set Initial adaptive refinement = 0
set Time steps between mesh refinement = 0
end
subsection Postprocess
set List of postprocessors = velocity statistics, visualization, temperature statistics, composition statistics, heat flux statistics
subsection Visualization
set List of output variables = density
set Time steps between graphical output = 1
end
end
```

Another problem I find when I try to edit the files "global melt " and “steinberger”. I decided to modify these two files because, among the various cookboks, they seem to me to have more realistic data to carry out the simulation.

Here I stumble into a couple of problems: as for the file "global melt ", I cannot insert a hyperdense and distinct layer as I did in the file “onset on convection”. Whatever I do, when I try to modify the function, mpirun does not allow me to do the simulation.

As for Steinberger, I modified the simulation trying to make the 2D section of a globe, and the simulation runs perfectly.

However, here I have absolutely no idea how I should add a hyperdense and circular layer around the CMB. There is no such example in the manual, and I really don’t know what I should write in the script to handle this, also because there is already a layer of material that should represent the density gradient of the planet.

Thank you in advance for your help.