texactTestK.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
---
texactTestK.c (8156B)
---
1 /*
2 Copyright (C) 2007, 2011, 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 #include <stdio.h>
22 #include <math.h>
23 #include <gsl/gsl_errno.h>
24 #include <gsl/gsl_math.h> /* M_PI */
25 #include <gsl/gsl_roots.h>
26 #include "exactTestK.h"
27
28 const double SperA = 31556926.0; /* seconds per year; 365.2422 days */
29
30 const double c_p_ice = 2009.0; /* J/(kg K) specific heat capacity of ice */
31 const double rho_ice = 910.0; /* kg/(m^3) density of ice */
32 const double k_ice = 2.10; /* J/(m K s) = W/(m K) thermal conductivity of ice */
33 const double c_p_BR = 1000.0; /* J/(kg K) specific heat capacity of bedrock */
34 const double rho_BR = 3300.0; /* kg/(m^3) density of bedrock */
35 const double k_BR = 3.0; /* J/(m K s) = W/(m K) thermal conductivity of bedrock */
36
37 const double H0 = 3000.0; /* m */
38 const double B0 = 1000.0; /* m */
39 const double Ts = 223.15; /* K */
40 const double G = 0.042; /* W/(m^2) */
41 const double phi = 0.0125; /* K/m */
42
43 /* number of terms in eigenfunction expansion; the exact
44 solution is deliberately chosen to have finite expansion */
45 #define Nsum 30
46
47 static int exactK_old(const double t, const double z, double *TT, double *FF, const int bedrockIsIce_p) {
48 int k;
49 int belowB0;
50 double ZZ, alpha, lambda, beta, my_gamma, XkSQR, Xk,
51 theta, dthetakdz, P, dPdz,
52 Ck, I1, I2, aH, bB, mI, mR;
53 double c_p_bedrock, rho_bedrock, k_bedrock;
54 /* following constants were produced by calling print_alpha_k(30) (below) */
55 double alf[Nsum] = {3.350087528822397e-04, 1.114576827617396e-03, 1.953590840303518e-03,
56 2.684088585781064e-03, 3.371114869333445e-03, 4.189442265117592e-03,
57 5.008367405382524e-03, 5.696044031764593e-03, 6.425563506942886e-03,
58 7.264372872913219e-03, 8.044853066396166e-03, 8.714877612414516e-03,
59 9.493529164160654e-03, 1.033273985210279e-02, 1.106421822502108e-02,
60 1.175060460132703e-02, 1.256832682090360e-02, 1.338784224692084e-02,
61 1.407617951778051e-02, 1.480472324161026e-02, 1.564331999062109e-02,
62 1.642470780103220e-02, 1.709475346624607e-02, 1.787248418996684e-02,
63 1.871188358061674e-02, 1.944434477688470e-02, 2.013010181370026e-02,
64 2.094721145334310e-02, 2.176730968036079e-02, 2.245631776169424e-02};
65
66 if (bedrockIsIce_p) {
67 c_p_bedrock = c_p_ice;
68 rho_bedrock = rho_ice;
69 k_bedrock = k_ice;
70 for (k = 0; k < Nsum; k++) { /* overwrite alpha_k with ice-meets-ice values; see preprint */
71 alf[k] = (2.0 * k + 1.0) * M_PI / (2.0 * (H0 + B0));
72 }
73 } else {
74 c_p_bedrock = c_p_BR;
75 rho_bedrock = rho_BR;
76 k_bedrock = k_BR;
77 }
78 if (z > H0) {
79 *TT = Ts;
80 return 0;
81 }
82 belowB0 = (z < -B0);
83
84 ZZ = sqrt((rho_bedrock * c_p_bedrock * k_ice) / (rho_ice * c_p_ice * k_bedrock));
85 mI = (G / k_ice) - phi; mR = (G / k_bedrock) - phi;
86 /* DEBUG: printf("ZZ = %10e, mI = %10e, mR = %10e\n", ZZ,mI,mR); */
87 *TT = 0.0;
88 *FF = 0.0;
89 for (k = Nsum-1; k >= 0; k--) {
90 /* constants only having to do with eigenfunctions; theta = theta_k(z) is the
91 normalized eigenfunction */
92 alpha = alf[k];
93 beta = ZZ * alpha;
94 my_gamma = sin(alpha * H0) / cos(beta * B0);
95 XkSQR = (rho_bedrock * c_p_bedrock * my_gamma * my_gamma * B0 + rho_ice * c_p_ice * H0) / 2.0;
96 Xk = sqrt(XkSQR);
97 /* theta = ((z > 0) ? sin(alpha * (H0 - z)) : my_gamma * cos(beta * (B0 + z))) / Xk; */
98 theta = (z > 0) ? sin(alpha * (H0 - z))
99 : my_gamma * cos(beta * (B0 + z));
100 theta /= Xk;
101 dthetakdz = (z > 0) ? - alpha * cos(alpha * (H0 - z))
102 : - beta * my_gamma * sin(beta * (B0 + z));
103 dthetakdz /= Xk;
104 lambda = (k_ice * alpha * alpha) / (rho_ice * c_p_ice);
105 /* DEBUG: printf("k = %3d: alpha = %10e, Xk = %10e, theta = %10e, dthetakdz = %10e, lambda = %10e,\n",
106 k,alpha,Xk,theta,dthetakdz,lambda); */
107 /* constants involved in computing the expansion coefficients */
108 aH = alpha * H0; bB = beta * B0;
109 I1 = - mI * (sin(aH) - aH * cos(aH)) / (alpha * alpha);
110 I2 = mR * (cos(bB) - 1.0 + bB * sin(bB)) / (beta * beta)
111 - (B0 * mR + H0 * mI) * sin(bB) / beta;
112 Ck = (rho_ice * c_p_ice * I1 + rho_bedrock * c_p_bedrock * my_gamma * I2) / Xk;
113 /* add the term to the expansion */
114 *TT += Ck * exp(- lambda * t) * theta;
115 *FF += - ((z > 0) ? k_ice : k_bedrock) * Ck * exp(- lambda * t) * dthetakdz;
116 /* DEBUG: printf(" I1 = %10e, I2 = %10e, Ck = %10e, term = %10f\n",
117 I1,I2,Ck, Ck * exp(- lambda * t) * theta); */
118 }
119 /* P = (z >= 0) ? (z / k_ice) - (H0 / k_ice) : (z / k_bedrock) - (H0 / k_ice); */
120 P = (z / ((z > 0) ? k_ice : k_bedrock)) - (H0 / k_ice);
121 dPdz = 1.0 / ((z > 0) ? k_ice : k_bedrock);
122 *TT += Ts - G * P;
123 *FF += ((z > 0) ? k_ice : k_bedrock) * G * dPdz;
124
125 return belowB0;
126
127 }
128
129 struct TestKParameters exactK(double t, double z, int bedrock_is_ice) {
130 struct TestKParameters result;
131
132 result.error_code = exactK_old(t, z, &result.T, &result.F, bedrock_is_ice);
133
134 return result;
135 }
136
137 #define ALPHA_RELTOL 1.0e-14
138 #define ITER_MAXED_OUT 999
139
140 /* parameters needed for root problem: */
141 struct coscross_params {
142 double Afrac, HZBsum, HZBdiff;
143 };
144
145 /* the root problem is to make this function zero: */
146 double coscross(double alpha, void *params) {
147 struct coscross_params *p = (struct coscross_params *) params;
148 return cos(p->HZBsum * alpha) - p->Afrac * cos(p->HZBdiff * alpha);
149 }
150
151 /* compute the first N roots alpha_k of the equation
152 ((A-1)/(A+1)) cos((H - Z B) alpha) = cos((H + Z B) alpha)
153 where H and B are heights and A, Z are defined in terms of material
154 constants */
155 int print_alpha_k(const int N) {
156 int status = 0, iter = 0, k = 0, max_iter = 200;
157 double Z = 0.0, A = 0.0;
158 double alpha = 0.0, alpha_lo = 0.0, alpha_hi = 0.0, temp_lo = 0.0;
159 const gsl_root_fsolver_type *solvT;
160 gsl_root_fsolver *solv;
161 gsl_function F;
162 struct coscross_params params;
163
164 Z = sqrt((rho_BR * c_p_BR * k_ice) / (rho_ice * c_p_ice * k_BR));
165 A = (k_BR / k_ice) * Z;
166 params.Afrac = (A - 1.0) / (A + 1.0);
167 params.HZBsum = H0 + Z * B0;
168 params.HZBdiff = H0 - Z * B0;
169
170 F.function = &coscross;
171 F.params = ¶ms;
172 solvT = gsl_root_fsolver_brent; /* faster than bisection but still bracketing */
173 solv = gsl_root_fsolver_alloc(solvT);
174
175 for (k = 0; k < N; k++) {
176 /* these numbers bracket exactly one solution */
177 alpha_lo = ((double)(k) * M_PI) / params.HZBsum;
178 alpha_hi = ((double)(k + 1) * M_PI) / params.HZBsum;
179 gsl_root_fsolver_set(solv, &F, alpha_lo, alpha_hi);
180
181 iter = 0;
182 do {
183 iter++;
184
185 status = gsl_root_fsolver_iterate(solv);
186 if (status != GSL_SUCCESS) {
187 goto cleanup;
188 }
189
190 alpha = gsl_root_fsolver_root(solv);
191 alpha_lo = gsl_root_fsolver_x_lower(solv);
192 alpha_hi = gsl_root_fsolver_x_upper(solv);
193 temp_lo = (alpha_lo > 0) ? alpha_lo : (alpha_hi/2.0);
194
195 status = gsl_root_test_interval(temp_lo, alpha_hi, 0, ALPHA_RELTOL);
196 } while ((status == GSL_CONTINUE) && (iter < max_iter));
197
198 if (iter >= max_iter) {
199 printf("!!!ERROR: root finding iteration reached maximum iterations; QUITTING!\n");
200 goto cleanup;
201 }
202
203 printf("%19.15e,\n",alpha);
204 }
205
206 cleanup:
207 gsl_root_fsolver_free(solv);
208 return status;
209 }