Plume_PhaseTransitions.dat

# This is 2D setup to simulate plume-lithosphere interaction with a heterogenous plume 
# impinging on a moving lithosphere, with the aim to understand how long heterogeneities 
# can stay stable underneath the moving lithosphere 


#===============================================================================
# Scaling
#===============================================================================
	units            = geo		# geological units 
	
	unit_temperature = 1000
	unit_length      = 100e3
	unit_viscosity   = 1e20
	unit_stress      = 1e9
	
#===============================================================================
# Time stepping parameters
#===============================================================================
	dt        = 0.01   # initial time step
	dt_min    = 1e-5  	# minimum time step (declare divergence if lower value is attempted)
	dt_max    = 1   	# maximum time step
	inc_dt    = 0.1   	# time step increment per time step (fraction of unit)
	CFL       = 0.5   	# CFL (Courant-Friedrichs-Lewy) criterion
	nstep_max = 70  	# maximum allowed number of steps (lower bound: time_end/dt_max)
	nstep_out = 10     	# save output every n steps
	nstep_rdb = 2000   	# save restart database every n steps


#===============================================================================
# Grid & discretization parameters
#===============================================================================

# Number of cells for all segments
	nel_x 	= 	64
	nel_y 	= 	1
	nel_z 	= 	64

# Coordinates of all segments (including start and end points)

	coord_x = -500  500
	coord_y = -10   10
	coord_z = -1000 50

#===============================================================================
# Free surface
#===============================================================================
	surf_use           = 0                # free surface activation flag
	surf_corr_phase    = 1                # air phase ratio correction flag (due to surface position)
	surf_level         = 0                # initial level
	surf_air_phase     = 0                # phase ID of sticky air layer
	surf_max_angle     = 10.0             # maximum angle with horizon (smoothed if larger)
	
#===============================================================================
# Boundary conditions
#===============================================================================
	
# temperature on the top & bottom boundaries
	temp_top  = 0.0
	temp_bot  = 1300.0;
	
# No-slip boundary flag mask (left right front back bottom top)
	noslip 					= 	0 0 0 0 1 0
	open_top_bound 			= 	1


#===============================================================================
# Solution parameters & controls
#===============================================================================

	gravity        = 0.0 0.0 -10    # gravity vector
	FSSA           = 1.0            # free surface stabilization parameter [0 - 1]
	act_p_shift    = 1              # pressure shift activation flag (zero pressure in the top cell layer)
	init_guess     = 0              # initial guess flag
	act_temp_diff  = 0  			# activate thermal diffusion
	shear_heat_eff = 0.0  
	eta_min        = 1e18           # viscosity upper bound
	eta_ref 	   = 1e20			# reference viscosity for initial guess	
	eta_max        = 1e24           # viscosity lower limit
	#p_shift        = 2.216e4        # pressure shift activation flag (zero pressure in the top cell layer)
    init_lith_pres = 1
	DII_ref        = 1e-15          # background (reference) strain-rate
	Phasetrans     = 1	
#===============================================================================
# Solver options
#===============================================================================
	SolverType 		=	direct 			# solver [direct or multigrid]
	DirectSolver 	=	superlu_dist			# mumps/superlu_dist/pastix	
	DirectPenalty 	=	1e3

#===============================================================================
# Model setup & advection
#===============================================================================
	msetup         = geom              # setup type
	nmark_x        = 3                 # markers per cell in x-direction
	nmark_y        = 3                 # ...                 y-direction
	nmark_z        = 3                 # ...                 z-direction
	bg_phase       = 0                 # background phase ID
	rand_noise     = 1                 # random noise flag
	advect         = rk2               # advection scheme
	interp         = stagp             # velocity interpolation scheme
	mark_ctrl      = subgrid           # marker control type
	nmark_lim      = 8 100             # min/max number per cell
	nmark_sub       = 1                 # max number of same phase markers per subcell (subgrid marker control)


# Geometric primitives to define the initial setup:
	# sticky air
	<LayerStart>
		phase  	= 	0
		top 	= 	300  
		bottom 	= 	0  
		
		Temperature = constant
		cstTemp     = 0
	<LayerEnd>
	
	# Define mantle and lithosphere (for Temperature structure)
	<LayerStart>
		phase  	= 	3
		top 	= 	0  
		bottom 	= 	-1000

		Temperature = halfspace
		thermalAge  = 100		# thermal age
	
		botTemp     = 1300		# bottom T for halfspace cooling
		topTemp     = 0			# top T for halfspace cooling
	<LayerEnd>
	
	# Define oceanic crust (for Phase)
	<LayerStart>
		phase  	= 	1
		top 	= 	0  
		bottom 	= 	-10
	<LayerEnd>
	
	# Define mantle lithosphere (for Phase)
	
	# Define plume (for Phase)
	<SphereStart>
		phase  		= 	4		# required; phase of sphere
		center 		= 	0.0 0.0 -550	
		radius 		= 	100		
		
		Temperature = 	constant
		cstTemp     = 	1400		# temperature
	<SphereEnd>


#===============================================================================
# Output
#===============================================================================

# Grid output options (output is always active)

	out_file_name       = Plume_PhaseTransitions # output file name
	out_pvd             = 1       	# activate writing .pvd file
	out_phase           = 1
	out_density         = 1
	out_visc_total      = 1
	out_visc_creep      = 1
	out_visc_plast      = 1
	out_velocity        = 1
	out_pressure        = 1
    out_tot_press       = 1  
    out_eff_press       = 1
	out_temperature     = 1
	out_dev_stress      = 1
	out_j2_dev_stress   = 1
	out_strain_rate     = 1
	out_j2_strain_rate  = 1
	out_yield           = 1
	out_plast_strain    = 1
	out_plast_dissip    = 1
	out_tot_displ       = 1
	out_moment_res      = 1
	out_cont_res        = 1
	
# AVD phase viewer output options (requires activation)

	out_avd     		= 1 # activate AVD phase output
	out_avd_pvd 		= 1 # activate writing .pvd file
	out_avd_ref 		= 3 # AVD grid refinement factor
	
#===============================================================================
# Material phase parameters
#===============================================================================

	# Define properties of sticky air
	<MaterialStart>
		ID  		= 	0 				# phase id
		rho 		= 	3300 			# density
		eta 		= 	1e22 			# viscosity
		alpha  	 	= 	3e-5 

		rho 		= 	0 				# density

		# Thermal properties
        k           =   100        		# conductivity
        Cp          =   1e6        		# heat capacity - should be artificially large for sticky air
	<MaterialEnd>
	
	# Define properties of oceanic crust
	<MaterialStart>
		ID  		= 1       # phase id

		rho        	= 	3300			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		eta  		= 	1e24

     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>
	
	# Define properties of oceanic mantle lithosphere
	<MaterialStart>
		ID  		= 	2       # phase id

		rho        	= 	3300			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		eta 	    =  	1e23
		
     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>
	
	# Define properties of upper mantle
	<MaterialStart>
		ID  		= 	3 				# phase id
		
		rho        	= 	3300			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		# Linear viscous	
		eta 		= 	1e20
		
     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>

	# Define properties of plume
	<MaterialStart>
		ID  		= 	4 	# phase id
		
		rho        	= 	3300			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		eta 		= 	1e20
		
     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>

	# Define properties of lower mantle
	<MaterialStart>
		ID  		= 	5 				# phase id
		
		rho        	= 	3300			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		# Linear viscous	
		eta 		= 	1e21
		
     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>

	# Define properties of plume #2
	<MaterialStart>
		ID  		= 	6 	# phase id
		
		rho        	= 	3350			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		eta 		= 	1e20
		
     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>


	# Define properties of upper mantle, different color
	<MaterialStart>
		ID  		= 	7 				# phase id
		
		rho        	= 	3300			# density [kg/m3]
		alpha  	 	= 	3e-5 			# coeff. of thermal expansion [1/K]
		
		# Linear viscous	
		eta 		= 	1e20
		
     	# Thermal properties
        k           =   3        		# conductivity
        Cp          =   1050        	# heat capacity
	<MaterialEnd>


#===============================================================================
# Define phase transitions
#===============================================================================
	
	# T-dependent:
	<PhaseTransitionStart>
	    ID   					= 	0  			# Phase_transition law ID
     	Type 					= 	Constant		# [Constant, Clapeyron, Box_type]
		Parameter_transition 	= 	T 			# T = Temperature, P = Pressure, Depth = z coordinate
		ConstantValue 			= 	1200			# The value @ which the phase transition occurs
		
	   	number_phases 			= 	1				# Number of phases involved			
		PhaseAbove 				= 	3 		        # Phases below the phase transition 
		PhaseBelow 				= 	2               # Phases below
        PhaseDirection          =   BothWays    # [BothWays=default; BelowToAbove; AboveToBelow]
	<PhaseTransitionEnd>

	# Depth-dependent:
	<PhaseTransitionStart>
	    ID   					= 	1   			# Phase_transition law ID
     	Type 					= 	Constant		# [Constant, Clapeyron, Box_type]
		Parameter_transition 	= 	Depth 			# T = Temperature, P = Pressure, Depth = z coordinate
		ConstantValue 			= 	-400			# The value @ which the phase transition occurs
		
	   	number_phases 			= 	1				# Number of phases involved			
		PhaseAbove 				= 	6 		        # Phases below the phase transition 
		PhaseBelow 				= 	4               # Phases below
        PhaseDirection          =   BelowToAbove    # [BothWays=default; BelowToAbove; AboveToBelow]
		
		ResetParam 				=	APS				# Set APS to zero if we are above the constant value

	<PhaseTransitionEnd>

	# Clapeyron slope
	<PhaseTransitionStart>
		ID                      =   2                           # Phase_transition law ID
		Type                    =   Clapeyron                   # Use the pressure retrieved by the clapeyron linear equation to trigger the phase transition
		Name_Clapeyron          =   Mantle_Transition_660km     # predefined profiles; see SetClapeyron_Eq in phase_transition.cpp for options 
        number_phases           =   1
		PhaseAbove              =   5
		PhaseBelow              =   3
        PhaseDirection          =   BothWays    # [BothWays=default; BelowToAbove; AboveToBelow]
	
	<PhaseTransitionEnd>

	# Box-like region with T-condition
	<PhaseTransitionStart>
    	ID                      =   3                               # Phase_transition law ID
		Type                    =   Box                             # A box-like region
	    PTBox_Bounds            =   200 400 -100 100 -1000 -500        # box bound coordinates: [left, right, front, back, bottom, top]
        number_phases           =   1
        PhaseInside             =   7							    # Phase within the box [use -1 if you don't want to change the phase inside the box]
        PhaseOutside            =   3							    # Phase outside the box
        PhaseDirection          =   BothWays                        # [BothWays=default; OutsideToInside; InsideToOutside]
	
        PTBox_TempType          =	linear 					        # Temperature condition witin the box [none, constant, linear, halfspace]
        PTBox_topTemp           =   20                              # Temp @ top of box [for linear & halfspace]
        PTBox_botTemp           =   1300                            # Temp @ bottom of box [for linear & halfspace]
        
        #PTBox_thermalAge        =   100                            # Thermal age 
        #PTBox_cstTemp          =   1200                            # Temp within box [for constant T]
       
	   	ResetParam 				=	APS

	<PhaseTransitionEnd>

#===============================================================================
# PETSc options
#===============================================================================

<PetscOptionsStart>
	# LINEAR & NONLINEAR SOLVER OPTIONS
#	-snes_ksp_ew
#	-snes_ksp_ew_rtolmax 1e-4
	-snes_rtol 1e-4					
	-snes_atol 1e-4					
	-snes_max_it 50					
	
	-snes_PicardSwitchToNewton_rtol 1e-4   # relative tolerance to switch to Newton (1e-2)
	-snes_NewtonSwitchToPicard_it  	20     # number of Newton iterations after which we switch back to Picard

	
	-js_ksp_type fgmres
	-js_ksp_monitor 				# display how the inner iterations converge
	-js_ksp_max_it 20				# inner
	-js_ksp_atol 1e-8
	-js_ksp_rtol 1e-4

	-snes_linesearch_type l2
	-snes_linesearch_monitor
	-snes_linesearch_maxstep 10
	
<PetscOptionsEnd>

#===============================================================================