texactTestsIJ.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
---
texactTestsIJ.c (4813B)
---
1 /*
2 Copyright (C) 2007-2008, 2011, 2014, 2015, 2016 Ed Bueler 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 <stdio.h>
22 #include <stdlib.h>
23 #include <math.h>
24 #include <gsl/gsl_math.h> /* M_PI */
25 #include "exactTestsIJ.h"
26
27 #define secpera 31556926.0 /* seconds per year; 365.2422 days */
28
29 struct TestIParameters exactI(const double m, const double x, const double y) {
30 /* see exact solution for an ice stream sliding over plastic till described
31 on pages 237 and 238 of C. Schoof 2006 "A variational approach to ice streams"
32 J Fluid Mech 556 pp 227--251 */
33 /* const double n = 3.0, p = 1.0 + 1.0/n; // p = 4/3 */
34
35 struct TestIParameters result;
36
37 const double L = 40e3; /* = 40km; y in [-3L,3L]; full width is 6L = 240 km */
38 const double aspect = 0.05;
39 const double h0 = aspect * L; /* if aspect = 0.05 and L = 40km then h0 = 2000 m */
40 const double theta = atan(0.001); /* a slope of 1/1000, a la Siple streams */
41 const double rho = 910, g = 9.81; /* kg/m^3 and m/s^2, resp. */
42 const double f = rho * g * h0 * tan(theta); /* about 18 kPa given above
43 values for rho,g,theta,aspect,L */
44 const double W = pow(m+1.0,1.0/m) * L; /* e.g. W = 1.2 * L if m = 10 */
45 const double B = 3.7e8; /* Pa s^{1/3}; hardness given on p. 239 of Schoof;
46 why so big? */
47
48 const double s = fabs(y/L);
49
50 double C0, C1, C2, C3, C4, z1, z2, z3, z4;
51
52 result.tauc = f * pow(s,m);
53
54 /* to compute bed, assume bed(x=0)=0 and bed is sloping down for increasing x;
55 if tan(theta) = 0.001 and -Lx < x < Lx with Lx = 240km then bed(x=Lx) = -240 m */
56 result.bed = - x * tan(theta);
57
58 /* formula (4.3) in Schoof; note u is indep of aspect because f/h0 ratio gives C0 */
59 if (fabs(y) < W) {
60 C0 = 2.0 * pow(f / (B * h0),3.0) * pow(L,4.0);
61 /* printf(" C0*secpera = %10.5e\n",C0*secpera); */
62 C1 = pow(m + 1.0, 4.0/m);
63 C2 = (m+1.0) * C1;
64 C3 = (m+1.0) * C2;
65 C4 = (m+1.0) * C3;
66 z1 = (pow(s,4.0) - C1) / 4.0;
67 z2 = (pow(s,m+4.0) - C2) / ((m+1.0) * (m+4.0));
68 z3 = (pow(s,2.0*m+4.0) - C3) / ((m+1.0)*(m+1.0) * (2.0*m+4.0));
69 z4 = (pow(s,3.0*m+4.0) - C4) / (pow((m+1.0),3.0) * (3.0*m+4.0));
70 /* printf(" u / C0 = %10.5e\n",- (z1 - 3.0 * z2 + 3.0 * z3 - z4)); */
71 result.u = - C0 * (z1 - 3.0 * z2 + 3.0 * z3 - z4); /* comes out positive */
72 } else {
73 result.u = 0.0;
74 }
75 result.v = 0.0; /* no transverse flow */
76 return result;
77 }
78
79
80 struct TestJParameters exactJ(const double x, const double y) {
81 /* documented only in preprint by Bueler */
82
83 struct TestJParameters result;
84
85 const double L = 300.0e3; /* 300 km half-width */
86 const double H0 = 500.0; /* 500 m typical thickness */
87 /* use Ritz et al (2001) value of 30 MPa year for typical
88 vertically-averaged viscosity */
89 const double nu0 = 30.0 * 1.0e6 * secpera; /* = 9.45e14 Pa s */
90 const double rho_ice = 910.0; /* kg/m^3 */
91 const double rho_sw = 1028.0; /* kg/m^3 */
92 const double g = 9.81; /* m/s^2 */
93 const double C = rho_ice * g * (1.0 - rho_ice/rho_sw) * H0 / (2.0 * nu0);
94 const double my_gamma[3][3] = {{1.0854, 0.108, 0.0027},
95 {0.402 , 0.04 , 0.001 },
96 {0.0402, 0.004, 0.0001}};
97 const double A = L / (4.0 * M_PI);
98 double uu = 0.0, vv = 0.0, denom, trig, kx, ly, B;
99 int k,l;
100
101 result.H = H0 * (1.0 + 0.4 * cos(M_PI * x / L)) * (1.0 + 0.1 * cos(M_PI * y / L));
102 result.nu = (nu0 * H0) / result.H; /* so \nu(x,y) H(x,y) = \nu_0 H_0 */
103 for (k=-2; k<=2; k++) {
104 for (l=-2; l<=2; l++) {
105 if ((k != 0) || (l != 0)) { /* note alpha_00 = beta_00 = 0 */
106 denom = (double)(k * k + l * l);
107 kx = (double)(k) * M_PI * x / L;
108 ly = (double)(l) * M_PI * y / L;
109 trig = cos(kx) * sin(ly) + sin(kx) * cos(ly);
110 B = (A / denom) * (C * my_gamma[abs(k)][abs(l)]) * trig;
111 uu += B * (double)(k);
112 vv += B * (double)(l);
113 }
114 }
115 }
116 result.u = uu;
117 result.v = vv;
118
119 return result;
120 }
121