Skip to content

Commit 6f4a021

Browse files
valeriabarravaleriabarra
committed
SWE: Define coordinate transformations for points expressed with local coordinate systems on different panels (faces) of the cube
1 parent 0136af5 commit 6f4a021

File tree

4 files changed

+251
-57
lines changed

4 files changed

+251
-57
lines changed

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

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -108,13 +108,13 @@ CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q,
108108
};
109109

110110
CeedScalar dxdX[3][2];
111-
for (int j=0; j<3; j++)
112-
for (int k=0; k<2; k++) {
111+
for (CeedInt j=0; j<3; j++) {
112+
for (CeedInt k=0; k<2; k++) {
113113
dxdX[j][k] = 0;
114-
for (int l=0; l<3; l++)
114+
for (CeedInt l=0; l<3; l++)
115115
dxdX[j][k] += dxdxx[j][l]*dxxdX[l][k];
116116
}
117-
117+
}
118118
// J is given by the cross product of the columns of dxdX
119119
const CeedScalar J[3] = {dxdX[1][0]*dxdX[2][1] - dxdX[2][0]*dxdX[1][1],
120120
dxdX[2][0]*dxdX[0][1] - dxdX[0][0]*dxdX[2][1],
@@ -129,13 +129,13 @@ CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q,
129129

130130
// dxdX_k,j * dxdX_j,k
131131
CeedScalar dxdXTdxdX[2][2];
132-
for (int j=0; j<2; j++)
133-
for (int k=0; k<2; k++) {
132+
for (CeedInt j=0; j<2; j++) {
133+
for (CeedInt k=0; k<2; k++) {
134134
dxdXTdxdX[j][k] = 0;
135-
for (int l=0; l<3; l++)
135+
for (CeedInt l=0; l<3; l++)
136136
dxdXTdxdX[j][k] += dxdX[l][j]*dxdX[l][k];
137137
}
138-
138+
}
139139
const CeedScalar detdxdXTdxdX = dxdXTdxdX[0][0] * dxdXTdxdX[1][1]
140140
-dxdXTdxdX[1][0] * dxdXTdxdX[0][1];
141141

@@ -151,12 +151,13 @@ CEED_QFUNCTION(SetupGeo)(void *ctx, CeedInt Q,
151151

152152
// Compute the pseudo inverse of dxdX
153153
CeedScalar pseudodXdx[2][3];
154-
for (int j=0; j<2; j++)
155-
for (int k=0; k<3; k++) {
154+
for (CeedInt j=0; j<2; j++) {
155+
for (CeedInt k=0; k<3; k++) {
156156
pseudodXdx[j][k] = 0;
157-
for (int l=0; l<2; l++)
157+
for (CeedInt l=0; l<2; l++)
158158
pseudodXdx[j][k] += dxdXTdxdXinv[j][l]*dxdX[k][l];
159159
}
160+
}
160161

161162
// Interp-to-Grad qdata
162163
// Pseudo inverse of dxdX: (x_i,j)+ = X_i,j

‎examples/fluids/shallow-water/shallowwater.c‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,6 @@ int main(int argc, char **argv) {
5656
Units units;
5757
problemType problemChoice;
5858
problemData *problem = NULL;
59-
EdgeNode edgenodes;
6059
StabilizationType stab;
6160
PetscBool implicit;
6261
PetscInt degree, qextra, outputfreq, steps, contsteps;
@@ -66,7 +65,7 @@ int main(int argc, char **argv) {
6665
const CeedInt ncompx = 3;
6766
PetscInt viz_refine = 0;
6867
PetscBool read_mesh, simplex, test;
69-
PetscInt topodim = 2, ncompq = 3, lnodes, nedgenodes;
68+
PetscInt topodim = 2, ncompq = 3, lnodes;
7069
// libCEED context
7170
char ceedresource[PETSC_MAX_PATH_LEN] = "/cpu/self",
7271
filename[PETSC_MAX_PATH_LEN];
@@ -369,14 +368,14 @@ int main(int argc, char **argv) {
369368
cEnd - cStart, gdofs/ncompq, odofs/ncompq, ncompq,
370369
gdofs, odofs); CHKERRQ(ierr);
371370
}
372-
373371
}
374372

375373
// Set up global mass vector
376374
ierr = VecDuplicate(Q, &user->M); CHKERRQ(ierr);
377375

378376
// Setup global lat-long vector for different panels (charts) of the cube
379-
ierr = FindPanelEdgeNodes(dm, ncompq, &nedgenodes, &edgenodes); CHKERRQ(ierr);
377+
Mat T;
378+
ierr = FindPanelEdgeNodes(dm, &phys_ctx, ncompq, &T); CHKERRQ(ierr);
380379

381380
// Setup libCEED's objects
382381
ierr = PetscMalloc1(1, &ceeddata); CHKERRQ(ierr);

‎examples/fluids/shallow-water/src/setup.c‎

Lines changed: 226 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -71,22 +71,39 @@ problemData problemOptions[] = {
7171
// Auxiliary function to determine if nodes belong to cube faces (panels)
7272
// -----------------------------------------------------------------------------
7373

74-
PetscErrorCode FindPanelEdgeNodes(DM dm, PetscInt ncomp, PetscInt *nedgenodes,
75-
EdgeNode *edgenodes) {
74+
PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx,
75+
PetscInt ncomp, Mat *T) {
7676

77-
PetscInt ierr;
78-
PetscInt cstart, cend, nstart, nend, nnodes, depth, edgenode = 0;
77+
PetscInt ierr, rankid;
78+
MPI_Comm comm;
79+
PetscInt cstart, cend, nstart, nend, lnodes, gdofs, depth, edgenodecnt = 0;
7980
PetscSection section;
81+
EdgeNode edgenodes;
82+
8083
PetscFunctionBeginUser;
8184
// Get Nelem
82-
ierr = DMGetSection(dm, &section); CHKERRQ(ierr);
85+
ierr = DMGetLocalSection(dm, &section); CHKERRQ(ierr);
8386
ierr = DMPlexGetHeightStratum(dm, 0, &cstart, &cend); CHKERRQ(ierr);
8487
ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr);
8588
ierr = DMPlexGetHeightStratum(dm, depth, &nstart, &nend); CHKERRQ(ierr);
86-
nnodes = nend - nstart;
87-
unsigned int bitmap[nnodes];
89+
lnodes = nend - nstart;
90+
91+
PetscSF sf;
92+
Vec bitmapVec;
93+
ierr = DMGetSectionSF(dm, &sf); CHKERRQ(ierr);
94+
ierr = DMCreateGlobalVector(dm, &bitmapVec); CHKERRQ(ierr);
95+
ierr = VecGetSize(bitmapVec, &gdofs); CHKERRQ(ierr);
96+
ierr = VecDestroy(&bitmapVec);
97+
98+
// Arrays for local and global bitmaps
99+
unsigned int bitmaploc[lnodes * ncomp];
100+
unsigned int bitmap[gdofs];
88101
ierr = PetscMemzero(bitmap, sizeof(bitmap)); CHKERRQ(ierr);
89102

103+
// Broadcast bitmap values to all ranks
104+
ierr = PetscSFBcastBegin(sf, MPI_UNSIGNED, bitmap, bitmaploc); CHKERRQ(ierr);
105+
ierr = PetscSFBcastEnd(sf, MPI_UNSIGNED, bitmap, bitmaploc); CHKERRQ(ierr);
106+
90107
// Get indices
91108
for (PetscInt c = cstart; c < cend; c++) { // Traverse elements
92109
PetscInt numindices, *indices, n, panel;
@@ -96,49 +113,220 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PetscInt ncomp, PetscInt *nedgenodes,
96113
&numindices, &indices, NULL, NULL);
97114
CHKERRQ(ierr);
98115
for (n = 0; n < numindices; n += ncomp) { // Traverse nodes per element
99-
PetscInt bitmapidx = indices[n] / ncomp;
100-
bitmap[bitmapidx] |= (1 << panel);
116+
PetscInt bitmapidx = indices[n];
117+
bitmaploc[bitmapidx] |= (1 << panel);
101118
}
102119
ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE,
103120
&numindices, &indices, NULL, NULL);
104121
CHKERRQ(ierr);
105122
}
106-
// Read the 1's in the resulting bitmap and extract edge nodes only
107-
ierr = PetscMalloc1(nnodes + 24, edgenodes); CHKERRQ(ierr);
108-
for (PetscInt i = 0; i < nnodes; i++) {
109-
PetscInt ones = 0, panels[3];
110-
for (PetscInt p = 0; p < 6; p++) {
111-
if (bitmap[i] & 1)
112-
panels[ones++] = p;
113-
bitmap[i] >>= 1;
114-
}
115-
if (ones == 2) {
116-
(*edgenodes)[edgenode].idx = i;
117-
(*edgenodes)[edgenode].panelA = panels[0];
118-
(*edgenodes)[edgenode].panelB = panels[1];
119-
edgenode++;
123+
124+
// Reduce on all ranks
125+
ierr = PetscSFReduceBegin(sf, MPI_UNSIGNED, bitmaploc, bitmap, MPI_BOR);
126+
CHKERRQ(ierr);
127+
ierr = PetscSFReduceEnd(sf, MPI_UNSIGNED, bitmaploc, bitmap, MPI_BOR);
128+
CHKERRQ(ierr);
129+
130+
// Rank 0 reads the 1's in the resulting bitmap and extracts edge nodes only
131+
PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr);
132+
MPI_Comm_rank(comm, &rankid);
133+
if (rankid == 0) {
134+
ierr = PetscMalloc1(gdofs + 24*ncomp, &edgenodes); CHKERRQ(ierr);
135+
for (PetscInt i = 0; i < gdofs; i += ncomp) {
136+
PetscInt ones = 0, panels[3];
137+
for (PetscInt p = 0; p < 6; p++) {
138+
if (bitmap[i] & 1)
139+
panels[ones++] = p;
140+
bitmap[i] >>= 1;
141+
}
142+
if (ones == 2) {
143+
edgenodes[edgenodecnt].idx = i;
144+
edgenodes[edgenodecnt].panelA = panels[0];
145+
edgenodes[edgenodecnt].panelB = panels[1];
146+
edgenodecnt++;
147+
}
148+
else if (ones == 3) {
149+
edgenodes[edgenodecnt].idx = i;
150+
edgenodes[edgenodecnt].panelA = panels[0];
151+
edgenodes[edgenodecnt].panelB = panels[1];
152+
edgenodecnt++;
153+
edgenodes[edgenodecnt].idx = i;
154+
edgenodes[edgenodecnt].panelA = panels[0];
155+
edgenodes[edgenodecnt].panelB = panels[2];
156+
edgenodecnt++;
157+
edgenodes[edgenodecnt].idx = i;
158+
edgenodes[edgenodecnt].panelA = panels[1];
159+
edgenodes[edgenodecnt].panelB = panels[2];
160+
edgenodecnt++;
161+
}
120162
}
121-
else if (ones == 3) {
122-
(*edgenodes)[edgenode].idx = i;
123-
(*edgenodes)[edgenode].panelA = panels[0];
124-
(*edgenodes)[edgenode].panelB = panels[1];
125-
edgenode++;
126-
(*edgenodes)[edgenode].idx = i;
127-
(*edgenodes)[edgenode].panelA = panels[0];
128-
(*edgenodes)[edgenode].panelB = panels[2];
129-
edgenode++;
130-
(*edgenodes)[edgenode].idx = i;
131-
(*edgenodes)[edgenode].panelA = panels[1];
132-
(*edgenodes)[edgenode].panelB = panels[2];
133-
edgenode++;
163+
ierr = SetupPanelCoordTransformations(dm, phys_ctx, ncomp, edgenodes,
164+
edgenodecnt, T); CHKERRQ(ierr);
165+
// Free edgenodes structure array
166+
ierr = PetscFree(edgenodes); CHKERRQ(ierr);
167+
}
168+
169+
PetscFunctionReturn(0);
170+
}
171+
172+
// -----------------------------------------------------------------------------
173+
// Auxiliary function that sets up all corrdinate transformations between panels
174+
// -----------------------------------------------------------------------------
175+
PetscErrorCode SetupPanelCoordTransformations(DM dm, PhysicsContext phys_ctx,
176+
PetscInt ncomp,
177+
EdgeNode edgenodes,
178+
PetscInt nedgenodes, Mat *T) {
179+
PetscInt ierr;
180+
MPI_Comm comm;
181+
Vec X;
182+
PetscInt gdofs;
183+
const PetscScalar *xarray;
184+
185+
PetscFunctionBeginUser;
186+
ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr);
187+
188+
ierr = DMGetCoordinates(dm, &X); CHKERRQ(ierr);
189+
ierr = VecGetSize(X, &gdofs); CHKERRQ(ierr);
190+
191+
// Preallocate sparse matrix
192+
ierr = MatCreateBAIJ(comm, 2, PETSC_DECIDE, PETSC_DECIDE, 4*gdofs/ncomp,
193+
4*gdofs/ncomp, 2, NULL, 0, NULL, T); CHKERRQ(ierr);
194+
for (PetscInt i=0; i < 4*gdofs/ncomp; i++) {
195+
ierr = MatSetValue(*T, i, i, 1., INSERT_VALUES); CHKERRQ(ierr);
196+
}
197+
198+
ierr = VecGetArrayRead(X, &xarray); CHKERRQ(ierr);
199+
PetscScalar R = phys_ctx->R;
200+
// Nodes loop
201+
for (PetscInt i=0; i < gdofs; i += ncomp) {
202+
// Read global Cartesian coordinates
203+
PetscScalar x[3] = {xarray[i*ncomp + 0],
204+
xarray[i*ncomp + 1],
205+
xarray[i*ncomp + 2]
206+
};
207+
// Normalize quadrature point coordinates to sphere
208+
PetscScalar rad = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
209+
x[0] *= R / rad;
210+
x[1] *= R / rad;
211+
x[2] *= R / rad;
212+
// Compute latitude and longitude
213+
const PetscScalar theta = asin(x[2] / R); // latitude
214+
const PetscScalar lambda = atan2(x[1], x[0]); // longitude
215+
216+
// For P_1 (east), P_3 (front), P_4 (west), P_5 (back):
217+
PetscScalar T00 = cos(theta)*cos(lambda) * cos(lambda);
218+
PetscScalar T01 = cos(theta)*cos(lambda) * 0.;
219+
PetscScalar T10 = cos(theta)*cos(lambda) * -sin(theta)*sin(lambda);
220+
PetscScalar T11 = cos(theta)*cos(lambda) * cos(theta);
221+
const PetscScalar T_lateral[2][2] = {{T00,
222+
T01},
223+
{T10,
224+
T11}
225+
};
226+
PetscScalar Tinv00 = 1./(cos(theta)*cos(lambda)) * (1./cos(lambda));
227+
PetscScalar Tinv01 = 1./(cos(theta)*cos(lambda)) * 0.;
228+
PetscScalar Tinv10 = 1./(cos(theta)*cos(lambda)) * tan(theta)*tan(lambda);
229+
PetscScalar Tinv11 = 1./(cos(theta)*cos(lambda)) * (1./cos(theta));
230+
const PetscScalar T_lateralinv[2][2] = {{Tinv00,
231+
Tinv01},
232+
{Tinv10,
233+
Tinv11}
234+
};
235+
// For P2 (north):
236+
T00 = sin(theta) * cos(lambda);
237+
T01 = sin(theta) * sin(lambda);
238+
T10 = sin(theta) * -sin(theta)*sin(lambda);
239+
T11 = sin(theta) * sin(theta)*cos(lambda);
240+
const PetscScalar T_top[2][2] = {{T00,
241+
T01},
242+
{T10,
243+
T11}
244+
};
245+
Tinv00 = 1./(sin(theta)*sin(theta)) * sin(theta)*cos(lambda);
246+
Tinv01 = 1./(sin(theta)*sin(theta)) * (-sin(lambda));
247+
Tinv10 = 1./(sin(theta)*sin(theta)) * sin(theta)*sin(lambda);
248+
Tinv11 = 1./(sin(theta)*sin(theta)) * cos(lambda);
249+
const PetscScalar T_topinv[2][2] = {{Tinv00,
250+
Tinv01},
251+
{Tinv10,
252+
Tinv11}
253+
};
254+
255+
// For P0 (south):
256+
T00 = sin(theta) * (-cos(theta));
257+
T01 = sin(theta) * sin(lambda);
258+
T10 = sin(theta) * sin(theta)*sin(lambda);
259+
T11 = sin(theta) * sin(theta)*cos(lambda);
260+
const PetscScalar T_bottom[2][2] = {{T00,
261+
T01},
262+
{T10,
263+
T11}
264+
};
265+
Tinv00 = 1./(sin(theta)*sin(theta)) * (-sin(theta)*cos(lambda));
266+
Tinv01 = 1./(sin(theta)*sin(theta)) * sin(lambda);
267+
Tinv10 = 1./(sin(theta)*sin(theta)) * sin(theta)*sin(lambda);
268+
Tinv11 = 1./(sin(theta)*sin(theta)) * cos(lambda);
269+
const PetscScalar T_bottominv[2][2] = {{Tinv00,
270+
Tinv01},
271+
{Tinv10,
272+
Tinv11}
273+
};
274+
275+
const PetscScalar (*transforms[6])[2][2] = {&T_bottom,
276+
&T_lateral,
277+
&T_top,
278+
&T_lateral,
279+
&T_lateral,
280+
&T_lateral
281+
};
282+
283+
const PetscScalar (*inv_transforms[6])[2][2] = {&T_bottominv,
284+
&T_lateralinv,
285+
&T_topinv,
286+
&T_lateralinv,
287+
&T_lateralinv,
288+
&T_lateralinv
289+
};
290+
291+
for (PetscInt e = 0; e < nedgenodes; e++) {
292+
if (edgenodes[e].idx == i) {
293+
const PetscScalar (*matrixA)[2][2] = inv_transforms[edgenodes[e].panelA];
294+
const PetscScalar (*matrixB)[2][2] = transforms[edgenodes[e].panelB];
295+
296+
// inv_transform * transform (A^{-1}*B)
297+
// This product represents the mapping from coordinate system A
298+
// to spherical coordinates and then to coordinate system B. Vice versa
299+
// for transform * inv_transform (B^{-1}*A)
300+
PetscScalar matrixAB[2][2], matrixBA[2][2];
301+
for (int j=0; j<2; j++) {
302+
for (int k=0; k<2; k++) {
303+
matrixAB[j][k] = matrixBA[j][k] = 0;
304+
for (int l=0; l<2; l++) {
305+
matrixAB[j][k] += (*matrixA)[j][l] * (*matrixB)[l][k];
306+
matrixBA[j][k] += (*matrixB)[j][l] * (*matrixA)[l][k];
307+
}
308+
}
309+
}
310+
PetscInt idxAB[2] = {4*i/ncomp + 0, 4*i/ncomp +1};
311+
PetscInt idxBA[2] = {4*i/ncomp + 2, 4*i/ncomp +3};
312+
ierr = MatSetValues(*T, 2, idxAB, 2, idxAB, (PetscScalar *)matrixAB,
313+
INSERT_VALUES); CHKERRQ(ierr);
314+
ierr = MatSetValues(*T, 2, idxBA, 2, idxBA, (PetscScalar *)matrixBA,
315+
INSERT_VALUES); CHKERRQ(ierr);
316+
}
134317
}
135318
}
136-
ierr = PetscRealloc(edgenode, edgenodes); CHKERRQ(ierr);
137-
*nedgenodes = edgenode;
319+
// Assemble matrix for all node transformations
320+
ierr = MatAssemblyBegin(*T, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
321+
ierr = MatAssemblyEnd(*T, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
322+
323+
// Restore array read
324+
ierr = VecRestoreArrayRead(X, &xarray); CHKERRQ(ierr);
138325

139326
PetscFunctionReturn(0);
140327
}
141328

329+
142330
// -----------------------------------------------------------------------------
143331
// Auxiliary function to create PETSc FE space for a given degree
144332
// -----------------------------------------------------------------------------

0 commit comments

Comments
 (0)