Skip to content

Commit fb02c59

Browse files
valeriabarravaleriabarra
committed
SWE: Initial implementation of generic restriction matrix
This generic restriction matrix does not only have 1's but also the coordinate transformations needed from the coord system of the nodes that are on the panel edges (lat, long) to the local ones of the nodes on the adjacent panels
1 parent db149ba commit fb02c59

File tree

4 files changed

+149
-164
lines changed

4 files changed

+149
-164
lines changed

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

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -374,8 +374,8 @@ int main(int argc, char **argv) {
374374
// Setup transformation matrix different panels of the cube
375375
Mat T;
376376
EdgeNode edgenodes;
377-
ierr = FindPanelEdgeNodes(dm, &physCtxData, ncompq, &nedgenodes, &edgenodes,
378-
&T);
377+
ierr = FindPanelEdgeNodes(dm, &physCtxData, ncompq, degree, topodim,
378+
&nedgenodes, &edgenodes, &T);
379379
CHKERRQ(ierr);
380380

381381
// Transform coordinate according to their panel systems

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

Lines changed: 128 additions & 143 deletions
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,8 @@ problemData problemOptions[] = {
7676
// -----------------------------------------------------------------------------
7777

7878
PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx,
79-
PetscInt ncomp, PetscInt *edgenodecnt,
79+
PetscInt ncomp, PetscInt degree,
80+
PetscInt topodim, PetscInt *edgenodecnt,
8081
EdgeNode *edgenodes, Mat *T) {
8182

8283
PetscInt ierr;
@@ -154,8 +155,8 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx,
154155
(*edgenodecnt)++;
155156
}
156157
}
157-
// ierr = SetupPanelCoordTransformations(dm, phys_ctx, ncomp, *edgenodes,
158-
// *edgenodecnt, T); CHKERRQ(ierr);
158+
ierr = SetupRestrictionMatrix(dm, phys_ctx, degree, ncomp, *edgenodes,
159+
*edgenodecnt, T); CHKERRQ(ierr);
159160

160161
// Free heap
161162
ierr = PetscFree2(bitmaploc, bitmap); CHKERRQ(ierr);
@@ -166,157 +167,148 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx,
166167
// -----------------------------------------------------------------------------
167168
// Auxiliary function that sets up all coordinate transformations between panels
168169
// -----------------------------------------------------------------------------
169-
PetscErrorCode SetupPanelCoordTransformations(DM dm, PhysicsContext phys_ctx,
170-
PetscInt ncomp,
171-
EdgeNode edgenodes,
172-
PetscInt nedgenodes, Mat *T) {
173-
PetscInt ierr;
170+
PetscErrorCode SetupRestrictionMatrix(DM dm, PhysicsContext phys_ctx,
171+
PetscInt degree, PetscInt ncomp,
172+
EdgeNode edgenodes, PetscInt nedgenodes,
173+
Mat *T) {
174+
PetscInt ierr, nnodes, c, cStart, cEnd, nelem, numP, gdofs;
174175
MPI_Comm comm;
175176
Vec X;
176-
PetscInt gdofs;
177+
PetscSection section;
177178
const PetscScalar *xarray;
178179

179180
PetscFunctionBeginUser;
180181
ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr);
182+
ierr = DMGetSection(dm, &section); CHKERRQ(ierr);
181183

182184
ierr = DMGetCoordinates(dm, &X); CHKERRQ(ierr);
183185
ierr = VecGetSize(X, &gdofs); CHKERRQ(ierr);
186+
nnodes = gdofs/ncomp;
187+
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
188+
nelem = cEnd - cStart;
189+
numP = degree + 1;
184190

185191
// Preallocate sparse matrix
186-
ierr = MatCreateBAIJ(comm, 2, PETSC_DECIDE, PETSC_DECIDE, 4*gdofs/ncomp,
187-
4*gdofs/ncomp, 2, NULL, 0, NULL, T); CHKERRQ(ierr);
188-
for (PetscInt i=0; i < 4*gdofs/ncomp; i++) {
189-
ierr = MatSetValue(*T, i, i, 1., INSERT_VALUES); CHKERRQ(ierr);
190-
}
192+
ierr = MatCreateBAIJ(comm, 3, nelem*numP*numP*ncomp, nnodes*ncomp,
193+
PETSC_DETERMINE, PETSC_DETERMINE, 3, NULL,
194+
nnodes*ncomp - 3, NULL, T); CHKERRQ(ierr);
191195

192-
ierr = VecGetArrayRead(X, &xarray); CHKERRQ(ierr);
193-
PetscScalar R = phys_ctx->R;
194-
// Nodes loop
195-
for (PetscInt i=0; i < gdofs; i += ncomp) {
196-
// Read global Cartesian coordinates
197-
PetscScalar x[3] = {xarray[i*ncomp + 0],
198-
xarray[i*ncomp + 1],
199-
xarray[i*ncomp + 2]
200-
};
201-
// Normalize quadrature point coordinates to sphere
202-
PetscScalar rad = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
203-
x[0] *= R / rad;
204-
x[1] *= R / rad;
205-
x[2] *= R / rad;
206-
// Compute latitude and longitude
207-
const PetscScalar theta = asin(x[2] / R); // latitude
208-
const PetscScalar lambda = atan2(x[1], x[0]); // longitude
209-
210-
// For P_0 (front), P_1 (east), P_2 (back), P_3 (west):
211-
PetscScalar T00 = cos(theta)*cos(lambda) * cos(lambda);
212-
PetscScalar T01 = cos(theta)*cos(lambda) * 0.;
213-
PetscScalar T10 = cos(theta)*cos(lambda) * -sin(theta)*sin(lambda);
214-
PetscScalar T11 = cos(theta)*cos(lambda) * cos(theta);
215-
const PetscScalar T_lateral[2][2] = {{T00,
216-
T01},
217-
{T10,
218-
T11}
219-
};
220-
PetscScalar Tinv00 = 1./(cos(theta)*cos(lambda)) * (1./cos(lambda));
221-
PetscScalar Tinv01 = 1./(cos(theta)*cos(lambda)) * 0.;
222-
PetscScalar Tinv10 = 1./(cos(theta)*cos(lambda)) * tan(theta)*tan(lambda);
223-
PetscScalar Tinv11 = 1./(cos(theta)*cos(lambda)) * (1./cos(theta));
224-
const PetscScalar T_lateralinv[2][2] = {{Tinv00,
225-
Tinv01},
226-
{Tinv10,
227-
Tinv11}
228-
};
229-
// For P4 (north):
230-
T00 = sin(theta) * cos(lambda);
231-
T01 = sin(theta) * sin(lambda);
232-
T10 = sin(theta) * -sin(theta)*sin(lambda);
233-
T11 = sin(theta) * sin(theta)*cos(lambda);
234-
const PetscScalar T_top[2][2] = {{T00,
235-
T01},
236-
{T10,
237-
T11}
238-
};
239-
Tinv00 = 1./(sin(theta)*sin(theta)) * sin(theta)*cos(lambda);
240-
Tinv01 = 1./(sin(theta)*sin(theta)) * (-sin(lambda));
241-
Tinv10 = 1./(sin(theta)*sin(theta)) * sin(theta)*sin(lambda);
242-
Tinv11 = 1./(sin(theta)*sin(theta)) * cos(lambda);
243-
const PetscScalar T_topinv[2][2] = {{Tinv00,
244-
Tinv01},
245-
{Tinv10,
246-
Tinv11}
247-
};
248-
249-
// For P5 (south):
250-
T00 = sin(theta) * (-cos(theta));
251-
T01 = sin(theta) * sin(lambda);
252-
T10 = sin(theta) * sin(theta)*sin(lambda);
253-
T11 = sin(theta) * sin(theta)*cos(lambda);
254-
const PetscScalar T_bottom[2][2] = {{T00,
255-
T01},
256-
{T10,
257-
T11}
258-
};
259-
Tinv00 = 1./(sin(theta)*sin(theta)) * (-sin(theta)*cos(lambda));
260-
Tinv01 = 1./(sin(theta)*sin(theta)) * sin(lambda);
261-
Tinv10 = 1./(sin(theta)*sin(theta)) * sin(theta)*sin(lambda);
262-
Tinv11 = 1./(sin(theta)*sin(theta)) * cos(lambda);
263-
const PetscScalar T_bottominv[2][2] = {{Tinv00,
264-
Tinv01},
265-
{Tinv10,
266-
Tinv11}
267-
};
268-
269-
const PetscScalar (*transforms[6])[2][2] = {&T_lateral,
270-
&T_lateral,
271-
&T_lateral,
272-
&T_lateral,
273-
&T_top,
274-
&T_bottom
275-
};
276-
277-
const PetscScalar (*inv_transforms[6])[2][2] = {&T_lateralinv,
278-
&T_lateralinv,
279-
&T_lateralinv,
280-
&T_lateralinv,
281-
&T_topinv,
282-
&T_bottominv
283-
};
196+
// Loop over elements
197+
for (c = cStart; c < cEnd; c++) {
198+
PetscInt numindices, *indices, i;
199+
ierr = DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE,
200+
&numindices, &indices, NULL, NULL);
201+
CHKERRQ(ierr);
202+
for (i = 0; i < numindices; i += ncomp) {
203+
for (PetscInt j = 0; j < ncomp; j++) {
204+
if (indices[i+j] != indices[i] + (PetscInt)(copysign(j, indices[i])))
205+
SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP,
206+
"Cell %D closure indices not interlaced", c);
207+
}
208+
// NO BC on closed surfaces
209+
PetscInt loc = indices[i];
284210

285-
for (PetscInt e = 0; e < nedgenodes; e++) {
286-
if (edgenodes[e].idx == i) {
287-
const PetscScalar (*matrixA)[2][2] = inv_transforms[edgenodes[e].panelA];
288-
const PetscScalar (*matrixB)[2][2] = transforms[edgenodes[e].panelB];
289-
290-
// inv_transform * transform (A^{-1}*B)
291-
// This product represents the mapping from coordinate system A
292-
// to spherical coordinates and then to coordinate system B. Vice versa
293-
// for transform * inv_transform (B^{-1}*A)
294-
PetscScalar matrixAB[2][2], matrixBA[2][2];
295-
for (int j=0; j<2; j++) {
296-
for (int k=0; k<2; k++) {
297-
matrixAB[j][k] = matrixBA[j][k] = 0;
298-
for (int l=0; l<2; l++) {
299-
matrixAB[j][k] += (*matrixA)[j][l] * (*matrixB)[l][k];
300-
matrixBA[j][k] += (*matrixB)[j][l] * (*matrixA)[l][k];
301-
}
302-
}
211+
// Query element panel
212+
PetscInt panel;
213+
ierr = DMGetLabelValue(dm, "panel", c, &panel); CHKERRQ(ierr);
214+
215+
// Query if edge node
216+
PetscBool isEdgeNode = PETSC_FALSE;
217+
218+
for (PetscInt e = 0; e < nedgenodes; e++) {
219+
if (edgenodes[e].idx == loc) {
220+
isEdgeNode = PETSC_TRUE;
221+
break;
303222
}
304-
PetscInt idxAB[2] = {4*i/ncomp + 0, 4*i/ncomp +1};
305-
PetscInt idxBA[2] = {4*i/ncomp + 2, 4*i/ncomp +3};
306-
ierr = MatSetValues(*T, 2, idxAB, 2, idxAB, (PetscScalar *)matrixAB,
307-
INSERT_VALUES); CHKERRQ(ierr);
308-
ierr = MatSetValues(*T, 2, idxBA, 2, idxBA, (PetscScalar *)matrixBA,
309-
INSERT_VALUES); CHKERRQ(ierr);
310223
}
224+
225+
PetscScalar entries[9];
226+
if (isEdgeNode) {
227+
ierr = VecGetArrayRead(X, &xarray); CHKERRQ(ierr);
228+
PetscScalar R = phys_ctx->R;
229+
// Read global Cartesian coordinates
230+
PetscScalar x[3] = {xarray[loc + 0],
231+
xarray[loc + 1],
232+
xarray[loc + 2]
233+
};
234+
// Restore array read
235+
ierr = VecRestoreArrayRead(X, &xarray); CHKERRQ(ierr);
236+
237+
// Normalize quadrature point coordinates to sphere
238+
PetscScalar rad = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
239+
x[0] *= R / rad;
240+
x[1] *= R / rad;
241+
x[2] *= R / rad;
242+
// Compute latitude and longitude
243+
const PetscScalar theta = asin(x[2] / R); // latitude
244+
const PetscScalar lambda = atan2(x[1], x[0]); // longitude
245+
246+
PetscScalar T00, T01, T10, T11;
247+
248+
switch (panel) {
249+
case 0:
250+
case 1:
251+
case 2:
252+
case 3:
253+
// For P_0 (front), P_1 (east), P_2 (back), P_3 (west):
254+
T00 = cos(theta)*cos(lambda) * cos(lambda);
255+
T01 = cos(theta)*cos(lambda) * 0.;
256+
T10 = cos(theta)*cos(lambda) * -sin(theta)*sin(lambda);
257+
T11 = cos(theta)*cos(lambda) * cos(theta);
258+
break;
259+
case 4:
260+
// For P4 (north):
261+
T00 = sin(theta) * cos(lambda);
262+
T01 = sin(theta) * sin(lambda);
263+
T10 = sin(theta) * -sin(theta)*sin(lambda);
264+
T11 = sin(theta) * sin(theta)*cos(lambda);
265+
break;
266+
case 5:
267+
// For P5 (south):
268+
T00 = sin(theta) * (-cos(theta));
269+
T01 = sin(theta) * sin(lambda);
270+
T10 = sin(theta) * sin(theta)*sin(lambda);
271+
T11 = sin(theta) * sin(theta)*cos(lambda);
272+
break;
273+
}
274+
275+
entries[0] = T00;
276+
entries[1] = T01;
277+
entries[2] = 0.;
278+
entries[3] = T10;
279+
entries[4] = T11;
280+
entries[5] = 0.;
281+
entries[6] = 0.;
282+
entries[7] = 0.;
283+
entries[8] = 1.;
284+
285+
} else {
286+
entries[0] = 1.;
287+
entries[1] = 0.;
288+
entries[2] = 0.;
289+
entries[3] = 0.;
290+
entries[4] = 1.;
291+
entries[5] = 0.;
292+
entries[6] = 0.;
293+
entries[7] = 0.;
294+
entries[8] = 1.;
295+
}
296+
PetscInt elem = c - cStart;
297+
PetscInt idxrow[3] = {elem*loc*ncomp + 0, elem*loc*ncomp + 1,
298+
elem*loc*ncomp + 2};
299+
PetscInt idxcol[3] = {loc*ncomp + 0, loc*ncomp + 1, loc*ncomp + 2};
300+
ierr = MatSetValues(*T, ncomp, idxrow, ncomp, idxcol, entries,
301+
INSERT_VALUES); CHKERRQ(ierr);
311302
}
303+
ierr = DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE,
304+
&numindices, &indices, NULL, NULL);
305+
CHKERRQ(ierr);
312306
}
307+
313308
// Assemble matrix for all node transformations
314309
ierr = MatAssemblyBegin(*T, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
315310
ierr = MatAssemblyEnd(*T, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
316311

317-
// Restore array read
318-
ierr = VecRestoreArrayRead(X, &xarray); CHKERRQ(ierr);
319-
320312
PetscFunctionReturn(0);
321313
}
322314

@@ -377,9 +369,7 @@ PetscErrorCode TransformCoords(DM dm, Vec Xloc, const PetscInt ncompx,
377369
xpanelslocarray[n*ncompx+1] = lambda;
378370
xpanelslocarray[n*ncompx+2] = -1;
379371
}
380-
381372
else {
382-
383373
switch (panel) {
384374
case 0:
385375
x = l * (Y / X);
@@ -680,9 +670,6 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra,
680670
CHKERRQ(ierr);
681671

682672
// CEED restrictions
683-
// TODO: After the resitrction matrix is implemented comment the following line
684-
ierr = CreateRestrictionPlex(ceed, dmcoord, numP, ncompq, &Erestrictq);
685-
CHKERRQ(ierr);
686673
ierr = CreateRestrictionPlex(ceed, dmcoord, 2, ncompx, &Erestrictx);
687674
CHKERRQ(ierr);
688675

@@ -693,11 +680,9 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra,
693680
CEED_STRIDES_BACKEND, &Erestrictqdi);
694681
// Solution vec restriction is a strided (identity) because we use a user
695682
// mat mult before and after operator apply
696-
// TODO: After the restriction matrix for the solution vector is implemented,
697-
// decomment the following line
698-
// CeedElemRestrictionCreateStrided(ceed, nelem, numP*numP, ncompq,
699-
// ncompq*nelem*numP*numP,
700-
// CEED_STRIDES_BACKEND, &Erestrictq);
683+
CeedElemRestrictionCreateStrided(ceed, nelem, numP*numP, ncompq,
684+
ncompq*nelem*numP*numP,
685+
CEED_STRIDES_BACKEND, &Erestrictq);
701686

702687
// Element coordinates
703688
ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr);

0 commit comments

Comments
 (0)