diff --git a/examples/fluids/shallow-water/shallowwater.c b/examples/fluids/shallow-water/shallowwater.c index 1e2a68e9a6..4a11f828fb 100644 --- a/examples/fluids/shallow-water/shallowwater.c +++ b/examples/fluids/shallow-water/shallowwater.c @@ -56,11 +56,11 @@ int main(int argc, char **argv) { problemData *problem = NULL; StabilizationType stab; PetscBool implicit; - PetscInt degree, qextra, outputfreq, steps, contsteps; + PetscInt degree, qextra, outputfreq, steps, contsteps, nedgenodes = 0; PetscMPIInt rank; PetscScalar ftime; - Vec Q, Qloc, Xloc; - const CeedInt ncompx = 3; + Vec Q, Qloc, Xloc, Xpanelsloc; + const PetscInt ncompx = 3; PetscInt viz_refine = 0; PetscBool read_mesh, simplex, test; PetscInt topodim = 2, ncompq = 3, lnodes; @@ -196,7 +196,7 @@ int main(int argc, char **argv) { g *= mpersquareds; // Set up the libCEED context - PhysicsContext_s phys_ctx = { + PhysicsContext_s physCtxData = { .u0 = 0., .v0 = 0., .h0 = .1, @@ -208,7 +208,7 @@ int main(int argc, char **argv) { .time = 0. }; - ProblemContext_s probl_ctx = { + ProblemContext_s problCtxData = { .g = g, .H0 = H0, .CtauS = CtauS, @@ -216,7 +216,7 @@ int main(int argc, char **argv) { .stabilization = stab }; - // Setup DM + // Create DM if (read_mesh) { ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, &dm); CHKERRQ(ierr); @@ -224,11 +224,11 @@ int main(int argc, char **argv) { // Create the mesh as a 0-refined sphere. This will create a cubic surface, not a box. PetscBool simplex = PETSC_FALSE; ierr = DMPlexCreateSphereMesh(PETSC_COMM_WORLD, topodim, simplex, - phys_ctx.R, &dm); + physCtxData.R, &dm); CHKERRQ(ierr); // Set the object name ierr = PetscObjectSetName((PetscObject)dm, "Sphere"); CHKERRQ(ierr); - // Define cube panels (charts) + // Define cube panels DMLabel label; PetscInt c, cStart, cEnd, npanel, permidx[6] = {5, 1, 4, 0, 3, 2}; ierr = DMCreateLabel(dm, "panel"); @@ -238,6 +238,7 @@ int main(int argc, char **argv) { for (c = cStart, npanel = 0; c < cEnd; c++) { ierr = DMLabelSetValue(label, c, permidx[npanel++]); CHKERRQ(ierr); } + ierr = DMPlexLabelComplete(dm, label); CHKERRQ(ierr); // Distribute mesh over processes { DM dmDist = NULL; @@ -259,11 +260,8 @@ int main(int argc, char **argv) { ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr); } - // Create DM - ierr = DMLocalizeCoordinates(dm); CHKERRQ(ierr); - ierr = DMSetFromOptions(dm); CHKERRQ(ierr); - ierr = SetupDM(dm, degree, ncompq, topodim); - CHKERRQ(ierr); + // Setup DM + ierr = SetupDMByDegree(dm, degree, ncompq, topodim); CHKERRQ(ierr); // Refine DM for high-order viz dmviz = NULL; @@ -283,7 +281,8 @@ int main(int argc, char **argv) { ierr = DMSetCoarseDM(dmhierarchy[i+1], dmhierarchy[i]); CHKERRQ(ierr); d = (d + 1) / 2; if (i + 1 == viz_refine) d = 1; - ierr = SetupDM(dmhierarchy[i+1], degree, ncompq, topodim); CHKERRQ(ierr); + ierr = SetupDMByDegree(dmhierarchy[i+1], degree, ncompq, topodim); + CHKERRQ(ierr); ierr = DMCreateInterpolation(dmhierarchy[i], dmhierarchy[i+1], &interp_next, NULL); CHKERRQ(ierr); if (!i) interpviz = interp_next; @@ -296,7 +295,7 @@ int main(int argc, char **argv) { interpviz = C; } } - for (PetscInt i=1; iM); CHKERRQ(ierr); - // Setup global lat-long vector for different panels (charts) of the cube + // Setup transformation matrix different panels of the cube Mat T; - ierr = FindPanelEdgeNodes(dm, &phys_ctx, ncompq, &T); CHKERRQ(ierr); + EdgeNode edgenodes; + ierr = FindPanelEdgeNodes(dm, &physCtxData, ncompq, &nedgenodes, &edgenodes, + &T); + CHKERRQ(ierr); + + // Pass transformation matrix to context for setup operator + problCtxData.T = T; + + ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); + ierr = VecDuplicate(Xloc, &Xpanelsloc); CHKERRQ(ierr); + // Transform coordinate according to their panel systems + ierr = TransformCoords(dm, Xloc, ncompx, edgenodes, nedgenodes, &physCtxData, + &Xpanelsloc); CHKERRQ(ierr); + + // Free edgenodes structure array + ierr = PetscFree(edgenodes); CHKERRQ(ierr); // Setup libCEED's objects ierr = PetscMalloc1(1, &ceeddata); CHKERRQ(ierr); - ierr = SetupLibceed(dm, ceed, degree, qextra, ncompx, ncompq, user, ceeddata, - problem, &phys_ctx, &probl_ctx); CHKERRQ(ierr); + ierr = SetupLibceed(dm, ceed, degree, qextra, ncompx, ncompq, user, + Xpanelsloc, ceeddata, problem, &physCtxData, &problCtxData); + CHKERRQ(ierr); // Set up PETSc context // Set up units structure @@ -402,9 +417,11 @@ int main(int argc, char **argv) { ierr = VecZeroEntries(Q); CHKERRQ(ierr); ierr = VectorPlacePetscVec(user->q0ceed, Qloc); CHKERRQ(ierr); + // Set up contex for QFunctions + CeedQFunctionSetContext(ceeddata->qf_setup, ceeddata->physCtx); + // Apply Setup Ceed Operators - ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); - ierr = VectorPlacePetscVec(ceeddata->xcorners, Xloc); CHKERRQ(ierr); + ierr = VectorPlacePetscVec(ceeddata->xcorners, Xpanelsloc); CHKERRQ(ierr); CeedOperatorApply(ceeddata->op_setup, ceeddata->xcorners, ceeddata->qdata, CEED_REQUEST_IMMEDIATE); ierr = ComputeLumpedMassMatrix(ceed, dm, ceeddata->Erestrictq, @@ -414,7 +431,7 @@ int main(int argc, char **argv) { // Apply IC operator and fix multiplicity of initial state vector ierr = ICs_FixMultiplicity(ceeddata->op_ics, ceeddata->xcorners, user->q0ceed, dm, Qloc, Q, ceeddata->Erestrictq, - &phys_ctx, 0.0); CHKERRQ(ierr); + &physCtxData, 0.0); CHKERRQ(ierr); MPI_Comm_rank(comm, &rank); if (!rank) { @@ -524,6 +541,8 @@ int main(int argc, char **argv) { CeedQFunctionDestroy(&ceeddata->qf_explicit); CeedQFunctionDestroy(&ceeddata->qf_implicit); CeedQFunctionDestroy(&ceeddata->qf_jacobian); + CeedQFunctionContextDestroy(&ceeddata->physCtx); + CeedQFunctionContextDestroy(&ceeddata->problCtx); CeedOperatorDestroy(&ceeddata->op_setup); CeedOperatorDestroy(&ceeddata->op_ics); CeedOperatorDestroy(&ceeddata->op_explicit); diff --git a/examples/fluids/shallow-water/src/setup.c b/examples/fluids/shallow-water/src/setup.c index 9dc30ae836..d59087931e 100644 --- a/examples/fluids/shallow-water/src/setup.c +++ b/examples/fluids/shallow-water/src/setup.c @@ -27,7 +27,7 @@ #include "../sw_headers.h" // Function prototytes #include "../qfunctions/setup_geo.h" // Geometric factors #include "../qfunctions/advection.h" // Physics point-wise functions -#include "../qfunctions/geostrophic.h" // Physics point-wise functions +#include "../qfunctions/geostrophic.h" // Physics point-wise functions #if PETSC_VERSION_LT(3,14,0) # define DMPlexGetClosureIndices(a,b,c,d,e,f,g,h,i) DMPlexGetClosureIndices(a,b,c,d,f,g,i) @@ -76,20 +76,19 @@ problemData problemOptions[] = { // ----------------------------------------------------------------------------- PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx, - PetscInt ncomp, Mat *T) { + PetscInt ncomp, PetscInt *edgenodecnt, + EdgeNode *edgenodes, Mat *T) { PetscInt ierr; MPI_Comm comm; - PetscInt cstart, cend, gdofs, depth, edgenodecnt = 0; + PetscInt cstart, cend, gdofs; PetscSection section, sectionloc; - EdgeNode edgenodes; PetscFunctionBeginUser; // Get Nelem ierr = DMGetGlobalSection(dm, §ion); CHKERRQ(ierr); ierr = DMGetLocalSection(dm, §ionloc); CHKERRQ(ierr); ierr = DMPlexGetHeightStratum(dm, 0, &cstart, &cend); CHKERRQ(ierr); - ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); PetscSF sf; Vec bitmapVec; @@ -126,7 +125,7 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx, MPI_Reduce(bitmaploc, bitmap, gdofs, MPI_UNSIGNED, MPI_BOR, 0, comm); // Read the resulting bitmap and extract edge nodes only - ierr = PetscMalloc1(gdofs + 24*ncomp, &edgenodes); CHKERRQ(ierr); + ierr = PetscMalloc1(gdofs + 24*ncomp, edgenodes); CHKERRQ(ierr); for (PetscInt i = 0; i < gdofs; i += ncomp) { PetscInt ones = 0, panels[3]; for (PetscInt p = 0; p < 6; p++) { @@ -135,31 +134,29 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx, bitmap[i] >>= 1; } if (ones == 2) { - edgenodes[edgenodecnt].idx = i; - edgenodes[edgenodecnt].panelA = panels[0]; - edgenodes[edgenodecnt].panelB = panels[1]; - edgenodecnt++; + (*edgenodes)[*edgenodecnt].idx = i; + (*edgenodes)[*edgenodecnt].panelA = panels[0]; + (*edgenodes)[*edgenodecnt].panelB = panels[1]; + (*edgenodecnt)++; } else if (ones == 3) { - edgenodes[edgenodecnt].idx = i; - edgenodes[edgenodecnt].panelA = panels[0]; - edgenodes[edgenodecnt].panelB = panels[1]; - edgenodecnt++; - edgenodes[edgenodecnt].idx = i; - edgenodes[edgenodecnt].panelA = panels[0]; - edgenodes[edgenodecnt].panelB = panels[2]; - edgenodecnt++; - edgenodes[edgenodecnt].idx = i; - edgenodes[edgenodecnt].panelA = panels[1]; - edgenodes[edgenodecnt].panelB = panels[2]; - edgenodecnt++; + (*edgenodes)[*edgenodecnt].idx = i; + (*edgenodes)[*edgenodecnt].panelA = panels[0]; + (*edgenodes)[*edgenodecnt].panelB = panels[1]; + (*edgenodecnt)++; + (*edgenodes)[*edgenodecnt].idx = i; + (*edgenodes)[*edgenodecnt].panelA = panels[0]; + (*edgenodes)[*edgenodecnt].panelB = panels[2]; + (*edgenodecnt)++; + (*edgenodes)[*edgenodecnt].idx = i; + (*edgenodes)[*edgenodecnt].panelA = panels[1]; + (*edgenodes)[*edgenodecnt].panelB = panels[2]; + (*edgenodecnt)++; } } - ierr = SetupPanelCoordTransformations(dm, phys_ctx, ncomp, edgenodes, - edgenodecnt, T); CHKERRQ(ierr); + ierr = SetupPanelCoordTransformations(dm, phys_ctx, ncomp, *edgenodes, + *edgenodecnt, T); CHKERRQ(ierr); - // Free edgenodes structure array - ierr = PetscFree(edgenodes); CHKERRQ(ierr); // Free heap ierr = PetscFree(bitmap); CHKERRQ(ierr); ierr = PetscFree(bitmaploc); CHKERRQ(ierr); @@ -168,7 +165,7 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx, } // ----------------------------------------------------------------------------- -// Auxiliary function that sets up all corrdinate transformations between panels +// Auxiliary function that sets up all coordinate transformations between panels // ----------------------------------------------------------------------------- PetscErrorCode SetupPanelCoordTransformations(DM dm, PhysicsContext phys_ctx, PetscInt ncomp, @@ -324,6 +321,103 @@ PetscErrorCode SetupPanelCoordTransformations(DM dm, PhysicsContext phys_ctx, PetscFunctionReturn(0); } +// ----------------------------------------------------------------------------- +// Auxiliary function that converts global 3D coors into local panel coords +// ----------------------------------------------------------------------------- + +PetscErrorCode TransformCoords(DM dm, Vec Xloc, const PetscInt ncompx, + EdgeNode edgenodes, const PetscInt nedgenodes, + PhysicsContext phys_ctx, Vec *Xpanelsloc) { + + PetscInt ierr; + PetscInt lsize, depth, nstart, nend; + const PetscScalar *xlocarray; + PetscScalar *xpanelslocarray; + + PetscFunctionBeginUser; + ierr = VecGetArrayRead(Xloc, &xlocarray); CHKERRQ(ierr); + ierr = VecGetSize(Xloc, &lsize); CHKERRQ(ierr); + ierr = VecGetArray(*Xpanelsloc, &xpanelslocarray); CHKERRQ(ierr); + + ierr = DMPlexGetDepth(dm, &depth); CHKERRQ(ierr); + ierr = DMPlexGetHeightStratum(dm, depth, &nstart, &nend); CHKERRQ(ierr); + + for (PetscInt n = 0; n < (nend - nstart); n++) { // Traverse nodes + // Query element panel + PetscInt panel; + PetscSection sectionloc; + ierr = DMGetLocalSection(dm, §ionloc); CHKERRQ(ierr); + ierr = DMGetLabelValue(dm, "panel", n + nstart, &panel); CHKERRQ(ierr); + + PetscScalar X = xlocarray[n*ncompx+0]; + PetscScalar Y = xlocarray[n*ncompx+1]; + PetscScalar Z = xlocarray[n*ncompx+2]; + + // Normalize quadrature point coordinates to sphere + CeedScalar rad = sqrt(X*X + Y*Y + Z*Z); + X *= phys_ctx->R / rad; + Y *= phys_ctx->R / rad; + Z *= phys_ctx->R / rad; + + PetscScalar x, y; + PetscScalar l = phys_ctx->R/sqrt(3.); + PetscBool isedgenode = PETSC_FALSE; + + // Determine if this node is an edge node + for (PetscInt e = 0; e < nedgenodes; e++) { + if (n * ncompx == edgenodes[e].idx) { + isedgenode = PETSC_TRUE; break; + } + } + + if (isedgenode) { + // Compute latitude and longitude + const CeedScalar theta = asin(Z / phys_ctx->R); // latitude + const CeedScalar lambda = atan2(Y, X); // longitude + xpanelslocarray[n*ncompx+0] = theta; + xpanelslocarray[n*ncompx+1] = lambda; + xpanelslocarray[n*ncompx+2] = -1; + } + + else { + + switch (panel) { + case 0: + x = l * (Y / X); + y = l * (Z / X); + break; + case 1: + x = l * (-X / Y); + y = l * (Z / Y); + break; + case 2: + x = l * (Y / X); + y = l * (-Z / X); + break; + case 3: + x = l * (-X / Y); + y = l * (-Z / Y); + break; + case 4: + x = l * (Y / Z); + y = l * (-X / Z); + break; + case 5: + x = l * (-Y / Z); + y = l * (-X / Z); + break; + } + xpanelslocarray[n*ncompx+0] = x; + xpanelslocarray[n*ncompx+1] = y; + xpanelslocarray[n*ncompx+2] = panel; + } + } // End of nodes for loop + + ierr = VecRestoreArrayRead(Xloc, &xlocarray); CHKERRQ(ierr); + ierr = VecRestoreArray(*Xpanelsloc, &xpanelslocarray); CHKERRQ(ierr); + + PetscFunctionReturn(0); +} // ----------------------------------------------------------------------------- // Auxiliary function to create PETSc FE space for a given degree @@ -400,8 +494,8 @@ PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, // Auxiliary function to setup DM FE space and info // ----------------------------------------------------------------------------- -PetscErrorCode SetupDM(DM dm, PetscInt degree, PetscInt ncompq, - PetscInt dim) { +PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompq, + PetscInt dim) { PetscErrorCode ierr; PetscFunctionBeginUser; @@ -558,11 +652,11 @@ PetscErrorCode VectorPlacePetscVec(CeedVector c, Vec p) { PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, PetscInt ncompx, PetscInt ncompq, User user, - CeedData data, problemData *problem, - PhysicsContext phys_ctx, ProblemContext probl_ctx) { + Vec Xloc, CeedData data, problemData *problem, + PhysicsContext physCtxData, + ProblemContext problCtxData) { int ierr; DM dmcoord; - Vec Xloc; CeedBasis basisx, basisxc, basisq; CeedElemRestriction Erestrictx, Erestrictq, Erestrictqdi; CeedQFunction qf_setup, qf_ics, qf_explicit, qf_implicit, @@ -729,10 +823,16 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, user->op_jacobian = op_jacobian; // Set up the libCEED context - CeedQFunctionSetContext(qf_ics, phys_ctx, sizeof *phys_ctx); - CeedQFunctionSetContext(qf_explicit, probl_ctx, sizeof *probl_ctx); - CeedQFunctionSetContext(qf_implicit, probl_ctx, sizeof *probl_ctx); - CeedQFunctionSetContext(qf_jacobian, probl_ctx, sizeof *probl_ctx); + CeedQFunctionContextCreate(ceed, &data->physCtx); + CeedQFunctionContextSetData(data->physCtx, CEED_MEM_HOST, CEED_USE_POINTER, + sizeof physCtxData, &physCtxData); + CeedQFunctionSetContext(qf_ics, data->physCtx); + CeedQFunctionContextCreate(ceed, &data->problCtx); + CeedQFunctionContextSetData(data->problCtx, CEED_MEM_HOST, CEED_USE_POINTER, + sizeof problCtxData, &problCtxData); + CeedQFunctionSetContext(qf_explicit, data->problCtx); + CeedQFunctionSetContext(qf_implicit, data->problCtx); + CeedQFunctionSetContext(qf_jacobian, data->problCtx); // Save libCEED data required for level // TODO: check how many of these are really needed outside data->basisx = basisx; diff --git a/examples/fluids/shallow-water/sw_headers.h b/examples/fluids/shallow-water/sw_headers.h index f0ad6489d5..58c4a00668 100644 --- a/examples/fluids/shallow-water/sw_headers.h +++ b/examples/fluids/shallow-water/sw_headers.h @@ -47,6 +47,7 @@ typedef struct { CeedScalar CtauS; CeedScalar strong_form; int stabilization; // See StabilizationType: 0=none, 1=SU, 2=SUPG + Mat T; } ProblemContext_s; typedef ProblemContext_s *ProblemContext; @@ -131,6 +132,7 @@ struct CeedData_ { CeedElemRestriction Erestrictx, Erestrictq, Erestrictqdi; CeedQFunction qf_setup, qf_mass, qf_ics, qf_explicit, qf_implicit, qf_jacobian; + CeedQFunctionContext physCtx, problCtx; CeedOperator op_setup, op_mass, op_ics, op_explicit, op_implicit, op_jacobian; CeedVector xcorners, xceed, qdata, q0ceed, mceed, hsceed, H0ceed; }; @@ -144,14 +146,20 @@ extern problemData problemOptions[]; // Auxiliary function to determine if nodes belong to cube face (panel) edges PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx, - PetscInt ncomp, Mat *T); + PetscInt ncomp, PetscInt *edgenodecnt, + EdgeNode *edgenodes, Mat *T); -// Auxiliary function that sets up all corrdinate transformations between panels +// Auxiliary function that sets up all coordinate transformations between panels PetscErrorCode SetupPanelCoordTransformations(DM dm, PhysicsContext phys_ctx, PetscInt ncomp, EdgeNode edgenodes, PetscInt nedgenodes, Mat *T); +// Auxiliary function that converts global 3D coors into local panel coords +PetscErrorCode TransformCoords(DM dm, Vec Xloc, const PetscInt ncompx, + EdgeNode edgenodes, const PetscInt nedgenodes, + PhysicsContext phys_ctx, Vec *Xpanelsloc); + // ----------------------------------------------------------------------------- // Setup DM functions // ----------------------------------------------------------------------------- @@ -162,7 +170,8 @@ PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, PetscInt order, PetscFE *fem); // Auxiliary function to setup DM FE space and info -PetscErrorCode SetupDM(DM dm, PetscInt degree, PetscInt ncompq, PetscInt dim); +PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompq, + PetscInt dim); // ----------------------------------------------------------------------------- // libCEED functions @@ -174,9 +183,10 @@ PetscErrorCode CreateRestrictionPlex(Ceed ceed, DM dm, CeedInt P, CeedInt ncomp, // Auxiliary function to set up libCEED objects for a given degree PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, - PetscInt ncompx, PetscInt ncompq, User user, - CeedData data, problemData *problem, - PhysicsContext phys_ctx, ProblemContext probl_ctx); + const PetscInt ncompx, PetscInt ncompq, User user, + Vec Xloc, CeedData data, problemData *problem, + PhysicsContext physCtxData, + ProblemContext problCtxData); // ----------------------------------------------------------------------------- // RHS (Explicit part in time-stepper) function setup