Skip to content

Commit 4afe695

Browse files
valeriabarravaleriabarra
committed
SWE: Put back geometric factors as they were
1 parent a41b89a commit 4afe695

File tree

1 file changed

+115
-103
lines changed

1 file changed

+115
-103
lines changed

‎examples/fluids/shallow-water/qfunctions/setup_geo.h‎

Lines changed: 115 additions & 103 deletions
Original file line numberDiff line numberDiff line change
@@ -26,44 +26,35 @@
2626

2727
// *****************************************************************************
2828
// This QFunction sets up the geometric factors required for integration and
29-
// coordinate transformations. See ref: "A Discontinuous Galerkin Transport
30-
// Scheme on the Cubed Sphere", by Nair et al. (2004).
29+
// coordinate transformations
3130
//
32-
// Reference (parent) 2D coordinates: X \in [-1, 1]^2.
31+
// Reference (parent) 2D coordinates: X \in [-1, 1]^2
3332
//
34-
// Local 2D physical coordinates on the 2D manifold: x \in [-l, l]^2
35-
// with l half edge of the cube inscribed in the sphere. These coordinate
36-
// systems vary locally on each face (or panel) of the cube.
33+
// Global 3D physical coordinates given by the mesh: xx \in [-R, R]^3
34+
// with R radius of the sphere
3735
//
38-
// (theta, lambda) represnt latitude and longitude coordinates.
36+
// Local 3D physical coordinates on the 2D manifold: x \in [-l, l]^3
37+
// with l half edge of the cube inscribed in the sphere
3938
//
40-
// Change of coordinates from x (on the 2D manifold) to xx (phyisical 3D on
41-
// the sphere), i.e., "cube-to-sphere" A, with equidistant central projection:
39+
// Change of coordinates matrix computed by the library:
40+
// (physical 3D coords relative to reference 2D coords)
41+
// dxx_j/dX_i (indicial notation) [3 * 2]
4242
//
43-
// For lateral panels (P0-P3):
44-
// A = R cos(theta)cos(lambda) / l * [cos(lambda) 0]
45-
// [-sin(theta)sin(lambda) cos(theta)]
43+
// Change of coordinates x (on the 2D manifold) relative to xx (phyisical 3D):
44+
// dx_i/dxx_j (indicial notation) [3 * 3]
4645
//
47-
// For top panel (P4):
48-
// A = R sin(theta) / l * [cos(lambda) sin(lambda)]
49-
// [-sin(theta)sin(lambda) sin(theta)cos(lambda)]
46+
// Change of coordinates x (on the 2D manifold) relative to X (reference 2D):
47+
// (by chain rule)
48+
// dx_i/dX_j [3 * 2] = dx_i/dxx_k [3 * 3] * dxx_k/dX_j [3 * 2]
5049
//
51-
// For bottom panel(P5):
52-
// A = R sin(theta) / l * [-cos(lambda) sin(lambda)]
53-
// [sin(theta)sin(lambda) sin(theta)cos(lambda)]
54-
//
55-
// The inverse of A, A^{-1}, is the "sphere-to-cube" change of coordinates.
56-
//
57-
// The metric tensor g_{ij} = A^TA, and its inverse,
58-
// g^{ij} = g_{ij}^{-1} = A^{-1}A^{-T}
59-
//
60-
// modJ represents the magnitude of the cross product of the columns of A, i.e.,
61-
// J = det(g_{ij})
50+
// modJ is given by the magnitude of the cross product of the columns of dx_i/dX_j
6251
//
6352
// The quadrature data is stored in the array qdata.
6453
//
6554
// We require the determinant of the Jacobian to properly compute integrals of
66-
// the form: int( u v ).
55+
// the form: int( u v )
56+
//
57+
// Qdata: modJ * w
6758
//
6859
// *****************************************************************************
6960
CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q,
@@ -77,92 +68,113 @@ CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q,
7768
CeedScalar (*qdata)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
7869
// *INDENT-ON*
7970

80-
CeedInt panel = -1;
81-
8271
CeedPragmaSIMD
83-
// Quadrature point loop to determine which panel the element belongs to
72+
// Quadrature Point Loop
8473
for (CeedInt i=0; i<Q; i++) {
85-
const CeedScalar pidx = X[2][i];
86-
if (pidx != -1)
87-
panel = pidx;
88-
break;
89-
}
90-
91-
// Check that we found nodes panel
92-
if (panel == -1)
93-
return -1;
94-
95-
CeedPragmaSIMD
96-
// Quadrature point loop for metric factors
97-
for (CeedInt i=0; i<Q; i++) {
98-
// Read local panel coordinates and relative panel index
99-
CeedScalar x[2] = {X[0][i],
100-
X[1][i]
101-
};
102-
const CeedScalar pidx = X[2][i];
103-
104-
if (pidx != panel) {
105-
const CeedScalar theta = X[0][i];
106-
const CeedScalar lambda = X[1][i];
107-
108-
CeedScalar T_inv[2][2];
109-
110-
switch (panel) {
111-
case 0:
112-
case 1:
113-
case 2:
114-
case 3:
115-
// For P_0 (front), P_1 (east), P_2 (back), P_3 (west):
116-
T_inv[0][0] = 1./(cos(theta)*cos(lambda)) * (1./cos(lambda));
117-
T_inv[0][1] = 1./(cos(theta)*cos(lambda)) * 0.;
118-
T_inv[1][0] = 1./(cos(theta)*cos(lambda)) * tan(theta)*tan(lambda);
119-
T_inv[1][1] = 1./(cos(theta)*cos(lambda)) * (1./cos(theta));
120-
break;
121-
case 4:
122-
// For P4 (north):
123-
T_inv[0][0] = 1./(sin(theta)*sin(theta)) * sin(theta)*cos(lambda);
124-
T_inv[0][1] = 1./(sin(theta)*sin(theta)) * (-sin(lambda));
125-
T_inv[1][0] = 1./(sin(theta)*sin(theta)) * sin(theta)*sin(lambda);
126-
T_inv[1][1] = 1./(sin(theta)*sin(theta)) * cos(lambda);
127-
break;
128-
case 5:
129-
// For P5 (south):
130-
T_inv[0][0] = 1./(sin(theta)*sin(theta)) * (-sin(theta)*cos(lambda));
131-
T_inv[0][1] = 1./(sin(theta)*sin(theta)) * sin(lambda);
132-
T_inv[1][0] = 1./(sin(theta)*sin(theta)) * sin(theta)*sin(lambda);
133-
T_inv[1][1] = 1./(sin(theta)*sin(theta)) * cos(lambda);
134-
break;
135-
}
136-
x[0] = T_inv[0][0]*theta + T_inv[0][1]*lambda;
137-
x[1] = T_inv[1][0]*theta + T_inv[1][1]*lambda;
138-
}
139-
140-
const CeedScalar xxT[2][2] = {{x[0]*x[0],
141-
x[0]*x[1]},
142-
{x[1]*x[0],
143-
x[1]*x[1]}
144-
};
145-
146-
// Read dxdX Jacobian entries, stored in columns
74+
// Read global Cartesian coordinates
75+
const CeedScalar xx[3] = {X[0][i],
76+
X[1][i],
77+
X[2][i]
78+
};
79+
// Read dxxdX Jacobian entries, stored in columns
14780
// J_00 J_10
14881
// J_01 J_11
149-
const CeedScalar dxdX[2][2] = {{J[0][0][i],
150-
J[1][0][i]},
151-
{J[0][1][i],
152-
J[1][1][i]}
82+
// J_02 J_12
83+
const CeedScalar dxxdX[3][2] = {{J[0][0][i],
84+
J[1][0][i]},
85+
{J[0][1][i],
86+
J[1][1][i]},
87+
{J[0][2][i],
88+
J[1][2][i]}
89+
};
90+
// Setup
91+
// x = xx (xx^T xx)^{-1/2}
92+
// dx/dxx = I (xx^T xx)^{-1/2} - xx xx^T (xx^T xx)^{-3/2}
93+
const CeedScalar modxxsq = xx[0]*xx[0]+xx[1]*xx[1]+xx[2]*xx[2];
94+
CeedScalar xxsq[3][3];
95+
for (int j=0; j<3; j++)
96+
for (int k=0; k<3; k++)
97+
xxsq[j][k] = xx[j]*xx[k] / (sqrt(modxxsq) * modxxsq);
98+
99+
const CeedScalar dxdxx[3][3] = {{1./sqrt(modxxsq) - xxsq[0][0],
100+
-xxsq[0][1],
101+
-xxsq[0][2]},
102+
{-xxsq[1][0],
103+
1./sqrt(modxxsq) - xxsq[1][1],
104+
-xxsq[1][2]},
105+
{-xxsq[2][0],
106+
-xxsq[2][1],
107+
1./sqrt(modxxsq) - xxsq[2][2]}
153108
};
154109

155-
CeedScalar dxxdX[2][2];
156-
for (int j=0; j<2; j++)
157-
for (int k=0; k<2; k++) {
158-
dxxdX[j][k] = 0;
159-
for (int l=0; l<2; l++)
160-
dxxdX[j][k] += xxT[j][l]*dxdX[l][k];
110+
CeedScalar dxdX[3][2];
111+
for (CeedInt j=0; j<3; j++) {
112+
for (CeedInt k=0; k<2; k++) {
113+
dxdX[j][k] = 0;
114+
for (CeedInt l=0; l<3; l++)
115+
dxdX[j][k] += dxdxx[j][l]*dxxdX[l][k];
161116
}
117+
}
118+
// J is given by the cross product of the columns of dxdX
119+
const CeedScalar J[3] = {dxdX[1][0]*dxdX[2][1] - dxdX[2][0]*dxdX[1][1],
120+
dxdX[2][0]*dxdX[0][1] - dxdX[0][0]*dxdX[2][1],
121+
dxdX[0][0]*dxdX[1][1] - dxdX[1][0]*dxdX[0][1]
122+
};
123+
124+
// Use the magnitude of J as our detJ (volume scaling factor)
125+
const CeedScalar modJ = sqrt(J[0]*J[0]+J[1]*J[1]+J[2]*J[2]);
162126

163127
// Interp-to-Interp qdata
164-
qdata[0][i] = (dxxdX[0][0]*dxxdX[1][1] - dxxdX[0][1]*dxxdX[1][0]) * w[i];
128+
qdata[0][i] = modJ * w[i];
129+
130+
// dxdX_k,j * dxdX_j,k
131+
CeedScalar dxdXTdxdX[2][2];
132+
for (CeedInt j=0; j<2; j++) {
133+
for (CeedInt k=0; k<2; k++) {
134+
dxdXTdxdX[j][k] = 0;
135+
for (CeedInt l=0; l<3; l++)
136+
dxdXTdxdX[j][k] += dxdX[l][j]*dxdX[l][k];
137+
}
138+
}
139+
const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1]
140+
-dxdXTdxdX[1][0] * dxdXTdxdX[0][1];
141+
142+
// Compute inverse of dxdXTdxdX, needed for the pseudoinverse. This is also
143+
// equivalent to the 2x2 metric tensor g^{ij}, needed for the
144+
// Grad-to-Grad qdata (pseudodXdx * pseudodXdxT, which simplifies to
145+
// dxdXTdxdXinv)
146+
CeedScalar dxdXTdxdXinv[2][2];
147+
dxdXTdxdXinv[0][0] = dxdXTdxdX[1][1] / detdxdXTdxdX;
148+
dxdXTdxdXinv[0][1] = -dxdXTdxdX[0][1] / detdxdXTdxdX;
149+
dxdXTdxdXinv[1][0] = -dxdXTdxdX[1][0] / detdxdXTdxdX;
150+
dxdXTdxdXinv[1][1] = dxdXTdxdX[0][0] / detdxdXTdxdX;
151+
152+
// Compute the pseudo inverse of dxdX
153+
CeedScalar pseudodXdx[2][3];
154+
for (CeedInt j=0; j<2; j++) {
155+
for (CeedInt k=0; k<3; k++) {
156+
pseudodXdx[j][k] = 0;
157+
for (CeedInt l=0; l<2; l++)
158+
pseudodXdx[j][k] += dxdXTdxdXinv[j][l]*dxdX[k][l];
159+
}
160+
}
165161

162+
// Interp-to-Grad qdata
163+
// Pseudo inverse of dxdX: (x_i,j)+ = X_i,j
164+
qdata[1][i] = pseudodXdx[0][0];
165+
qdata[2][i] = pseudodXdx[0][1];
166+
qdata[3][i] = pseudodXdx[0][2];
167+
qdata[4][i] = pseudodXdx[1][0];
168+
qdata[5][i] = pseudodXdx[1][1];
169+
qdata[6][i] = pseudodXdx[1][2];
170+
171+
// Grad-to-Grad qdata stored in Voigt convention
172+
qdata[7][i] = dxdXTdxdXinv[0][0];
173+
qdata[8][i] = dxdXTdxdXinv[1][1];
174+
qdata[9][i] = dxdXTdxdXinv[0][1];
175+
176+
// Terrain topography, hs
177+
qdata[10][i] = sin(xx[0]) + cos(xx[1]); // put 0 for constant flat topography
166178
} // End of Quadrature Point Loop
167179

168180
// Return

0 commit comments

Comments
 (0)