erosion_model=2 (rescribed rate with given level) does not work

[wenrongcao]

Dear Boris and LaMEM.jl developers,

The erosion_model=2 (rescribed rate with given level) does not work in LaMEM.jl (v.0.4.1). Below is the minimal script to reproduce the issue.

Basically, the script generates an initial surface level at 10 km and sets a constant erosion rate at 1 km/Myr above the 0 km baseline. Yet, after the code runs for 5 Myr, nothing happens.

erosion_model=1 (infinitely fast) does work.

Wenrong

using LaMEM, GeophysicalModelGenerator,  Plots

model = Model(Grid(nel=(32,1,32), x=[-50,50], z=[-50,20], y=[-1,1] ), 
                Scaling(GEO_units(stress=1000MPa, viscosity=1e20Pa*s)),
                 Time(dt=1e-2, dt_min=1e-5, dt_max=1e-1, nstep_out=5, nstep_max=200, time_end=5))


model.FreeSurface = FreeSurface(    surf_use        = 1,                # free surface activation flag
                                    surf_corr_phase = 1,                # air phase ratio correction flag (due to surface position)
                                    surf_level      = 10.0,              # initial level
                                    surf_air_phase  = 0,                # phase ID of sticky air layer
                                    surf_max_angle  = 40.0,             # maximum angle with horizon (smoothed if larger))
                                    erosion_model   = 2,                # 2-prescribed rate with given level
                                    er_num_phases   = 2,                # number of erosion phases
                                    er_time_delims  = [1],         # erosion time delimiters (one less than number)
                                    er_rates        = [0.1, 0.1],       # constant erosion rates in different time periods (0.1 cm/yr = 1mm/yr=1km/Myr)
                                    er_levels       = [0,0]        # levels above which we apply constant erosion rates in different time period
                                )
      
model.Output = Output(  out_surf            = 1, 	
                        out_surf_pvd        = 1,
                        out_surf_topography = 1,
                                 )                                

# add an air phase, phase =0 
add_box!(model;  xlim    = (-50, 50), 
                ylim    = (model.Grid.coord_y...,), 
                zlim    = (0, 20),
                phase   = ConstantPhase(0),
                T       = nothing )    
                
# add a crust phase, phase =1
add_box!(model;  xlim    = (-50, 50), 
                ylim    = (model.Grid.coord_y...,), 
                zlim    = (-50.0, 10.0),
                phase   = ConstantPhase(1),
                T       = nothing )  

air    = Phase(ID=0, Name="Air",    eta=1e19, rho=50, ch=10e6, fr=0);
crust  =  Phase(ID=1, Name="crust",  eta=1e21, rho=2700, ch=30e6, fr=20);

rm_phase!(model)
add_phase!(model, air, crust)

run_lamem(model,1)

[wenrongcao]

The problem is caused by the incorrect translation from Julia code to the .bat file: the keyword erosion_model = 2 is missing in the resulted output.dat file.

After manually adding " erosion_model = 2" to the output.dat file. The erosion_model=2 seems to work.

I will see if I can fix that and make a pull request.