PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
exactTestK.c
Go to the documentation of this file.
1 /*
2  Copyright (C) 2007, 2011, 2014, 2016, 2023 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 "pism/verification/tests/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: */
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 = &params;
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 }
const double k_ice
Definition: exactTestK.c:32
const double B0
Definition: exactTestK.c:38
#define Nsum
Definition: exactTestK.c:45
const double c_p_ice
Definition: exactTestK.c:30
#define ALPHA_RELTOL
Definition: exactTestK.c:137
const double H0
Definition: exactTestK.c:37
static int exactK_old(const double t, const double z, double *TT, double *FF, const int bedrockIsIce_p)
Definition: exactTestK.c:47
const double SperA
Definition: exactTestK.c:28
const double rho_ice
Definition: exactTestK.c:31
struct TestKParameters exactK(double t, double z, int bedrock_is_ice)
Definition: exactTestK.c:129
const double G
Definition: exactTestK.c:40
double coscross(double alpha, void *params)
Definition: exactTestK.c:146
int print_alpha_k(const int N)
Definition: exactTestK.c:155
const double c_p_BR
Definition: exactTestK.c:33
const double Ts
Definition: exactTestK.c:39
const double k_BR
Definition: exactTestK.c:35
const double rho_BR
Definition: exactTestK.c:34
const double phi
Definition: exactTestK.c:41
static double F(double SL, double B, double H, double alpha)
std::string printf(const char *format,...)
static const double k
Definition: exactTestP.cc:42