diff --git a/examples/fluids/shallow-water/shallowwater.c b/examples/fluids/shallow-water/shallowwater.c index 4a11f828fb..bfb661f720 100644 --- a/examples/fluids/shallow-water/shallowwater.c +++ b/examples/fluids/shallow-water/shallowwater.c @@ -59,11 +59,11 @@ int main(int argc, char **argv) { PetscInt degree, qextra, outputfreq, steps, contsteps, nedgenodes = 0; PetscMPIInt rank; PetscScalar ftime; - Vec Q, Qloc, Xloc, Xpanelsloc; + Vec Q, Qloc, Xloc; const PetscInt ncompx = 3; PetscInt viz_refine = 0; PetscBool read_mesh, simplex, test; - PetscInt topodim = 2, ncompq = 3, lnodes; + PetscInt topodim = 2, dim = 3, ncompq = 3, lnodes; // libCEED context char ceedresource[PETSC_MAX_PATH_LEN] = "/cpu/self", filename[PETSC_MAX_PATH_LEN]; @@ -196,7 +196,7 @@ int main(int argc, char **argv) { g *= mpersquareds; // Set up the libCEED context - PhysicsContext_s physCtxData = { + PhysicsContext_s physCtxData = { .u0 = 0., .v0 = 0., .h0 = .1, @@ -208,7 +208,7 @@ int main(int argc, char **argv) { .time = 0. }; - ProblemContext_s problCtxData = { + ProblemContext_s problCtxData = { .g = g, .H0 = H0, .CtauS = CtauS, @@ -378,14 +378,9 @@ int main(int argc, char **argv) { &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); +// ierr = TransformCoords(dm, Xloc, ncompx, edgenodes, nedgenodes, &physCtxData, +// &Xpanelsloc); CHKERRQ(ierr); // Free edgenodes structure array ierr = PetscFree(edgenodes); CHKERRQ(ierr); @@ -393,7 +388,7 @@ int main(int argc, char **argv) { // Setup libCEED's objects ierr = PetscMalloc1(1, &ceeddata); CHKERRQ(ierr); ierr = SetupLibceed(dm, ceed, degree, qextra, ncompx, ncompq, user, - Xpanelsloc, ceeddata, problem, &physCtxData, &problCtxData); + ceeddata, problem, &physCtxData, &problCtxData); CHKERRQ(ierr); // Set up PETSc context @@ -411,6 +406,7 @@ int main(int argc, char **argv) { user->dmviz = dmviz; user->interpviz = interpviz; user->ceed = ceed; + user->T = T; // Calculate qdata and ICs // Set up state global and local vectors @@ -421,7 +417,8 @@ int main(int argc, char **argv) { CeedQFunctionSetContext(ceeddata->qf_setup, ceeddata->physCtx); // Apply Setup Ceed Operators - ierr = VectorPlacePetscVec(ceeddata->xcorners, Xpanelsloc); CHKERRQ(ierr); + ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); + ierr = VectorPlacePetscVec(ceeddata->xcorners, Xloc); CHKERRQ(ierr); CeedOperatorApply(ceeddata->op_setup, ceeddata->xcorners, ceeddata->qdata, CEED_REQUEST_IMMEDIATE); ierr = ComputeLumpedMassMatrix(ceed, dm, ceeddata->Erestrictq, diff --git a/examples/fluids/shallow-water/src/setup.c b/examples/fluids/shallow-water/src/setup.c index d59087931e..0d542de45c 100644 --- a/examples/fluids/shallow-water/src/setup.c +++ b/examples/fluids/shallow-water/src/setup.c @@ -154,8 +154,8 @@ PetscErrorCode FindPanelEdgeNodes(DM dm, PhysicsContext phys_ctx, (*edgenodecnt)++; } } - ierr = SetupPanelCoordTransformations(dm, phys_ctx, ncomp, *edgenodes, - *edgenodecnt, T); CHKERRQ(ierr); +// ierr = SetupPanelCoordTransformations(dm, phys_ctx, ncomp, *edgenodes, +// *edgenodecnt, T); CHKERRQ(ierr); // Free heap ierr = PetscFree(bitmap); CHKERRQ(ierr); @@ -494,7 +494,7 @@ PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, // Auxiliary function to setup DM FE space and info // ----------------------------------------------------------------------------- -PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompq, +PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncomp, PetscInt dim) { PetscErrorCode ierr; @@ -502,10 +502,10 @@ PetscErrorCode SetupDMByDegree(DM dm, PetscInt degree, PetscInt ncompq, { // Configure the finite element space PetscFE fe; - ierr = PetscFECreateByDegree(dm, dim, ncompq, PETSC_FALSE, NULL, degree, + ierr = PetscFECreateByDegree(dm, dim, ncomp, PETSC_FALSE, NULL, degree, &fe); ierr = PetscObjectSetName((PetscObject)fe, "Q"); CHKERRQ(ierr); - ierr = DMAddField(dm,NULL,(PetscObject)fe); CHKERRQ(ierr); + ierr = DMAddField(dm, NULL, (PetscObject)fe); CHKERRQ(ierr); ierr = DMCreateDS(dm); CHKERRQ(ierr); ierr = DMPlexSetClosurePermutationTensor(dm, PETSC_DETERMINE, NULL); @@ -652,11 +652,12 @@ PetscErrorCode VectorPlacePetscVec(CeedVector c, Vec p) { PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, PetscInt ncompx, PetscInt ncompq, User user, - Vec Xloc, CeedData data, problemData *problem, + 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, @@ -682,7 +683,8 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, CHKERRQ(ierr); // CEED restrictions - ierr = CreateRestrictionPlex(ceed, dm, numP, ncompq, &Erestrictq); + // TODO: After the resitrction matrix is implemented comment the following line + ierr = CreateRestrictionPlex(ceed, dmcoord, numP, ncompq, &Erestrictq); CHKERRQ(ierr); ierr = CreateRestrictionPlex(ceed, dmcoord, 2, ncompx, &Erestrictx); CHKERRQ(ierr); @@ -692,6 +694,13 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, CeedInt degree, CeedInt qextra, CeedElemRestrictionCreateStrided(ceed, nelem, numQ*numQ, qdatasize, qdatasize*nelem*numQ*numQ, CEED_STRIDES_BACKEND, &Erestrictqdi); + // Solution vec restriction is a strided (identity) because we use a user + // mat mult before and after operator apply + // TODO: After the restriction matrix for the solution vector is implemented, + // decomment the following line +// CeedElemRestrictionCreateStrided(ceed, nelem, numP*numP, ncompq, +// ncompq*nelem*numP*numP, +// CEED_STRIDES_BACKEND, &Erestrictq); // Element coordinates ierr = DMGetCoordinatesLocal(dm, &Xloc); CHKERRQ(ierr); diff --git a/examples/fluids/shallow-water/src/wrappers.c b/examples/fluids/shallow-water/src/wrappers.c index 8104bd10c4..d5a283dd48 100644 --- a/examples/fluids/shallow-water/src/wrappers.c +++ b/examples/fluids/shallow-water/src/wrappers.c @@ -46,10 +46,18 @@ PetscErrorCode FormRHSFunction_SW(TS ts, PetscReal t, Vec Q, Vec G, CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); CeedVectorSetArray(user->gceed, CEED_MEM_HOST, CEED_USE_POINTER, g); + // TODO: + // L-vector to E-vector: + // Apply user-defined restriction with appropriate coordinate transformations + // Apply CEED operator CeedOperatorApply(user->op_explicit, user->qceed, user->gceed, CEED_REQUEST_IMMEDIATE); + // TODO: + // E-vector to L-vector: + // Apply transpose of user-defined restriction + // Restore vectors ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); ierr = VecRestoreArray(Gloc, &g); CHKERRQ(ierr); @@ -101,10 +109,19 @@ PetscErrorCode FormIFunction_SW(TS ts, PetscReal t, Vec Q, Vec Qdot, (PetscScalar *)qdot); CeedVectorSetArray(user->fceed, CEED_MEM_HOST, CEED_USE_POINTER, f); + // TODO: + // L-vector to E-vector: + // Apply user-defined restriction with appropriate coordinate transformations + + // Apply CEED operator CeedOperatorApply(user->op_implicit, user->qceed, user->fceed, CEED_REQUEST_IMMEDIATE); + // TODO: + // E-vector to L-vector: + // Apply transpose of user-defined restriction + // Restore vectors ierr = VecRestoreArrayRead(Qloc, &q); CHKERRQ(ierr); ierr = VecRestoreArrayRead(Qdotloc, &qdot); CHKERRQ(ierr); @@ -179,11 +196,19 @@ PetscErrorCode ApplyJacobian_SW(Mat mat, Vec Q, Vec JVec) { CeedVectorSetArray(user->qceed, CEED_MEM_HOST, CEED_USE_POINTER, q); CeedVectorSetArray(user->jceed, CEED_MEM_HOST, CEED_USE_POINTER, j); + // TODO: + // L-vector to E-vector: + // Apply user-defined restriction with appropriate coordinate transformations + // Apply the CEED operator for the dF/dQ terms CeedOperatorApply(user->op_jacobian, user->qceed, user->jceed, CEED_REQUEST_IMMEDIATE); CeedVectorSyncArray(user->jceed, CEED_MEM_HOST); + // TODO: + // E-vector to L-vector: + // Apply transpose of user-defined restriction + // Restore PETSc vectors ierr = VecRestoreArrayRead(Qloc, (const PetscScalar **)&q); CHKERRQ(ierr); diff --git a/examples/fluids/shallow-water/sw_headers.h b/examples/fluids/shallow-water/sw_headers.h index 58c4a00668..412fac62ee 100644 --- a/examples/fluids/shallow-water/sw_headers.h +++ b/examples/fluids/shallow-water/sw_headers.h @@ -47,7 +47,6 @@ 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; @@ -101,7 +100,7 @@ struct User_ { PetscInt outputfreq; DM dm; DM dmviz; - Mat interpviz; + Mat interpviz, T; Ceed ceed; Units units; CeedVector qceed, q0ceed, qdotceed, fceed, gceed, jceed; @@ -184,7 +183,7 @@ 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, const PetscInt ncompx, PetscInt ncompq, User user, - Vec Xloc, CeedData data, problemData *problem, + CeedData data, problemData *problem, PhysicsContext physCtxData, ProblemContext problCtxData);