Sedimentation not active when run LaMEM using Julia #62
Description
Hi @boriskaus,
I moved the thread from the LaMEM main directory (UniMainzGeo/LaMEM#20) to this one.
Please find attached the .jl example where I find this issue. I also attached the log file where it should appear "Applying sedimentation model (1) to internal free surface." at the end of each timestep but it doesn't. If I run the simulation for 20 Myr for example, there is no sedimentation visible either.
Marc
using LaMEM, GeophysicalModelGenerator, Plots, LaMEM_jll
model = Model(
# Define the grid
Grid(coord_x=[-15.0, 0.0], bias_x=[1.0], nel_x=[512],
coord_y=[0.0, 0.1], bias_y=[1.0], nel_y=[1],
coord_z=[-1.5, 0, 3.5], bias_z=[1.0,1.0], nel_z=[128,16],
nmark_x = 3,
nmark_y = 3,
nmark_z = 3),
# 1 period strainrate, de 1 a 50 Myr and -1e-15 v 1/s compression
BoundaryConditions(exx_num_periods = 2,
exx_time_delims = [5],
exx_strain_rates = [9e-16, 0],
temp_top = 20,
temp_bot = 370,
open_top_bound = 1),
# We use a multigrid solver with 4 levels:
Solver(SolverType = "direct",
DirectSolver = "mumps",
DirectPenalty = 1e3,
MGCoarseSolver = "mumps",
PETSc_options=[
## SNES
"-snes_ksp_ew",
#"snes_ksp_ew_rtol0 1e-4",
#"-snes_ksp_ew_rtolmax 1e-1", #1e-5
#"-snes_max_linear_solve_fail 100",
#"-snes_ksp_ew_gamma 1e-4",
#"-snes_ksp_ew_alpha 0.1",
#"-snes_ksp_ew_alpha2 0.5",
#"-snes_ksp_ew_threshold",
"-snes_max_it 60",
#Newton/Picard options
"-snes_PicardSwitchToNewton_rtol 1e-4",
"-snes_npicard 2",
"-snes_monitor",
"-snes_atol 1e-4",
"-snes_rtol 1e-4",
"-snes_stol 1e-50",
"-snes_linesearch_type l2",
"-res_log",
#"-snes_NewtonSwitchToPicard_it 30", # number of Newton iterations after which we switch back to Picard
# Linesearch options
#"-snes_linesearch_type l2", #Linesearch type (one of) shell basic l2 bt cp (SNESLineSearchSetType)
#"-snes_linesearch_maxstep 1.0", # very important to prevent the code from "blowing up"
#Jacobian solver
"-js_ksp_type fgmres",
"-js_ksp_max_it 60",
#"-js_ksp_min_it 6",
"-js_ksp_converged_reason",
"-js_ksp_monitor",
"-js_ksp_rtol 1e-6",
#Preconditioner
#"-pcmat_type block",
#"-pcmat_pgamma 1e3",
#"-jp_type bf",
#"-bf_vs_type user",
#"-vs_ksp_type preonly",
#"-vs_pc_type lu",
#"-vs_pc_factor_mat_solver_package mumps",
#"-pcmat_type block",
#"-pcmat_no_dev_proj",
#"-jp_type bf",
#"-bf_vs_type mg",
#"-vs_ksp_type preonly",
#"-pcmat_type mono",
#"-jp_type mg",
#"-gmg_pc_view",
#"gmg_dump"
#"-gmg_pc_type mg",
# "-gmg_pc_mg_levels 3",
# "-gmg_pc_mg_galerkin",
# "-gmg_pc_mg_type multiplicative",
# "-gmg_pc_mg_cycle_type v",
# "-gmg_mg_levels_ksp_type richardson",
# "-gmg_mg_levels_ksp_richardson_scale 0.5",
# "-gmg_mg_levels_ksp_max_it 5",
# "-gmg_mg_levels_pc_type jacobi",
#"-crs_ksp_type preonly",
#"-crs_pc_type lu",
#"-crs_pc_factor_mat_solver_package superlu_dist",
#"-crs_pc_type redundant",
#"-crs_pc_redundant_number 2",
#"-crs_redundant_pc_factor_mat_solver_type mumps",
"-crs_ksp_type preonly",
"-crs_pc_type lu",
"-crs_pc_factor_mat_solver_package mumps",
"-objects_dump"]),
ModelSetup(bg_phase = 0,
rand_noise = 0,
advect = "rk2",
interp = "minmod",
stagp_a = 0.7,
mark_ctrl = "subgrid",
nmark_sub = 3),
# Output filename
Output(out_file_name ="Sedimentation_test", # output file name
out_dir ="Sedimentation_test", # output directory name
write_VTK_setup = true,
out_pvd = 1, # activate writing .pvd file
out_phase = 1,
out_density = 0,
out_visc_total = 0,
out_visc_creep = 0,
out_velocity = 0,
out_pressure = 0,
out_tot_press = 0,
out_eff_press = 0,
out_over_press = 0,
out_litho_press = 0,
out_pore_press = 0,
out_temperature = 0,
out_dev_stress = 0,
out_j2_dev_stress = 0,
out_strain_rate = 0,
out_j2_strain_rate = 0,
out_shmax = 0,
out_ehmax = 0,
out_yield = 0,
out_rel_dif_rate = 0,
out_rel_dis_rate = 0,
out_rel_prl_rate = 0,
out_rel_pl_rate = 0,
out_plast_strain = 0,
out_plast_dissip = 0,
out_tot_displ = 0,
out_moment_res = 0,
out_cont_res = 0,
out_energ_res = 0,
out_melt_fraction = 0,
out_fluid_density = 0,
out_conductivity = 0,
out_vel_gr_tensor = 0,
out_surf = 0,
out_surf_pvd = 0,
out_surf_velocity = 0,
out_surf_topography = 0,
out_surf_amplitude = 0,
out_mark = 0,
out_mark_pvd = 0,
out_avd = 1, # activate AVD phase output
out_avd_pvd = 1, # activate writing .pvd file
out_avd_ref = 3, # AVD grid refinement factor
out_ptr = 0,
out_ptr_ID = 0,
out_ptr_phase = 0,
out_ptr_Pressure = 0,
out_ptr_Temperature = 0,
out_ptr_MeltFraction = 0,
out_ptr_Active = 0),
#out_ptr_Grid_Mf = 0),
# Timestepping etc. If the values are in blue in the terminal when you check the script, it means that the've been changed from the default values.
Time(time_end = 30,
dt = 0.001,
dt_min = 1e-10,
dt_max = 0.01,
inc_dt = 0.1,
CFL = 0.5,
CFLMAX = 0.8,
nstep_max = 50,
nstep_out = 5,
nstep_rdb = 100),
SolutionParams(gravity = [0.0, 0.0, -9.81],
FSSA = 1.0,
shear_heat_eff = 0.0,
act_temp_diff = 0,
init_guess = 1,
p_litho_visc = 1,
p_litho_plast = 1,
eta_min = 1e18,
eta_max = 1e23,
eta_ref = 1e22,
act_p_shift = 1),
#quasi_harm_avg = 1),
#PassiveTracers(Passive_Tracer = 1,
# PassiveTracer_Resolution = [400, 1, 3],
# PassiveTracer_Box = [-7,0, 0,10, -7.85,-6.85]),
# #Passive_Tracer_ActiveType = "Always",
# #Passive_Tracer_ActiveValue = 0.1),
FreeSurface(surf_use = 1,
surf_corr_phase = 1,
surf_level = 0.0,
surf_air_phase = 0,
surf_max_angle = 10,
#erosion_model = 1,
#er_num_phases = 10,
#er_time_delims = [0, 3, 6, 9, 12, 15, 18, 21, 24, 30, 35],
#er_rates = [2000, 200, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004], # constant erosion rates in different time periods
#er_levels = [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1], # levels above which we apply constant erosion rates in different time periods
sediment_model = 1, # sedimentation model [0-none (dafault), 1-prescribed rate, 2-cont. margin]
sed_num_layers = 10, # number of sediment layers
sed_time_delims = [0, 3, 6, 9, 12, 15, 18, 21, 24], # sediment layers time delimiters (one less than number)
sed_rates = [0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004, 0.004], # sedimentation rates
sed_phases = [5, 6, 5, 6, 5, 6, 5, 6, 5, 6],
sed_levels = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]), # sediment layers phase numbers
Scaling(GEO_units(temperature = 100,
length = 1km,
viscosity = 1e18,
stress = 1e6Pa)),
)
# Define Phases
################################
# Geometical primitives
################################
model.Grid.Phases .= 0; #Define background phase. The dot is for something that I dont remember
X = model.Grid.Grid.X;
Y = model.Grid.Grid.Y;
Z = model.Grid.Grid.Z; # FreeSurface Geometry
#Variables
FreeSurf = 0;
Angle_Faults = 60;
Left_Boundary = -15;
Right_Boundary = 0;
Bottom_Boundary = -1.5;
Top_Boundary = 3.5;
Thickness_cover_1 = 0.15;
Thickness_cover_2 = 0.15;
Cover_1_deep = FreeSurf - Thickness_cover_1;
Cover_2_deep = FreeSurf - Thickness_cover_1 - Thickness_cover_2;
Salt_Thickness = 0.5;
Basement_Height = FreeSurf - Thickness_cover_1 - Thickness_cover_2 - Salt_Thickness;
Basement_Fault_Position = -8;
#Background phase
model.Grid.Phases[Z.>0.0] .= 0; #
#sticky air
add_box!(model,
xlim = [-15, 0], # los vectores tienen que ser entre corchetes
ylim = [0, 0.1],
zlim = [0, 3.5],
phase = ConstantPhase(0));
#salt
add_box!(model,
xlim = [-15, 0], # los vectores tienen que ser entre corchetes
ylim = [0, 0.1],
zlim = [-1.5, 0],
phase = ConstantPhase(2));
#Cover_1
add_box!(model,
xlim = [-15, 0], # los vectores tienen que ser entre corchetes
ylim = [0, 0.1],
zlim = [Cover_1_deep, FreeSurf],
phase = ConstantPhase(3));
#Cover_2
add_box!(model,
xlim = [-15, 0], # los vectores tienen que ser entre corchetes
ylim = [0, 0.1],
zlim = [Cover_2_deep, Cover_1_deep],
phase = ConstantPhase(4));
model.Grid.Phases[Z.< -Basement_Fault_Position .+ X.*(tand(Angle_Faults)) .&& Z.< Basement_Height] .= 1;
#################################
# Phases
#################################
air = Phase( Name = "Air",
ID = 0,
rho = 1200,
G = 5e10,
ch = 200e6,
fr = 30,
eta = 1e16);
basement = Phase( Name = "basement",
ID = 1,
rho = 2800,
G = 1e10,
ch = 20e6,
fr = 30,
eta = 1e24);
salt = Phase( Name = "salt",
ID = 2,
rho = 2200,
G = 1e10,
ch = 1e6,
fr = 5,
eta = 1e18);
cover_1 = Phase( Name = "cover_1",
ID = 3,
rho = 2700,
G = 1e10,
ch = 10e6,
fr = 30,
eta = 1e22);
cover_2 = Phase( Name = "cover_2",
ID = 4,
rho = 2700,
G = 1e10,
ch = 10e6,
fr = 30,
eta = 1e22);
sed_1 = Phase( Name = "sed_1",
ID = 5,
rho = 2500,
#G = 1e10,
#ch = 10e6,
#fr = 30,
eta = 1e21);
sed_2 = Phase( Name = "sed_2",
ID = 6,
rho = 2500,
#G = 1e10,
#ch = 10e6,
#fr = 30,
eta = 1e21);
add_phase!(model, air, basement, salt, cover_1, cover_2, sed_1, sed_2)
plot_cross_section(model, y=0.05, field=:phase)
savefig("Sedimentation_test.png")
run_lamem(model, 6)