PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
ssafd_pik_code.cc
Go to the documentation of this file.
1 const double dx2 = dx*dx, dy2 = dy*dy, d4 = 4*dx*dy, d2 = 2*dx*dy;
2 
3 /* Coefficients of the discretization of the first equation; u first, then v. */
4 double eq1[] = {
5  0, -c_n*bPP/dy2, 0,
6  -4*c_w*aMM/dx2, (c_n*bPP+c_s*bMM)/dy2+(4*c_e*aPP+4*c_w*aMM)/dx2, -4*c_e*aPP/dx2,
7  0, -c_s*bMM/dy2, 0,
8  c_w*aMM*bPw/d2+c_n*aMn*bPP/d4, (c_n*aPn*bPP-c_n*aMn*bPP)/d4+(c_w*aMM*bPP-c_e*aPP*bPP)/d2, -c_e*aPP*bPe/d2-c_n*aPn*bPP/d4,
9  (c_w*aMM*bMw-c_w*aMM*bPw)/d2+(c_n*aMM*bPP-c_s*aMM*bMM)/d4, (c_n*aPP*bPP-c_n*aMM*bPP-c_s*aPP*bMM+c_s*aMM*bMM)/d4+(c_e*aPP*bPP-c_w*aMM*bPP-c_e*aPP*bMM+c_w*aMM*bMM)/d2, (c_e*aPP*bPe-c_e*aPP*bMe)/d2+(c_s*aPP*bMM-c_n*aPP*bPP)/d4,
10  -c_w*aMM*bMw/d2-c_s*aMs*bMM/d4, (c_s*aMs*bMM-c_s*aPs*bMM)/d4+(c_e*aPP*bMM-c_w*aMM*bMM)/d2, c_e*aPP*bMe/d2+c_s*aPs*bMM/d4,
11 };
12 
13 /* Coefficients of the discretization of the second equation; u first, then v. */
14 double eq2[] = {
15  c_w*aMM*bPw/d4+c_n*aMn*bPP/d2, (c_n*aPn*bPP-c_n*aMn*bPP)/d2+(c_w*aMM*bPP-c_e*aPP*bPP)/d4, -c_e*aPP*bPe/d4-c_n*aPn*bPP/d2,
16  (c_w*aMM*bMw-c_w*aMM*bPw)/d4+(c_n*aMM*bPP-c_s*aMM*bMM)/d2, (c_n*aPP*bPP-c_n*aMM*bPP-c_s*aPP*bMM+c_s*aMM*bMM)/d2+(c_e*aPP*bPP-c_w*aMM*bPP-c_e*aPP*bMM+c_w*aMM*bMM)/d4, (c_e*aPP*bPe-c_e*aPP*bMe)/d4+(c_s*aPP*bMM-c_n*aPP*bPP)/d2,
17  -c_w*aMM*bMw/d4-c_s*aMs*bMM/d2, (c_s*aMs*bMM-c_s*aPs*bMM)/d2+(c_e*aPP*bMM-c_w*aMM*bMM)/d4, c_e*aPP*bMe/d4+c_s*aPs*bMM/d2,
18  0, -4*c_n*bPP/dy2, 0,
19  -c_w*aMM/dx2, (4*c_n*bPP+4*c_s*bMM)/dy2+(c_e*aPP+c_w*aMM)/dx2, -c_e*aPP/dx2,
20  0, -4*c_s*bMM/dy2, 0,
21 };
22 
23 /* i indices */
24 const int I[] = {
25  i-1, i, i+1,
26  i-1, i, i+1,
27  i-1, i, i+1,
28  i-1, i, i+1,
29  i-1, i, i+1,
30  i-1, i, i+1,
31 };
32 
33 /* j indices */
34 const int J[] = {
35  j+1, j+1, j+1,
36  j, j, j,
37  j-1, j-1, j-1,
38  j+1, j+1, j+1,
39  j, j, j,
40  j-1, j-1, j-1,
41 };
42 
43 /* component indices */
44 const int C[] = {
45  0, 0, 0,
46  0, 0, 0,
47  0, 0, 0,
48  1, 1, 1,
49  1, 1, 1,
50  1, 1, 1,
51 };
52 
const int J[]
double eq1[]
const int I[]
const double d2
const int C[]
const double d4
const double dy2
double eq2[]
const double dx2