Sedimentation not active when run LaMEM using Julia

Hi @boriskaus,
I moved the thread from the LaMEM main directory (https://github.com/UniMainzGeo/LaMEM/issues/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
sed_test_log.txt
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)
JuliaGeodynamics / LaMEM.jl

Sedimentation not active when run LaMEM using Julia #62
