1 | import numpy as np
|
---|
2 | from SetIceSheetBC import SetIceSheetBC
|
---|
3 | from frictionshakti import frictionshakti
|
---|
4 | from verbose import verbose
|
---|
5 |
|
---|
6 | #Start defining model parameters here
|
---|
7 | # Set up bed topography and ice geometry for a tilted 500m thick slab
|
---|
8 | md.geometry.base = .02 * md.mesh.x
|
---|
9 | md.geometry.bed = md.geometry.base
|
---|
10 | md.geometry.surface = .02 * md.mesh.x + 500
|
---|
11 | md.geometry.thickness = md.geometry.surface - md.geometry.bed
|
---|
12 |
|
---|
13 | # Define ice sliding velocity (m / yr)
|
---|
14 | md.initialization.vx = 1.0e-6 * md.constants.yts * np.ones(md.mesh.numberofvertices)
|
---|
15 | md.initialization.vy = np.zeros(md.mesh.numberofvertices)
|
---|
16 | md.initialization.pressure = np.zeros(md.mesh.numberofvertices)
|
---|
17 |
|
---|
18 | # Materials
|
---|
19 | # Ice flow law parameter (note that the standard parameter A = B^(- 3))
|
---|
20 | md.materials.rheology_B = 5e-25**(-1 / 3) * np.ones(md.mesh.numberofvertices)
|
---|
21 | md.initialization.temperature = 273. * np.ones(md.mesh.numberofvertices)
|
---|
22 | md.materials.rheology_n = 3. * np.ones(md.mesh.numberofelements)
|
---|
23 |
|
---|
24 | #Calving
|
---|
25 | md.calving.calvingrate = np.zeros(md.mesh.numberofvertices)
|
---|
26 | #md.calving.spclevelset = NaN(md.mesh.numberofvertices, 1)
|
---|
27 |
|
---|
28 | # Friction law and coefficient
|
---|
29 | md.friction = frictionshakti(md)
|
---|
30 | md.friction.coefficient = 20. * np.ones(md.mesh.numberofvertices)
|
---|
31 |
|
---|
32 | #Numerical parameters
|
---|
33 | #md.stressbalance.viscosity_overshoot = 0.0
|
---|
34 | md.masstransport.stabilization = 1.
|
---|
35 | md.thermal.stabilization = 1.
|
---|
36 | md.verbose = verbose(0)
|
---|
37 | md.settings.waitonlock = 30
|
---|
38 | md.stressbalance.restol = 0.05
|
---|
39 | md.steadystate.reltol = 0.05
|
---|
40 | md.stressbalance.reltol = 0.05
|
---|
41 | md.stressbalance.abstol = np.nan
|
---|
42 | md.timestepping.time_step = 1.
|
---|
43 | md.timestepping.final_time = 3.
|
---|
44 |
|
---|
45 | #GIA:
|
---|
46 | md.gia.lithosphere_thickness = 100. * np.ones(md.mesh.numberofvertices) # in km
|
---|
47 | md.gia.mantle_viscosity = 1.0e21 * np.ones(md.mesh.numberofvertices) # in Pa.s
|
---|
48 | md.materials.lithosphere_shear_modulus = 6.7e10 # in Pa
|
---|
49 | md.materials.lithosphere_density = 3.32 # in g / cm^ - 3
|
---|
50 | md.materials.mantle_shear_modulus = 1.45e11 # in Pa
|
---|
51 | md.materials.mantle_density = 3.34 # in g / cm^ - 3
|
---|
52 |
|
---|
53 | #Boundary conditions:
|
---|
54 | md = SetIceSheetBC(md)
|
---|
55 |
|
---|
56 | #Change name so that no test have the same name
|
---|
57 | md.private.runtimename = True
|
---|