texactTestL.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
---
texactTestL.cc (7040B)
---
1 /* Copyright (C) 2016, 2017 PISM Authors
2 *
3 * This file is part of PISM.
4 *
5 * PISM is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 3 of the License, or (at your option) any later
8 * version.
9 *
10 * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 * details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with PISM; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20 #include "exactTestL.hh"
21
22 #include <cstdlib>
23 #include <cmath>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_matrix.h>
26 #include <gsl/gsl_math.h> /* M_PI */
27
28 #include <gsl/gsl_version.h>
29 #if (defined GSL_MAJOR_VERSION) && (defined GSL_MINOR_VERSION) && \
30 ((GSL_MAJOR_VERSION >= 1 && GSL_MINOR_VERSION >= 15) || (GSL_MAJOR_VERSION >= 2))
31 #define PISM_USE_ODEIV2 1
32 #include <gsl/gsl_odeiv2.h>
33 #endif
34
35 #include "pism/util/error_handling.hh"
36 using pism::RuntimeError;
37
38 static const double SperA = 31556926.0; /* seconds per year; 365.2422 days */
39
40 static const double L = 750.0e3; /* m; i.e. 750 km */
41 static const double b0 = 500.0; /* m */
42 static const double z0 = 1.2;
43 static const double g = 9.81;
44 static const double rho = 910.0;
45 static const double n = 3.0; /* Glen power */
46
47 int funcL(double r, const double u[], double f[], void *params) {
48 /*
49 RHS for differential equation:
50 du 5/8 / a_0 r (L^2 - r^2) \ 1/3
51 -- = - (8/3) b'(r) u - |---------------------|
52 dr \ 2 L^2 \tilde\Gamma /
53 */
54
55 const double Lsqr = L * L;
56 const double a0 = 0.3 / SperA; /* m/s; i.e. 0.3 m/a */
57 const double A = 1.0e-16 / SperA; /* = 3.17e-24 1/(Pa^3 s); EISMINT I flow law */
58 const double Gamma = 2 * pow(rho * g,n) * A / (n+2);
59 const double tilGamma = Gamma * pow(n,n) / (pow(2.0 * n + 2.0, n));
60 const double C = a0 / (2.0 * Lsqr * tilGamma);
61
62 if (params == NULL) {} /* quash warning "unused parameters" */
63
64 if ((r >= 0.0) && (r <= L)) {
65 const double freq = z0 * M_PI / L;
66 const double bprime = b0 * freq * sin(freq * r);
67 f[0] = - (8.0/3.0) * bprime * pow(u[0], 5.0/8.0)
68 - pow(C * r * (Lsqr - r * r), 1.0/3.0);
69 } else {
70 f[0] = 0.0; /* no changes outside of defined interval */
71 }
72 return GSL_SUCCESS;
73 }
74
75 /* exit status could be one of these; returned zero indicates success */
76 #define TESTL_NOT_DONE 8966
77 #define TESTL_NOT_DECREASING 8967
78 #define TESTL_INVALID_METHOD 8968
79 #define TESTL_NO_LIST 8969
80
81 #ifdef PISM_USE_ODEIV2
82
83 /* combination EPS_ABS = 1e-12, EPS_REL=0.0, method = 1 = RK Cash-Karp
84 is believed to be predictable and accurate */
85 int getU(const double *r, int N, double *u,
86 const double EPS_ABS, const double EPS_REL, const int ode_method) {
87 /* solves ODE for u(r)=H(r)^{8/3}, 0 <= r <= L, for test L
88 r and u must be allocated vectors of length N; r[] must be decreasing */
89 int k, count;
90 int status = TESTL_NOT_DONE;
91 double rr, hstart;
92 const gsl_odeiv2_step_type* Tpossible[4];
93 const gsl_odeiv2_step_type *T;
94 gsl_odeiv2_system sys = {funcL, NULL, 1, NULL}; /* Jac-free method and no params */
95 gsl_odeiv2_driver *d;
96
97 /* setup for GSL ODE solver; these choices don't need Jacobian */
98 Tpossible[0] = gsl_odeiv2_step_rk8pd;
99 Tpossible[1] = gsl_odeiv2_step_rk2;
100 Tpossible[2] = gsl_odeiv2_step_rkf45;
101 Tpossible[3] = gsl_odeiv2_step_rkck;
102 if ((ode_method > 0) && (ode_method < 5))
103 T = Tpossible[ode_method-1];
104 else {
105 return TESTL_INVALID_METHOD;
106 }
107
108 /* check first: we have a list, and r is decreasing */
109 if (N < 1) return TESTL_NO_LIST;
110 for (k = 1; k<N; k++) { if (r[k] > r[k-1]) return TESTL_NOT_DECREASING; }
111
112 /* outside of ice cap, u = 0 */
113 k = 0;
114 while (r[k] >= L) {
115 u[k] = 0.0;
116 k++;
117 if (k == N) return GSL_SUCCESS;
118 }
119
120 /* initialize GSL ODE solver */
121 hstart = -10000.0;
122 d = gsl_odeiv2_driver_alloc_y_new(&sys, T, hstart, EPS_ABS, EPS_REL);
123
124 /* initial conditions: (r,u) = (L,0); r decreases from L */
125 rr = L;
126 for (count = k; count < N; count++) {
127 /* except at start, use value at end of last interval as initial guess */
128 u[count] = (count == 0) ? 0.0 : u[count-1];
129 while (rr > r[count]) {
130 status = gsl_odeiv2_driver_apply(d, &rr, r[count], &(u[count]));
131 if (status != GSL_SUCCESS){
132 break;
133 }
134 }
135 }
136
137 gsl_odeiv2_driver_free(d);
138 return status;
139 }
140
141 #else /* old GSL (< 1.15) */
142 int getU(const double *r, int N, double *u,
143 const double EPS_ABS, const double EPS_REL, const int ode_method) {
144 (void) r;
145 (void) N;
146 (void) u;
147 (void) EPS_ABS;
148 (void) EPS_REL;
149 (void) ode_method;
150 throw RuntimeError(PISM_ERROR_LOCATION, "Test L requires GSL version 1.15 or later.");
151 return 0;
152 }
153 #endif
154
155 int exactL_list(const double *r, int N, double *H, double *b, double *a) {
156 /* N values r[0] > r[1] > ... > r[N-1] (decreasing)
157 assumes r, H, b, a are allocated length N arrays */
158
159 const double Lsqr = L * L;
160 const double a0 = 0.3 / SperA; /* m/s; i.e. 0.3 m/a */
161 double *u;
162 int stat, i;
163
164 u = (double *) malloc((size_t)N * sizeof(double)); /* temporary arrays */
165
166 /* combination EPS_ABS = 1e-12, EPS_REL=0.0, method = 1 = RK Cash-Karp
167 believed to be predictable and accurate */
168 stat = getU(r,N,u,1.0e-12,0.0,1);
169 if (stat != GSL_SUCCESS) {
170 free(u);
171 return stat;
172 }
173
174 for (i = 0; i < N; i++) {
175 H[i] = pow(u[i],3.0/8.0);
176 b[i] = - b0 * cos(z0 * M_PI * r[i] / L);
177 a[i] = a0 * (1.0 - (2.0 * r[i] * r[i] / Lsqr));
178 }
179
180 free(u);
181 return 0;
182 }
183
184 ExactLParameters::ExactLParameters(size_t size)
185 : H(size), a(size), b(size) {
186 }
187
188 ExactLParameters exactL(const std::vector<double> &r) {
189 ExactLParameters result(r.size());
190
191 int error_code = exactL_list(&r[0], r.size(), &result.H[0], &result.b[0], &result.a[0]);
192
193 switch (error_code) {
194 case TESTL_NOT_DONE:
195 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Test L ERROR: exactL_list() returns 'NOT_DONE' (%d) ...",
196 error_code);
197 break;
198 case TESTL_NOT_DECREASING:
199 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Test L ERROR: exactL_list() returns 'NOT_DECREASING' (%d) ...",
200 error_code);
201 break;
202 case TESTL_INVALID_METHOD:
203 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Test L ERROR: exactL_list() returns 'INVALID_METHOD' (%d) ...",
204 error_code);
205 break;
206 case TESTL_NO_LIST:
207 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Test L ERROR: exactL_list() returns 'NO_LIST' (%d) ...",
208 error_code);
209 break;
210 default:
211 break;
212 }
213
214 return result;
215 }