Dear all . I want to make a subduction model with the bottom boundary open . But “boundary traction” can’t achieve this . Let me know if anyone has any ideas. Thanks .

Hi @Xiao,

Using the traction boundary condition is in fact the only mechanism to allow dynamically driven flow through the bottom boundary.

To achieve this, you will need to specify the lithostatic pressure as the bottom traction value.

Can you let us know what when wrong in your simulation using a bottom traction boundary condition? Images and relevant parts (or all) of the prm file will help.

Cheers,

John

Hi John,

Thank you for your suggestion regarding the traction boundary condition to allow dynamically driven flow through the bottom boundary. I understand that specifying the lithostatic pressure as the bottom traction value is crucial for this.

However, I’ve encountered some issues with this approach. The main challenge is that the lithostatic pressure can only be specified at a single point. This makes it difficult to accurately define the bottom boundary traction. Additionally, as the subduction process progresses, the bottom boundary traction can become unrealistic. The discrepancy between the initial lithostatic pressure and the evolving conditions results in significant deviations, making the boundary condition less reliable.

Best regards,

Xiao

Xiao:

You are of course right: Prescribing “lithostatic” pressure along the entire bottom boundary is unrealistic. If you simulated the entire earth, the pressure along a fixed depth would have lateral variations. At the core of your problem is that you are not modeling the entire earth, but only a part that you want to truncate at a certain depth. Mathematically, you have to prescribe some boundary condition there, and the question then is: What is the right boundary condition in your situation? The problem is that there is no answer to this: To the best of my knowledge, nobody has come up with a better idea than prescribing a constant pressure at the bottom, and so nobody has implemented a better idea in ASPECT either. Among the reasons why this is a difficult question is because in reality, what that pressure *should* be does not just depend on the state of the material *above* the interface, but also *below* the interface (think about the case of a rising plume that has not yet reached the depth at which you truncate your domain: it will surely affect the pressure at the bottom of your domain). In other words, the “correct” boundary condition needs to make assumptions about what happens outside the domain you are simulating, but since we are not simulating the “outside”, we cannot know what “correct” is supposed to mean.

Does that make sense?

Best

W.

Hi , Wolf :

Thank you for your explanation . I understand !

Best

Xiao

Dear all ,

I try to make a subduction model with bottom boundary open . This is my setup:

**subsection Boundary velocity model**

** set Prescribed velocity boundary indicators = left :function , right x:function **

** subsection Function**

** set Variable names = x,y**

** set Function constants = cm=0.01 , yr = 1 , VelContinental = 3 , Ay = 300e3 , HFore = 100e3 , XLong = 1200e3 , VelOcean = 3 , VelConOut = 0.942 , VelOceOut = 0.942 , Htrans = 100e3 , HOcean = 90e3**

** set Function expression = if( y>=(Ay - HOcean) & x<1e3 , VelOcean*cm/yr , **

** if( y>=(Ay-HFore) & x<1e3 , -VelContinental*cm/yr , 0 ) );\ **

** 0**

** end**

**end**

**####### Boundary traction model #######**

**####### All the bottom are open #######**

**subsection Boundary traction model**

** set Prescribed traction boundary indicators = bottom :initial lithostatic pressure **

** subsection Initial lithostatic pressure**

** set Number of integration points = 1800**

** set Representative point = 100000, 300000**

** end**

**end**

**####### Boundary composition model #######**

**subsection Boundary composition model**

** set List of model names = function **

** set Fixed composition boundary indicators = top , left , right ,bottom**

** #set Allow fixed composition on outflow boundaries = true**

** subsection Function**

** set Coordinate system = cartesian**

** set Variable names = x,y,t**

** set Function constants = tsed=2e6 # duration of sedimentation switch time intervall **

** set Function expression =**

** if( (y<220e3 & x<=1e3) | (y<200e3 & x>1e3 ) , 1 , 0 );**

** if( (y>=275e3 && t==0 & x > 500e3 & x <= 1199e3) | (y>=275e3 & x>1199e3) , 1,0); **

** if( (y<275e3 & y>=265e3 & x>1199e3 ) , 1 , 0 );**

** if( (y<265e3 & y>=200e3 & x>1199e3 ), 1 , 0);**

** if( (y>=292e3 && t==0 & x >=1e3 & x <= 500e3) | ( x<1e3 & y>=292e3 ) , 1,0); **

** if( (y<292e3 & y>=220e3 & x<10e3 ) , 1 , 0);**

** if( y>292e3 && t>0 &x>=1e3 & x<=1199e3 ,1,0); **

** if( y>292e3 && t>0 &x>=1e3 & x<=1199e3 , t/1e6, 0) ;**

** 0;**

** 0;**

** 0**

** end**

**end**

**####### Boundary temperature model #######**

**subsection Boundary temperature model**

** set Fixed temperature boundary indicators = top , bottom , left , right**

** set List of model names = initial temperature**

**end**

**set Dimension = 2**

**set Start time = 0**

**set End time = 2e6**

**set CFL number = 0.5**

**set Use years in output instead of seconds = true**

**set Maximum time step = 5e4**

**set Pressure normalization = no**

**set Surface pressure = 0**

**set Output directory = subduction**

**set Timing output frequency = 10**

**# Nonlinear solver **

**set Nonlinear solver scheme = single Advection, iterated Stokes**

**set Max nonlinear iterations = 100**

**set Nonlinear solver tolerance = 1e-3**

**set Max nonlinear iterations in pre-refinement = 5 **

**# Linear solver parameters**

**subsection Solver parameters**

** subsection Stokes solver parameters**

** set Number of cheap Stokes solver steps = 200**

** end**

**end**

But the results don’t converge .

**Currently Loaded Modulefiles:**

*** 1) compiler/gcc/9.3.0***

*** 2) mpi/openmpi/openmpi-4.1.5-gcc9.3.0***

*** 3) compiler/cmake/3.23.3***

**-----------------------------------------------------------------------------**

**– This is ASPECT, the Advanced Solver for Problems in Earth’s ConvecTion.**

**– . version 2.3.0-pre (fault_analysis, ead281c)**

**– . using deal.II 9.2.0**

**– . with 32 bit indices and vectorization level 1 (128 bits)**

**– . using Trilinos 12.17.0**

**– . using p4est 2.8.5**

**– . running in OPTIMIZED mode**

**– . running with 64 MPI processes**

**-----------------------------------------------------------------------------**

*** Loading Ascii data initial file /work/home/hanx23/aspect/aspect_Fastscape/aspect_2.3_rift/build/workspace/test/down_open/Temperature.txt.***

*** Loading Ascii data initial file /work/home/hanx23/aspect/aspect_Fastscape/aspect_2.3_rift/build/workspace/test/down_open/Composition.txt.***

**-----------------------------------------------------------------------------**

**– For information on how to cite ASPECT, see:**

**– The ASPECT mantle convection code: How to cite?**

**-----------------------------------------------------------------------------**

**Number of active cells: 1,024 (on 6 levels)**

**Number of degrees of freedom: 60,239 (8,450+1,089+4,225+4,225+4,225+4,225+4,225+4,225+4,225+4,225+4,225+4,225+4,225+4,225)**

**Number of mesh deformation degrees of freedom: 2178**

****** Timestep 0: t=0 years, dt=0 years***

*** Solving mesh velocity system… 0 iterations.***

*** Solving temperature system… 0 iterations.***

*** Solving Mantle system … 0 iterations.***

*** Solving ContinentalUpperCrust system … 0 iterations.***

*** Solving ContinentalLowerCrust system … 0 iterations.***

*** Solving ContinentalMantle system … 0 iterations.***

*** Solving OceanicCrust system … 0 iterations.***

*** Solving OceanicMantle system … 0 iterations.***

*** Skipping sediment composition solve because RHS is zero.***

*** Skipping sediment_age composition solve because RHS is zero.***

*** Skipping noninitial_plastic_strain composition solve because RHS is zero.***

*** Skipping plastic_strain composition solve because RHS is zero.***

*** Skipping viscous_strain composition solve because RHS is zero.***

*** Rebuilding Stokes preconditioner…***

*** Solving Stokes system… 200+30 iterations.***

*** Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1***

*** Rebuilding Stokes preconditioner…***

*** Solving Stokes system… ---------------------------------------------------------***

**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 (“The " + solver_name + " did not converge. It reported the following error:\n\n” + exc.what() + "\n The required residual for convergence is: " + std::to_string(solver_controls.front().tolerance()) + “.\n See " + output_filename + " for convergence history.”)’ on rank 0 on processing: ***

**--------------------------------------------------------**

**An error occurred in line <398> of file </work/home/hanx23/aspect/aspect_Fastscape/aspect_2.3_rift/source/simulator/solver.cc> in function**

*** void aspect::{anonymous}::linear_solver_failed(const string&, const string&, const std::vectordealii::SolverControl&, const std::exception&)***

**The violated condition was: ***
*** false**

**Additional information: *****

*** The iterative Stokes solver did not converge. It reported the following error:*** The iterative Stokes solver did not converge. It reported the following error:

**--------------------------------------------------------**

**An error occurred in line <356> of file </work/home/hanx23/aspect/aspect_Fastscape/aspect_2.3_rift/source/simulator/solver.cc> in function**

*** void aspect::internal::BlockSchurPreconditioner<PreconditionerA, PreconditionerMp>::vmult(aspect::LinearAlgebra::BlockVector&, const BlockVector&) const [with PreconditionerA = dealii::TrilinosWrappers::PreconditionAMG; PreconditionerMp = dealii::TrilinosWrappers::PreconditionBase; aspect::LinearAlgebra::BlockVector = dealii::TrilinosWrappers::MPI::BlockVector]***

**The violated condition was: ***
*** false**

**Additional information: *****

*** The iterative (top left) solver in BlockSchurPreconditioner::vmult did not converge to a tolerance of -nan. It reported the following error:*** The iterative (top left) solver in BlockSchurPreconditioner::vmult did not converge to a tolerance of -nan. It reported the following error:

**--------------------------------------------------------**

**An error occurred in line <497> of file </work/share/ac2pv79iif/soft/dealii-9.2.0/source/lac/trilinos_solver.cc> in function**

*** void dealii::TrilinosWrappers::SolverBase::do_solve(const Preconditioner&) [with Preconditioner = dealii::TrilinosWrappers::PreconditionBase]***

**The violated condition was: ***
*** false**

**Additional information: *****

*** AztecOO::Iterate error code -2: numerical breakdown*** AztecOO::Iterate error code -2: numerical breakdown

**--------------------------------------------------------****--------------------------------------------------------**

*** The required residual for convergence is: -nan.***

*** See subduction/solver_history.txt for convergence history.***

**--------------------------------------------------------**

**Aborting!**

**----------------------------------------------------**

**--------------------------------------------------------------------------**

**Primary job terminated normally, but 1 process returned**

**a non-zero exit code. Per user-direction, the job has been aborted.**

**--------------------------------------------------------------------------**

**--------------------------------------------------------------------------**

**mpirun detected that one or more processes exited with non-zero status, thus causing**

**the job to be terminated. The first process to do so was:**

*** Process name: [[47071,1],0]***

*** Exit code: 1***

Has anybody got any good ideas or ways for them? Thanks!

Add : Top boundary is free surface

Dear Xiao,

I have not looked at your input file, but the error message tells you that something in your model becomes not a number. Can you rerun the model in DEBUG mode (rather than optimized mode)? This might give you additional error messages that could help you figure out what exactly is going wrong.

Best,

Juliane

original.prm (19.0 KB)

Dear all

This is my input file . I make the bottom boundary open , but it can not run . I don’t know why . After running the program, it never iterates and has no output . It stopped at the first step.

Hi @Xiao,

Thank you for sending over the parameter file. As Juliane noted, it would be really helpful to see the results of running this model in debug mode. However, to first order the result of the crash is a failure of the Stokes solver.

To diagnose all of this, I think you may want to take a step back and begin with simpler models run in debug mode and slowly introduce all the complexity in your PRM file. This would likely take a few days of work, but will go a long way to helping diagnose future issues that may arise.

Such a workflow may look this:

- Start with constant viscosity (say 1e22 Pa s), but have the current density and temperature variations throughout the model domain. Run this model with free slip boundary conditions everywhere. The end result should be flow patterns driven by the buoyancy variations.
- Modify step 1 to include the boundary conditions on the left and right boundaries.
- Modify step 2 to include the lithostatic pressure boundary conditions (see a note on this below)
- Modify step 3 to include the more complex viscosity structure (you may want to do this in sub-steps by systemically increasing the min-max viscosity contrast)
- Modify step 4 to include a free surface
- Modify step 5 to included coupling with Fastscape (current model setup).

Regarding the lithostatic pressure boundary condition, I think the bottom boundary should be deeper (600 or 1200 km depth) to make the simulation a bit more stable. At 300 km depth in combination with a free surface, I would not be surprised to see odd dynamic behavior arise.

Cheers,

John