Loading [MathJax]/jax/input/TeX/config.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
6namespace pism {
7
8Vector2d 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
18Vector2d 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
35Vector2d 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
43Vector2d 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
57Vector2d 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
72Vector2d 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
87Vector2d 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
95Vector2d 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
104Vector2d 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
112Vector2d 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
120Vector2d 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
132Vector2d 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
151Vector2d 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
166Vector2d 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
181Vector2d 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
189double 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
194Vector2d 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
203Vector2d 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
212Vector2d 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
221double 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
Generated by doxygen 1.9.8