texactTestM.c - 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
---
texactTestM.c (6330B)
---
1 /*
2 Copyright (C) 2008, 2014, 2016 Ed Bueler
3
4 This file is part of PISM.
5
6 PISM is free software; you can redistribute it and/or modify it under the
7 terms of the GNU General Public License as published by the Free Software
8 Foundation; either version 3 of the License, or (at your option) any later
9 version.
10
11 PISM is distributed in the hope that it will be useful, but WITHOUT ANY
12 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
14 details.
15
16 You should have received a copy of the GNU General Public License
17 along with PISM; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <math.h>
25 #include <gsl/gsl_errno.h>
26 #include <gsl/gsl_matrix.h>
27 #include <gsl/gsl_odeiv.h>
28 #include "exactTestM.h"
29
30 #define MAX(a, b) (a > b ? a : b)
31 #define MIN(a, b) (a > b ? b : a)
32
33 #define SperA 31556926.0 /* seconds per year; 365.2422 days */
34 #define g 9.81
35 #define rho 910.0 /* ice density; kg/m^3 */
36 #define rhow 1028.0 /* sea water density; kg/m^3 */
37 #define n 3.0 /* Glen power */
38 #define barB 3.7e8 /* strength of shelf; Pa s^(1/3); as in Schoof 2006;
39 compare 1.9e8 from MacAyeal et al 1996 */
40 #define H0 500.0 /* m */
41 #define Rg 300.0e3 /* m; 300 km */
42 #define Rc 600.0e3 /* m; 600 km */
43
44
45 double F_M(double x, double alpha, double r, double Q) {
46 const double
47 aor = alpha / r,
48 DD = x * x + x * aor + pow(aor, 2.0);
49 return Q * pow(DD, 1./3.) - 2.0 * r * x - alpha;
50 }
51
52
53 double dF_M(double x, double alpha, double r, double Q) {
54 const double
55 aor = alpha / r,
56 DD = x * x + x * aor + pow(aor, 2.0);
57 return (1. / 3.) * Q * pow(DD, - 2./3.) * (2.0 * x + aor) - 2.0 * r;
58 }
59
60
61 int funcM_ode_G(double r, const double alpha[], double f[], void* params) {
62 (void) params;
63 /* RHS G for differential equation:
64 alpha' = G(alpha,r)
65 but where we solve this equation to find alpha':
66 F(alpha',alpha,r) = 0
67 heuristic: guess is about 1/7 th of solution to a nearby problem;
68 no range checking on r, so use away from zero */
69
70 const double Q = (1.0 - rho / rhow) * rho * g * Rc * H0 / (2.0 * barB),
71 guess = 0.15 * (pow(Q/r,n) - alpha[0]/r);
72 /* in Python (exactM.py): f[0] = fsolve(F_M,guess,args=(alpha[0],r));
73 we could call GSL to find root, but hand-coding Newton's is easier */
74 double Old = guess, New; /* capitalized to avoid the C++ keyword name
75 clash */
76 int i;
77 for (i = 1; i < 100; i++) {
78 New = Old - F_M(Old,alpha[0],r,Q) / dF_M(Old,alpha[0],r,Q);
79 if (fabs((New-Old)/Old) < 1.0e-12) break;
80 Old = New;
81 }
82 if (i >= 90)
83 printf("exactTestM WARNING: Newton iteration not converged in funcM_ode_G!\n");
84 f[0] = New;
85 return GSL_SUCCESS;
86 }
87
88
89 #define NOT_DONE 8966
90 #define INVALID_METHOD 8968
91 #define NEGATIVE_R 8969
92
93 /* combination EPS_ABS = 1e-12, EPS_REL=0.0, method = 1 = RK Cash-Karp
94 is believed to be predictable and accurate; returns GSL_SUCCESS=0 if success */
95 int exactM_old(double r,
96 double *alpha, double *Drr,
97 const double EPS_ABS, const double EPS_REL, const int ode_method) {
98
99 double ug = 100.0 / SperA; /* velocity across grounding line is 100 m/a */
100 double DrrRg, xx, xA, nu, aa, rr, myalf, step;
101 const gsl_odeiv_step_type* T;
102 int status = NOT_DONE;
103 gsl_odeiv_step* s;
104 gsl_odeiv_control* c;
105 gsl_odeiv_evolve* e;
106 gsl_odeiv_system sys = {funcM_ode_G, NULL, 1, NULL}; /* Jac-free method and no params */
107
108 if (r < 0) {
109 return NEGATIVE_R; /* only nonnegative radial coord allowed */
110 } else if (r <= Rg/4.0) {
111 *alpha = 0.0; /* zero velocity near center */
112 *Drr = 0.0;
113 return GSL_SUCCESS;
114 } else if (r <= Rg) {
115 /* power law from alpha=0 to alpha=ug in Rg/4 < r <= Rg;
116 f(r) w: f(Rg/4)=f'(Rg/4)=0 and f(Rg)=ug and f(Rg) = DrrRg */
117 funcM_ode_G(Rg, &ug, &DrrRg, NULL); /* first get Drr = alpha' at Rg where alpha=ug */
118 /* printf("DrrRg=%e (1/a)\n",DrrRg*SperA); */
119 xx = r - 0.25 * Rg;
120 xA = 0.75 * Rg;
121 nu = DrrRg * xA / ug;
122 aa = ug / pow(xA, nu);
123 /* printf("power nu=%e\n",nu); */
124 *alpha = aa * pow(xx, nu);
125 *Drr = aa * nu * pow(xx, nu - 1);
126 return GSL_SUCCESS;
127 } else if (r >= Rc + 1.0) {
128 *alpha = 0.0; /* zero velocity beyond calving front */
129 *Drr = 0.0;
130 return GSL_SUCCESS;
131 }
132
133 /* need to solve ODE to find alpha, so setup for GSL ODE solver */
134 switch (ode_method) {
135 case 1:
136 T = gsl_odeiv_step_rkck; /* RK Cash-Karp */
137 break;
138 case 2:
139 T = gsl_odeiv_step_rk2;
140 break;
141 case 3:
142 T = gsl_odeiv_step_rk4;
143 break;
144 case 4:
145 T = gsl_odeiv_step_rk8pd;
146 break;
147 default:
148 printf("INVALID ode_method in exactM(): must be 1,2,3,4\n");
149 return INVALID_METHOD;
150 }
151 s = gsl_odeiv_step_alloc(T, (size_t)1); /* one scalar ode */
152 c = gsl_odeiv_control_y_new(EPS_ABS,EPS_REL);
153 e = gsl_odeiv_evolve_alloc((size_t)1); /* one scalar ode */
154
155 /* initial conditions: (r,alf) = (Rg,ug); r increases */
156 rr = Rg;
157 myalf = ug;
158 /* printf (" r (km) alpha (m/a)\n");
159 printf (" %11.5e %11.5e\n", rr/1000.0, myalf * SperA); */
160 while (rr < r) {
161 /* step = r - rr; try to get to solution in one step; trust stepping algorithm */
162 step = MIN(r-rr,20.0e3);
163 status = gsl_odeiv_evolve_apply(e, c, s, &sys, &rr, r, &step, &myalf);
164 if (status != GSL_SUCCESS) break;
165 /* printf (" %11.5e %11.5e\n", rr/1000.0, myalf * SperA); */
166 }
167
168 gsl_odeiv_evolve_free(e);
169 gsl_odeiv_control_free(c);
170 gsl_odeiv_step_free(s);
171
172 *alpha = myalf;
173 funcM_ode_G(r, alpha, Drr, NULL);
174 return status;
175 }
176
177 struct TestMParameters exactM(double r,
178 double EPS_ABS, double EPS_REL, int ode_method) {
179 struct TestMParameters result;
180 result.error_code = exactM_old(r, &result.alpha, &result.Drr,
181 EPS_ABS, EPS_REL, ode_method);
182 return result;
183 }