texactTestsABCD.c - 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
---
texactTestsABCD.c (6985B)
---
1 /*
2 Copyright (C) 2004-2006, 2014, 2016 Jed Brown and 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 #include <stdio.h>
22 #include <math.h>
23 #include <gsl/gsl_math.h> /* M_PI */
24 #include "exactTestsABCD.h"
25
26 #define SperA 31556926.0 /* seconds per year; 365.2422 days */
27
28 int exactA_old(const double r, double *H, double *M) {
29 /* NOTE: t is in seconds */
30 const double L = 750000.0; /* m; distance of margin from center */
31 const double M0 = 0.3 / SperA; /* 30 cm year-1 constant accumulation */
32 const double g = 9.81; /* m/s^2; accel of gravity */
33 const double rho = 910.0; /* kg/m^3; density */
34 const double n = 3.0; /* Glen exponent */
35 const double A = 1.0E-16/SperA; /* = 3.17e-24 1/(Pa^3 s); */
36 /* (EISMINT value) flow law parameter */
37 const double Gamma = 2 * pow(rho * g,n) * A / (n+2);
38
39 double m, p, C;
40
41 if (r < L) {
42 m = 2.0 * n + 2.0;
43 p = 1.0 + 1.0 / n;
44 C = pow(pow(2.0, n - 1) * M0 / Gamma, 1.0 / m);
45 *H = C * pow(pow(L, p) - pow(r, p), n / m);
46 } else {
47 *H = 0.0;
48 }
49
50 *M = M0;
51
52 return 0;
53 }
54
55 struct TestABCDParameters exactA(double r) {
56 struct TestABCDParameters result;
57
58 result.error_code = exactA_old(r, &result.H, &result.M);
59
60 return result;
61 }
62
63 int exactB_old(const double t, const double r, double *H, double *M) {
64 /* NOTE: t and t0 are in seconds */
65 double alpha, beta, t0, Rmargin;
66 const double n = 3.0, H0 = 3600.0, R0=750000.0;
67
68 /* lambda=0.0 case of Bueler et al (2005) family of similarity solns;
69 is Halfar (1983) soln */
70 alpha=1.0/9.0; /* alpha=(2-(n+1)*lambda)/(5*n+3)=1/9 */
71 beta=1.0/18.0; /* beta=(1+(2*n+1)*lambda)/(5*n+3)=1/18 */
72 t0=422.45*SperA; /* t0 = (beta/Gamma)
73 * pow((2n+1)/(n+1),n)*(pow(R0,n+1)/pow(H0,2n+1)) */
74
75 Rmargin=R0*pow(t/t0,beta);
76 if (r<Rmargin)
77 *H = H0 * pow(t/t0,-alpha)
78 * pow(1.0-pow(pow(t/t0,-beta)*(r/R0), (n+1)/n), n/(2*n+1));
79 else
80 *H = 0.0;
81
82 *M=0.0;
83 return 0;
84 }
85
86 struct TestABCDParameters exactB(const double t, const double r) {
87 struct TestABCDParameters result;
88
89 result.error_code = exactB_old(t, r, &result.H, &result.M);
90
91 return result;
92 }
93
94 int exactC_old(const double t, const double r, double *H, double *M) {
95 double lambda, alpha, beta, t0, Rmargin;
96 const double n = 3.0, H0 = 3600.0, R0=750000.0;
97
98 lambda=5.0;
99 alpha=-1.0; /* alpha=(2-(n+1)*lambda)/(5*n+3) */
100 beta=2.0; /* beta=(1+(2*n+1)*lambda)/(5*n+3) */
101 t0=15208.0*SperA; /* t0 = (beta/Gamma)
102 * pow((2n+1)/(n+1),n)*(pow(R0,n+1)/pow(H0,2n+1)) */
103
104 Rmargin=R0*pow(t/t0,beta);
105 if (r<Rmargin)
106 *H = H0 * pow(t/t0,-alpha)
107 * pow(1.0-pow(pow(t/t0,-beta)*(r/R0), (n+1)/n), n/(2*n+1));
108 else
109 *H = 0.0;
110
111 if (t>0.1*SperA)
112 *M = (lambda/t)* (*H);
113 else { /* when less than 0.1 year, avoid division by time */
114 Rmargin=R0*pow(0.1*SperA/t0,beta);
115 if (r<Rmargin)
116 *M=5*H0/t0; /* constant value in disc of Rmargin radius */
117 else
118 *M=0.0;
119 }
120 return 0;
121 }
122
123 struct TestABCDParameters exactC(const double t, const double r) {
124 struct TestABCDParameters result;
125
126 result.error_code = exactC_old(t, r, &result.H, &result.M);
127
128 return result;
129 }
130
131 int exactD_old(const double t, const double rin, double *H, double *M) {
132
133 /* parameters describing extent of sheet: */
134 const double H0=3600.0; /* m */
135 const double L=750000.0; /* m */
136 /* parameters for perturbation: */
137 const double Tp=5000.0*SperA; /* s */
138 const double Cp=200.0; /* m */
139 /* fundamental physical constants */
140 const double g=9.81; /* m/s^2; accel of gravity */
141 /* ice properties; parameters which appear in constitutive relation: */
142 const double rho=910.0; /* kg/m^3; density */
143 const double n=3.0; /* Glen exponent */
144 const double A=1.0E-16/SperA; /* = 3.17e-24 1/(Pa^3 s); */
145 /* (EISMINT value) flow law parameter */
146
147 double r=rin;
148 double Gamma, power, s, lamhat, f, Hs, temp, C, Ms, df, ddf, chi, dchi,
149 ddchi, c1, dHs, ddHs, dH, ddH, divterms, Mc;
150
151 if (r < 0.0) r=-r;
152
153 if (r >= L - 0.01) {
154 *H = 0.0;
155 *M = -0.1 / SperA;
156 } else {
157 if (r < 0.01) r = 0.01; /* avoid r=0 singularity in formulas */
158
159 /* important derived quantities */
160 Gamma = 2 * pow(rho * g,n) * A / (n+2);
161 power = n / (2*n+2);
162 s = r/L;
163
164 /* compute H from analytical steady state Hs plus perturbation */
165 lamhat = (1+1/n)*s - (1/n) + pow(1-s,1+1/n) - pow(s,1+1/n);
166 if ((r>0.3*L) && (r<0.9*L)) f = pow(cos(M_PI*(r-0.6*L)/(0.6*L)) ,2.0);
167 else f = 0.0;
168 Hs = (H0 / pow(1-1/n,power)) * pow(lamhat,power);
169 *H = Hs + Cp * sin(2.0*M_PI*t/Tp) * f;
170
171 /* compute steady part of accumulation */
172 temp = pow(s,1/n) + pow(1-s,1/n) - 1;
173 C = Gamma * pow(H0,2*n+2) / pow(2.0 * L * (1-1/n),n);
174 Ms = (C/r) * pow(temp,n-1)
175 * (2*pow(s,1/n) + pow(1-s,(1/n)-1) * (1 - 2*s) - 1);
176
177 /* derivs of H */
178 if ((r>0.3*L) && (r<0.9*L)) {
179 df = -(M_PI/(0.6*L)) * sin(M_PI*(r-0.6*L)/(0.3*L));
180 ddf = -(M_PI*M_PI/(0.18*L*L)) * cos(M_PI*(r-0.6*L)/(0.3*L));
181 chi = (4.0/3.0)*s - 1.0/3.0 + pow(1-s,4.0/3.0) - pow(s,4.0/3.0);
182 dchi = -(4.0/(3.0*L)) * (pow(s,1.0/3.0) + pow(1-s,1.0/3.0) - 1);
183 ddchi = -(4.0/(9.0*L*L)) * (pow(s,-2.0/3.0) - pow(1-s,-2.0/3.0));
184 c1 = (3.0*H0) / (8.0*pow(2.0/3.0,3.0/8.0));
185 dHs = c1 * pow(chi,-5.0/8.0) * dchi;
186 ddHs = c1 * ((-5.0/8.0) * pow(chi,-13.0/8.0) * dchi * dchi
187 + pow(chi,-5.0/8.0) * ddchi);
188 dH = dHs + Cp * sin(2.0*M_PI*t/Tp) * df;
189 ddH = ddHs + Cp * sin(2.0*M_PI*t/Tp) * ddf;
190 divterms = Gamma * pow(*H,4.) * dH * dH
191 * ((1.0/r) * (*H) * dH + 5.0 * dH * dH + 3.0 * (*H) * ddH);
192 Mc = (2.0*M_PI*Cp/Tp) * cos(2.0*M_PI*t/Tp) * f - Ms - divterms;
193 } else {
194 Mc = 0.0;
195 }
196
197 /* actual calculate total M */
198 *M = Ms + Mc;
199 }
200 return 0;
201 }
202
203 struct TestABCDParameters exactD(const double t, const double r) {
204 struct TestABCDParameters result;
205
206 result.error_code = exactD_old(t, r, &result.H, &result.M);
207
208 return result;
209 }