PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
exactTestM.c
Go to the documentation of this file.
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 }
#define g
Definition: exactTestM.c:34
#define H0
Definition: exactTestM.c:39
#define Rc
Definition: exactTestM.c:41
#define SperA
Definition: exactTestM.c:33
#define INVALID_METHOD
Definition: exactTestM.c:89
#define MIN(a, b)
Definition: exactTestM.c:31
#define NEGATIVE_R
Definition: exactTestM.c:90
double F_M(double x, double alpha, double r, double Q)
Definition: exactTestM.c:44
struct TestMParameters exactM(double r, double EPS_ABS, double EPS_REL, int ode_method)
Definition: exactTestM.c:176
#define Rg
Definition: exactTestM.c:40
#define rhow
Definition: exactTestM.c:36
double dF_M(double x, double alpha, double r, double Q)
Definition: exactTestM.c:52
int exactM_old(double r, double *alpha, double *Drr, const double EPS_ABS, const double EPS_REL, const int ode_method)
Definition: exactTestM.c:94
#define barB
Definition: exactTestM.c:38
#define NOT_DONE
Definition: exactTestM.c:88
#define n
Definition: exactTestM.c:37
int funcM_ode_G(double r, const double alpha[], double f[], void *params)
Definition: exactTestM.c:60
#define rho
Definition: exactTestM.c:35
std::string printf(const char *format,...)