diff --git a/examples/petsc/multigrid.c b/examples/petsc/multigrid.c index db54e0449f..1318872d40 100644 --- a/examples/petsc/multigrid.c +++ b/examples/petsc/multigrid.c @@ -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; @@ -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, <og_row, <og_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); @@ -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); diff --git a/examples/solids/elasticity.c b/examples/solids/elasticity.c index 8bca063215..84a0aaf93b 100644 --- a/examples/solids/elasticity.c +++ b/examples/solids/elasticity.c @@ -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; @@ -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, @@ -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, <og_row, <og_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; } } @@ -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); @@ -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); diff --git a/examples/solids/include/structs.h b/examples/solids/include/structs.h index e4aa5bb408..e0cb273183 100644 --- a/examples/solids/include/structs.h +++ b/examples/solids/include/structs.h @@ -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 diff --git a/examples/solids/src/misc.c b/examples/solids/src/misc.c index 9989a43168..a92e3f7849 100644 --- a/examples/solids/src/misc.c +++ b/examples/solids/src/misc.c @@ -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);