Problem using a free surface with sudden change in topography

ASPECT
1

Dear all,

I am running a model using a free surface, but I get an error that I do not understand when the simulation reaches ~4 million years. I guess the problem is due to the sudden change in topography, but I would like to know if anyone knows how to solve this problem. The error message is as follows:

ERROR MESSAGE:

*** Timestep 2837: t=4.64329e+06 years, dt=276.765 years
Solving mesh velocity system… 6 iterations.
Solving temperature system… 1 iterations.
Solving oceanic_crust system … 0 iterations.
Solving continental_uppercrust system … 1 iterations.
Solving continental_lowercrust system … 1 iterations.
Solving lith system … 1 iterations.
Rebuilding Stokes preconditioner…
Solving Stokes system… 0+23 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.011231

Rebuilding Stokes preconditioner…
Solving Stokes system… 0+9 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 2: 0.00020052

Rebuilding Stokes preconditioner…
Solving Stokes system… 0+5 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 3: 3.76123e-05

Postprocessing:
RMS, max velocity: 0.019 m/year, 0.464 m/year
Temperature min/avg/max: 293.2 K, 1522 K, 1643 K
Topography min/max: -1.147e+04 m, 2092 m

*** Timestep 2838: t=4.64356e+06 years, dt=271.978 years
Solving mesh velocity system… 6 iterations.
mpirun(20952,0x16ba7f000) malloc: can’t allocate region
:*** mach_vm_map(size=1125899906859008, flags: 100) failed (error code=3)
mpirun(20952,0x16ba7f000) malloc: *** set a breakpoint in malloc_error_break to debug
[cvi067114:20952] *** Process received signal ***
[cvi067114:20952] Signal: Segmentation fault: 11 (11)
[cvi067114:20952] Signal code: Invalid permissions (2)
[cvi067114:20952] Failing at address: 0x0
[cvi067114:20952] [ 0] 0 libsystem_platform.dylib 0x000000019bb32c44 _sigtramp + 56
[cvi067114:20952] [ 1] 0 mca_rml_oob.so 0x0000000104fe7d18 orte_rml_oob_send_buffer_nb + 892
[cvi067114:20952] [ 2] 0 libopen-rte.40.dylib 0x000000010469a2fc pmix_server_log_fn + 492
[cvi067114:20952] [ 3] 0 mca_pmix_pmix3x.so 0x0000000104d4cd24 server_log + 804
[cvi067114:20952] [ 4] 0 mca_plog_default.so 0x0000000104f6b8cc mylog + 472
[cvi067114:20952] [ 5] 0 mca_pmix_pmix3x.so 0x0000000104dc8bc8 pmix_plog_base_log + 1140
[cvi067114:20952] [ 6] 0 mca_pmix_pmix3x.so 0x0000000104d91490 pmix_server_log + 1788
[cvi067114:20952] [ 7] 0 mca_pmix_pmix3x.so 0x0000000104d794f4 pmix_server_message_handler + 5092
[cvi067114:20952] [ 8] 0 mca_pmix_pmix3x.so 0x0000000104dd0c64 pmix_ptl_base_process_msg + 632
[cvi067114:20952] [ 9] 0 libevent_core-2.1.7.dylib 0x0000000104845250 event_process_active_single_queue + 520
[cvi067114:20952] [10] 0 libevent_core-2.1.7.dylib 0x00000001048422d4 event_base_loop + 948
[cvi067114:20952] [11] 0 mca_pmix_pmix3x.so 0x0000000104d9bd38 progress_engine + 36
[cvi067114:20952] [12] 0 libsystem_pthread.dylib 0x000000019bae7878 _pthread_start + 320
[cvi067114:20952] [13] 0 libsystem_pthread.dylib 0x000000019bae25e0 thread_start + 8
[cvi067114:20952] *** End of error message ***
zsh: segmentation fault mpirun -np 8 ~/aspect/aspect-release Model_Gibraltar.prm
(base) pedro_gea@cvi067114 Velocidades_Gibraltar %

I copy a part of my input file in case it is helpful:

A PART OF THE INPUT FILE

set Dimension = 2
set World builder file = Slab_Geometry.wb
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 10
set CFL number = 0.1
set Pressure normalization = no

subsection Solver parameters
set Composition solver tolerance = 1e-5
set Temperature solver tolerance = 1e-5
subsection Stokes solver parameters
set Linear solver tolerance = 1e-5
set Linear solver A block tolerance = 5e-1
set Linear solver S block tolerance = 1e-6
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 0
set GMRES solver restart length = 200
set Maximum number of expensive Stokes solver steps = 1000
set Use full A block as preconditioner = false
end
end

subsection Geometry model
set Model name = box
subsection Box
set X extent = 1320e3
set Y extent = 660e3
set X repetitions = 2
set Y repetitions = 1
end
end

subsection Compositional fields
set Number of fields = 4
set Names of fields = oceanic_crust, continental_uppercrust, continental_lowercrust, lith
end

subsection Initial composition model
set Model name = world builder
end

subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end

subsection Material model
set Material averaging = harmonic average
set Model name = visco plastic
subsection Visco Plastic
set Reference temperature = 293
set Reference viscosity = 1e21
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-15
set Minimum viscosity = 1e19
set Maximum viscosity = 1e23
set Viscosity averaging scheme = harmonic
set Grain size = 1.0e-2
set Thermal diffusivities = 0.8e-6
set Heat capacities = 1250.
set Densities = 3300, 3000, 2750, 2900, 3300
set Thermal expansivities = 3.5e-5
set Viscous flow law = composite
set Prefactors for dislocation creep = 2.28e-18, 2.5e-23, 2.5e-23, 2.5e-50, 2.28e-18
set Stress exponents for dislocation creep = 3.5, 1, 1, 1, 3.5
set Activation energies for dislocation creep = 480e3, 0, 0, 0, 480e3
set Activation volumes for dislocation creep = 1.1e-5, 0, 0, 0, 1.1e-5
set Prefactors for diffusion creep = 4.7e-16, 2.5e-23, 2.5e-23, 5e-19, 4.7e-16
set Activation energies for diffusion creep = 335e3, 0, 0, 170e3, 335e3
set Activation volumes for diffusion creep = 4e-6, 0, 0, 0, 4e-6
set Grain size exponents for diffusion creep = 3, 0, 0, 7, 3
set Angles of internal friction = 0
set Cohesions = 1e30
end
end

subsection Mesh refinement
set Initial global refinement = 4
set Initial adaptive refinement = 4
set Minimum refinement level = 0
set Normalize individual refinement criteria = true
set Refinement criteria merge operation = max
set Skip solvers on initial refinement = false
set Strategy = viscosity, composition, temperature, density, slope
set Time steps between mesh refinement = 10
set Coarsening fraction = 0.1
set Refinement fraction = 0.9
end

subsection Discretization
subsection Stabilization parameters
set Use limiter for discontinuous composition solution = false
set Use limiter for discontinuous temperature solution = false
set Stabilization method = entropy viscosity
set alpha = 2
set beta = 0.052
set cR = 0.11
set gamma = 0.0
end
end

subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface & diffusion
subsection Diffusion
set Time steps between diffusion = 5
end

subsection Free surface
set Free surface stabilization theta = 0.8
set Surface velocity projection = vertical
end
end

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

END OF THE INPUT FILE

I have been playing with the parameter Free surface stabilization theta, giving it values from 0.5 to 1 but it has not worked. Also, I tested other CL values and tried other values for cR and beta in the subsection Discretization (specifically cR=0.5 and beta=0.078).

Could someone give me some advice on how to improve the behavior of the simulation?

I attach a picture of the simulation before it crashes.

Thank you very much!
Pedro

Free surface
Free surface1920×689 204 KB

2

Pedro:
The important part of the error message is this one:

In other words, your system has run out of memory. This may be because you are running too large a problem, or because there is a bug in the code – it’s hard to tell without access to the actual system. Can you say how many unknowns your problem has, and how many MPI processes you are using?
Best
Wolfgang

3

Thanks for your reply!

At the start of the simulation (after the initial refinement steps) there are 650000 degrees of freedom. Once the simulation starts and the mesh is refined every 10 time steps, the problem reaches ~1200000 degrees of freedom. I’m using 8 MPI processors.

The simulation works well if I use a free slip boundary condition at the top of the model (in that case the number of DOF is also ~1200000). The problem seems to appear when I add the free surface but I don’t know where the error is.

Pedro

4

Interesting. Does the problem happen reproducibly in the same place? Does it depend on the number of processes you use? If so: Do you think you could provide all of the input files for others to reproduce the issue?

Best
Wolfgang

5

Yes, the problem happens after 4 million years and does not seem to depend on the number of processors used.
I attach the .prm and .wb files used in the simulation in case you see any bugs in the code or if others can reproduce the same issue.

Thank you!

Best
Pedro

Gibraltar model.zip (5.9 KB)

6

Pedro – thanks for the input files. I’ve tried to run them, but I run into an error in time step 3:

*** Timestep 3:  t=76.8314 years, dt=30.2566 years
   Solving mesh velocity system... 7 iterations.

--------------------------------------------------------
An error occurred in line <103> of file <../source/material_model/rheology/diffusion_creep.cc> in function
    double aspect::MaterialModel::Rheology::DiffusionCreep<dim>::compute_viscosity(double, double, unsigned int, const std::vector<double, std::allocator<double> >&, const std::vector<unsigned int>&) const [with int dim = 2]
The violated condition was: 
    viscosity_diffusion > 0.0
Additional information: 
    Negative diffusion viscosity detected. This is unphysical and should
    not happen. Check for negative parameters. Temperature and pressure
    are -42.604121044089531 K, -242761935.67231774 Pa.

That raises the question whether you are running in release or debug mode when you encounter the bug you show above?

Best
W.

7

First of all, thanks for your time and your patient.

I ran the model in release mode directly because the simulation seemed to work well without the free surface. However, I did not see the error you have shown. It is weird. I guess the error must be in the input files but I am not able to find it. Any ideas as to the cause of the problem?

Thanks again!

Pedro

8

Well, you end up with negative temperatures and pressures. I would start by
looking at where that happens and why.

Best
W.