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