texactM.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
texactM.py (2776B)
---
1 #!/usr/bin/env python3
2
3 # Computes and plots exact solution "test M", in preparation for implementing C
4 # version for PISM verification. Preliminary tests are possible with gridded
5 # version of this output, even before a C implementation.
6
7 from pylab import *
8 from scipy.optimize import fsolve
9 from scipy.integrate import odeint
10
11 SperA = 31556926.0
12 g = 9.81
13 rho = 910.0 # density of ice; kg/m^3
14 rhow = 1028.0 # density of ocean water; kg/m^3
15 barB = 3.7e8 # strength of shelf; Pa s^(1/3); as in Schoof 2006; compare 1.9e8
16 # from MacAyeal et al 1996
17 Rc = 600.0e3 # calving front at 600 km
18 Rg = 300.0e3 # grounding line at 300 km
19 H0 = 500.0 # uniform 500 m thickness
20 ug = 100.0 / SperA # velocity across grounding line is 100 m/year
21
22 # compute physical constant in ODE
23 Q = (1.0 - rho / rhow) * rho * g * Rc * H0 / (2.0 * barB)
24 print('physical constant: Q = %f' % Q)
25
26 # solve FF(x)=0 to get strain rate alpha'; x is RHS of ODE
27
28
29 def FF(x, alpha, r):
30 DD = x * x + x * (alpha / r) + (alpha / r) ** 2
31 return Q * DD ** (1. / 3.) - 2.0 * r * x - alpha
32
33
34 # ODE we are solving is alpha' = GG(alpha,r)
35 qual = [] # quality information goes here
36
37
38 def GG(alpha, r):
39 # heuristic: guess is about 1/7 th of solution to a nearby problem
40 guess = 0.15 * ((Q / r) ** 3 - alpha[0] / r)
41 result = fsolve(FF, guess, args=(alpha[0], r))
42 qual.append([r, guess, result, abs(guess - result) / abs(result)])
43 return [result]
44
45
46 # build vector of locations to solve (and plot)
47 dr = 1000.0 # plot on 1 km grid
48 r = linspace(Rg, Rc, ((Rc - Rg) / dr) + 1)
49
50 # solve:
51 print('solving with odeint from scipy.integrate; it reports: ', end=' ')
52 alpha = odeint(GG, [ug], r, printmessg=1,
53 rtol=0.0, # ask for no change in digits
54 atol=0.00001 / SperA) # ask for abs tol of 0.01 mm/year
55
56 # info to evaluate solution
57 qual = array(qual)
58 rused = qual[:, 0]
59 strainrate = qual[:, 2]
60 guesserr = qual[:, 3]
61 print('maximum relative error in initial guess for fsolve() is %f'
62 % max(guesserr))
63
64 print("at calving front: alpha(Rc) = %f (m/year)"
65 % (alpha[-1] * SperA))
66 # print " (last r used = %f km; strain rate at last r = %f (1/a))" \
67 # % (rused[-1] / 1000.0, strainrate[-1] * SperA)
68
69 # plot velocity solution alpha(r), and also strain rate alpha'(r)
70 figure(1)
71 subplot(211)
72 plot(r / 1000.0, alpha * SperA, 'k', linewidth=3)
73 ylabel(r'velocity (m/year)', size=14)
74 subplot(212)
75 plot(rused / 1000.0, strainrate * SperA, 'ko-')
76 axis([Rg / 1000.0, Rc / 1000.0, 0, 1.1 * max(strainrate * SperA)])
77 xlabel(r'r (km)', size=14)
78 ylabel(r'strain rate (1/a)', size=14)
79 # print "saving figure 'combinedM.png'" # optional: save as PNG
80 #savefig('combinedM.png', dpi=300, facecolor='w', edgecolor='w')
81 print('close figure to end')
82 show()