Skip to content

PETSc examples: assemble matrix analytically instead of coloring #807

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jeremylt merged 5 commits into from
Feb 1, 2022
Merged

PETSc examples: assemble matrix analytically instead of coloring #807

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
53b04fa
examples/petsc/multigrid: assemble matrix analytically instead of col…
jedbrown Sep 9, 2021
17f0843
solids - add libCEED coarse level assembly
jeremylt Sep 14, 2021
caae649
build/solids: fix order/uninitialized variable
jedbrown Sep 14, 2021
05b9c82
examples/petsc/multigrid: add event for assembly
jedbrown Nov 4, 2021
cffe6a5
petsc - drop mat assembly by coloring
jeremylt Feb 1, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 30 additions & 18 deletions examples/petsc/multigrid.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,8 @@ int main(int argc, char **argv) {
PetscScalar eps = 1.0;
PetscBool test_mode, benchmark_mode, read_mesh, write_solution;
PetscLogStage solve_stage;
PetscLogEvent assemble_event;
DM *dm, dm_orig;
SNES snes_dummy;
KSP ksp;
PC pc;
Mat *mat_O, *mat_pr, mat_coarse;
Expand Down Expand Up @@ -430,24 +430,37 @@ int main(int argc, char **argv) {
}
}

// Setup dummy SNES for AMG coarse solve
ierr = SNESCreate(comm, &snes_dummy); CHKERRQ(ierr);
ierr = SNESSetDM(snes_dummy, dm[0]); CHKERRQ(ierr);
ierr = SNESSetSolution(snes_dummy, X[0]); CHKERRQ(ierr);

// -- Jacobian matrix
ierr = DMSetMatType(dm[0], MATAIJ); CHKERRQ(ierr);
// Assemble coarse grid Jacobian for AMG (or other sparse matrix) solve
ierr = DMCreateMatrix(dm[0], &mat_coarse); CHKERRQ(ierr);
ierr = SNESSetJacobian(snes_dummy, mat_coarse, mat_coarse, NULL,
NULL); CHKERRQ(ierr);

// -- Residual evaluation function
ierr = SNESSetFunction(snes_dummy, X[0], FormResidual_Ceed,
user_O[0]); CHKERRQ(ierr);

// -- Form Jacobian
ierr = SNESComputeJacobianDefaultColor(snes_dummy, X[0], mat_O[0],
mat_coarse, NULL); CHKERRQ(ierr);
ierr = PetscLogEventRegister("AssembleMatrix", MAT_CLASSID, &assemble_event);
CHKERRQ(ierr);
{
// Assemble matrix analytically
CeedInt num_entries, *rows, *cols;
CeedVector coo_values;
CeedOperatorLinearAssembleSymbolic(user_O[0]->op, &num_entries, &rows, &cols);
ISLocalToGlobalMapping ltog_row, ltog_col;
ierr = MatGetLocalToGlobalMapping(mat_coarse, &ltog_row, &ltog_col);
CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingApply(ltog_row, num_entries, rows, rows);
CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingApply(ltog_col, num_entries, cols, cols);
CHKERRQ(ierr);
ierr = MatSetPreallocationCOO(mat_coarse, num_entries, rows, cols);
CHKERRQ(ierr);
free(rows);
free(cols);
CeedVectorCreate(ceed, num_entries, &coo_values);
ierr = PetscLogEventBegin(assemble_event, mat_coarse, 0, 0, 0); CHKERRQ(ierr);
CeedOperatorLinearAssemble(user_O[0]->op, coo_values);
const CeedScalar *values;
CeedVectorGetArrayRead(coo_values, CEED_MEM_HOST, &values);
ierr = MatSetValuesCOO(mat_coarse, values, ADD_VALUES); CHKERRQ(ierr);
CeedVectorRestoreArrayRead(coo_values, &values);
ierr = PetscLogEventEnd(assemble_event, mat_coarse, 0, 0, 0); CHKERRQ(ierr);
CeedVectorDestroy(&coo_values);
}

// Set up KSP
ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr);
Expand Down Expand Up @@ -650,7 +663,6 @@ int main(int argc, char **argv) {
ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr);
ierr = MatDestroy(&mat_coarse); CHKERRQ(ierr);
ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
ierr = SNESDestroy(&snes_dummy); CHKERRQ(ierr);
ierr = DMDestroy(&dm_orig); CHKERRQ(ierr);
CeedVectorDestroy(&target);
CeedQFunctionDestroy(&qf_error);
Expand Down
46 changes: 24 additions & 22 deletions examples/solids/elasticity.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ int main(int argc, char **argv) {
Vec U, *U_g, *U_loc; // U: solution, R: residual, F: forcing
Vec R, R_loc, F, F_loc; // g: global, loc: local
Vec neumann_bcs = NULL, bcs_loc = NULL;
SNES snes, snes_coarse = NULL;
SNES snes;
Mat *jacob_mat, jacob_mat_coarse, *prolong_restr_mat;
// PETSc data
UserMult res_ctx, jacob_coarse_ctx = NULL, *jacob_ctx;
Expand Down Expand Up @@ -314,12 +314,12 @@ int main(int argc, char **argv) {
// Note: use high order ceed, if specified and degree > 4
PetscErrorCode (*SetupLibceedLevel)(DM, Ceed, AppCtx, PetscInt,
PetscInt, PetscInt, PetscInt, CeedVector, CeedData *);
if (!SetupLibceedLevel)
SETERRQ1(PETSC_COMM_SELF, 1, "Coarse grid setup for '%s' not found",
app_ctx->name);
ierr = PetscFunctionListFind(problem_functions->setupLibceedLevel,
app_ctx->name, &SetupLibceedLevel);
CHKERRQ(ierr);
if (!SetupLibceedLevel)
SETERRQ1(PETSC_COMM_SELF, 1, "Coarse grid setup for '%s' not found",
app_ctx->name);
ierr = (*SetupLibceedLevel)(level_dms[level], ceed, app_ctx,
level, num_comp_u, U_g_size[level],
U_loc_size[level], ceed_data[level+1]->x_ceed,
Expand Down Expand Up @@ -510,32 +510,34 @@ int main(int argc, char **argv) {
}

// ---------------------------------------------------------------------------
// Setup dummy SNES for AMG coarse solve
// Setup for AMG coarse solve
// ---------------------------------------------------------------------------
if (app_ctx->multigrid_choice != MULTIGRID_NONE) {
// -- Jacobian Matrix
ierr = DMSetMatType(level_dms[0], MATAIJ); CHKERRQ(ierr);
ierr = DMCreateMatrix(level_dms[0], &jacob_mat_coarse); CHKERRQ(ierr);

if (app_ctx->degree > 1) {
ierr = SNESCreate(comm, &snes_coarse); CHKERRQ(ierr);
ierr = SNESSetDM(snes_coarse, level_dms[0]); CHKERRQ(ierr);
ierr = SNESSetSolution(snes_coarse, U_g[0]); CHKERRQ(ierr);

// -- Jacobian function
ierr = SNESSetJacobian(snes_coarse, jacob_mat_coarse, jacob_mat_coarse, NULL,
NULL); CHKERRQ(ierr);

// -- Residual evaluation function
ierr = PetscMalloc1(1, &jacob_coarse_ctx); CHKERRQ(ierr);
ierr = PetscMemcpy(jacob_coarse_ctx, jacob_ctx[0], sizeof(*jacob_ctx[0]));
// -- Assemble sparsity pattern
CeedInt num_entries, *rows, *cols;
CeedVector coo_values;
CeedOperatorLinearAssembleSymbolic(ceed_data[0]->op_jacobian, &num_entries,
&rows, &cols);
ISLocalToGlobalMapping ltog_row, ltog_col;
ierr = MatGetLocalToGlobalMapping(jacob_mat_coarse, &ltog_row, &ltog_col);
CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingApply(ltog_row, num_entries, rows, rows);
CHKERRQ(ierr);
ierr = ISLocalToGlobalMappingApply(ltog_col, num_entries, cols, cols);
CHKERRQ(ierr);
ierr = MatSetPreallocationCOO(jacob_mat_coarse, num_entries, rows, cols);
CHKERRQ(ierr);
ierr = SNESSetFunction(snes_coarse, U_g[0], ApplyJacobianCoarse_Ceed,
jacob_coarse_ctx); CHKERRQ(ierr);
free(rows);
free(cols);
CeedVectorCreate(ceed, num_entries, &coo_values);

// -- Update form_jacob_ctx
form_jacob_ctx->u_coarse = U_g[0];
form_jacob_ctx->snes_coarse = snes_coarse;
form_jacob_ctx->coo_values = coo_values;
form_jacob_ctx->op_coarse = ceed_data[0]->op_jacobian;
form_jacob_ctx->jacob_mat_coarse = jacob_mat_coarse;
}
}
Expand Down Expand Up @@ -957,6 +959,7 @@ int main(int argc, char **argv) {
ierr = PetscFree(ceed_data); CHKERRQ(ierr);

// libCEED objects
CeedVectorDestroy(&form_jacob_ctx->coo_values);
CeedQFunctionContextDestroy(&ctx_phys);
CeedQFunctionContextDestroy(&ctx_phys_smoother);
CeedDestroy(&ceed);
Expand All @@ -971,7 +974,6 @@ int main(int argc, char **argv) {
ierr = VecDestroy(&bcs_loc); CHKERRQ(ierr);
ierr = MatDestroy(&jacob_mat_coarse); CHKERRQ(ierr);
ierr = SNESDestroy(&snes); CHKERRQ(ierr);
ierr = SNESDestroy(&snes_coarse); CHKERRQ(ierr);
ierr = DMDestroy(&dm_orig); CHKERRQ(ierr);
ierr = DMDestroy(&dm_energy); CHKERRQ(ierr);
ierr = DMDestroy(&dm_diagnostic); CHKERRQ(ierr);
Expand Down
4 changes: 2 additions & 2 deletions examples/solids/include/structs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,9 +109,9 @@ typedef struct FormJacobCtx_ *FormJacobCtx;
struct FormJacobCtx_ {
UserMult *jacob_ctx;
PetscInt num_levels;
SNES snes_coarse;
Mat *jacob_mat, jacob_mat_coarse;
Vec u_coarse;
CeedVector coo_values;
CeedOperator op_coarse;
};

// Data for PETSc Prolongation/Restriction Matshell
Expand Down
11 changes: 6 additions & 5 deletions examples/solids/src/misc.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,12 +87,13 @@ PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx) {
}

// Form coarse assembled matrix
ierr = VecZeroEntries(form_jacob_ctx->u_coarse); CHKERRQ(ierr);
ierr = SNESComputeJacobianDefaultColor(form_jacob_ctx->snes_coarse,
form_jacob_ctx->u_coarse,
form_jacob_ctx->jacob_mat[0],
form_jacob_ctx->jacob_mat_coarse, NULL);
CeedOperatorLinearAssemble(form_jacob_ctx->op_coarse,
form_jacob_ctx->coo_values);
const CeedScalar *values;
CeedVectorGetArrayRead(form_jacob_ctx->coo_values, CEED_MEM_HOST, &values);
ierr = MatSetValuesCOO(form_jacob_ctx->jacob_mat_coarse, values, ADD_VALUES);
CHKERRQ(ierr);
CeedVectorRestoreArrayRead(form_jacob_ctx->coo_values, &values);

// J_pre might be AIJ (e.g., when using coloring), so we need to assemble it
ierr = MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
Expand Down