texactTestP.cc - 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
---
texactTestP.cc (9005B)
---
1 /*
2 Copyright (C) 2012-2017 Ed Bueler and Constantine Khroulev
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 #include <cassert>
22 #include <cstdlib>
23 #include <cmath>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_matrix.h>
26
27 #include <gsl/gsl_version.h>
28 #if (defined GSL_MAJOR_VERSION) && (defined GSL_MINOR_VERSION) && \
29 ((GSL_MAJOR_VERSION >= 1 && GSL_MINOR_VERSION >= 15) || (GSL_MAJOR_VERSION >= 2))
30 #define PISM_USE_ODEIV2 1
31 #include <gsl/gsl_odeiv2.h>
32 #endif
33
34 #include "exactTestP.hh"
35
36 namespace pism {
37
38 static const double SperA = 31556926.0; // seconds per year; 365.2422 days
39 static const double g = 9.81; // m s-2
40 static const double rhoi = 910.0; // kg m-3
41 static const double rhow = 1000.0; // kg m-3
42
43 // major model parameters:
44 static const double Aglen = 3.1689e-24; // Pa-3 s-1
45 static const double k = (0.01 / (rhow * g)); // FIXME: this is extremely low but it matches what we were using
46 static const double Wr = 1.0; // m
47 static const double c1 = 0.500; // m-1
48 static const double c2 = 0.040; // [pure]
49
50 // specific to exact solution
51 static const double m0 = ((0.20 / SperA) * rhow); // kg m-2 s-1; = 20 cm year-1
52 static const double h0 = 500.0; // m
53 static const double v0 = (100.0 / SperA); // m s-1
54 static const double R1 = 5000.0; // m
55
56 int getsb(double r, double *sb, double *dsbdr) {
57 double CC;
58 if (r < R1) {
59 *sb = 0.0;
60 *dsbdr = 0.0;
61 } else {
62 CC = pow((c1 * v0) / (c2 * Aglen * pow((TESTP_L - R1), 5.0)), (1.0 / 3.0));
63 *sb = CC * pow(r - R1, (5.0 / 3.0));
64 *dsbdr = (5.0 / 3.0) * CC * pow(r - R1, (2.0 / 3.0));
65 }
66 return 0;
67 }
68
69
70 double criticalW(double r) {
71 double
72 h = h0 * (1.0 - (r / TESTP_R0) * (r / TESTP_R0)),
73 Po = rhoi * g * h;
74 double sb, dsb;
75 getsb(r, &sb, &dsb);
76
77 double sbcube = sb * sb * sb;
78 double Pocube = Po * Po * Po;
79
80 return (sbcube / (sbcube + Pocube)) * Wr;
81 }
82
83
84 int funcP(double r, const double W[], double f[], void *params) {
85 /* Computes RHS f(r,W) for differential equation as given in dampnotes.pdf
86 at https://github.com/bueler/hydrolakes:
87 dW
88 -- = f(r,W)
89 dr
90 Compare doublediff.m. Assumes Glen power n=3.
91 */
92
93 double sb, dsb, dPo, tmp1, omega0, numer, denom;
94
95 (void)params; /* quash warning "unused parameters" */
96
97 if (r < 0.0) {
98 f[0] = 0.0; /* place-holder */
99 return TESTP_R_NEGATIVE;
100 } else if (r > TESTP_L) {
101 f[0] = 0.0;
102 return GSL_SUCCESS;
103 } else {
104 getsb(r, &sb, &dsb);
105 omega0 = m0 / (2.0 * rhow * k);
106 dPo = -(2.0 * rhoi * g * h0 / (TESTP_R0 * TESTP_R0)) * r;
107 tmp1 = pow(W[0], 4.0 / 3.0) * pow(Wr - W[0], 2.0 / 3.0);
108 numer = dsb * W[0] * (Wr - W[0]);
109 numer = numer - (omega0 * r / W[0] + dPo) * tmp1;
110 denom = (1.0 / 3.0) * sb * Wr + rhow * g * tmp1;
111 f[0] = numer / denom;
112 return GSL_SUCCESS;
113 }
114 }
115
116
117 /* Computes initial condition W(r=L) = W_c(L^-). */
118 double initialconditionW() {
119 double hL, PoL, sbL;
120 hL = h0 * (1.0 - (TESTP_L / TESTP_R0) * (TESTP_L / TESTP_R0));
121 PoL = rhoi * g * hL;
122 sbL = pow(c1 * v0 / (c2 * Aglen), 1.0 / 3.0);
123 return (pow(sbL, 3.0) / (pow(sbL, 3.0) + pow(PoL, 3.0))) * Wr;
124 }
125
126
127 double psteady(double W, double magvb, double Po) {
128 double sbcube, frac, P;
129 sbcube = c1 * fabs(magvb) / (c2 * Aglen);
130 frac = (W < Wr) ? (Wr - W) / W : 0.0;
131 P = Po - pow(sbcube * frac, 1.0 / 3.0);
132 if (P < 0.0) {
133 P = 0.0;
134 }
135 return P;
136 }
137
138 #ifdef PISM_USE_ODEIV2
139
140 /* Solves ODE for W(r), the exact solution. Input r[] and output W[] must be
141 allocated vectors of length N. Input r[] must be decreasing. The combination
142 EPS_ABS = 1e-12, EPS_REL=0.0, method = RK Dormand-Prince O(8)/O(9)
143 is believed for now to be predictable and accurate. Note hstart is negative
144 so that the ODE solver does negative steps. Assumes
145 0 <= r[N-1] <= r[N-2] <= ... <= r[1] <= r[0] <= L. */
146 int getW(const double *r, int N, double *W, double EPS_ABS, double EPS_REL, int ode_method) {
147 int count;
148 int status = TESTP_NOT_DONE;
149 double rr, hstart;
150 const gsl_odeiv2_step_type *Tpossible[4];
151 const gsl_odeiv2_step_type *T;
152 gsl_odeiv2_system sys = { funcP, NULL, 1, NULL }; /* Jac-free method and no params */
153 gsl_odeiv2_driver *d;
154
155 /* setup for GSL ODE solver; these choices don't need Jacobian */
156 Tpossible[0] = gsl_odeiv2_step_rk8pd;
157 Tpossible[1] = gsl_odeiv2_step_rk2;
158 Tpossible[2] = gsl_odeiv2_step_rkf45;
159 Tpossible[3] = gsl_odeiv2_step_rkck;
160 if ((ode_method > 0) && (ode_method < 5)) {
161 T = Tpossible[ode_method - 1];
162 } else {
163 return TESTP_INVALID_METHOD;
164 }
165
166 hstart = -1000.0;
167 d = gsl_odeiv2_driver_alloc_y_new(&sys, T, hstart, EPS_ABS, EPS_REL);
168
169 /* initial conditions: (r,W) = (L,W_c(L^-)); r decreases from L toward 0 */
170 rr = TESTP_L;
171 for (count = 0; count < N; count++) {
172 /* except at start, use value at end of last interval as initial value for subinterval */
173 W[count] = (count == 0) ? initialconditionW() : W[count - 1];
174 while (rr > r[count]) {
175 status = gsl_odeiv2_driver_apply(d, &rr, r[count], &(W[count]));
176 if (status != GSL_SUCCESS) {
177 break;
178 }
179 if (W[count] > Wr) {
180 return TESTP_W_EXCEEDS_WR;
181 } else if (W[count] < criticalW(r[count])) {
182 return TESTP_W_BELOW_WCRIT;
183 }
184 }
185 }
186
187 gsl_odeiv2_driver_free(d);
188 return status;
189 }
190
191 #else
192 int getW(const double *r, int N, double *W,
193 double EPS_ABS, double EPS_REL, int ode_method) {
194 (void) r;
195 (void) EPS_ABS;
196 (void) EPS_REL;
197 (void) ode_method;
198
199 for (int j = 0; j < N; ++j) {
200 W[j] = 0;
201 }
202 return TESTP_OLD_GSL;
203 }
204 #endif
205
206 int exactP_list(const double *r, int N, double *h, double *magvb, double *Wcrit, double *W, double *P,
207 double EPS_ABS, double EPS_REL, int ode_method) {
208
209 int i, M, status;
210 /* check first: we have a list, r is decreasing, r is in range [0,L) */
211 if (N < 1) {
212 return TESTP_NO_LIST;
213 }
214
215 for (i = 0; i < N; i++) {
216 if ((i > 0) && (r[i] > r[i - 1])) {
217 return TESTP_LIST_NOT_DECREASING;
218 }
219 if (r[i] < 0.0) {
220 return TESTP_R_NEGATIVE;
221 }
222 }
223
224 M = 0;
225 while (r[M] > TESTP_L) {
226 h[M] = 0.0;
227 magvb[M] = 0.0;
228 Wcrit[M] = 0.0;
229 W[M] = 0.0;
230 P[M] = 0.0;
231 M++;
232 }
233
234 for (i = M; i < N; i++) {
235 h[i] = h0 * (1.0 - (r[i] / TESTP_R0) * (r[i] / TESTP_R0));
236
237 if (r[i] > R1) {
238 magvb[i] = v0 * pow((r[i] - R1) / (TESTP_L - R1), 5.0);
239 } else {
240 magvb[i] = 0.0;
241 }
242
243 Wcrit[i] = criticalW(r[i]);
244 }
245
246 status = getW(&(r[M]), N - M, &(W[M]), EPS_ABS, EPS_REL, ode_method);
247
248 if (status) {
249 for (i = M; i < N; i++) {
250 P[i] = 0.0;
251 }
252 return status;
253 } else {
254 for (i = M; i < N; i++) {
255 P[i] = psteady(W[i], magvb[i], rhoi * g * h[i]);
256 }
257 return 0;
258 }
259 }
260
261 TestPParameters exactP(const std::vector<double> &r,
262 double EPS_ABS, double EPS_REL, int ode_method) {
263 TestPParameters result(r.size());
264 result.r = r;
265
266 result.error_code = exactP_list(&r[0], r.size(),
267 &result.h[0],
268 &result.magvb[0],
269 &result.Wcrit[0],
270 &result.W[0],
271 &result.P[0],
272 EPS_ABS, EPS_REL, ode_method);
273
274 switch (result.error_code) {
275 case 0:
276 result.error_message = "success";
277 break;
278 case TESTP_R_NEGATIVE:
279 result.error_message = "error: r < 0";
280 break;
281 case TESTP_W_EXCEEDS_WR:
282 result.error_message = "error: W > W_r";
283 break;
284 case TESTP_W_BELOW_WCRIT:
285 result.error_message = "error: W < W_crit";
286 break;
287 case TESTP_INVALID_METHOD:
288 result.error_message = "error: invalid choice for ODE method";
289 break;
290 case TESTP_NOT_DONE:
291 result.error_message = "error: ODE integrator not done";
292 break;
293 case TESTP_NO_LIST:
294 result.error_message = "error: no list of r values at input to exactP_list()";
295 break;
296 case TESTP_LIST_NOT_DECREASING:
297 result.error_message = "error: input list of r values to exactP_list() is not decreasing";
298 break;
299 default:
300 result.error_message = "unknown error code";
301 }
302
303 return result;
304 }
305
306 } // end of namespace pism