-
Notifications
You must be signed in to change notification settings - Fork 0
/
ssaErebus.py
92 lines (79 loc) · 2.7 KB
/
ssaErebus.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
"""
Simplified script for running and exploring the Fiber-Bundle enabled shallow
shelf model.
Recommended use is to run this script in an ipython session or notebook:
>> %run ssaRun.py
Can plot the calving events
>> plt.plot(ssaModel.obslist[0].ts, ssaModel.obslist[0].xs)
"""
# Import the models
from ssaModel import *
# Setup some fenics log stuff to output diagnostic information
set_log_level(00)
set_log_active(False)
# Set simulation parameters
DT = 8640000
Nx = 2**8 # Number of points in the x-direction
Lx = 200e3/1 # Length of domain in the x-direction
DX = Lx/Nx
NT = 1000
# Set ice shelf parameters
accum = -2./time_factor # m/s
H0 = 434. # Ice thickness at the grounding line (m)
U0 = 95./time_factor # Velocity of ice at the grounding line (m/a)
B = (2.54e-17/time_factor)**(-1./3.)
Lx = -H0*U0/accum - 1e-6
# Set ice shel parameters
#accum = 0.5/time_factor
#H0 = 500.
#U0 = 50./time_factor
#B = 0.5e8
print('CFL cond numb: {}'.format(U0*DT/DX))
# Initialize model
mesh = IntervalMesh(Nx, 0.0, Lx)
fbmkwargs={'Lx':Lx,
'N0':1000,
'Nf':10,
'xsep':200,
'fbm_type':'full'}
# In the large-Nf limit, the maximum stress per fiber of a bundle of fibers
# with strictly increasing strengths (with maximum xmax) is Fmax/Nf = xmax^2/4.
# For Fmax=2.8, we need xmax=sqrt(4*Fmax/Nf)=1.058
# Accounting for fluctuations due to finite N (Fmax_actual = Fmax + d N**1/3),
# with d~0.3138 for the strict distribution over [0,1], get xmax=0.923
#fbmkwargs['dist'] = strict_dist(0,0.923)
fbmkwargs['dist'] = strict_dist(0,0.179)
#fbmkwargs['dist'] = uni_dist(0,0.179)
results = []
for i in range(1):
mesh = IntervalMesh(Nx, 0.0, Lx)
ssaModel = ssa1D(mesh,order=1,U0=U0,H0=H0,B=B,
advect_front=True, calve_flag=True,
fbm_type='full', fbmkwargs=fbmkwargs,
Lmax=-H0*U0/accum) ;
del mesh
x,H,U = ssaModel.steady_state(accum)
H,U=ssaModel.init_shelf(accum)
ssaModel.H = H
ssaModel.U = U
# Run the model in time
#H,U = ssaModel.integrate(H,U,dt=DT,Nt=NT,accum=accum);
results.append(ssaModel.frontobs.data)
#plt.plot(*ssaModel.fbmobs.data, marker='o')
#plt.xlabel('Position of particles (m)')
#plt.ylabel('Damage of particles (dimless)')
#plt.show()
#
#for entry in results:
# plt.plot(entry[0]/time_factor, entry[1])
#plt.xlabel('Time (yrs)')
#plt.ylabel('Front pos (m)')
#plt.show()
# get latest damage of fibers with ssaModel.fbmobs.data
# get front position over time with ssaModel.frontobs.data
# get calving events with ssaModel.calveobs.data
#fbmkwargs={'Lx':Lx,
# 'N0':1000,
# 'xsep':200,
# 'fbm_type':'max',
# 'dist':retconst(2.8)}