PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
manufactured_solutions.cc
Go to the documentation of this file.
1 // DO NOT EDIT. This code was generated by a Python script.
2 #include <cmath>
3 
4 #include "pism/stressbalance/blatter/verification/manufactured_solutions.hh"
5 
6 namespace pism {
7 
8 Vector2d blatter_xy_exact(double x, double y)
9 {
10  double x0 = exp(x);
11  double x1 = 2*M_PI*y;
12  return {
13  x0*sin(x1),
14  x0*cos(x1)
15  };
16 }
17 
18 Vector2d blatter_xy_source(double x, double y, double B)
19 {
20  double x0 = 2*M_PI*y;
21  double x1 = cos(x0);
22  double x2 = pow(x1, 2);
23  double x3 = pow(M_PI, 2);
24  double x4 = pow(M_PI, 3);
25  double x5 = pow(M_PI, 4);
26  double x6 = sin(x0);
27  double x7 = pow(x6, 2);
28  double x8 = pow(2, 2.0/3.0)*B*exp((1.0/3.0)*x)/pow(x2*(4*x3 + 1 + 4*M_PI) + x7*(16*x3 - 8*M_PI + 4), 4.0/3.0);
29  return {
30  -1.0/3.0*x6*x8*(x2*(6*x3 + 20*x4 + 72*x5 - 2 + 3*M_PI) + x7*(-48*x3 + 32*x4 + 96*x5 - 8 + 36*M_PI)),
31  -1.0/6.0*x1*x8*(x2*(-12*x3 + 136*x4 + 192*x5 - 18*M_PI - 1) + x7*(96*x3 - 128*x4 + 384*x5 - 24*M_PI - 4))
32  };
33 }
34 
35 Vector2d blatter_xz_exact(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta) {
36 
37  return {
38  -4*A*pow(alpha, 3)*pow(g, 3)*pow(rho, 3)*pow(x, 3)*(-pow(H, 4) + pow(-alpha*pow(x, 2) + s_0 - z, 4)) + 2*H*alpha*g*rho*x/beta,
39  0.0
40  };
41 }
42 
43 Vector2d blatter_xz_source(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta) {
44  double phi_1 = alpha*pow(x, 2) - s_0 + z;
45  double phi_2 = 4*A*pow(alpha, 3)*pow(g, 3)*pow(rho, 3)*x;
46  double phi_3 = 4*pow(phi_1, 5)*pow(phi_2, 2)*pow(x, 3);
47  double phi_4 = -2*H*alpha*g*rho/beta + 8*alpha*pow(phi_1, 3)*phi_2*pow(x, 3) + 3*phi_2*x*(-pow(H, 4) + pow(phi_1, 4));
48  double phi_5 = 48*pow(alpha, 2)*pow(phi_1, 2)*phi_2*pow(x, 4) + 56*alpha*pow(phi_1, 3)*phi_2*pow(x, 2) + 6*phi_2*(-pow(H, 4) + pow(phi_1, 4));
49  double mu = (1.0/2.0)/cbrt(A*phi_1*phi_3*x + A*pow(phi_4, 2));
50 
51  return {
52  (16.0/3.0)*A*pow(mu, 4)*(-6*pow(phi_1, 3)*phi_2*phi_3*pow(x, 3) - 18*pow(phi_1, 2)*phi_2*pow(phi_4, 2)*pow(x, 2) - 6*phi_1*phi_3*phi_5*x + 24*phi_3*phi_4*(2*alpha*pow(x, 2) + phi_1) - 2*pow(phi_4, 2)*phi_5),
53  0.0
54  };
55 }
56 
57 Vector2d blatter_xz_source_bed(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta) {
58  double phi_1 = alpha*pow(x, 2) - s_0 + z;
59  double phi_2 = 4*A*pow(alpha, 3)*pow(g, 3)*pow(rho, 3)*x;
60  double phi_3 = 4*pow(phi_1, 5)*pow(phi_2, 2)*pow(x, 3);
61  double phi_4 = -2*H*alpha*g*rho/beta + 8*alpha*pow(phi_1, 3)*phi_2*pow(x, 3) + 3*phi_2*x*(-pow(H, 4) + pow(phi_1, 4));
62  double mu = (1.0/2.0)/cbrt(A*phi_1*phi_3*x + A*pow(phi_4, 2));
63  double n_x = -2*alpha*x/sqrt(4*pow(alpha, 2)*pow(x, 2) + 1);
64  double n_z = -1/sqrt(4*pow(alpha, 2)*pow(x, 2) + 1);
65 
66  return {
67  2*H*alpha*g*rho*x - beta*phi_2*pow(x, 2)*(-pow(H, 4) + pow(phi_1, 4)) - 4*mu*n_x*phi_4 - 4*mu*n_z*pow(phi_1, 3)*phi_2*pow(x, 2),
68  0.0
69  };
70 }
71 
72 Vector2d blatter_xz_source_surface(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta) {
73  double phi_1 = alpha*pow(x, 2) - s_0 + z;
74  double phi_2 = 4*A*pow(alpha, 3)*pow(g, 3)*pow(rho, 3)*x;
75  double phi_3 = 4*pow(phi_1, 5)*pow(phi_2, 2)*pow(x, 3);
76  double phi_4 = -2*H*alpha*g*rho/beta + 8*alpha*pow(phi_1, 3)*phi_2*pow(x, 3) + 3*phi_2*x*(-pow(H, 4) + pow(phi_1, 4));
77  double mu = (1.0/2.0)/cbrt(A*phi_1*phi_3*x + A*pow(phi_4, 2));
78  double n_x = 2*alpha*x/sqrt(4*pow(alpha, 2)*pow(x, 2) + 1);
79  double n_z = pow(4*pow(alpha, 2)*pow(x, 2) + 1, -1.0/2.0);
80 
81  return {
82  -4*mu*n_x*phi_4 - 4*mu*n_z*pow(phi_1, 3)*phi_2*pow(x, 2),
83  0.0
84  };
85 }
86 
87 Vector2d blatter_xz_cfbc_exact(double x, double z, double B, double L, double rho_i, double rho_w, double g)
88 {
89  return {
90  (1.0/2.0)*L*g*z*(rho_i - rho_w)*sin(M_PI*x/L)/(M_PI*B),
91  0
92  };
93 }
94 
95 Vector2d blatter_xz_cfbc_source(double x, double z, double L, double rho_i, double rho_w, double g)
96 {
97  double x0 = M_PI/L;
98  return {
99  -g*x0*z*(rho_i - rho_w)*sin(x*x0),
100  0
101  };
102 }
103 
104 Vector2d blatter_xz_cfbc_surface(double x, double L, double rho_i, double rho_w, double g)
105 {
106  return {
107  (1.0/4.0)*L*g*(rho_i - rho_w)*sin(M_PI*x/L)/M_PI,
108  0
109  };
110 }
111 
112 Vector2d blatter_xz_cfbc_base(double x, double L, double rho_i, double rho_w, double g)
113 {
114  return {
115  -1.0/4.0*L*g*(rho_i - rho_w)*sin(M_PI*x/L)/M_PI,
116  0
117  };
118 }
119 
120 Vector2d blatter_xz_halfar_exact(double x, double z, double H_0, double R_0, double rho_i, double g, double B) {
121  double C_0 = H_0;
122  double C_1 = 1.0/R_0;
123  double C_2 = (1.0/2.0)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
124  double h0 = C_0*pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 3.0/7.0);
125  double h_x = -4.0/7.0*C_0*pow(C_1, 4.0/3.0)*cbrt(x)/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
126  return {
127  -C_2*pow(h_x, 3)*(pow(h0, 4) - pow(h0 - z, 4)),
128  0
129  };
130 }
131 
132 Vector2d blatter_xz_halfar_source(double x, double z, double H_0, double R_0, double rho_i, double g, double B) {
133  double C_0 = H_0;
134  double C_1 = 1.0/R_0;
135  double C_2 = (1.0/2.0)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
136  double h0 = C_0*pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 3.0/7.0);
137  double h_x = -4.0/7.0*C_0*pow(C_1, 4.0/3.0)*cbrt(x)/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
138  double h_xx = (4.0/147.0)*C_0*pow(C_1, 4.0/3.0)*(16*pow(C_1, 4.0/3.0)*pow(x, 2.0/3.0)/(pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) - 1) - 7/pow(x, 2.0/3.0))/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
139  double h_xxx = -8.0/3087.0*C_0*(-168*pow(C_1, 8.0/3.0)/(cbrt(x)*(pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) - 1)) - 49*pow(C_1, 4.0/3.0)/pow(x, 5.0/3.0) + 352*pow(C_1, 4)*x/pow(pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) - 1, 2))/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
140  double u_x = -C_2*pow(h_x, 3)*(4*pow(h0, 3)*h_x - 4*h_x*pow(h0 - z, 3)) - 3*C_2*pow(h_x, 2)*h_xx*(pow(h0, 4) - pow(h0 - z, 4));
141  double u_z = -4*C_2*pow(h_x, 3)*pow(h0 - z, 3);
142  double u_xx = -C_2*pow(h_x, 3)*(4*pow(h0, 3)*h_xx + 12*pow(h0, 2)*pow(h_x, 2) - 12*pow(h_x, 2)*pow(h0 - z, 2) - 4*h_xx*pow(h0 - z, 3)) - 6*C_2*pow(h_x, 2)*h_xx*(4*pow(h0, 3)*h_x - 4*h_x*pow(h0 - z, 3)) - 3*C_2*pow(h_x, 2)*h_xxx*(pow(h0, 4) - pow(h0 - z, 4)) - 6*C_2*h_x*pow(h_xx, 2)*(pow(h0, 4) - pow(h0 - z, 4));
143  double u_xz = -12*C_2*pow(h_x, 4)*pow(h0 - z, 2) - 12*C_2*pow(h_x, 2)*h_xx*pow(h0 - z, 3);
144  double u_zz = 12*C_2*pow(h_x, 3)*pow(h0 - z, 2);
145  return {
146  2*B*u_x*(-2.0/3.0*u_x*u_xx - 1.0/6.0*u_xz*u_z)/pow(pow(u_x, 2) + (1.0/4.0)*pow(u_z, 2), 4.0/3.0) + 2*B*u_xx/cbrt(pow(u_x, 2) + (1.0/4.0)*pow(u_z, 2)) + (1.0/2.0)*B*u_z*(-2.0/3.0*u_x*u_xz - 1.0/6.0*u_z*u_zz)/pow(pow(u_x, 2) + (1.0/4.0)*pow(u_z, 2), 4.0/3.0) + (1.0/2.0)*B*u_zz/cbrt(pow(u_x, 2) + (1.0/4.0)*pow(u_z, 2)),
147  0.0
148  };
149 }
150 
151 Vector2d blatter_xz_halfar_source_lateral(double x, double z, double H_0, double R_0, double rho_i, double g, double B) {
152  double C_0 = H_0;
153  double C_1 = 1.0/R_0;
154  double C_2 = (1.0/2.0)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
155  double h0 = C_0*pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 3.0/7.0);
156  double h_x = -4.0/7.0*C_0*pow(C_1, 4.0/3.0)*cbrt(x)/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
157  double h_xx = (4.0/147.0)*C_0*pow(C_1, 4.0/3.0)*(16*pow(C_1, 4.0/3.0)*pow(x, 2.0/3.0)/(pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) - 1) - 7/pow(x, 2.0/3.0))/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
158  double u_x = -C_2*pow(h_x, 3)*(4*pow(h0, 3)*h_x - 4*h_x*pow(h0 - z, 3)) - 3*C_2*pow(h_x, 2)*h_xx*(pow(h0, 4) - pow(h0 - z, 4));
159  double u_z = -4*C_2*pow(h_x, 3)*pow(h0 - z, 3);
160  return {
161  2*pow(2, 2.0/3.0)*B*u_x/cbrt(4*pow(u_x, 2) + pow(u_z, 2)),
162  0.0
163  };
164 }
165 
166 Vector2d blatter_xz_halfar_source_surface(double x, double H_0, double R_0, double rho_i, double g, double B) {
167  double C_0 = H_0;
168  double C_1 = 1.0/R_0;
169  double C_2 = (1.0/2.0)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
170  double h0 = C_0*pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 3.0/7.0);
171  double z = h0;
172  double h_x = -4.0/7.0*C_0*pow(C_1, 4.0/3.0)*cbrt(x)/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
173  double h_xx = (4.0/147.0)*C_0*pow(C_1, 4.0/3.0)*(16*pow(C_1, 4.0/3.0)*pow(x, 2.0/3.0)/(pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) - 1) - 7/pow(x, 2.0/3.0))/pow(-pow(C_1, 4.0/3.0)*pow(x, 4.0/3.0) + 1, 4.0/7.0);
174  double u_x = -C_2*pow(h_x, 3)*(4*pow(h0, 3)*h_x - 4*h_x*pow(h0 - z, 3)) - 3*C_2*pow(h_x, 2)*h_xx*(pow(h0, 4) - pow(h0 - z, 4));
175  return {
176  -2*B*h_x*u_x/(sqrt(pow(h_x, 2) + 1)*cbrt(pow(u_x, 2))),
177  0.0
178  };
179 }
180 
181 Vector2d blatter_xz_halfar_source_base(double x, double H_0, double R_0, double rho_i, double g, double B) {
182  (void) B;
183  return {
184  -4.0/7.0*pow(H_0, 2)*g*rho_i*cbrt(x)/(pow(R_0, 8.0/7.0)*pow(pow(R_0, 4.0/3.0) - pow(x, 4.0/3.0), 1.0/7.0)),
185  0.0
186  };
187 }
188 
189 double blatter_xz_vanderveen_thickness(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B) {
190  double C = (1.0/8.0)*pow(alpha, 3)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
191  return pow(4*C*x/Q_0 + pow(H_0, -4), -1.0/4.0);
192 }
193 
194 Vector2d blatter_xz_vanderveen_exact(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B) {
195  double C = (1.0/8.0)*pow(alpha, 3)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
196  double thickness = pow(4*C*x/Q_0 + pow(H_0, -4), -1.0/4.0);
197  return {
198  Q_0/thickness,
199  0
200  };
201 }
202 
203 Vector2d blatter_xz_vanderveen_source_lateral(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B) {
204  double C = (1.0/8.0)*pow(alpha, 3)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
205  double thickness = pow(4*C*x/Q_0 + pow(H_0, -4), -1.0/4.0);
206  return {
207  2*B*cbrt(C)*thickness,
208  0.0
209  };
210 }
211 
212 Vector2d blatter_xz_vanderveen_source_surface(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B) {
213  double C = (1.0/8.0)*pow(alpha, 3)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
214  double thickness = pow(4*C*x/Q_0 + pow(H_0, -4), -1.0/4.0);
215  return {
216  2*B*pow(C, 4.0/3.0)*alpha*pow(thickness, 6)/sqrt(pow(C, 2)*pow(alpha, 2)*pow(thickness, 10) + pow(Q_0, 2)),
217  0.0
218  };
219 }
220 
221 double blatter_xz_vanderveen_beta(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B) {
222  double C = (1.0/8.0)*pow(alpha, 3)*pow(g, 3)*pow(rho_i, 3)/pow(B, 3);
223  double thickness = pow(4*C*x/Q_0 + pow(H_0, -4), -1.0/4.0);
224  return 2*B*pow(C, 4.0/3.0)*pow(thickness, 7)*(alpha - 1)/(Q_0*sqrt(pow(C, 2)*pow(thickness, 10)*pow(alpha - 1, 2) + pow(Q_0, 2)));
225 }
226 
227 } // end of namespace pism
pism::Vector2d
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2d.hh:29
L
static const double L
Definition: exactTestL.cc:40
rho
#define rho
Definition: exactTestM.c:35
pism::blatter_xz_vanderveen_beta
double blatter_xz_vanderveen_beta(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:221
pism::blatter_xy_source
Vector2d blatter_xy_source(double x, double y, double B)
Definition: manufactured_solutions.cc:18
pism::blatter_xz_cfbc_source
Vector2d blatter_xz_cfbc_source(double x, double z, double L, double rho_i, double rho_w, double g)
Definition: manufactured_solutions.cc:95
pism::blatter_xz_vanderveen_thickness
double blatter_xz_vanderveen_thickness(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:189
pism::blatter_xz_cfbc_base
Vector2d blatter_xz_cfbc_base(double x, double L, double rho_i, double rho_w, double g)
Definition: manufactured_solutions.cc:112
pism::blatter_xz_vanderveen_source_lateral
Vector2d blatter_xz_vanderveen_source_lateral(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:203
pism::blatter_xz_source_surface
Vector2d blatter_xz_source_surface(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta)
Definition: manufactured_solutions.cc:72
pism::blatter_xz_vanderveen_source_surface
Vector2d blatter_xz_vanderveen_source_surface(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:212
pism::blatter_xz_halfar_source_lateral
Vector2d blatter_xz_halfar_source_lateral(double x, double z, double H_0, double R_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:151
pism::g
static const double g
Definition: exactTestP.cc:36
pism::blatter_xz_halfar_source
Vector2d blatter_xz_halfar_source(double x, double z, double H_0, double R_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:132
pism::blatter_xz_vanderveen_exact
Vector2d blatter_xz_vanderveen_exact(double x, double alpha, double H_0, double Q_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:194
pism::blatter_xz_halfar_exact
Vector2d blatter_xz_halfar_exact(double x, double z, double H_0, double R_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:120
pism::blatter_xz_halfar_source_surface
Vector2d blatter_xz_halfar_source_surface(double x, double H_0, double R_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:166
pism::blatter_xz_cfbc_exact
Vector2d blatter_xz_cfbc_exact(double x, double z, double B, double L, double rho_i, double rho_w, double g)
Definition: manufactured_solutions.cc:87
pism::blatter_xz_source
Vector2d blatter_xz_source(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta)
Definition: manufactured_solutions.cc:43
pism::blatter_xz_exact
Vector2d blatter_xz_exact(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta)
Definition: manufactured_solutions.cc:35
pism::blatter_xz_source_bed
Vector2d blatter_xz_source_bed(double x, double z, double A, double rho, double g, double s_0, double alpha, double H, double beta)
Definition: manufactured_solutions.cc:57
pism::h0
static const double h0
Definition: exactTestP.cc:49
pism::blatter_xz_halfar_source_base
Vector2d blatter_xz_halfar_source_base(double x, double H_0, double R_0, double rho_i, double g, double B)
Definition: manufactured_solutions.cc:181
pism::blatter_xy_exact
Vector2d blatter_xy_exact(double x, double y)
Definition: manufactured_solutions.cc:8
pism::blatter_xz_cfbc_surface
Vector2d blatter_xz_cfbc_surface(double x, double L, double rho_i, double rho_w, double g)
Definition: manufactured_solutions.cc:104
C
const int C[]
Definition: ssafd_code.cc:44
Generated by doxygen 1.9.1