PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
exactTestsFG.cc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2004-2008, 2014, 2015, 2016, 2023 Ed Bueler and Jed Brown 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/exactTestsFG.hh"
22 #include <cmath> // for pow, exp, cos, sin, M_PI, fabs, sqrt
23 #include <cstddef> // for size_t
24 #include <stdexcept> // for runtime_error
25 
26 namespace pism {
27 
28 static double p3(double x) {
29  // p_3=x^3-3*x^2+6*x-6, using Horner's
30  return -6.0 + x*(6.0 + x*(-3.0 + x));
31 }
32 
33 static double p4(double x) {
34  // p_4=x^4-4*x^3+12*x^2-24*x+24, using Horner's
35  return 24.0 + x*(-24.0 + x*(12.0 + x*(-4.0 + x)));
36 }
37 
38 TestFGParameters exactFG(double t, double r, const std::vector<double> &z, double Cp) {
39 
40  const double SperA = 31556926.0; // seconds per year; 365.2422 days
41 
42  // parameters describing extent of sheet:
43  const double H0 = 3000.0; // m
44  const double L = 750000.0; // m
45 
46  // period of perturbation; inactive in Test F:
47  const double Tp = 2000.0*SperA; // s
48 
49  // fundamental physical constants
50  const double g = 9.81; // m / s^2; accel of gravity
51  const double Rgas = 8.314; // J / (mol K)
52 
53  // ice properties; parameters which appear in constitutive relation:
54  const double rho = 910.0; // kg / m^3; density
55  const double k = 2.1; // J / m K s; thermal conductivity
56  const double cpheat = 2009.0; // J / kg K; specific heat capacity
57  const double n = 3.0; // Glen exponent
58 
59  // next two are EISMINT II values; Paterson - Budd for T < 263
60  const double A = 3.615E-13; // Pa^-3 s^-1
61  const double Q = 6.0E4; // J / mol
62 
63  // EISMINT II temperature boundary condition (Experiment F):
64  const double Ggeo = 0.042; // J / m^2 s; geo. heat flux
65  const double ST = 1.67E-5; // K m^-1
66  const double Tmin = 223.15; // K
67  const double Kcond = k / (rho*cpheat); // constant in temp equation
68 
69  const size_t Mz = z.size();
70  TestFGParameters result(Mz);
71  double
72  &H = result.H,
73  &M = result.M;
74  std::vector<double>
75  &T = result.T,
76  &U = result.U,
77  &w = result.w,
78  &Sig = result.Sig,
79  &Sigc = result.Sigc;
80 
81  // temporary storage
82  std::vector<double> I3(Mz);
83 
84  if ((r <= 0) or (r >= L)) {
85  throw std::runtime_error("exactFG(): code and derivation assume 0 < r < L !");
86  }
87 
88  // compute H from analytical steady state Hs (Test D) plus perturbation
89  const double power = n / (2*n + 2);
90  const double Hconst = H0 / pow(1 - 1 / n, power);
91  const double s = r / L;
92  const double lamhat = (1 + 1 / n)*s - (1 / n) + pow(1 - s, 1 + 1 / n) - pow(s, 1 + 1 / n);
93 
94  double f = 0.0;
95  if ((r > 0.3*L) and (r < 0.9*L)) {
96  f = pow(cos(M_PI*(r - 0.6*L) / (0.6*L)) , 2.0);
97  } else {
98  f = 0.0;
99  }
100 
101  const double goft = Cp*sin(2.0*M_PI*t / Tp);
102 
103  H = Hconst*pow(lamhat, power) + goft*f;
104 
105  // compute T = temperature
106  const double Ts = Tmin + ST*r;
107  const double nusqrt = sqrt(1 + (4.0*H*Ggeo) / (k*Ts));
108  const double nu = (k*Ts / (2.0*Ggeo))*(1 + nusqrt);
109  for (size_t i = 0; i < Mz; i++) {
110  if (z[i] < H) {
111  T[i] = Ts * (nu + H) / (nu + z[i]);
112  } else { // surface value above ice surface; matches numerical way
113  T[i] = Ts;
114  }
115  // old way: extend formula above surface: T[i] = Ts * (nu + H) / (nu + z[i]);
116  }
117 
118  // compute surface slope and horizontal velocity
119  const double lamhatr = ((1 + 1 / n) / L)*(1 - pow(1 - s, 1 / n) - pow(s, 1 / n));
120 
121  double fr = 0.0;
122  if ((r > 0.3*L) and (r < 0.9*L)) {
123  fr = - (M_PI / (0.6*L)) * sin(2.0*M_PI*(r - 0.6*L) / (0.6*L));
124  } else {
125  fr = 0.0;
126  }
127 
128  const double Hr = Hconst * power * pow(lamhat, power - 1) * lamhatr + goft*fr; // chain rule
129  if (Hr > 0) {
130  throw std::runtime_error("exactFG(): assumes H_r negative for all 0 < r < L !");
131  }
132 
133  const double mu = Q / (Rgas*Ts*(nu + H));
134  const double surfArr = exp(-Q / (Rgas*Ts));
135  const double Uconst = 2.0 * pow(rho*g, n) * A;
136  const double omega = Uconst * pow( - Hr, n) * surfArr * pow(mu, -n - 1);
137 
138  for (size_t i = 0; i < Mz; i++) {
139  if (z[i] < H) {
140  I3[i] = p3(mu*H) * exp(mu*H) - p3(mu*(H - z[i])) * exp(mu*(H - z[i]));
141  U[i] = omega * I3[i];
142  } else { // surface value above ice surface; matches numerical way
143  I3[i] = p3(mu*H) * exp(mu*H) - p3(0.0); // z[i] = H case in above
144  U[i] = omega * I3[i];
145  }
146  }
147 
148  // compute strain heating
149  for (size_t i = 0; i < Mz; i++) {
150  if (z[i] < H) {
151  const double Sigmu = - (Q*(nu + z[i])) / (Rgas*Ts*(nu + H));
152  Sig[i] = (Uconst*g / cpheat) * exp(Sigmu) * pow(fabs(Hr)*(H - z[i]) , n + 1);
153  } else {
154  Sig[i] = 0.0;
155  }
156  }
157 
158  // compute vertical velocity
159  const double lamhatrr = ((1 + 1 / n) / (n*L*L)) * (pow(1 - s, (1 / n) - 1) - pow(s, (1 / n) - 1));
160 
161  double frr = 0.0;
162  if ((r > 0.3*L) and (r < 0.9*L)) {
163  frr = - (2.0*M_PI*M_PI / (0.36*L*L)) * cos(2.0*M_PI*(r - 0.6*L) / (0.6*L));
164  } else {
165  frr = 0.0;
166  }
167 
168  const double Hrr = Hconst*power*(power - 1)*pow(lamhat, power - 2.0) * pow(lamhatr, 2.0) +
169  Hconst*power*pow(lamhat, power - 1)*lamhatrr + goft*frr;
170  const double Tsr = ST;
171  const double nur = (k*Tsr / (2.0*Ggeo)) * (1 + nusqrt) + (1 / Ts) * (Hr*Ts - H*Tsr) / nusqrt;
172  const double mur = ( - Q / (Rgas*Ts*Ts*pow(nu + H, 2.0))) * (Tsr*(nu + H) + Ts*(nur + Hr));
173  const double phi = 1 / r + n*Hrr / Hr + Q*Tsr / (Rgas*Ts*Ts) - (n + 1)*mur / mu; // division by r
174  const double gam = pow(mu, n) * exp(mu*H) * (mur*H + mu*Hr) * pow(H, n);
175  for (size_t i = 0; i < Mz; i++) {
176  const double I4 = p4(mu*H) * exp(mu*H) - p4(mu*(H - z[i])) * exp(mu*(H - z[i]));
177  w[i] = omega * ((mur / mu - phi)*I4 / mu + (phi*(H - z[i]) + Hr)*I3[i] - gam*z[i]);
178  }
179 
180  // compute compensatory accumulation M
181  const double I4H = p4(mu*H) * exp(mu*H) - 24.0;
182  const double divQ = - omega * (mur / mu - phi) * I4H / mu + omega * gam * H;
183  const double Ht = (Cp*2.0*M_PI / Tp) * cos(2.0*M_PI*t / Tp) * f;
184  M = Ht + divQ;
185 
186  // compute compensatory heating
187  const double nut = Ht / nusqrt;
188  for (size_t i = 0; i < Mz; i++) {
189  if (z[i] < H) {
190  const double dTt = Ts * ((nut + Ht)*(nu + z[i]) - (nu + H)*nut) * pow(nu + z[i], - 2.0);
191  const double Tr = Tsr*(nu + H) / (nu + z[i])
192  + Ts * ((nur + Hr)*(nu + z[i]) - (nu + H)*nur) * pow(nu + z[i], - 2.0);
193  const double Tz = - Ts * (nu + H) * pow(nu + z[i], - 2.0);
194  const double Tzz = 2.0 * Ts * (nu + H) * pow(nu + z[i], - 3.0);
195  Sigc[i] = dTt + U[i]*Tr + w[i]*Tz - Kcond*Tzz - Sig[i];
196  } else {
197  Sigc[i] = 0.0;
198  }
199  }
200 
201  return result;
202 }
203 
204 } // end of namespace pism
const double H0
Definition: exactTestK.c:37
const double Ts
Definition: exactTestK.c:39
const double phi
Definition: exactTestK.c:41
static const double L
Definition: exactTestL.cc:40
#define n
Definition: exactTestM.c:37
#define rho
Definition: exactTestM.c:35
static const double ST
static const double Tmin
static double p4(double x)
Definition: exactTestsFG.cc:33
static const double g
Definition: exactTestP.cc:36
static double p3(double x)
Definition: exactTestsFG.cc:28
static const double SperA
Definition: exactTestP.cc:35
static const double k
Definition: exactTestP.cc:42
TestFGParameters exactFG(double t, double r, const std::vector< double > &z, double Cp)
Definition: exactTestsFG.cc:38
std::vector< double > U
Definition: exactTestsFG.hh:51
std::vector< double > Sig
Definition: exactTestsFG.hh:51
std::vector< double > Sigc
Definition: exactTestsFG.hh:51
std::vector< double > w
Definition: exactTestsFG.hh:51
std::vector< double > T
Definition: exactTestsFG.hh:51