Merge remote-tracking branch 'origin/master'

README.md

@@ -1,6 +1,7 @@
 # LaMEM
 [![Build Status](https://github.com/UniMainzGeo/LaMEM/workflows/CI/badge.svg)](https://github.com/UniMainzGeo/LaMEM/actions)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://unimainzgeo.github.io/LaMEM/dev)
+[![DOI](https://zenodo.org/badge/641438500.svg)](https://zenodo.org/doi/10.5281/zenodo.10211634)
 
 LaMEM (Lithosphere and Mantle Evolution Model) is a parallel 3D numerical code that can be used to simulate various thermo-mechanical 
 geodynamical processes such as mantle-lithosphere interaction for rocks 

input_models/ScalingTests/Volcano.dat

@@ -8,23 +8,23 @@
 units                       =	geo
 
 # Always in SI units!!
-unit_temperature = 500
+unit_temperature = 1000
 unit_length      = 1000
-unit_viscosity   = 1e19
+unit_viscosity   = 1e20
 unit_stress      = 1e9
 
 #===============================================================================
 # Time stepping parameters
 #===============================================================================
 
-	dt        = 5e-4	# time step
-	dt_min    = 4.5e-4	# minimum time step (declare divergence if lower value is attempted)
+	dt        = 1e-7	# time step
+	dt_min    = 1e-9	# minimum time step (declare divergence if lower value is attempted)
 	dt_max    = 1e-3	# maximum time step
-	inc_dt    = 1		# time step increment per time step (fraction of unit)
+	inc_dt    = 0.5		# time step increment per time step (fraction of unit)
 	CFL       = 0.5		# CFL (Courant-Friedrichs-Lewy) criterion
 	CFLMAX    = 0.8		# CFL criterion for elasticity
 	nstep_ini = 0		# save output for n initial steps
-	nstep_max = 5		# maximum allowed number of steps (lower bound: time_end/dt_max)
+	nstep_max = 1		# maximum allowed number of steps (lower bound: time_end/dt_max)
 	nstep_out = 1		# save output every n steps
 
 #===============================================================================
@@ -64,8 +64,8 @@ unit_stress      = 1e9
 
 #	noslip   = 0 0 0 0 0 0
 
-	temp_top	    =   0			# Temperature @ top
-	temp_bot		=   1350     	# Temperature @ bottom; side BC's are flux-free
+	temp_top	=   0			# Temperature @ top
+	temp_bot	=   800     	# Temperature @ bottom; side BC's are flux-free
 
 	# Background strain rate parameters
 #	exx_num_periods  = 1      # number intervals of constant strain rate (x-axis)
@@ -80,27 +80,30 @@ unit_stress      = 1e9
 
 	gravity         = 0.0 0.0 -9.81  # gravity vector
 	act_temp_diff   = 1		 		 # temperature diffusion activation flag
-	#act_steady_temp = 0              # steady-state temperature initial guess activation flag
-	#steady_temp_t   = 0.1            # time for (quasi-)steady-state temperature initial guess
-	#nstep_steady    = 50             # number of steps for (quasi-)steady-state temperature initial guess (default = 1)
+	act_steady_temp = 1              # steady-state temperature initial guess activation flag
+	steady_temp_t   = 0.1            # time for (quasi-)steady-state temperature initial guess
+	nstep_steady    = 150             # number of steps for (quasi-)steady-state temperature initial guess (default = 1)
 	act_heat_rech   = 0              # heat recharge activation flag
 	init_lith_pres  = 1              # initial pressure with lithostatic pressure
 	init_guess      = 1              # initial guess flag
-	eta_min         = 1e19           # viscosity upper bound
+	eta_min         = 1e17           # viscosity upper bound
 	eta_max         = 1e23           # viscosity lower limit
 	eta_ref         = 1e22           # reference viscosity (initial guess)
 	T_ref           = 20             # reference temperature
 	p_lim_plast		= 1
-
+     	min_cohes       = 1e6            # cohesion lower bound  [Pa]
+   	min_fric        = 1.0            # friction lower bound  [degree]
+    	tau_ult         = 1e9            # ultimate yield stress [Pa]
+
 #===============================================================================
 # Solver options
 #===============================================================================
 
 	SolverType 	 =	multigrid	# solver [direct or multigrid]
-	MGLevels 	 =	6			# number of MG levels [default=3]
-    MGSweeps 	 =	5			# number of MG smoothening steps per level [default=10]
-    MGSmoother 	 =	jacobi 		# type of smoothener used [chebyshev or jacobi]
-    MGJacobiDamp =	0.5			# Dampening parameter [only employed for Jacobi smoothener; default=0.6]
+	MGLevels 	 =	3			# number of MG levels
+	MGSweeps 	 =	25			# number of MG smoothening steps per level
+	MGSmoother 	 =	jacobi 		# type of smoothener used [chebyshev or jacobi]
+	MGJacobiDamp     =	0.5			# Dampening parameter [only employed for Jacobi smoothener; default=0.6]
 
 
 #===============================================================================
@@ -125,6 +128,13 @@ unit_stress      = 1e9
     <EllipsoidEnd>
 
 
+    <EllipsoidStart>
+        phase       = 1
+        axes        = 14.0 14.0 1.0  # semi-axes of ellipsoid in x, y and z
+        center      = 0.0 0.0 -4.0
+        Temperature = constant     # optional: Temperature of the sphere. possibilities: [constant]
+        cstTemp     = 1000         # required in case of [constant]: temperature value [in Celcius in case of GEO units]
+    <EllipsoidEnd>
 #===============================================================================
 # Output
 #===============================================================================
@@ -177,7 +187,7 @@ unit_stress      = 1e9
 	<MaterialStart>
 		Name = air
 		ID	= 0
-		rho	= 100
+		rho	= 1
 		alpha   = 3e-5 
 
 	# Linear Viscosity	
@@ -192,7 +202,7 @@ unit_stress      = 1e9
 		Cp 	= 1000
 
 	# Plastic parameters
-		ch 	= 10e6
+		ch 	= 20e6
 		fr 	= 30
 	<MaterialEnd>
 
@@ -216,9 +226,9 @@ unit_stress      = 1e9
 		Cp 	= 1000
 
 	# Plastic parameters
-		ch 		   =  10e6
-		fr 		   =  30
-	  	eta_st     =  5e20   # stabilization viscosity (default is eta_min)
+		ch 		   =  5e6
+		fr 		   =  10
+	  	eta_st     	   =  5e20   # stabilization viscosity (default is eta_min)
 
 	<MaterialEnd>
 
@@ -244,7 +254,7 @@ unit_stress      = 1e9
 
 	# Plastic parameters
 		ch 		   =  1e6
-		fr 		   =  10
+		fr 		   =  5
 	  	eta_st     =  5e20   # stabilization viscosity (default is eta_min)
 
 	<MaterialEnd>
@@ -276,7 +286,7 @@ unit_stress      = 1e9
 	-snes_max_it 200
 	-snes_rtol 1e-5
 	-snes_atol 1e-6
-	-snes_PicardSwitchToNewton_rtol 5e-2 #-7
+	-snes_PicardSwitchToNewton_rtol 1e-3 #-7
 	#-snes_NewtonSwitchToPicard_it   1       # number of Newton iterations after which we switch back to Picard
 #	-snes_monitor
 	-snes_linesearch_monitor