Problems during the setup of a thermochemical layer at the CMB to study the origin and evolution of Superplumes (LLSVP)

[Francyrad]

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.

[Francyrad]

Here are the other snapshot that i could insert because i’m a new user
Internal structure of ultralow-velocity zones consistent with origin from a basal magma ocean. Pachhai et al. 2022

[Francyrad]

Globe on stenberger simulation

[bangerth]

Francesco:
let’s start with the very first error you show and work from there. The screenshot you show says

Line <24> of file <input string>: No entry with name <Phi periodic> was declared in the current subsection.

but at least the version of that file in the current repository does not actually have a line in which “Phi periodic” is set. The question then is: Where did you get that input file from, and does it match the version of ASPECT you have?

As for the later errors, when you say “Whatever I do, when I try to modify the function, mpirun does not allow me to do the simulation.” You will have to be more descriptive in what you did and what happened as a consequence.

Best
W.

[Francyrad]

Dear Bangheret and other readers
I’m so sorry for my late response. I did a correction in my question deleting the line 24, which is not present in the actual repository.
That error is present because i edited the file to add a periodic boundary, so i added that line:

  subsection Box
    set X extent = 9.424777e6 # pi * 3e6
    set Y extent = 2.9e6
    set X repetitions = 3
    set Y repetitions = 1
    set Phi periodic = true

As there is written in the page 166 of the manual. but i’ve been unsuccessful to make that work, maybe because it only works in a quarter shell and not in a box.
In truth this was not the principal problem that i encountered, if you could help me to add a periodic boundary to the box as i saw in some Citcom’s simulation it would be fantastic! Although the main problem of my simulation-box is that it crashed, as described in the post.

Regarding the later errors, i’m sorry to do not have been more specific. I attach the code that i edited to have 2D steinberger sphere that actually works:

set Dimension                              = 2
set End time                               = 3e8
set Adiabatic surface temperature          = 1613.0
set Output directory                       = output-steinberger


set Nonlinear solver scheme                = iterated Advection and Stokes
set Nonlinear solver tolerance             = 1e-3

# spherical shell

subsection Geometry model
  set Model name = spherical shell
  subsection Spherical shell
    set Inner radius  = 3480000
    set Outer radius  = 6370000
  end
end

  set Tangential velocity boundary indicators = top, bottom
end

subsection Nullspace removal
  set Remove nullspace = net rotation
end

# We use the gravity model from PREM.
subsection Gravity model
  set Model name = ascii data
end

# Set temperature
  set Fixed temperature boundary indicators = top, bottom
  set List of model names = spherical constant

  subsection Spherical constant
    set Inner temperature = 3773
    set Outer temperature = 273
  end
end

# The initial temperature model consists of an adiabatic profile,
# thermal boundary layers at the surface and the core-mantle boundary,
# and a small harmonic perturbation to initiate convection.
subsection Initial temperature model
  set List of model names = adiabatic, function

  subsection Function
    set Coordinate system   = spherical
    set Variable names      = r, phi
    set Function constants  = r0 = 3481000, r1 = 6371000
    set Function expression = 30 * sin((r-r0)/(r1-r0)*pi)*sin(16*phi) + 20 * sin(2*(r-r0)/(r1-r0)*pi)*sin(12*phi)
  end

  subsection Adiabatic
    set Age top boundary layer = 1e8
    set Age bottom boundary layer = 1e8
  end
end

subsection Material model
  set Model name = Steinberger
    

  subsection Steinberger model

    set Data directory                                = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
    set Lateral viscosity file name                   = temp-viscosity-prefactor.txt
    set Material file names                           = pyr-ringwood88.txt
    set Radial viscosity file name                    = radial-visc-simple.txt
    
    set Maximum lateral viscosity variation           = 1e3
    set Maximum viscosity                             = 5e23
    set Minimum viscosity                             = 1e20

    set Use lateral average temperature for viscosity = false

    set Latent heat                                   = false
  end
end


subsection Heating model
  set List of model names = adiabatic heating, shear heating
end

subsection Formulation
  set Mass conservation = projected density field
  set Temperature equation = real density
end

subsection Compositional fields
  set Number of fields                     = 1
  set Names of fields                     = density_field
  set Types of fields                     = density
  set Compositional field methods        = prescribed field
end

subsection Mesh refinement
  set Initial adaptive refinement        = 2
  set Initial global refinement          = 5
  set Refinement fraction                = 0.95
  set Coarsening fraction                = 0.05
  set Strategy                           = temperature, boundary, minimum refinement function
  set Time steps between mesh refinement = 500

  subsection Boundary
    set Boundary refinement indicators = top, bottom
  end

  subsection Minimum refinement function
    set Coordinate system   = depth
    set Variable names      = depth, phi
    set Function expression = if(depth<700000, 6, 5)
  end
end


subsection Postprocess
  set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics

  subsection Visualization
    set Output format                 = vtu
    set List of output variables      = material properties, nonadiabatic temperature, named additional outputs
    set Time between graphical output = 1e6
  end
end


subsection Solver parameters
  subsection Stokes solver parameters
    set GMRES solver restart length = 200
    set Number of cheap Stokes solver steps = 500
  end
end

subsection Checkpointing
  set Steps between checkpoint = 10
end

The difficulties that now i encounter is that i don’t know how to add a circular hyperdense compositional field layer around the CMB like i did in the upper simulation (where, instead of a circle, i added a rectangular layer because i was simulating a box).
In the user manual (pg 92) there is written that, to add an active compositional field, i’ve to add this to my script:

# This is the new part: We declare that there will
# be two compositional fields that will be
# advected along. Their initial conditions are given by
# a function that is one for the lowermost 0.2 height
# units of the domain and zero otherwise in the first case,
# and one in the top most 0.2 height units in the latter.
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<0.2, 1, 0) ; if(y>0.8, 1, 0)
#obviously the number changes if i use a box with real kilometers.
end
end

Anyway, in the steinberger simulation that i posted there are already compositional fields lines inside:

subsection Compositional fields
  set Number of fields                     = 1
  set Names of fields                     = density_field
  set Types of fields                     = density
  set Compositional field methods        = prescribed field
end

So i don’t get where i should put the mine. Even if i try to do that:

subsection Compositional fields
  set Number of fields                     = 2
  set Names of fields                     = density_field
  set Types of fields                     = density
  set Compositional field methods        = prescribed field
end

subsection Material model
  set Model name = Steinberger
  
    

  subsection Steinberger model

    set Data directory                                = $ASPECT_SOURCE_DIR/data/material-model/steinberger/
    set Lateral viscosity file name                   = temp-viscosity-prefactor.txt
    set Material file names                           = pyr-ringwood88.txt
    set Radial viscosity file name                    = radial-visc-simple.txt
    
    set Maximum lateral viscosity variation           = 1e3
    set Maximum viscosity                             = 5e23
    set Minimum viscosity                             = 1e20

    set Use lateral average temperature for viscosity = false

    set Latent heat                                   = false
    set Density differential for compositional field 2 = 1.5
  end
end

subsection Initial composition model
  set Model name = function
  
  subsection Function
    set Variable names = x,y
    set Function expression = if((sqrt((x-3600000)^2+(y-3600000)^2)>807400) , 0 , 6370000 )
  end
end

The simulation doesn’t even start

mpirun -np 8 aspect steinberger_edit.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
-----------------------------------------------------------------------------

    Line <21> of file <input string>: No entry with name <Tangential velocity boundary indicators> was declared in the current subsection.



----------------------------------------------------
Exception 'dealii::ExcMessage ("Invalid input parameter file.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <462> of file </Users/francyrad/Documents/aspect_new/source/main.cc> in function
    void parse_parameters(const std::string &, dealii::ParameterHandler &)
The violated condition was: 
    false
Additional information: 
    Invalid input parameter file.

Stacktrace:
-----------
#0  2   aspect                              0x00000001030704e4 _Z16parse_parametersRKNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEERN6dealii16ParameterHandlerE + 472: 2   aspect                              0x00000001030704e4 _Z16parse_parametersRKNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEERN6dealii16ParameterHandlerE 
#1  3   aspect                              0x000000010307113c _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb + 176: 3   aspect                              0x000000010307113c _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb 
#2  4   aspect                              0x0000000103070ac0 main + 968: 4   aspect                              0x0000000103070ac0 main 
#3  5   dyld                                0x000000018e403e50 start + 2544: 5   dyld                                0x000000018e403e50 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.

And also it is not clear what function should i add (and where), because there is not any example in the manual regarding this particular step, so i don’t know if my function is correct.

Another problem that i encountered is to add a basal hyperdense compositional field to the “global melt” script. Here i suppose that the solution to the problem would be easier. Anyway, i don’t know where to add the compositional field because there are already 2:

subsection Compositional fields
  set Number of fields = 2
  set Names of fields = porosity, peridotite
end

So i tried to do this and add a compositional field like in the classic “onset of convection” simulation. I release here all the script:

set Dimension                              = 2
set Adiabatic surface temperature          = 1600
set Maximum time step                      = 1e6
set Output directory                       = output-global_melt
set Use years in output instead of seconds = true


set Nonlinear solver scheme                = iterated Advection and Stokes
set Max nonlinear iterations               = 20
set Nonlinear solver tolerance             = 1e-5


set Use operator splitting                     = true
subsection Solver parameters
  subsection Operator splitting parameters
    set Reaction time step                     = 1e4
    set Reaction time steps per advection step = 10
  end
  subsection Stokes solver parameters
    set Linear solver tolerance = 1e-8
    set Number of cheap Stokes solver steps = 100
  end
end


set End time                               = 1.4e8


subsection Boundary velocity model
  set Tangential velocity boundary indicators = left, right, top, bottom
end

subsection Melt settings
  
  set Include melt transport                  = true
end

#I EDITED HERE
subsection Compositional fields
  set Number of fields = 3
  set Names of fields = porosity, peridotite, superplumes
end
#I EDITED HERE


subsection Discretization
  subsection Stabilization parameters
    set beta  = 0.5
    set cR    = 1
  end
end


subsection Initial temperature model
  set List of model names = adiabatic, function
  subsection Adiabatic
    set Age bottom boundary layer = 5e8
    set Age top boundary layer    = 3e8

    subsection Function
      set Function expression       = 0;0
    end
  end
  

  subsection Function
    set Function constants        = r=350000, amplitude=50
    set Function expression       = if((x-2900000)*(x-2900000)+y*y<r,amplitude,0)
  end
end


#I EDITED HERE
subsection Initial composition model
  set Model name = function
  subsection Function
    set Function expression = (y<341176, 2.9e6, 0) ; (y>2558824, 2.9e6, 0)
    set Variable names      = x,y
  end
end
#I EDITED HERE

subsection Boundary temperature model
  set Fixed temperature boundary indicators = top, bottom
  set List of model names = initial temperature

  subsection Initial temperature
    set Minimal temperature = 293
    set Maximal temperature = 3700
  end
end

# We again choose the initial composition as boundary condition
# for all compositional fields.
subsection Boundary composition model
  set List of model names = initial composition
end


subsection Boundary fluid pressure model
  set Plugin name = density
  subsection Density
    set Density formulation = fluid density
  end
end

subsection Geometry model
  set Model name = box

  subsection Box
    set X extent = 8700000
    set Y extent = 2900000
    set X repetitions = 3
    
  end

end


subsection Gravity model
  set Model name = vertical

  subsection Vertical
    set Magnitude = 9.81
  end
end


subsection Material model

  set Model name = melt global
  subsection Melt global

    
    set Thermal conductivity              = 4.7
    set Reference solid density           = 3400
    set Thermal expansion coefficient     = 2e-5
    set Reference shear viscosity         = 5e21
    set Thermal viscosity exponent        = 7
    set Reference temperature             = 1600
    set Solid compressibility             = 4.2e-12
    set Density differential for compositional field 3 = 1.5

    
    set Reference melt density            = 3000


    set Reference permeability            = 1e-8

    set Reference bulk viscosity          = 1e19


    set Exponential melt weakening factor = 10

    
    set Thermal bulk viscosity exponent   = 7

    
    set Melt compressibility              = 1.25e-11
    set Melt bulk modulus derivative      = 0.0


    set Reference melt viscosity          = 1


    set Depletion density change          = -200.0

    
    set Depletion solidus change          = 200

    
    set Melting time scale for operator splitting = 1e4
  end
end


subsection Mesh refinement
  set Coarsening fraction                      = 0.05
  set Refinement fraction                      = 0.8

  set Initial adaptive refinement              = 2
  set Initial global refinement                = 4
  set Strategy                                 = composition threshold, minimum refinement function
  set Time steps between mesh refinement       = 5

  
  subsection Minimum refinement function
    set Function expression = 4
  end

  
  subsection Composition threshold
    set Compositional field thresholds = 1e-8,1.0
  end
end

subsection Heating model
  set List of model names = adiabatic heating
end


subsection Postprocess

  set List of postprocessors = visualization, composition statistics, velocity statistics, temperature statistics, depth average


  subsection Visualization
    set List of output variables      = material properties, nonadiabatic temperature, strain rate, melt material properties

    subsection Material properties
      set List of material properties = density, viscosity, thermal expansivity
    end

    subsection Melt material properties
      set List of properties = fluid density, permeability, fluid viscosity, compaction viscosity
    end

    set Number of grouped files       = 0
    set Output format                 = vtu
    set Time between graphical output = 6e5
  end

  subsection Depth average
    set Number of zones = 12
    set Time between graphical output = 6e5
  end

end


subsection Checkpointing
  set Time between checkpoint = 1700
end

But my simulation doesn’t start

mpirun aspect global_melt.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 10 MPI processes
-----------------------------------------------------------------------------

    Line <144> of file <input string>: No entry with name <Density differential for compositional field 3> was declared in the current subsection.



----------------------------------------------------
Exception 'dealii::ExcMessage ("Invalid input parameter file.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <462> of file </Users/francyrad/Documents/aspect_new/source/main.cc> in function
    void parse_parameters(const std::string &, dealii::ParameterHandler &)
The violated condition was: 
    false
Additional information: 
    Invalid input parameter file.

Stacktrace:
-----------
#0  2   aspect                              0x00000001033084e4 _Z16parse_parametersRKNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEERN6dealii16ParameterHandlerE + 472: 2   aspect                              0x00000001033084e4 _Z16parse_parametersRKNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEERN6dealii16ParameterHandlerE 
#1  3   aspect                              0x000000010330913c _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb + 176: 3   aspect                              0x000000010330913c _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb 
#2  4   aspect                              0x0000000103308ac0 main + 968: 4   aspect                              0x0000000103308ac0 main 
#3  5   dyld                                0x000000018e403e50 start + 2544: 5   dyld                                0x000000018e403e50 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.
--------------------------------------------------------------------------

I’m so sorry for my late and long response. But I have no idea how to solve these problem. The solution is surely easy and the problem is surely given by the fact that I’m new to the scripts, but i have no idea how to solve them and i’d like to have some help, so i can learn.

Thank you for your availability.

Francesco