PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
exactTestP.cc
Go to the documentation of this file.
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^-). */
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])) {
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;
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;
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
#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:56
static const double rhoi
Definition: exactTestP.cc:40
static const double rhow
Definition: exactTestP.cc:41
double initialconditionW()
Definition: exactTestP.cc:118
static const double g
Definition: exactTestP.cc:39
int funcP(double r, const double W[], double f[], void *params)
Definition: exactTestP.cc:84
static const double c2
Definition: exactTestP.cc:48
static const double SperA
Definition: exactTestP.cc:38
double criticalW(double r)
Definition: exactTestP.cc:70
static const double v0
Definition: exactTestP.cc:53
static const double k
Definition: exactTestP.cc:45
static const double Aglen
Definition: exactTestP.cc:44
static const double Wr
Definition: exactTestP.cc:46
TestPParameters exactP(const std::vector< double > &r, double EPS_ABS, double EPS_REL, int ode_method)
Definition: exactTestP.cc:261
static const double c1
Definition: exactTestP.cc:47
static const double h0
Definition: exactTestP.cc:52
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)
Definition: exactTestP.cc:206
static const double R1
Definition: exactTestP.cc:54
static const double m0
Definition: exactTestP.cc:51
int getW(const double *r, int N, double *W, double EPS_ABS, double EPS_REL, int ode_method)
Definition: exactTestP.cc:192
double psteady(double W, double magvb, double Po)
Definition: exactTestP.cc:127
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