Hi,
As someone who is not very familiar with ASPECT software, I aim to create a visco plastic model in a 2-dimensional and 2-layer model where I can see buckling due to compression with newtonian behavior depending on the geothermal gradient (power law) and I want to evaluate the behavior of the first thin layer. Is there a starting point or cook book you can recommend for this?
Theodor:
I don’t recall whether there is a particularly good example to start with, but I am pretty sure you can piece together what you want to do by looking at the cookbooks that explain how to prescribe extensional or compressional motion at the edges of a model, and the cookbooks that show how to select visco-plastic models such as the continental extension model.
I will add that the term “buckling” describes a phenomenon that is specific to elastic models. Buckling cannot happen in a (viscous) fluid. You can have dynamic topography in a fluid, however, and perhaps that is what you mean.
Best
W.
Dear @bangerth,
Thank you very much for your message and suggestion. Taking this into consideration, I decided to start with the model below.
cookbooks/continental_extension/continental_extension.prm
However, MPI gives an error for the ASPECT 2.6-pre version after the first timestep (attached) in the ASPECT 2.4 version it works.
log.txt (55.2 KB)
In addition to this, I have a problem with the following parameter file
compres4buckle_test.prm (13.8 KB)
Despite adjusting the boundary velocity model and mesh deformation, for each time step, there is an increase in y (increase in the overall topography). This increase is quite obvious from the edges of the box model. According your previous explanation, I would like to obtain the intermittent folding of the upper crust layer on this model.
Thedor:
Let’s keep it to one issue per thread – feel free to open a different thread for the boundary compression. In an ideal world, you might also want to show us a picture that illustrates what you get. If you only give us the input file, you require us to run and visualize it ourselves, when you’ve already done that step.
As for the error you get with the continental extension model: I can’t reproduce this. The only difference between your run and mine is that you are using deal.II 9.4.2, and I am using the current development sources. Could you try again with deal.II 9.5?
Best
W.
Dear @bangerth thanks for your reply.
Regarding the error, when I download and run the deal.II virtual machine image, deal version 9.4.0,
and aspect 2.5 work fine. In my local installation, I was able to run the parameter file in deal.II version 9.5.1 and aspect 2.5, but I get the same error when I use deal.II version 9.4.2 or deal.II version 9.5.1 and aspect 2.6.
Regarding the parameter file, I was having this problem because I did not specify influx due to my mistake. But I would be very grateful if you could help me with another issue. If you can look at the attached GIF animation file, I can get 2 bucklings in the upper crust during the compression process. There are two parts of the parameter file that I don’t understand, the first one is its relationship with the weak seed. By increasing the amount of this, I want to get more than one successive buckle fold with the relationship of viscosity and thermal boundary condition. My other question is about the hillslope diffusion effect, I wonder about the contribution of the softening effect here.
compres4buckle_test.prm (12.6 KB)
Best
W.
Hi Wolfgang,
Thanks for your reply. Unfortunately, since I don’t have enough experience, I thought it would be more useful to reduce the model to 2 layers I removed the lower crust and it seemed simpler to follow the buiter_et_al_2008_jgr benchmark file in the parameter file. This seems to be more useful for the buckle folding I want to achieve.
# Global parameters
set Dimension = 2
set Start time = 0
set End time = 150e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-3
set Max nonlinear iterations = 120
set CFL number = 0.7
set Output directory = out_test
set Pressure normalization = no
set Timing output frequency = 1
# Model geometry (200x60 km)
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 100 # for test
set Y repetitions = 20 # for test
set X extent = 200e3
set Y extent = 60e3
end
end
# No mesh refinement, but the coarse mesh is 200x100
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end
# The nonlinear solver typically converges more quickly
# when no cheap Stokes solver steps are used for
# problems with Drucker-Prager plasticity.
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 0
end
end
# Free surface boundary classifications.
# Advecting the free surface vertically rather than
# in the surface normal direction can result in a
# more stable mesh when the deformation is large
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Surface velocity projection = vertical
end
end
# The top boundary is open (zero traction), which allows the sticky air to
# flow freely through it as topography develops along the wedge. Additional
# testing revealed that using a true free surface leads to significant mesh
# distortion and associated numerical instabilities.
subsection Boundary traction model
set Prescribed traction boundary indicators = top: zero traction
end
# Velocity on boundaries characterized by functions
# Total extension rate is 1 cm/yr (0.5 cm/yr on each side)
subsection Boundary velocity model
set Prescribed velocity boundary indicators = left x: function, right x:function
set Tangential velocity boundary indicators = bottom
subsection Function
set Variable names = x,y
set Function constants = m=-0.05, year=1
set Function expression = if (x<100e3 , -1*m/year, 1*m/year); m*2*100e3/200e3 # Include influx
end
end
# Number and name of compositional fields
# The field plastic_strain is used for tracking the plastic finite strain invariant
# upper: brittle upper crust; seed: 'weaker' brittle region
# lower: viscous lower crust
subsection Compositional fields
set Number of fields = 4
set Names of fields = plastic_strain, upper, seed, lower
end
# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0; \
if((y>=50e3 && x<=98.0e3) || (y>=50e3 && x>=102.0e3) || \
(y>=51e3 && x>=96.0e3 && x<=104.0e3), 1, 0); \
if(y>=50e3 && y<51e3 && x>95.0e3 && x<100.0e3, 1, 0); \
if(y<50e3, 1, 0);
end
end
# Composition boundary conditions
subsection Boundary composition model
set List of model names = initial composition
end
# Use discontinuous composition bound preserving limiter
subsection Discretization
set Composition polynomial degree = 2
set Stokes velocity polynomial degree = 2
set Temperature polynomial degree = 2
set Use discontinuous composition discretization = true
subsection Stabilization parameters
set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions
set Global composition maximum = 100.0, 1.0, 1.0, 1.0
set Global composition minimum = 0.0, 0.0, 0.0, 1.0
end
end
# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1400
set Top temperature = 550
end
end
# Temperature initial conditions (linear)
subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 550 + (1400 - 550) * y / 60e3
end
end
# Material model
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Reference temperature = 273
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-16
set Minimum viscosity = 1e18
set Maximum viscosity = 1e26
set Thermal diffusivities = 1.8e-6
set Heat capacities = 750.
set Densities = 2800
set Thermal expansivities = 0.
set Viscosity averaging scheme = harmonic
set Viscous flow law = dislocation
set Prefactors for dislocation creep = 5.e-25, 5.e-25, 5.e-26, 5.e-25, 5.e-11
set Stress exponents for dislocation creep = 1.0
set Activation energies for dislocation creep = 0.
set Activation volumes for dislocation creep = 0.
set Angles of internal friction = 30., 30., 30., 4., 30.
set Cohesions = 20.e6, 20.e6, 20.e6, 2.e6, 1.e11
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.5, 0.5, 0.5, 0.5, 0.5
set End plasticity strain weakening intervals = 1.5, 1.5, 1.5, 1.5, 1.5
set Cohesion strain weakening factors = 0.1, 0.1, 0.1, 1.0, 1.0
set Friction strain weakening factors = 0.13333, 0.1333, 0.1333, 1.0, 1.0
end
end
# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 8.87
end
end
# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator, partition, surface stress, shear stress, surface dynamic topography
set Time between graphical output = 1e5
set Interpolate output = false
end
end
I changed the model geometry to a thicker lower layer and a thinner upper layer. Regarding the seeds, I placed seeds to get decoupling from the viscosity differences between the top and bottom layers and organized them as in the parameter file.
log.txt (537.0 KB)
-- Total wallclock time elapsed including restarts: 9374s
*** Timestep 56: t=178130 years, dt=1061.51 years
Solving mesh displacement system... 9 iterations.
Solving temperature system... 10 iterations.
Solving plastic_strain system ... 4 iterations.
Solving upper system ... 4 iterations.
Skipping seed composition solve because RHS is zero.
Solving lower system ... 4 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 0+--------------------------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 0 on node Computer exited on signal 8 (Floating point exception).
The process stops after the above timestep. I would appreciate it if you could help with this issue.
@Thedor Is there an error file in the directory you ran this in? Something like log.err or similar?
If not, can you run the model with only one processor and see whether you get a better error message?
Best
W.
Hi all,
Apologies for jumping into the conversation midway through.
@Thedor - Regarding the first issue, can you update to the most recent ASPECT development version? The issue should have been resolved with this PR, and it was related to how we were using deal.II functionality to access the “old” solution in the strain dependent rheology.
The lower resolution tests did not catch the issue, and this is a good reminder to update those tests and think about doing an update release since the bug is in ASPECT 2.5.0 
If updating to a more recent version of ASPECT in combination with deal.II 9.5 does not resolve the issue, can you report back.
Regarding the buckling models, the buiter et al. 2008 models are indeed the same place I would start as well or similar recently published geodynamic models. In the end, you may need to combine elements from the continental extension cookbook (continental geotherm, temperature and pressure-dependent viscous flow) into the buiter et al. cookbook.
The process stops after the above timestep. I would appreciate it if you could help with this issue.
Adding to @bangerth suggestions above, it looks like the Stokes solver is not converging. You ran in debug mode, so hopefully there is an additional error file with more detailed information.
Regardless, here are some suggestions for improving the Stokes solver behavior:
Try reducing the min/max viscosity contrast. Currently, it is quite high (1e18 - 1e26 Pa s). Maybe reduce it to six orders of magnitude (1e19 - 1e25) and see if that enables convergence for the same model?
Use a stricter CFL value - I would systematically try 0.5, 0.2, and 0.1. For this class of problems, I often end up using something like 0.1 or 0.2.
Use a stricter nonlinear solver tolerance (1e-4 or 1-e5) in combination with additional adjustments to the Stokes solver section. The class of model you are doing is highly nonlinear, so quite a bit of testing at the onset might be needed here. Here are some suggestions for the Stokes solver section:
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 500
set Linear solver tolerance = 1e-7
set GMRES solver restart length = 100
set Use full A block as preconditioner = true
end
end
If the above settings don’t result in convergence for the current or adjusted model setup, try using set GMRES solver restart length = 150
It is on my list to update this and other similar cookbooks with many of the above suggestions.
Let’s keep the discussion going here and with enough adjustments and testing we will arrive at a stable base model for your investigation 
Cheers,
John
Hi Wolfgang,
I ran deal 9.5.1 and ASPECT 2.5.0 using a single core with debug mode. The error message is as follows,
-- Total wallclock time elapsed including restarts: 76257s
*** Timestep 55: t=176548 years, dt=1090.52 years
Solving mesh displacement system... 8 iterations.
Solving temperature system... 6 iterations.
Solving plastic_strain system ... 3 iterations.
Solving upper system ... 3 iterations.
Skipping seed composition solve because RHS is zero.
Solving lower system ... 3 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 0+6 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.000496133
Postprocessing:
RMS, max velocity: 0.0341 m/year, 0.497 m/year
Temperature min/avg/max: 550 K, 975 K, 1401 K
Pressure at top/bottom of domain: 1.519e+06 Pa, 1.617e+09 Pa
Computing dynamic topography
+----------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 7.63e+04s | |
| | | |
| Section | no. calls | wall time | % of total |
+----------------------------------+-----------+------------+------------+
| Assemble Stokes system | 2379 | 3.05e+03s | 4% |
| Assemble composition system | 224 | 1.19e+03s | 1.6% |
| Assemble temperature system | 56 | 103s | 0.13% |
| Build Stokes preconditioner | 2379 | 1.39e+03s | 1.8% |
| Build composition preconditioner | 167 | 7.69s | 0% |
| Build temperature preconditioner | 56 | 0.492s | 0% |
| Initialization | 1 | 0.753s | 0% |
| Mesh deformation | 56 | 16.6s | 0% |
| Mesh deformation initialize | 1 | 0.268s | 0% |
| Postprocessing | 56 | 24s | 0% |
| Setup dof systems | 1 | 0.437s | 0% |
| Setup initial conditions | 1 | 1.54s | 0% |
| Setup matrices | 1 | 8.23s | 0% |
| Solve Stokes system | 2379 | 7.05e+04s | 92% |
| Solve composition system | 167 | 24.6s | 0% |
| Solve temperature system | 56 | 0.219s | 0% |
+----------------------------------+-----------+------------+------------+
-- Total wallclock time elapsed including restarts: 76297s
*** Timestep 56: t=177646 years, dt=1097.71 years
Solving mesh displacement system... 8 iterations.
Solving temperature system... 6 iterations.
Solving plastic_strain system ... 3 iterations.
Solving upper system ... 3 iterations.
Skipping seed composition solve because RHS is zero.
Solving lower system ... 3 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 0+---------------------------------------------------------
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 <2943> of file </Users/okaragoz/aspect/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 Stokes solver in Simulator::solve_stokes did not
converge.
The initial residual was: 2.880029e+12
The final residual is: 2.880029e+12
The required residual for convergence is: 1.128283e+08
See new_test_output/solver_history.txt for the full convergence
history.
The solver reported the following error:
--------------------------------------------------------
An error occurred in line <2943> of file
</Users/okaragoz/aspect/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 (top left) solver in BlockSchurPreconditioner::vmult did
not converge.
The initial residual was: nan
The final residual is: 1.141692e+01
The required residual for convergence is: 1.559773e-02
The solver reported the following error:
--------------------------------------------------------
An error occurred in line <558> of file
</Users/okaragoz/deal.II-candi/tmp/unpack/deal.II-v9.5.1/source/lac/trilinos_solver.cc>
in function
void dealii::TrilinosWrappers::SolverBase::do_solve(const
Preconditioner &) [Preconditioner =
dealii::TrilinosWrappers::PreconditionBase]
The violated condition was:
false
Additional information:
Iterative method reported convergence failure in step 10000. The
residual in the last step was 11.4169.
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 libdeal_II.g.9.5.1.dylib 0x00000001180cd606
_ZN6dealii16TrilinosWrappers10SolverBase8do_solveINS0_16PreconditionBaseEEEvRKT_
+ 1430: 2 libdeal_II.g.9.5.1.dylib 0x00000001180cd606
_ZN6dealii16TrilinosWrappers10SolverBase8do_solveINS0_16PreconditionBaseEEEvRKT_
#1 3 libdeal_II.g.9.5.1.dylib 0x00000001180ce416
_ZN6dealii16TrilinosWrappers10SolverBase5solveERKNS0_12SparseMatrixERNS0_3MPI6VectorERKS6_RKNS0_16PreconditionBaseE
+ 118: 3 libdeal_II.g.9.5.1.dylib 0x00000001180ce416
_ZN6dealii16TrilinosWrappers10SolverBase5solveERKNS0_12SparseMatrixERNS0_3MPI6VectorERKS6_RKNS0_16PreconditionBaseE
#2 4 aspect 0x0000000103b38e5e
_ZNK6aspect8internal24BlockSchurPreconditionerIN6dealii16TrilinosWrappers15PreconditionAMGENS3_16PreconditionBaseEE5vmultERNS3_3MPI11BlockVectorERKS8_
+ 638: 4 aspect 0x0000000103b38e5e
_ZNK6aspect8internal24BlockSchurPreconditionerIN6dealii16TrilinosWrappers15PreconditionAMGENS3_16PreconditionBaseEE5vmultERNS3_3MPI11BlockVectorERKS8_
#3 5 aspect 0x0000000103a67de6
_ZN6dealii12SolverFGMRESINS_16TrilinosWrappers3MPI11BlockVectorEE5solveIN6aspect8internal11StokesBlockENS7_24BlockSchurPreconditionerINS1_15PreconditionAMGENS1_16PreconditionBaseEEEEEvRKT_RS3_RKS3_RKT0_
+ 854: 5 aspect 0x0000000103a67de6
_ZN6dealii12SolverFGMRESINS_16TrilinosWrappers3MPI11BlockVectorEE5solveIN6aspect8internal11StokesBlockENS7_24BlockSchurPreconditionerINS1_15PreconditionAMGENS1_16PreconditionBaseEEEEEvRKT_RS3_RKS3_RKT0_
#4 6 aspect 0x0000000103a670dc
_ZN6aspect9SimulatorILi2EE12solve_stokesEv + 4988: 6 aspect
0x0000000103a670dc _ZN6aspect9SimulatorILi2EE12solve_stokesEv
#5 7 aspect 0x0000000103a6b751
_ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd + 193: 7
aspect 0x0000000103a6b751
_ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd
#6 8 aspect 0x0000000103a6c249
_ZN6aspect9SimulatorILi2EE38solve_single_advection_iterated_stokesEv +
217: 8 aspect 0x0000000103a6c249
_ZN6aspect9SimulatorILi2EE38solve_single_advection_iterated_stokesEv
#7 9 aspect 0x00000001038ea28d
_ZN6aspect9SimulatorILi2EE14solve_timestepEv + 237: 9 aspect
0x00000001038ea28d _ZN6aspect9SimulatorILi2EE14solve_timestepEv
#8 10 aspect 0x00000001038e91ad
_ZN6aspect9SimulatorILi2EE3runEv + 1165: 10 aspect
0x00000001038e91ad _ZN6aspect9SimulatorILi2EE3runEv
#9 11 aspect 0x000000010473b0b9
_Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
+ 569: 11 aspect 0x000000010473b0b9
_Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
#10 12 aspect 0x000000010473a8d5 main +
1093: 12 aspect 0x000000010473a8d5 main
#11 13 libdyld.dylib 0x00007fff205eef3d start
+ 1: 13 libdyld.dylib 0x00007fff205eef3d start
#12 14 ??? 0x0000000000000002 0x0 +
2: 14 ??? 0x0000000000000002 0x0
--------------------------------------------------------
Stacktrace:
-----------
#0 2 aspect 0x00000001039f0e9d
_ZN6aspect9Utilities20linear_solver_failedERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_RKNS1_6vectorIN6dealii13SolverControlENS5_ISC_EEEERKSt9exceptionRKiS9_
+ 1085: 2 aspect 0x00000001039f0e9d
_ZN6aspect9Utilities20linear_solver_failedERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_RKNS1_6vectorIN6dealii13SolverControlENS5_ISC_EEEERKSt9exceptionRKiS9_
#1 3 aspect 0x0000000103b3911c
_ZNK6aspect8internal24BlockSchurPreconditionerIN6dealii16TrilinosWrappers15PreconditionAMGENS3_16PreconditionBaseEE5vmultERNS3_3MPI11BlockVectorERKS8_
+ 1340: 3 aspect 0x0000000103b3911c
_ZNK6aspect8internal24BlockSchurPreconditionerIN6dealii16TrilinosWrappers15PreconditionAMGENS3_16PreconditionBaseEE5vmultERNS3_3MPI11BlockVectorERKS8_
#2 4 aspect 0x0000000103a67de6
_ZN6dealii12SolverFGMRESINS_16TrilinosWrappers3MPI11BlockVectorEE5solveIN6aspect8internal11StokesBlockENS7_24BlockSchurPreconditionerINS1_15PreconditionAMGENS1_16PreconditionBaseEEEEEvRKT_RS3_RKS3_RKT0_
+ 854: 4 aspect 0x0000000103a67de6
_ZN6dealii12SolverFGMRESINS_16TrilinosWrappers3MPI11BlockVectorEE5solveIN6aspect8internal11StokesBlockENS7_24BlockSchurPreconditionerINS1_15PreconditionAMGENS1_16PreconditionBaseEEEEEvRKT_RS3_RKS3_RKT0_
#3 5 aspect 0x0000000103a670dc
_ZN6aspect9SimulatorILi2EE12solve_stokesEv + 4988: 5 aspect
0x0000000103a670dc _ZN6aspect9SimulatorILi2EE12solve_stokesEv
#4 6 aspect 0x0000000103a6b751
_ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd + 193: 6
aspect 0x0000000103a6b751
_ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd
#5 7 aspect 0x0000000103a6c249
_ZN6aspect9SimulatorILi2EE38solve_single_advection_iterated_stokesEv +
217: 7 aspect 0x0000000103a6c249
_ZN6aspect9SimulatorILi2EE38solve_single_advection_iterated_stokesEv
#6 8 aspect 0x00000001038ea28d
_ZN6aspect9SimulatorILi2EE14solve_timestepEv + 237: 8 aspect
0x00000001038ea28d _ZN6aspect9SimulatorILi2EE14solve_timestepEv
#7 9 aspect 0x00000001038e91ad
_ZN6aspect9SimulatorILi2EE3runEv + 1165: 9 aspect
0x00000001038e91ad _ZN6aspect9SimulatorILi2EE3runEv
#8 10 aspect 0x000000010473b0b9
_Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
+ 569: 10 aspect 0x000000010473b0b9
_Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
#9 11 aspect 0x000000010473a8d5 main +
1093: 11 aspect 0x000000010473a8d5 main
#10 12 libdyld.dylib 0x00007fff205eef3d start
+ 1: 12 libdyld.dylib 0x00007fff205eef3d start
#11 13 ??? 0x0000000000000002 0x0 +
2: 13 ??? 0x0000000000000002 0x0
--------------------------------------------------------
Stacktrace:
-----------
#0 2 aspect 0x00000001039f0e9d _ZN6aspect9Utilities20linear_solver_failedERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_RKNS1_6vectorIN6dealii13SolverControlENS5_ISC_EEEERKSt9exceptionRKiS9_ + 1085: 2 aspect 0x00000001039f0e9d _ZN6aspect9Utilities20linear_solver_failedERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_RKNS1_6vectorIN6dealii13SolverControlENS5_ISC_EEEERKSt9exceptionRKiS9_
#1 3 aspect 0x0000000103a67330 _ZN6aspect9SimulatorILi2EE12solve_stokesEv + 5584: 3 aspect 0x0000000103a67330 _ZN6aspect9SimulatorILi2EE12solve_stokesEv
#2 4 aspect 0x0000000103a6b751 _ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd + 193: 4 aspect 0x0000000103a6b751 _ZN6aspect9SimulatorILi2EE25assemble_and_solve_stokesEbPd
#3 5 aspect 0x0000000103a6c249 _ZN6aspect9SimulatorILi2EE38solve_single_advection_iterated_stokesEv + 217: 5 aspect 0x0000000103a6c249 _ZN6aspect9SimulatorILi2EE38solve_single_advection_iterated_stokesEv
#4 6 aspect 0x00000001038ea28d _ZN6aspect9SimulatorILi2EE14solve_timestepEv + 237: 6 aspect 0x00000001038ea28d _ZN6aspect9SimulatorILi2EE14solve_timestepEv
#5 7 aspect 0x00000001038e91ad _ZN6aspect9SimulatorILi2EE3runEv + 1165: 7 aspect 0x00000001038e91ad _ZN6aspect9SimulatorILi2EE3runEv
#6 8 aspect 0x000000010473b0b9 _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb + 569: 8 aspect 0x000000010473b0b9 _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbb
#7 9 aspect 0x000000010473a8d5 main + 1093: 9 aspect 0x000000010473a8d5 main
#8 10 libdyld.dylib 0x00007fff205eef3d start + 1: 10 libdyld.dylib 0x00007fff205eef3d start
#9 11 ??? 0x0000000000000002 0x0 + 2: 11 ??? 0x0000000000000002 0x0
--------------------------------------------------------
Aborting!
----------------------------------------------------
Abort(1) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 0
Hi John,
@jbnaliboff - I made couple of test, ASPECT 2.5 and 2.6 really on the same output (please see my answer to Wolfgang). Thank you very much for your detailed explanation and help regarding parameters. After making the edits you suggested, the code looks like below,
# Global parameters
set Dimension = 2
set Start time = 0
set End time = 150e6
set Use years in output instead of seconds = true
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 120
set CFL number = 0.1
set Output directory = new_aspect26_output
set Pressure normalization = no
set Timing output frequency = 1
# Model geometry (200x60 km)
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 100 # for test
set Y repetitions = 20 # for test
set X extent = 200e3
set Y extent = 60e3
end
end
# No mesh refinement, but the coarse mesh is 200x100
subsection Mesh refinement
set Initial adaptive refinement = 0
set Initial global refinement = 0
set Time steps between mesh refinement = 0
end
# The nonlinear solver typically converges more quickly
# when no cheap Stokes solver steps are used for
# problems with Drucker-Prager plasticity.
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 500
set Linear solver tolerance = 1e-7
set GMRES solver restart length = 100
set Use full A block as preconditioner = true
end
end
# Free surface boundary classifications.
# Advecting the free surface vertically rather than
# in the surface normal direction can result in a
# more stable mesh when the deformation is large
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
subsection Free surface
set Surface velocity projection = vertical
end
end
# The top boundary is open (zero traction), which allows the sticky air to
# flow freely through it as topography develops along the wedge. Additional
# testing revealed that using a true free surface leads to significant mesh
# distortion and associated numerical instabilities.
subsection Boundary traction model
set Prescribed traction boundary indicators = top: zero traction
end
# Velocity on boundaries characterized by functions
# Total extension rate is 1 cm/yr (0.5 cm/yr on each side)
subsection Boundary velocity model
set Prescribed velocity boundary indicators = left x: function, right x:function
set Tangential velocity boundary indicators = bottom
subsection Function
set Variable names = x,y
set Function constants = m=-0.05, year=1
set Function expression = if (x<100e3 , -1*m/year, 1*m/year); m*2*100e3/200e3 # Include influx
end
end
# Number and name of compositional fields
# The field plastic_strain is used for tracking the plastic finite strain invariant
# upper: brittle upper crust; seed: 'weaker' brittle region
# lower: viscous lower crust
subsection Compositional fields
set Number of fields = 4
set Names of fields = plastic_strain, upper, seed, lower
end
# Spatial domain of different compositional fields
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0; \
if((y>=50e3 && x<=98.0e3) || (y>=50e3 && x>=102.0e3) || \
(y>=51e3 && x>=96.0e3 && x<=104.0e3), 1, 0); \
if(y>=50e3 && y<51e3 && x>95.0e3 && x<100.0e3, 1, 0); \
if(y<50e3, 1, 0);
end
end
# Composition boundary conditions
subsection Boundary composition model
set List of model names = initial composition
end
# Use discontinuous composition bound preserving limiter
subsection Discretization
set Composition polynomial degree = 2
set Stokes velocity polynomial degree = 2
set Temperature polynomial degree = 2
set Use discontinuous composition discretization = true
subsection Stabilization parameters
set Use limiter for discontinuous composition solution = true # apply the limiter to the DG solutions
set Global composition maximum = 100.0, 1.0, 1.0, 1.0
set Global composition minimum = 0.0, 0.0, 0.0, 1.0
end
end
# Temperature boundary conditions
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1400
set Top temperature = 550
end
end
# Temperature initial conditions (linear)
subsection Initial temperature model
set Model name = function
subsection Function
set Function expression = 550 + (1400 - 550) * y / 60e3
end
end
# Material model
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Reference temperature = 273
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-16
set Minimum viscosity = 1e18
set Maximum viscosity = 1e26
set Thermal diffusivities = 1.8e-6
set Heat capacities = 750.
set Densities = 2800
set Thermal expansivities = 0.
set Viscosity averaging scheme = harmonic
set Viscous flow law = dislocation
set Prefactors for dislocation creep = 5.e-25, 5.e-25, 5.e-26, 5.e-25, 5.e-11
set Stress exponents for dislocation creep = 1.0
set Activation energies for dislocation creep = 0.
set Activation volumes for dislocation creep = 0.
set Angles of internal friction = 30., 30., 30., 4., 30.
set Cohesions = 20.e6, 20.e6, 20.e6, 2.e6, 1.e11
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.5, 0.5, 0.5, 0.5, 0.5
set End plasticity strain weakening intervals = 1.5, 1.5, 1.5, 1.5, 1.5
set Cohesion strain weakening factors = 0.1, 0.1, 0.1, 1.0, 1.0
set Friction strain weakening factors = 0.13333, 0.1333, 0.1333, 1.0, 1.0
end
end
# Gravity model
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 8.87
end
end
# Post processing
subsection Postprocess
set List of postprocessors = velocity statistics, basic statistics, temperature statistics, visualization
subsection Visualization
set List of output variables = density, viscosity, strain rate, error indicator, partition, surface stress, shear stress, surface dynamic topography
set Time between graphical output = 1e5
set Interpolate output = false
end
end
and I can run it without getting an error. However, as you can see in the attached gif file, from the first timestep, the mesh deforms from the edges and starts to rise.
I hope you can help me with this.
GMRES solver restart length = 150, works better.
In post-processing, when I wanted to output the particles for the material (which is not set above), I saw an incredible increase in the calculation process. I wonder if this is normal.
Thanks,
T.
and I can run it without getting an error.
Great news!
However, as you can see in the attached gif file, from the first timestep , the mesh deforms from the edges and starts to rise.
That is happening due to the cells not being allowed to deform on the left and right side. Allowing the mesh to deform on the sides (see below) would likely result in the whole upper boundary moving up, as it appears deformation is not localizing well in the center of your model.
Can you try using a different size weak zone or type of weak seed to initialize deformation? I think that is best path forward. Happy to help with strategies for how to achieve this!
Code to enable deformation on the left and right boundaries:
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface
set Additional tangential mesh velocity boundary indicators = left,right
end
In post-processing, when I wanted to output the particles for the material (which is not set above), I saw an incredible increase in the calculation process. I wonder if this is normal.
This can certainly happen depending on the number of particles being used and other settings. Can you post the particle-related section of the PRM file?
Cheers,
John