Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
exactTestP.cc
Go to the documentation of this file.
1/*
2 Copyright (C) 2012-2017, 2023 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 "pism/verification/tests/exactTestP.hh"
22#include <gsl/gsl_errno.h> // for GSL_SUCCESS
23#include <cmath> // for pow, fabs
24#include <cstdlib> // for NULL
25
26#include <gsl/gsl_version.h>
27#if (defined GSL_MAJOR_VERSION) && (defined GSL_MINOR_VERSION) && \
28 ((GSL_MAJOR_VERSION >= 1 && GSL_MINOR_VERSION >= 15) || (GSL_MAJOR_VERSION >= 2))
29#define PISM_USE_ODEIV2 1
30#include <gsl/gsl_odeiv2.h>
31#endif
32
33namespace pism {
34
35static const double SperA = 31556926.0; // seconds per year; 365.2422 days
36static const double g = 9.81; // m s-2
37static const double rhoi = 910.0; // kg m-3
38static const double rhow = 1000.0; // kg m-3
39
40// major model parameters:
41static const double Aglen = 3.1689e-24; // Pa-3 s-1
42static const double k = (0.01 / (rhow * g)); // FIXME: this is extremely low but it matches what we were using
43static const double Wr = 1.0; // m
44static const double c1 = 0.500; // m-1
45static const double c2 = 0.040; // [pure]
46
47// specific to exact solution
48static const double m0 = ((0.20 / SperA) * rhow); // kg m-2 s-1; = 20 cm year-1
49static const double h0 = 500.0; // m
50static const double v0 = (100.0 / SperA); // m s-1
51static const double R1 = 5000.0; // m
52
53int getsb(double r, double *sb, double *dsbdr) {
54 double CC;
55 if (r < R1) {
56 *sb = 0.0;
57 *dsbdr = 0.0;
58 } else {
59 CC = pow((c1 * v0) / (c2 * Aglen * pow((TESTP_L - R1), 5.0)), (1.0 / 3.0));
60 *sb = CC * pow(r - R1, (5.0 / 3.0));
61 *dsbdr = (5.0 / 3.0) * CC * pow(r - R1, (2.0 / 3.0));
62 }
63 return 0;
64}
65
66
67double criticalW(double r) {
68 double
69 h = h0 * (1.0 - (r / TESTP_R0) * (r / TESTP_R0)),
70 Po = rhoi * g * h;
71 double sb, dsb;
72 getsb(r, &sb, &dsb);
73
74 double sbcube = sb * sb * sb;
75 double Pocube = Po * Po * Po;
76
77 return (sbcube / (sbcube + Pocube)) * Wr;
78}
79
80
81int funcP(double r, const double W[], double f[], void *params) {
82 /* Computes RHS f(r,W) for differential equation as given in dampnotes.pdf
83 at https://github.com/bueler/hydrolakes:
84 dW
85 -- = f(r,W)
86 dr
87 Compare doublediff.m. Assumes Glen power n=3.
88 */
89
90 double sb, dsb, dPo, tmp1, omega0, numer, denom;
91
92 (void)params; /* quash warning "unused parameters" */
93
94 if (r < 0.0) {
95 f[0] = 0.0; /* place-holder */
96 return TESTP_R_NEGATIVE;
97 }
98
99 if (r > TESTP_L) {
100 f[0] = 0.0;
101 return GSL_SUCCESS;
102 }
103
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/* Computes initial condition W(r=L) = W_c(L^-). */
118 double hL, PoL, sbL;
119 hL = h0 * (1.0 - (TESTP_L / TESTP_R0) * (TESTP_L / TESTP_R0));
120 PoL = rhoi * g * hL;
121 sbL = pow(c1 * v0 / (c2 * Aglen), 1.0 / 3.0);
122 return (pow(sbL, 3.0) / (pow(sbL, 3.0) + pow(PoL, 3.0))) * Wr;
123}
124
125
126double psteady(double W, double magvb, double Po) {
127 double sbcube, frac, P;
128 sbcube = c1 * fabs(magvb) / (c2 * Aglen);
129 frac = (W < Wr) ? (Wr - W) / W : 0.0;
130 P = Po - pow(sbcube * frac, 1.0 / 3.0);
131 if (P < 0.0) {
132 P = 0.0;
133 }
134 return P;
135}
136
137#ifdef PISM_USE_ODEIV2
138
139/* Solves ODE for W(r), the exact solution. Input r[] and output W[] must be
140allocated vectors of length N. Input r[] must be decreasing. The combination
141EPS_ABS = 1e-12, EPS_REL=0.0, method = RK Dormand-Prince O(8)/O(9)
142is believed for now to be predictable and accurate. Note hstart is negative
143so that the ODE solver does negative steps. Assumes
144 0 <= r[N-1] <= r[N-2] <= ... <= r[1] <= r[0] <= L. */
145int getW(const double *r, int N, double *W, double EPS_ABS, double EPS_REL, int ode_method) {
146 int count;
147 int status = TESTP_NOT_DONE;
148 double rr, hstart;
149 const gsl_odeiv2_step_type *Tpossible[4];
150 const gsl_odeiv2_step_type *T;
151 gsl_odeiv2_system sys = { funcP, NULL, 1, NULL }; /* Jac-free method and no params */
152 gsl_odeiv2_driver *d;
153
154 /* setup for GSL ODE solver; these choices don't need Jacobian */
155 Tpossible[0] = gsl_odeiv2_step_rk8pd;
156 Tpossible[1] = gsl_odeiv2_step_rk2;
157 Tpossible[2] = gsl_odeiv2_step_rkf45;
158 Tpossible[3] = gsl_odeiv2_step_rkck;
159 if ((ode_method > 0) && (ode_method < 5)) {
160 T = Tpossible[ode_method - 1];
161 } else {
163 }
164
165 hstart = -1000.0;
166 d = gsl_odeiv2_driver_alloc_y_new(&sys, T, hstart, EPS_ABS, EPS_REL);
167
168 /* initial conditions: (r,W) = (L,W_c(L^-)); r decreases from L toward 0 */
169 rr = TESTP_L;
170 for (count = 0; count < N; count++) {
171 /* except at start, use value at end of last interval as initial value for subinterval */
172 W[count] = (count == 0) ? initialconditionW() : W[count - 1];
173 while (rr > r[count]) {
174 status = gsl_odeiv2_driver_apply(d, &rr, r[count], &(W[count]));
175 if (status != GSL_SUCCESS) {
176 break;
177 }
178 if (W[count] > Wr) {
179 return TESTP_W_EXCEEDS_WR;
180 }
181 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
192int 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
206int 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])) {
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 != 0) {
249 for (i = M; i < N; i++) {
250 P[i] = 0.0;
251 }
252 return status;
253 }
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
260TestPParameters exactP(const std::vector<double> &r,
261 double EPS_ABS, double EPS_REL, int ode_method) {
262 TestPParameters result(r.size());
263 result.r = r;
264
265 result.error_code = exactP_list(r.data(), r.size(),
266 (result.h).data(),
267 (result.magvb).data(),
268 (result.Wcrit).data(),
269 (result.W).data(),
270 (result.P).data(),
271 EPS_ABS, EPS_REL, ode_method);
272
273 switch (result.error_code) {
274 case 0:
275 result.error_message = "success";
276 break;
277 case TESTP_R_NEGATIVE:
278 result.error_message = "error: r < 0";
279 break;
281 result.error_message = "error: W > W_r";
282 break;
284 result.error_message = "error: W < W_crit";
285 break;
287 result.error_message = "error: invalid choice for ODE method";
288 break;
289 case TESTP_NOT_DONE:
290 result.error_message = "error: ODE integrator not done";
291 break;
292 case TESTP_NO_LIST:
293 result.error_message = "error: no list of r values at input to exactP_list()";
294 break;
296 result.error_message = "error: input list of r values to exactP_list() is not decreasing";
297 break;
298 default:
299 result.error_message = "unknown error code";
300 }
301
302 return result;
303}
304
305} // end of namespace pism
#define TESTP_NOT_DONE
Definition exactTestP.hh:52
#define TESTP_LIST_NOT_DECREASING
Definition exactTestP.hh:54
#define TESTP_NO_LIST
Definition exactTestP.hh:53
#define TESTP_OLD_GSL
Definition exactTestP.hh:55
#define TESTP_R_NEGATIVE
Definition exactTestP.hh:48
#define TESTP_INVALID_METHOD
Definition exactTestP.hh:51
#define TESTP_L
Definition exactTestP.hh:45
#define TESTP_W_EXCEEDS_WR
Definition exactTestP.hh:49
#define TESTP_R0
Definition exactTestP.hh:44
#define TESTP_W_BELOW_WCRIT
Definition exactTestP.hh:50
int getsb(double r, double *sb, double *dsbdr)
Definition exactTestP.cc:53
static const double rhoi
Definition exactTestP.cc:37
static const double rhow
Definition exactTestP.cc:38
double initialconditionW()
static const double g
Definition exactTestP.cc:36
int funcP(double r, const double W[], double f[], void *params)
Definition exactTestP.cc:81
static const double c2
Definition exactTestP.cc:45
static const double SperA
Definition exactTestP.cc:35
double criticalW(double r)
Definition exactTestP.cc:67
static const double v0
Definition exactTestP.cc:50
static const double k
Definition exactTestP.cc:42
static const double Aglen
Definition exactTestP.cc:41
static const double Wr
Definition exactTestP.cc:43
TestPParameters exactP(const std::vector< double > &r, double EPS_ABS, double EPS_REL, int ode_method)
static const double c1
Definition exactTestP.cc:44
static const double h0
Definition exactTestP.cc:49
int exactP_list(const double *r, int N, double *h, double *magvb, double *Wcrit, double *W, double *P, double EPS_ABS, double EPS_REL, int ode_method)
static const double R1
Definition exactTestP.cc:51
static const double m0
Definition exactTestP.cc:48
int getW(const double *r, int N, double *W, double EPS_ABS, double EPS_REL, int ode_method)
double psteady(double W, double magvb, double Po)
std::vector< double > h
Definition exactTestP.hh:65
std::vector< double > Wcrit
Definition exactTestP.hh:65
std::vector< double > magvb
Definition exactTestP.hh:65
std::vector< double > P
Definition exactTestP.hh:65
std::string error_message
Definition exactTestP.hh:64
std::vector< double > W
Definition exactTestP.hh:65
std::vector< double > r
Definition exactTestP.hh:65
int count
Definition test_cube.c:16