Skip to content

Commit

Permalink
WIP SWE: pass sparse matrix to user context to be used in wrapper fun…
Browse files Browse the repository at this point in the history 
…ctions
  • Loading branch information
@valeriabarra
valeriabarra committed Sep 3, 2020
1 parent 4afe695 commit 123df56
Show file tree
Hide file tree
Showing 4 changed files with 53 additions and 23 deletions.
23 changes: 10 additions & 13 deletions examples/fluids/shallow-water/shallowwater.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -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,
Expand All @@ -208,7 +208,7 @@ int main(int argc, char **argv) {
.time = 0.
};

ProblemContext_s problCtxData = {
ProblemContext_s problCtxData = {
.g = g,
.H0 = H0,
.CtauS = CtauS,
Expand Down Expand Up @@ -378,22 +378,17 @@ 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);

// 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
Expand All @@ -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
Expand All @@ -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,
Expand Down
23 changes: 16 additions & 7 deletions examples/fluids/shallow-water/src/setup.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -494,18 +494,18 @@ 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;

PetscFunctionBeginUser;
{
// 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);
Expand Down Expand Up @@ -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,
Expand All @@ -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);
Expand All @@ -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);
Expand Down
25 changes: 25 additions & 0 deletions examples/fluids/shallow-water/src/wrappers.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down
5 changes: 2 additions & 3 deletions examples/fluids/shallow-water/sw_headers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);

Expand Down

0 comments on commit 123df56

Please sign in to comment.