time-dependent phase transition not working

[wenrongcao]

I tested a simple time-dependent phase transition in which phase # 1 transforms to phase # 0 when simulation time = 1 Myr. I am following the example given in the lamem_input.dat file. However, the phase transition was immediately performed at step 1 (time=0.011 Myr) but not at the required time at 1 Myr.

A minimal code to reproduce this using LaMEM.jl v0.4.8:

using LaMEM
using GeophysicalModelGenerator

# Main model setup
model = Model(Grid(nel=(32,1,32), x=[-50,50], z=[-50,0], 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),
        )
# add a crust_upper phase, phase =0
add_box!(model;  xlim    = (-50, 50), 
        ylim    = (model.Grid.coord_y...,), 
        zlim    = (-25.0, 0.0),
        phase   = ConstantPhase(0),
        T       = nothing )  
        
# add a crust_lower phase, phase =1
add_box!(model;  xlim    = (-50, 50), 
        ylim    = (model.Grid.coord_y...,), 
        zlim    = (-50.0, -25.0),
        phase   = ConstantPhase(1),
        T       = nothing )           

crust_upper  =  Phase(ID=0, Name="crust_upper",  eta=1e21, rho=2700, ch=30e6, fr=20);
crust_lower  =  Phase(ID=1, Name="crust_lower",  eta=1e21, rho=3000, ch=30e6, fr=20);

# Define Phase Transition laws.
# time dependent phase transition: at 1 Myr, crust_lower (phase #1) transforms to crust_upper (phase #0)

PT0 = PhaseTransition(ID=0, Type="Constant", Parameter_transition="t", PhaseBelow=[0], PhaseAbove=[1], PhaseDirection="AboveToBelow", ConstantValue=1, ResetParam="APS")

model.Materials.PhaseTransitions = [PT0]

add_phase!(model, crust_upper, crust_lower)
run_lamem(model,1);

Note, the depth-dependent transition works well. If I only replace the phase transition line with the following line, it works properly.

PT0 = PhaseTransition(ID=0, Type="Constant", Parameter_transition="Depth",  PhaseBelow=[0], PhaseAbove=[1], PhaseDirection="AboveToBelow", ConstantValue=-40, ResetParam="APS")

I also tested the same code using LaMEM v2.2.0, same result. The time-dependent phase transition was performed at step 1. The phase transition section in the output.dat file is:

#===============================================================================
# Define phase transitions
#===============================================================================

   # Define Phase Transition laws (maximum 10)
   <PhaseTransitionStart>
        ID               =  0     # Phase_transition law ID
        Type             =  Constant     # [Constant, Clapeyron, Box]: Constant - the phase transition occurs only at a fixed value of the parameter; Clapeyron - clapeyron slope
        Parameter_transition  =  t     # [T = Temperature, P = Pressure, Depth = z-coord, X=x-coord, Y=y-coord, APS = accumulated plastic strain, MeltFraction, t = time] parameter that triggers the phase transition
        ConstantValue    =  1.0     # Value of the parameter [unit of T,P,z, APS] 
        number_phases    =  1     # The number of involved phases [default=1]
        PhaseAbove       =  1     # Above the chosen value the phase is 1, below it, the value is PhaseBelow
        PhaseBelow       =  0     # Below the chosen value the phase is PhaseBelow, above it, the value is 1
        PhaseDirection   =  AboveToBelow     # [BothWays=default; BelowToAbove; AboveToBelow] Direction in which transition works
        ResetParam       =  APS     # [APS] Parameter to reset on particles below PT or within box
   <PhaseTransitionEnd>