From c54385ef8afc8639cc4f1a093cf7ded388cfccaa Mon Sep 17 00:00:00 2001 From: rezgarshakeri Date: Fri, 16 Dec 2022 09:56:18 -0700 Subject: [PATCH] some clean-up in Hdiv-mixed folder --- examples/Hdiv-mixed/include/post-processing.h | 6 +- examples/Hdiv-mixed/include/setup-dm.h | 2 +- examples/Hdiv-mixed/include/setup-fe.h | 7 + examples/Hdiv-mixed/include/setup-libceed.h | 12 +- examples/Hdiv-mixed/include/setup-solvers.h | 2 +- examples/Hdiv-mixed/include/setup-ts.h | 2 +- examples/Hdiv-mixed/main.c | 116 ++++++------- examples/Hdiv-mixed/src/cl-options.c | 2 +- examples/Hdiv-mixed/src/post-processing.c | 10 +- examples/Hdiv-mixed/src/setup-dm.c | 26 ++- examples/Hdiv-mixed/src/setup-fe.c | 158 +++++++++++++++++- examples/Hdiv-mixed/src/setup-libceed.c | 137 --------------- examples/Hdiv-mixed/src/setup-solvers.c | 46 ++--- examples/Hdiv-mixed/src/setup-ts.c | 4 +- 14 files changed, 277 insertions(+), 253 deletions(-) diff --git a/examples/Hdiv-mixed/include/post-processing.h b/examples/Hdiv-mixed/include/post-processing.h index 9c3506cfcf..5bd940be39 100644 --- a/examples/Hdiv-mixed/include/post-processing.h +++ b/examples/Hdiv-mixed/include/post-processing.h @@ -6,10 +6,10 @@ #include "../include/setup-libceed.h" #include "structs.h" -PetscErrorCode PrintOutput(DM dm, Ceed ceed, AppCtx app_ctx, PetscBool has_ts, CeedMemType mem_type_backend, TS ts, SNES snes, KSP ksp, Vec U, - CeedScalar l2_error_u, CeedScalar l2_error_p); +PetscErrorCode PrintOutput(DM dm, Ceed ceed, AppCtx app_ctx, PetscBool has_ts, TS ts, SNES snes, KSP ksp, Vec U, CeedScalar l2_error_u, + CeedScalar l2_error_p); PetscErrorCode SetupProjectVelocityCtx_Hdiv(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_Hdiv); -PetscErrorCode SetupProjectVelocityCtx_H1(MPI_Comm comm, DM dm_H1, Ceed ceed, CeedData ceed_data, VecType vec_type, OperatorApplyContext ctx_H1); +PetscErrorCode SetupProjectVelocityCtx_H1(MPI_Comm comm, DM dm_H1, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_H1); PetscErrorCode ProjectVelocity(AppCtx app_ctx, Vec U, Vec *U_H1); PetscErrorCode CtxVecDestroy(AppCtx app_ctx); #endif // post_processing_h diff --git a/examples/Hdiv-mixed/include/setup-dm.h b/examples/Hdiv-mixed/include/setup-dm.h index 71469cbd37..775b685b92 100644 --- a/examples/Hdiv-mixed/include/setup-dm.h +++ b/examples/Hdiv-mixed/include/setup-dm.h @@ -11,7 +11,7 @@ // --------------------------------------------------------------------------- // Setup DM // --------------------------------------------------------------------------- -PetscErrorCode CreateDM(MPI_Comm comm, MatType mat_type, VecType vec_type, DM *dm); +PetscErrorCode CreateDM(MPI_Comm comm, Ceed ceed, DM *dm); PetscErrorCode PerturbVerticesSmooth(DM dm); PetscErrorCode PerturbVerticesRandom(DM dm); diff --git a/examples/Hdiv-mixed/include/setup-fe.h b/examples/Hdiv-mixed/include/setup-fe.h index 382b047740..8836489f8b 100644 --- a/examples/Hdiv-mixed/include/setup-fe.h +++ b/examples/Hdiv-mixed/include/setup-fe.h @@ -11,6 +11,13 @@ // --------------------------------------------------------------------------- // Setup FE // --------------------------------------------------------------------------- +CeedMemType MemTypeP2C(PetscMemType mtype); PetscErrorCode SetupFEHdiv(MPI_Comm comm, DM dm, DM dm_u0, DM dm_p0); PetscErrorCode SetupFEH1(ProblemData problem_data, AppCtx app_ctx, DM dm_H1); +PetscInt Involute(PetscInt i); +PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr); +// Utility function to create local CEED Oriented restriction from DMPlex +PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, DM dm_u0, DM dm_p0, CeedInt P, CeedElemRestriction *elem_restr_u, + CeedElemRestriction *elem_restr_p, CeedElemRestriction *elem_restr_u0, + CeedElemRestriction *elem_restr_p0); #endif // setupfe_h diff --git a/examples/Hdiv-mixed/include/setup-libceed.h b/examples/Hdiv-mixed/include/setup-libceed.h index 3b509b4a65..d04c3b1fa5 100644 --- a/examples/Hdiv-mixed/include/setup-libceed.h +++ b/examples/Hdiv-mixed/include/setup-libceed.h @@ -1,20 +1,10 @@ #ifndef setuplibceed_h #define setuplibceed_h +#include "setup-fe.h" #include "structs.h" -// Convert PETSc MemType to libCEED MemType -CeedMemType MemTypeP2C(PetscMemType mtype); // Destroy libCEED objects PetscErrorCode CeedDataDestroy(CeedData ceed_data, ProblemData problem_data); -// Utility function - essential BC dofs are encoded in closure indices as -(i+1) -PetscInt Involute(PetscInt i); -// Utility function to create local CEED restriction from DMPlex -PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr); -// Utility function to create local CEED Oriented restriction from DMPlex -PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, DM dm_u0, DM dm_p0, CeedInt P, CeedElemRestriction *elem_restr_u, - CeedElemRestriction *elem_restr_p, CeedElemRestriction *elem_restr_u0, - CeedElemRestriction *elem_restr_p0); -// Set up libCEED for a given degree PetscErrorCode SetupLibceed(DM dm, DM dm_u0, DM dm_p0, DM dm_H1, Ceed ceed, AppCtx app_ctx, ProblemData problem_data, CeedData ceed_data); #endif // setuplibceed_h diff --git a/examples/Hdiv-mixed/include/setup-solvers.h b/examples/Hdiv-mixed/include/setup-solvers.h index 5502172408..36d50a210d 100644 --- a/examples/Hdiv-mixed/include/setup-solvers.h +++ b/examples/Hdiv-mixed/include/setup-solvers.h @@ -7,7 +7,7 @@ #include "petscvec.h" #include "structs.h" -PetscErrorCode SetupJacobianOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data, VecType vec_type, OperatorApplyContext ctx_jacobian); +PetscErrorCode SetupJacobianOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_jacobian); PetscErrorCode SetupResidualOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_residual); PetscErrorCode SetupErrorOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_error); PetscErrorCode ApplyMatOp(Mat A, Vec X, Vec Y); diff --git a/examples/Hdiv-mixed/include/setup-ts.h b/examples/Hdiv-mixed/include/setup-ts.h index d490ba9bab..9db7f31203 100644 --- a/examples/Hdiv-mixed/include/setup-ts.h +++ b/examples/Hdiv-mixed/include/setup-ts.h @@ -6,7 +6,7 @@ #include "structs.h" -PetscErrorCode CreateInitialConditions(CeedData ceed_data, AppCtx app_ctx, VecType vec_type, Vec U); +PetscErrorCode CreateInitialConditions(CeedData ceed_data, AppCtx app_ctx, Vec U); PetscErrorCode SetupResidualOperatorCtx_Ut(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_residual_ut); PetscErrorCode SetupResidualOperatorCtx_U0(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_initial_u0); PetscErrorCode SetupResidualOperatorCtx_P0(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_initial_p0); diff --git a/examples/Hdiv-mixed/main.c b/examples/Hdiv-mixed/main.c index a0df25b77e..c990c08e8f 100644 --- a/examples/Hdiv-mixed/main.c +++ b/examples/Hdiv-mixed/main.c @@ -41,6 +41,7 @@ int main(int argc, char **argv) { // Initialize PETSc // --------------------------------------------------------------------------- PetscCall(PetscInitialize(&argc, &argv, NULL, help)); + MPI_Comm comm = PETSC_COMM_WORLD; // --------------------------------------------------------------------------- // Create structs @@ -77,6 +78,12 @@ int main(int argc, char **argv) { // Context for post-processing app_ctx->ctx_Hdiv = ctx_Hdiv; app_ctx->ctx_H1 = ctx_H1; + app_ctx->comm = comm; + + // --------------------------------------------------------------------------- + // Process command line options + // --------------------------------------------------------------------------- + PetscCall(ProcessCommandLineOptions(app_ctx)); // --------------------------------------------------------------------------- // Initialize libCEED @@ -85,41 +92,20 @@ int main(int argc, char **argv) { Ceed ceed; CeedInit("/cpu/self/ref/serial", &ceed); // CeedInit(app_ctx->ceed_resource, &ceed); - CeedMemType mem_type_backend; - CeedGetPreferredMemType(ceed, &mem_type_backend); - - VecType vec_type = NULL; - MatType mat_type = NULL; - switch (mem_type_backend) { - case CEED_MEM_HOST: - vec_type = VECSTANDARD; - break; - case CEED_MEM_DEVICE: { - const char *resolved; - CeedGetResource(ceed, &resolved); - if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA; - else if (strstr(resolved, "/gpu/hip")) vec_type = VECKOKKOS; - else vec_type = VECSTANDARD; - } - } - if (strstr(vec_type, VECCUDA)) mat_type = MATAIJCUSPARSE; - else if (strstr(vec_type, VECKOKKOS)) mat_type = MATAIJKOKKOS; - else mat_type = MATAIJ; // -- Process general command line options - MPI_Comm comm = PETSC_COMM_WORLD; // --------------------------------------------------------------------------- // Create DM // --------------------------------------------------------------------------- DM dm, dm_u0, dm_p0, dm_H1; // DM for mixed problem - PetscCall(CreateDM(comm, mat_type, vec_type, &dm)); + PetscCall(CreateDM(app_ctx->comm, ceed, &dm)); // DM for projecting initial velocity to Hdiv space - PetscCall(CreateDM(comm, mat_type, vec_type, &dm_u0)); + PetscCall(CreateDM(app_ctx->comm, ceed, &dm_u0)); // DM for projecting initial pressure in L2 - PetscCall(CreateDM(comm, mat_type, vec_type, &dm_p0)); + PetscCall(CreateDM(app_ctx->comm, ceed, &dm_p0)); // DM for projecting solution U into H1 space for PetscViewer - PetscCall(CreateDM(comm, mat_type, vec_type, &dm_H1)); + PetscCall(CreateDM(app_ctx->comm, ceed, &dm_H1)); // TODO: add mesh option // perturb dm to have smooth random mesh // PetscCall( PerturbVerticesSmooth(dm) ); @@ -129,18 +115,10 @@ int main(int argc, char **argv) { // PetscCall(PerturbVerticesRandom(dm) ); // PetscCall(PerturbVerticesRandom(dm_H1) ); - // --------------------------------------------------------------------------- - // Process command line options - // --------------------------------------------------------------------------- - // -- Register problems to be available on the command line - PetscCall(RegisterProblems_Hdiv(app_ctx)); - - app_ctx->comm = comm; - PetscCall(ProcessCommandLineOptions(app_ctx)); - // --------------------------------------------------------------------------- // Choose the problem from the list of registered problems // --------------------------------------------------------------------------- + PetscCall(RegisterProblems_Hdiv(app_ctx)); { PetscErrorCode (*p)(Ceed, ProblemData, DM, void *); PetscCall(PetscFunctionListFind(app_ctx->problems, app_ctx->problem_name, &p)); @@ -151,7 +129,7 @@ int main(int argc, char **argv) { // --------------------------------------------------------------------------- // Setup FE for H(div) mixed-problem and H1 projection in post-processing.c // --------------------------------------------------------------------------- - PetscCall(SetupFEHdiv(comm, dm, dm_u0, dm_p0)); + PetscCall(SetupFEHdiv(app_ctx->comm, dm, dm_u0, dm_p0)); PetscCall(SetupFEH1(problem_data, app_ctx, dm_H1)); // --------------------------------------------------------------------------- @@ -167,36 +145,36 @@ int main(int argc, char **argv) { PetscCall(SetupLibceed(dm, dm_u0, dm_p0, dm_H1, ceed, app_ctx, problem_data, ceed_data)); // --------------------------------------------------------------------------- - // Setup pressure boundary conditions + // Setup pressure boundary conditions (not working) // --------------------------------------------------------------------------- // --Create empty local vector for libCEED - Vec P_loc; - PetscInt P_loc_size; - CeedScalar *p0; - CeedVector P_ceed; - PetscMemType pressure_mem_type; - PetscCall(DMCreateLocalVector(dm, &P_loc)); - PetscCall(VecGetSize(P_loc, &P_loc_size)); - PetscCall(VecZeroEntries(P_loc)); - PetscCall(VecGetArrayAndMemType(P_loc, &p0, &pressure_mem_type)); - CeedVectorCreate(ceed, P_loc_size, &P_ceed); - CeedVectorSetArray(P_ceed, MemTypeP2C(pressure_mem_type), CEED_USE_POINTER, p0); - // -- Apply operator to create local pressure vector on boundary - PetscCall(DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm, P_ceed)); - // CeedVectorView(P_ceed, "%12.8f", stdout); - // -- Map local to global - Vec P; - CeedVectorTakeArray(P_ceed, MemTypeP2C(pressure_mem_type), NULL); - PetscCall(VecRestoreArrayAndMemType(P_loc, &p0)); - PetscCall(DMCreateGlobalVector(dm, &P)); - PetscCall(VecZeroEntries(P)); - PetscCall(DMLocalToGlobal(dm, P_loc, ADD_VALUES, P)); + // Vec P_loc; + // PetscInt P_loc_size; + // CeedScalar *p0; + // CeedVector P_ceed; + // PetscMemType pressure_mem_type; + // PetscCall(DMCreateLocalVector(dm, &P_loc)); + // PetscCall(VecGetSize(P_loc, &P_loc_size)); + // PetscCall(VecZeroEntries(P_loc)); + // PetscCall(VecGetArrayAndMemType(P_loc, &p0, &pressure_mem_type)); + // CeedVectorCreate(ceed, P_loc_size, &P_ceed); + // CeedVectorSetArray(P_ceed, MemTypeP2C(pressure_mem_type), CEED_USE_POINTER, p0); + //// -- Apply operator to create local pressure vector on boundary + // PetscCall(DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm, P_ceed)); + //// CeedVectorView(P_ceed, "%12.8f", stdout); + //// -- Map local to global + // Vec P; + // CeedVectorTakeArray(P_ceed, MemTypeP2C(pressure_mem_type), NULL); + // PetscCall(VecRestoreArrayAndMemType(P_loc, &p0)); + // PetscCall(DMCreateGlobalVector(dm, &P)); + // PetscCall(VecZeroEntries(P)); + // PetscCall(DMLocalToGlobal(dm, P_loc, ADD_VALUES, P)); // --------------------------------------------------------------------------- // Setup context for projection problem; post-processing.c // --------------------------------------------------------------------------- - PetscCall(SetupProjectVelocityCtx_Hdiv(comm, dm, ceed, ceed_data, app_ctx->ctx_Hdiv)); - PetscCall(SetupProjectVelocityCtx_H1(comm, dm_H1, ceed, ceed_data, vec_type, app_ctx->ctx_H1)); + PetscCall(SetupProjectVelocityCtx_Hdiv(app_ctx->comm, dm, ceed, ceed_data, app_ctx->ctx_Hdiv)); + PetscCall(SetupProjectVelocityCtx_H1(app_ctx->comm, dm_H1, ceed, ceed_data, app_ctx->ctx_H1)); // --------------------------------------------------------------------------- // Setup TSSolve for Richard problem @@ -206,13 +184,13 @@ int main(int argc, char **argv) { // --------------------------------------------------------------------------- // Setup context for initial conditions // --------------------------------------------------------------------------- - PetscCall(SetupResidualOperatorCtx_U0(comm, dm_u0, ceed, ceed_data, app_ctx->ctx_initial_u0)); - PetscCall(SetupResidualOperatorCtx_P0(comm, dm_p0, ceed, ceed_data, app_ctx->ctx_initial_p0)); - PetscCall(SetupResidualOperatorCtx_Ut(comm, dm, ceed, ceed_data, app_ctx->ctx_residual_ut)); - PetscCall(CreateInitialConditions(ceed_data, app_ctx, vec_type, U)); + PetscCall(SetupResidualOperatorCtx_U0(app_ctx->comm, dm_u0, ceed, ceed_data, app_ctx->ctx_initial_u0)); + PetscCall(SetupResidualOperatorCtx_P0(app_ctx->comm, dm_p0, ceed, ceed_data, app_ctx->ctx_initial_p0)); + PetscCall(SetupResidualOperatorCtx_Ut(app_ctx->comm, dm, ceed, ceed_data, app_ctx->ctx_residual_ut)); + PetscCall(CreateInitialConditions(ceed_data, app_ctx, U)); // VecView(U, PETSC_VIEWER_STDOUT_WORLD); // Solve Richards problem - PetscCall(TSCreate(comm, &ts)); + PetscCall(TSCreate(app_ctx->comm, &ts)); PetscCall(VecZeroEntries(app_ctx->ctx_residual_ut->X_loc)); PetscCall(VecZeroEntries(app_ctx->ctx_residual_ut->X_t_loc)); PetscCall(TSSolveRichard(ceed_data, app_ctx, ts, &U)); @@ -225,10 +203,10 @@ int main(int argc, char **argv) { SNES snes; KSP ksp; if (!problem_data->has_ts) { - PetscCall(SetupJacobianOperatorCtx(dm, ceed, ceed_data, vec_type, app_ctx->ctx_jacobian)); + PetscCall(SetupJacobianOperatorCtx(dm, ceed, ceed_data, app_ctx->ctx_jacobian)); PetscCall(SetupResidualOperatorCtx(dm, ceed, ceed_data, app_ctx->ctx_residual)); // Create SNES - PetscCall(SNESCreate(comm, &snes)); + PetscCall(SNESCreate(app_ctx->comm, &snes)); PetscCall(SNESGetKSP(snes, &ksp)); PetscCall(PDESolver(ceed_data, app_ctx, snes, ksp, &U)); // VecView(U, PETSC_VIEWER_STDOUT_WORLD); @@ -244,14 +222,14 @@ int main(int argc, char **argv) { // --------------------------------------------------------------------------- // Print solver iterations and final norms // --------------------------------------------------------------------------- - PetscCall(PrintOutput(dm, ceed, app_ctx, problem_data->has_ts, mem_type_backend, ts, snes, ksp, U, l2_error_u, l2_error_p)); + PetscCall(PrintOutput(dm, ceed, app_ctx, problem_data->has_ts, ts, snes, ksp, U, l2_error_u, l2_error_p)); // --------------------------------------------------------------------------- // Save solution (paraview) // --------------------------------------------------------------------------- if (app_ctx->view_solution) { PetscViewer viewer_p; - PetscCall(PetscViewerVTKOpen(comm, "darcy_pressure.vtu", FILE_MODE_WRITE, &viewer_p)); + PetscCall(PetscViewerVTKOpen(app_ctx->comm, "darcy_pressure.vtu", FILE_MODE_WRITE, &viewer_p)); PetscCall(VecView(U, viewer_p)); PetscCall(PetscViewerDestroy(&viewer_p)); @@ -260,7 +238,7 @@ int main(int argc, char **argv) { PetscCall(ProjectVelocity(app_ctx, U, &U_H1)); PetscViewer viewer_u; - PetscCall(PetscViewerVTKOpen(comm, "darcy_velocity.vtu", FILE_MODE_WRITE, &viewer_u)); + PetscCall(PetscViewerVTKOpen(app_ctx->comm, "darcy_velocity.vtu", FILE_MODE_WRITE, &viewer_u)); PetscCall(VecView(U_H1, viewer_u)); PetscCall(PetscViewerDestroy(&viewer_u)); PetscCall(VecDestroy(&U_H1)); diff --git a/examples/Hdiv-mixed/src/cl-options.c b/examples/Hdiv-mixed/src/cl-options.c index 0d67fedb59..0b8bf4e387 100644 --- a/examples/Hdiv-mixed/src/cl-options.c +++ b/examples/Hdiv-mixed/src/cl-options.c @@ -25,7 +25,7 @@ PetscErrorCode ProcessCommandLineOptions(AppCtx app_ctx) { PetscBool ceed_flag = PETSC_FALSE; PetscFunctionBeginUser; - PetscOptionsBegin(app_ctx->comm, NULL, "H(div) examples in PETSc with libCEED", NULL); + PetscOptionsBegin(app_ctx->comm, NULL, "H(div) mixed-problem in PETSc with libCEED", NULL); PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, app_ctx->ceed_resource, app_ctx->ceed_resource, sizeof(app_ctx->ceed_resource), &ceed_flag)); diff --git a/examples/Hdiv-mixed/src/post-processing.c b/examples/Hdiv-mixed/src/post-processing.c index 2a99a61cf3..8bc99d0a04 100644 --- a/examples/Hdiv-mixed/src/post-processing.c +++ b/examples/Hdiv-mixed/src/post-processing.c @@ -5,12 +5,14 @@ // ----------------------------------------------------------------------------- // This function print the output // ----------------------------------------------------------------------------- -PetscErrorCode PrintOutput(DM dm, Ceed ceed, AppCtx app_ctx, PetscBool has_ts, CeedMemType mem_type_backend, TS ts, SNES snes, KSP ksp, Vec U, - CeedScalar l2_error_u, CeedScalar l2_error_p) { +PetscErrorCode PrintOutput(DM dm, Ceed ceed, AppCtx app_ctx, PetscBool has_ts, TS ts, SNES snes, KSP ksp, Vec U, CeedScalar l2_error_u, + CeedScalar l2_error_p) { PetscFunctionBeginUser; const char *used_resource; + CeedMemType mem_type_backend; CeedGetResource(ceed, &used_resource); + CeedGetPreferredMemType(ceed, &mem_type_backend); char hostname[PETSC_MAX_PATH_LEN]; PetscCall(PetscGetHostName(hostname, sizeof hostname)); PetscInt comm_size; @@ -128,9 +130,11 @@ PetscErrorCode SetupProjectVelocityCtx_Hdiv(MPI_Comm comm, DM dm, Ceed ceed, Cee PetscFunctionReturn(0); } -PetscErrorCode SetupProjectVelocityCtx_H1(MPI_Comm comm, DM dm_H1, Ceed ceed, CeedData ceed_data, VecType vec_type, OperatorApplyContext ctx_H1) { +PetscErrorCode SetupProjectVelocityCtx_H1(MPI_Comm comm, DM dm_H1, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_H1) { PetscFunctionBeginUser; + VecType vec_type; + PetscCall(DMGetVecType(dm_H1, &vec_type)); ctx_H1->comm = comm; ctx_H1->dm = dm_H1; PetscCall(DMCreateLocalVector(dm_H1, &ctx_H1->X_loc)); diff --git a/examples/Hdiv-mixed/src/setup-dm.c b/examples/Hdiv-mixed/src/setup-dm.c index 3e39ebd7e6..186dd90eca 100644 --- a/examples/Hdiv-mixed/src/setup-dm.c +++ b/examples/Hdiv-mixed/src/setup-dm.c @@ -3,11 +3,33 @@ #include "petscerror.h" // --------------------------------------------------------------------------- -// Setup DM +// Create DM // --------------------------------------------------------------------------- -PetscErrorCode CreateDM(MPI_Comm comm, MatType mat_type, VecType vec_type, DM *dm) { +PetscErrorCode CreateDM(MPI_Comm comm, Ceed ceed, DM *dm) { PetscFunctionBeginUser; + CeedMemType mem_type_backend; + CeedGetPreferredMemType(ceed, &mem_type_backend); + + VecType vec_type = NULL; + MatType mat_type = NULL; + switch (mem_type_backend) { + case CEED_MEM_HOST: + vec_type = VECSTANDARD; + break; + case CEED_MEM_DEVICE: { + const char *resolved; + CeedGetResource(ceed, &resolved); + if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA; + else if (strstr(resolved, "/gpu/hip/occa")) vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678 + else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP; + else vec_type = VECSTANDARD; + } + } + if (strstr(vec_type, VECCUDA)) mat_type = MATAIJCUSPARSE; + else if (strstr(vec_type, VECKOKKOS)) mat_type = MATAIJKOKKOS; + else mat_type = MATAIJ; + // Create DMPLEX PetscCall(DMCreate(comm, dm)); PetscCall(DMSetType(*dm, DMPLEX)); diff --git a/examples/Hdiv-mixed/src/setup-fe.c b/examples/Hdiv-mixed/src/setup-fe.c index e4d059bb03..da4b000055 100644 --- a/examples/Hdiv-mixed/src/setup-fe.c +++ b/examples/Hdiv-mixed/src/setup-fe.c @@ -2,6 +2,11 @@ #include "petscerror.h" +// ----------------------------------------------------------------------------- +// Convert PETSc MemType to libCEED MemType +// ----------------------------------------------------------------------------- +CeedMemType MemTypeP2C(PetscMemType mem_type) { return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; } + // --------------------------------------------------------------------------- // Setup FE // --------------------------------------------------------------------------- @@ -122,4 +127,155 @@ PetscErrorCode SetupFEH1(ProblemData problem_data, AppCtx app_ctx, DM dm_H1) { PetscCall(PetscSectionSetComponentName(section, 0, 2, "Velocity_Z")); } PetscFunctionReturn(0); -}; \ No newline at end of file +}; + +// ----------------------------------------------------------------------------- +// Utility function - convert from DMPolytopeType to CeedElemTopology +// ----------------------------------------------------------------------------- +CeedElemTopology ElemTopologyP2C(DMPolytopeType cell_type) { + switch (cell_type) { + case DM_POLYTOPE_TRIANGLE: + return CEED_TOPOLOGY_TRIANGLE; + case DM_POLYTOPE_QUADRILATERAL: + return CEED_TOPOLOGY_QUAD; + case DM_POLYTOPE_TETRAHEDRON: + return CEED_TOPOLOGY_TET; + case DM_POLYTOPE_HEXAHEDRON: + return CEED_TOPOLOGY_HEX; + default: + return 0; + } +}; + +// ----------------------------------------------------------------------------- +// Utility function - essential BC dofs are encoded in closure indices as -(i+1) +// ----------------------------------------------------------------------------- +PetscInt Involute(PetscInt i) { return i >= 0 ? i : -(i + 1); }; + +// ----------------------------------------------------------------------------- +// Get CEED restriction data from DMPlex +// ----------------------------------------------------------------------------- +PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { + PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; + + PetscFunctionBeginUser; + + PetscCall(DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); + + CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restr_offsets, elem_restr); + PetscCall(PetscFree(elem_restr_offsets)); + + PetscFunctionReturn(0); +}; + +// ----------------------------------------------------------------------------- +// Get Oriented CEED restriction data from DMPlex +// ----------------------------------------------------------------------------- +PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, DM dm_u0, DM dm_p0, CeedInt P, CeedElemRestriction *elem_restr_u, + CeedElemRestriction *elem_restr_p, CeedElemRestriction *elem_restr_u0, + CeedElemRestriction *elem_restr_p0) { + PetscSection section, section_u0, section_p0; + PetscInt p, num_elem, num_dof, num_dof_u0, num_dof_p0, *restr_indices_u, *restr_indices_p, *restr_indices_u0, *restr_indices_p0, elem_offset, + num_fields, num_fields_u0, num_fields_p0, dim, c_start, c_end; + Vec U_loc; + const PetscInt *ornt; // this is for orientation of dof + PetscFunctionBeginUser; + // Section for mixed problem + PetscCall(DMGetDimension(dm, &dim)); + PetscCall(DMGetLocalSection(dm, §ion)); + PetscCall(PetscSectionGetNumFields(section, &num_fields)); + PetscInt num_comp[num_fields], field_offsets[num_fields + 1]; + field_offsets[0] = 0; + for (PetscInt f = 0; f < num_fields; f++) { + PetscCall(PetscSectionGetFieldComponents(section, f, &num_comp[f])); + field_offsets[f + 1] = field_offsets[f] + num_comp[f]; + } + // Section for initial conditions u0 + PetscCall(DMGetLocalSection(dm_u0, §ion_u0)); + PetscCall(PetscSectionGetNumFields(section_u0, &num_fields_u0)); + PetscInt num_comp_u0[num_fields_u0], field_offsets_u0[num_fields_u0 + 1]; + field_offsets_u0[0] = 0; + for (PetscInt f = 0; f < num_fields_u0; f++) { + PetscCall(PetscSectionGetFieldComponents(section_u0, f, &num_comp_u0[f])); + field_offsets_u0[f + 1] = field_offsets_u0[f] + num_comp_u0[f]; + } + // Section for initial conditions p0 + PetscCall(DMGetLocalSection(dm_p0, §ion_p0)); + PetscCall(PetscSectionGetNumFields(section_p0, &num_fields_p0)); + PetscInt num_comp_p0[num_fields_p0], field_offsets_p0[num_fields_p0 + 1]; + field_offsets_p0[0] = 0; + for (PetscInt f = 0; f < num_fields_p0; f++) { + PetscCall(PetscSectionGetFieldComponents(section_p0, f, &num_comp_p0[f])); + field_offsets_p0[f + 1] = field_offsets_p0[f] + num_comp_p0[f]; + } + + PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); + num_elem = c_end - c_start; + PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &restr_indices_u)); + PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &restr_indices_u0)); + PetscCall(PetscMalloc1(num_elem, &restr_indices_p)); + PetscCall(PetscMalloc1(num_elem, &restr_indices_p0)); + bool *orient_indices_u, *orient_indices_u0; // to flip the dof + PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &orient_indices_u)); + PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &orient_indices_u0)); + for (p = 0, elem_offset = 0; p < num_elem; p++) { + PetscInt num_indices, *indices, faces_per_elem, dofs_per_face, num_indices_u0, *indices_u0, num_indices_p0, *indices_p0; + PetscCall(DMPlexGetClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL)); + PetscCall(DMPlexGetClosureIndices(dm_u0, section_u0, section_u0, p, PETSC_TRUE, &num_indices_u0, &indices_u0, NULL, NULL)); + PetscCall(DMPlexGetClosureIndices(dm_p0, section_p0, section_p0, p, PETSC_TRUE, &num_indices_p0, &indices_p0, NULL, NULL)); + restr_indices_p[p] = indices[num_indices - 1]; + restr_indices_p0[p] = indices_p0[0]; + PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); + // Get number of faces per element + PetscCall(DMPlexGetConeSize(dm, p, &faces_per_elem)); + dofs_per_face = faces_per_elem - 2; + for (PetscInt f = 0; f < faces_per_elem; f++) { + for (PetscInt i = 0; i < dofs_per_face; i++) { + PetscInt ii = dofs_per_face * f + i; + // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. + PetscInt loc = Involute(indices[ii * num_comp[0]]); + restr_indices_u[elem_offset] = loc; + // Set orientation + orient_indices_u[elem_offset] = ornt[f] < 0; + PetscInt loc_u0 = Involute(indices_u0[ii * num_comp_u0[0]]); + restr_indices_u0[elem_offset] = loc_u0; + // Set orientation + orient_indices_u0[elem_offset] = ornt[f] < 0; + elem_offset++; + } + } + PetscCall(DMPlexRestoreClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL)); + PetscCall(DMPlexRestoreClosureIndices(dm_u0, section_u0, section_u0, p, PETSC_TRUE, &num_indices_u0, &indices_u0, NULL, NULL)); + PetscCall(DMPlexRestoreClosureIndices(dm_p0, section_p0, section_p0, p, PETSC_TRUE, &num_indices_p0, &indices_p0, NULL, NULL)); + } + // if (elem_offset != num_elem*dim*PetscPowInt(P, dim)) + // SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, + // "ElemRestriction of size (%" PetscInt_FMT ", %" PetscInt_FMT" ) + // initialized %" PetscInt_FMT " nodes", num_elem, + // dim*PetscPowInt(P, dim),elem_offset); + + PetscCall(DMGetLocalVector(dm, &U_loc)); + PetscCall(VecGetLocalSize(U_loc, &num_dof)); + PetscCall(DMRestoreLocalVector(dm, &U_loc)); + // dof per element in Hdiv is dim*P^dim, for linear element P=2 + CeedElemRestrictionCreateOriented(ceed, num_elem, dim * PetscPowInt(P, dim), 1, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_u, + orient_indices_u, elem_restr_u); + CeedElemRestrictionCreate(ceed, num_elem, 1, 1, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_p, elem_restr_p); + PetscCall(DMGetLocalVector(dm_u0, &U_loc)); + PetscCall(VecGetLocalSize(U_loc, &num_dof_u0)); + PetscCall(DMRestoreLocalVector(dm_u0, &U_loc)); + // dof per element in Hdiv is dim*P^dim, for linear element P=2 + CeedElemRestrictionCreateOriented(ceed, num_elem, dim * PetscPowInt(P, dim), 1, 1, num_dof_u0, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_u0, + orient_indices_u0, elem_restr_u0); + PetscCall(DMGetLocalVector(dm_p0, &U_loc)); + PetscCall(VecGetLocalSize(U_loc, &num_dof_p0)); + PetscCall(DMRestoreLocalVector(dm_p0, &U_loc)); + CeedElemRestrictionCreate(ceed, num_elem, 1, 1, 1, num_dof_p0, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_p0, elem_restr_p0); + PetscCall(PetscFree(restr_indices_p)); + PetscCall(PetscFree(restr_indices_u)); + PetscCall(PetscFree(orient_indices_u)); + PetscCall(PetscFree(restr_indices_u0)); + PetscCall(PetscFree(orient_indices_u0)); + PetscCall(PetscFree(restr_indices_p0)); + PetscFunctionReturn(0); +}; diff --git a/examples/Hdiv-mixed/src/setup-libceed.c b/examples/Hdiv-mixed/src/setup-libceed.c index dc091f35e0..e147300890 100644 --- a/examples/Hdiv-mixed/src/setup-libceed.c +++ b/examples/Hdiv-mixed/src/setup-libceed.c @@ -9,10 +9,6 @@ #include "../include/setup-boundary.h" #include "ceed/ceed.h" -// ----------------------------------------------------------------------------- -// Convert PETSc MemType to libCEED MemType -// ----------------------------------------------------------------------------- -CeedMemType MemTypeP2C(PetscMemType mem_type) { return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; } // ----------------------------------------------------------------------------- // Destroy libCEED objects // ----------------------------------------------------------------------------- @@ -82,139 +78,6 @@ PetscErrorCode CeedDataDestroy(CeedData ceed_data, ProblemData problem_data) { PetscFunctionReturn(0); }; -// ----------------------------------------------------------------------------- -// Utility function - essential BC dofs are encoded in closure indices as -(i+1) -// ----------------------------------------------------------------------------- -PetscInt Involute(PetscInt i) { return i >= 0 ? i : -(i + 1); }; - -// ----------------------------------------------------------------------------- -// Get CEED restriction data from DMPlex -// ----------------------------------------------------------------------------- -PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt value, CeedElemRestriction *elem_restr) { - PetscInt num_elem, elem_size, num_dof, num_comp, *elem_restr_offsets; - - PetscFunctionBeginUser; - - PetscCall(DMPlexGetLocalOffsets(dm, domain_label, value, height, 0, &num_elem, &elem_size, &num_comp, &num_dof, &elem_restr_offsets)); - - CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, elem_restr_offsets, elem_restr); - PetscCall(PetscFree(elem_restr_offsets)); - - PetscFunctionReturn(0); -}; - -// ----------------------------------------------------------------------------- -// Get Oriented CEED restriction data from DMPlex -// ----------------------------------------------------------------------------- -PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, DM dm_u0, DM dm_p0, CeedInt P, CeedElemRestriction *elem_restr_u, - CeedElemRestriction *elem_restr_p, CeedElemRestriction *elem_restr_u0, - CeedElemRestriction *elem_restr_p0) { - PetscSection section, section_u0, section_p0; - PetscInt p, num_elem, num_dof, num_dof_u0, num_dof_p0, *restr_indices_u, *restr_indices_p, *restr_indices_u0, *restr_indices_p0, elem_offset, - num_fields, num_fields_u0, num_fields_p0, dim, c_start, c_end; - Vec U_loc; - const PetscInt *ornt; // this is for orientation of dof - PetscFunctionBeginUser; - // Section for mixed problem - PetscCall(DMGetDimension(dm, &dim)); - PetscCall(DMGetLocalSection(dm, §ion)); - PetscCall(PetscSectionGetNumFields(section, &num_fields)); - PetscInt num_comp[num_fields], field_offsets[num_fields + 1]; - field_offsets[0] = 0; - for (PetscInt f = 0; f < num_fields; f++) { - PetscCall(PetscSectionGetFieldComponents(section, f, &num_comp[f])); - field_offsets[f + 1] = field_offsets[f] + num_comp[f]; - } - // Section for initial conditions u0 - PetscCall(DMGetLocalSection(dm_u0, §ion_u0)); - PetscCall(PetscSectionGetNumFields(section_u0, &num_fields_u0)); - PetscInt num_comp_u0[num_fields_u0], field_offsets_u0[num_fields_u0 + 1]; - field_offsets_u0[0] = 0; - for (PetscInt f = 0; f < num_fields_u0; f++) { - PetscCall(PetscSectionGetFieldComponents(section_u0, f, &num_comp_u0[f])); - field_offsets_u0[f + 1] = field_offsets_u0[f] + num_comp_u0[f]; - } - // Section for initial conditions p0 - PetscCall(DMGetLocalSection(dm_p0, §ion_p0)); - PetscCall(PetscSectionGetNumFields(section_p0, &num_fields_p0)); - PetscInt num_comp_p0[num_fields_p0], field_offsets_p0[num_fields_p0 + 1]; - field_offsets_p0[0] = 0; - for (PetscInt f = 0; f < num_fields_p0; f++) { - PetscCall(PetscSectionGetFieldComponents(section_p0, f, &num_comp_p0[f])); - field_offsets_p0[f + 1] = field_offsets_p0[f] + num_comp_p0[f]; - } - - PetscCall(DMPlexGetHeightStratum(dm, 0, &c_start, &c_end)); - num_elem = c_end - c_start; - PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &restr_indices_u)); - PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &restr_indices_u0)); - PetscCall(PetscMalloc1(num_elem, &restr_indices_p)); - PetscCall(PetscMalloc1(num_elem, &restr_indices_p0)); - bool *orient_indices_u, *orient_indices_u0; // to flip the dof - PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &orient_indices_u)); - PetscCall(PetscMalloc1(num_elem * dim * PetscPowInt(P, dim), &orient_indices_u0)); - for (p = 0, elem_offset = 0; p < num_elem; p++) { - PetscInt num_indices, *indices, faces_per_elem, dofs_per_face, num_indices_u0, *indices_u0, num_indices_p0, *indices_p0; - PetscCall(DMPlexGetClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL)); - PetscCall(DMPlexGetClosureIndices(dm_u0, section_u0, section_u0, p, PETSC_TRUE, &num_indices_u0, &indices_u0, NULL, NULL)); - PetscCall(DMPlexGetClosureIndices(dm_p0, section_p0, section_p0, p, PETSC_TRUE, &num_indices_p0, &indices_p0, NULL, NULL)); - restr_indices_p[p] = indices[num_indices - 1]; - restr_indices_p0[p] = indices_p0[0]; - PetscCall(DMPlexGetConeOrientation(dm, p, &ornt)); - // Get number of faces per element - PetscCall(DMPlexGetConeSize(dm, p, &faces_per_elem)); - dofs_per_face = faces_per_elem - 2; - for (PetscInt f = 0; f < faces_per_elem; f++) { - for (PetscInt i = 0; i < dofs_per_face; i++) { - PetscInt ii = dofs_per_face * f + i; - // Essential boundary conditions are encoded as -(loc+1), but we don't care so we decode. - PetscInt loc = Involute(indices[ii * num_comp[0]]); - restr_indices_u[elem_offset] = loc; - // Set orientation - orient_indices_u[elem_offset] = ornt[f] < 0; - PetscInt loc_u0 = Involute(indices_u0[ii * num_comp_u0[0]]); - restr_indices_u0[elem_offset] = loc_u0; - // Set orientation - orient_indices_u0[elem_offset] = ornt[f] < 0; - elem_offset++; - } - } - PetscCall(DMPlexRestoreClosureIndices(dm, section, section, p, PETSC_TRUE, &num_indices, &indices, NULL, NULL)); - PetscCall(DMPlexRestoreClosureIndices(dm_u0, section_u0, section_u0, p, PETSC_TRUE, &num_indices_u0, &indices_u0, NULL, NULL)); - PetscCall(DMPlexRestoreClosureIndices(dm_p0, section_p0, section_p0, p, PETSC_TRUE, &num_indices_p0, &indices_p0, NULL, NULL)); - } - // if (elem_offset != num_elem*dim*PetscPowInt(P, dim)) - // SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, - // "ElemRestriction of size (%" PetscInt_FMT ", %" PetscInt_FMT" ) - // initialized %" PetscInt_FMT " nodes", num_elem, - // dim*PetscPowInt(P, dim),elem_offset); - - PetscCall(DMGetLocalVector(dm, &U_loc)); - PetscCall(VecGetLocalSize(U_loc, &num_dof)); - PetscCall(DMRestoreLocalVector(dm, &U_loc)); - // dof per element in Hdiv is dim*P^dim, for linear element P=2 - CeedElemRestrictionCreateOriented(ceed, num_elem, dim * PetscPowInt(P, dim), 1, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_u, - orient_indices_u, elem_restr_u); - CeedElemRestrictionCreate(ceed, num_elem, 1, 1, 1, num_dof, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_p, elem_restr_p); - PetscCall(DMGetLocalVector(dm_u0, &U_loc)); - PetscCall(VecGetLocalSize(U_loc, &num_dof_u0)); - PetscCall(DMRestoreLocalVector(dm_u0, &U_loc)); - // dof per element in Hdiv is dim*P^dim, for linear element P=2 - CeedElemRestrictionCreateOriented(ceed, num_elem, dim * PetscPowInt(P, dim), 1, 1, num_dof_u0, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_u0, - orient_indices_u0, elem_restr_u0); - PetscCall(DMGetLocalVector(dm_p0, &U_loc)); - PetscCall(VecGetLocalSize(U_loc, &num_dof_p0)); - PetscCall(DMRestoreLocalVector(dm_p0, &U_loc)); - CeedElemRestrictionCreate(ceed, num_elem, 1, 1, 1, num_dof_p0, CEED_MEM_HOST, CEED_COPY_VALUES, restr_indices_p0, elem_restr_p0); - PetscCall(PetscFree(restr_indices_p)); - PetscCall(PetscFree(restr_indices_u)); - PetscCall(PetscFree(orient_indices_u)); - PetscCall(PetscFree(restr_indices_u0)); - PetscCall(PetscFree(orient_indices_u0)); - PetscCall(PetscFree(restr_indices_p0)); - PetscFunctionReturn(0); -}; - // ----------------------------------------------------------------------------- // Set up libCEED on the fine grid for a given degree // ----------------------------------------------------------------------------- diff --git a/examples/Hdiv-mixed/src/setup-solvers.c b/examples/Hdiv-mixed/src/setup-solvers.c index 3fef27c924..72038d53a8 100644 --- a/examples/Hdiv-mixed/src/setup-solvers.c +++ b/examples/Hdiv-mixed/src/setup-solvers.c @@ -7,9 +7,11 @@ // ----------------------------------------------------------------------------- // Setup operator context data // ----------------------------------------------------------------------------- -PetscErrorCode SetupJacobianOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data, VecType vec_type, OperatorApplyContext ctx_jacobian) { +PetscErrorCode SetupJacobianOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data, OperatorApplyContext ctx_jacobian) { PetscFunctionBeginUser; + VecType vec_type; + PetscCall(DMGetVecType(dm, &vec_type)); ctx_jacobian->dm = dm; PetscCall(DMCreateLocalVector(dm, &ctx_jacobian->X_loc)); PetscCall(VecDuplicate(ctx_jacobian->X_loc, &ctx_jacobian->Y_loc)); @@ -88,29 +90,29 @@ PetscErrorCode SNESFormResidual(SNES snes, Vec X, Vec Y, void *ctx_residual) { // Jacobian setup // ----------------------------------------------------------------------------- PetscErrorCode SNESFormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx_jacobian) { - OperatorApplyContext ctx = (OperatorApplyContext)ctx_jacobian; + // OperatorApplyContext ctx = (OperatorApplyContext)ctx_jacobian; PetscFunctionBeginUser; - Mat A; - PetscCall(DMCreateMatrix(ctx->dm, &A)); - // Assemble matrix analytically - PetscCount num_entries; - CeedInt *rows, *cols; - CeedVector coo_values; - CeedOperatorLinearAssembleSymbolic(ctx->op_apply, &num_entries, &rows, &cols); - PetscCall(MatSetPreallocationCOO(A, num_entries, rows, cols)); - free(rows); - free(cols); - CeedVectorCreate(ctx->ceed, num_entries, &coo_values); - CeedOperatorLinearAssemble(ctx->op_apply, coo_values); - const CeedScalar *values; - CeedVectorGetArrayRead(coo_values, CEED_MEM_HOST, &values); - PetscCall(MatSetValuesCOO(A, values, ADD_VALUES)); - CeedVectorRestoreArrayRead(coo_values, &values); - MatView(A, PETSC_VIEWER_STDOUT_WORLD); - // CeedVectorView(coo_values, "%12.8f", stdout); - CeedVectorDestroy(&coo_values); - PetscCall(MatDestroy(&A)); + // Mat A; + // PetscCall(DMCreateMatrix(ctx->dm, &A)); + //// Assemble matrix analytically + // PetscCount num_entries; + // CeedInt *rows, *cols; + // CeedVector coo_values; + // CeedOperatorLinearAssembleSymbolic(ctx->op_apply, &num_entries, &rows, &cols); + // PetscCall(MatSetPreallocationCOO(A, num_entries, rows, cols)); + // free(rows); + // free(cols); + // CeedVectorCreate(ctx->ceed, num_entries, &coo_values); + // CeedOperatorLinearAssemble(ctx->op_apply, coo_values); + // const CeedScalar *values; + // CeedVectorGetArrayRead(coo_values, CEED_MEM_HOST, &values); + // PetscCall(MatSetValuesCOO(A, values, ADD_VALUES)); + // CeedVectorRestoreArrayRead(coo_values, &values); + // MatView(A, PETSC_VIEWER_STDOUT_WORLD); + //// CeedVectorView(coo_values, "%12.8f", stdout); + // CeedVectorDestroy(&coo_values); + // PetscCall(MatDestroy(&A)); // J_pre might be AIJ (e.g., when using coloring), so we need to assemble it PetscCall(MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY)); diff --git a/examples/Hdiv-mixed/src/setup-ts.c b/examples/Hdiv-mixed/src/setup-ts.c index e93bca2495..8600967c00 100644 --- a/examples/Hdiv-mixed/src/setup-ts.c +++ b/examples/Hdiv-mixed/src/setup-ts.c @@ -69,7 +69,7 @@ PetscErrorCode SetupResidualOperatorCtx_Ut(MPI_Comm comm, DM dm, Ceed ceed, Ceed // ----------------------------------------------------------------------------- // Create global initial conditions vector // ----------------------------------------------------------------------------- -PetscErrorCode CreateInitialConditions(CeedData ceed_data, AppCtx app_ctx, VecType vec_type, Vec U) { +PetscErrorCode CreateInitialConditions(CeedData ceed_data, AppCtx app_ctx, Vec U) { PetscFunctionBeginUser; // ---------------------------------------------- // Create local rhs for u field @@ -77,6 +77,8 @@ PetscErrorCode CreateInitialConditions(CeedData ceed_data, AppCtx app_ctx, VecTy Vec rhs_u_loc; PetscScalar *ru; PetscMemType ru_mem_type; + VecType vec_type; + PetscCall(DMGetVecType(app_ctx->ctx_initial_u0->dm, &vec_type)); PetscCall(DMCreateLocalVector(app_ctx->ctx_initial_u0->dm, &rhs_u_loc)); PetscCall(VecZeroEntries(rhs_u_loc)); PetscCall(VecGetArrayAndMemType(rhs_u_loc, &ru, &ru_mem_type));