From cf2dd0e2e102a37ac49a242289a4c93242ec244b Mon Sep 17 00:00:00 2001 From: jeremylt Date: Tue, 13 Apr 2021 14:59:38 -0600 Subject: [PATCH] wip --- examples/petsc/bddc.c | 182 +++++++++++++++++++++--------- examples/petsc/include/matops.h | 3 +- examples/petsc/include/structs.h | 18 +-- examples/petsc/multigrid.c | 17 +++ examples/petsc/src/libceedsetup.c | 35 +----- examples/petsc/src/matops.c | 43 +++++++ 6 files changed, 200 insertions(+), 98 deletions(-) diff --git a/examples/petsc/bddc.c b/examples/petsc/bddc.c index fc06948005..c1505b1be7 100644 --- a/examples/petsc/bddc.c +++ b/examples/petsc/bddc.c @@ -46,24 +46,25 @@ int main(int argc, char **argv) { char filename[PETSC_MAX_PATH_LEN], ceed_resource[PETSC_MAX_PATH_LEN] = "/cpu/self"; double my_rt_start, my_rt, rt_min, rt_max; - PetscInt degree = 3, q_extra, l_size, xl_size, g_size, l_vertex_size, - xl_vertex_size, g_vertex_size, dim = 3, m_elem[3] = {3, 3, 3}, - num_comp_u = 1; + PetscInt degree = 3, q_extra, l_size, xl_size, g_size, l_Pi_size, + xl_Pi_size, g_Pi_size, dim = 3, m_elem[3] = {3, 3, 3}, num_comp_u = 1; PetscScalar *r; PetscScalar eps = 1.0; PetscBool test_mode, benchmark_mode, read_mesh, write_solution; PetscLogStage solve_stage; - DM dm, dm_vertex; + DM dm, dm_Pi; KSP ksp; PC pc; - Mat mat_O, mat_vertex, mat_pr; - Vec X, X_loc, X_vertex, X_vertex_loc, mult, rhs, rhs_loc; + Mat mat_O, mat_Pi, mat_pr; + Vec X, X_loc, X_Pi, X_Pi_loc, rhs, rhs_loc, mult, mask_r, mask_Gamma, + mask_I; PetscMemType mem_type; - UserO user_O, user_vertex; - UserProlongRestr user_pr; + UserO user_O, user_Pi; + UserBDDC user_bddc; Ceed ceed; - CeedData ceed_data, ceed_data_vertex; - CeedVector rhs_ceed, target; + CeedData ceed_data + CeedDataBDDC ceed_data_bddc; + CeedVector rhs_ceed, mult_ceed, target; CeedQFunction qf_error, qf_restrict, qf_prolong; CeedOperator op_error; BPType bp_choice; @@ -153,7 +154,7 @@ int main(int argc, char **argv) { dm = dm_dist; } } - ierr = DMClone(dm, &dm_vertex); CHKERRQ(ierr); + ierr = DMClone(dm, &dm_Pi); CHKERRQ(ierr); // Apply Kershaw mesh transformation ierr = Kershaw(dm, eps); CHKERRQ(ierr); @@ -179,10 +180,10 @@ int main(int argc, char **argv) { bp_options[bp_choice].enforce_bc, bp_options[bp_choice].bc_func); CHKERRQ(ierr); - // Set up vertex DM - ierr = DMClone(dm, &dm_vertex); CHKERRQ(ierr); - ierr = DMSetVecType(dm_vertex, vec_type); CHKERRQ(ierr); - ierr = SetupVertexDMFromDM(dm, dm_vertex, num_comp_u, + // Set up subdomain vertex DM + ierr = DMClone(dm, &dm_Pi); CHKERRQ(ierr); + ierr = DMSetVecType(dm_Pi, vec_type); CHKERRQ(ierr); + ierr = SetupVertexDMFromDM(dm, dm_Pi, num_comp_u, bp_options[bp_choice].enforce_bc, bp_options[bp_choice].bc_func); CHKERRQ(ierr); @@ -193,11 +194,11 @@ int main(int argc, char **argv) { ierr = VecGetSize(X, &g_size); CHKERRQ(ierr); ierr = DMCreateLocalVector(dm, &X_loc); CHKERRQ(ierr); ierr = VecGetSize(X_loc &xl_size); CHKERRQ(ierr); - ierr = DMCreateGlobalVector(dm_vertex, &X); CHKERRQ(ierr); - ierr = VecGetLocalSize(X_vertex, &l_vertex_size); CHKERRQ(ierr); - ierr = VecGetSize(X_vertex, &g_vertex_size); CHKERRQ(ierr); - ierr = DMCreateLocalVector(dm_vertex, &X_vertex_loc); CHKERRQ(ierr); - ierr = VecGetSize(X_vertex_loc &xl_vertex_size); CHKERRQ(ierr); + ierr = DMCreateGlobalVector(dm_Pi, &X); CHKERRQ(ierr); + ierr = VecGetLocalSize(X_Pi, &l_Pi_size); CHKERRQ(ierr); + ierr = VecGetSize(X_Pi, &g_Pi_size); CHKERRQ(ierr); + ierr = DMCreateLocalVector(dm_Pi, &X_Pi_loc); CHKERRQ(ierr); + ierr = VecGetSize(X_Pi_loc &xl_Pi_size); CHKERRQ(ierr); // Operator ierr = PetscMalloc1(1, &user_O); CHKERRQ(ierr); @@ -210,18 +211,18 @@ int main(int argc, char **argv) { ierr = MatShellSetVecType(mat_O, vec_type); CHKERRQ(ierr); // Interface vertex operator - ierr = PetscMalloc1(1, &user_vertex); CHKERRQ(ierr); - ierr = MatCreateShell(comm, l_vertex_size, l_vertex_size, g_vertex_size, - g_vertex_size, user_vertex, &mat_vertex); CHKERRQ(ierr); + ierr = PetscMalloc1(1, &user_Pi); CHKERRQ(ierr); + ierr = MatCreateShell(comm, l_Pi_size, l_Pi_size, g_Pi_size, + g_Pi_size, user_Pi, &mat_Pi); CHKERRQ(ierr); ierr = MatShellSetOperation(mat_O, MATOP_MULT, (void(*)(void))MatMult_Ceed); CHKERRQ(ierr); - ierr = MatShellSetOperation(mat_vertex, MATOP_GET_DIAGONAL, + ierr = MatShellSetOperation(mat_Pi, MATOP_GET_DIAGONAL, (void(*)(void))MatGetDiag); CHKERRQ(ierr); - ierr = MatShellSetVecType(mat_vertex, vec_type); CHKERRQ(ierr); + ierr = MatShellSetVecType(mat_Pi, vec_type); CHKERRQ(ierr); // Injection operator ierr = PetscMalloc1(1, &user_pr); CHKERRQ(ierr); - ierr = MatCreateShell(comm, l_size, l_vertex_size, g_size, g_vertex_size, + ierr = MatCreateShell(comm, l_size, l_Pi_size, g_size, g_Pi_size, user_pr, &mat_pr); CHKERRQ(ierr); ierr = MatShellSetOperation(mat_pr, MATOP_MULT, (void(*)(void))MatMult_Inject); CHKERRQ(ierr); @@ -273,9 +274,9 @@ int main(int argc, char **argv) { CHKERRQ(ierr); // Set up libCEED operator on interface vertices - ierr = PetscMalloc1(1, &ceed_data_vertex); CHKERRQ(ierr); - ierr = SetupLibceedBDDC(dm, ceed_data, ceed_data_vertex, g_vertex_size, - xl_vertex_size, bp_options[bp_choice]); + ierr = PetscMalloc1(1, &ceed_data_bddc); CHKERRQ(ierr); + ierr = SetupLibceedBDDC(dm, ceed_data, ceed_data_bddc, g_Pi_size, + xl_Pi_size, bp_options[bp_choice]); CHKERRQ(ierr); // Gather RHS @@ -291,12 +292,6 @@ int main(int argc, char **argv) { CeedQFunctionCreateIdentity(ceed, num_comp_u, CEED_EVAL_INTERP, CEED_EVAL_NONE, &qf_prolong); - // Set up libCEED BDDC injection operators -// TODO - ierr = CeedInjectionSetup(ceed, num_comp_u, ceed_data, qf_restrict, qf_prolong); - CHKERRQ(ierr); -// TODO - // Create the error QFunction CeedQFunctionCreateInterior(ceed, 1, bp_options[bp_choice].error, bp_options[bp_choice].error_loc, &qf_error); @@ -315,8 +310,96 @@ int main(int argc, char **argv) { CeedOperatorSetField(op_error, "error", ceed_data[fine_level]->elem_restr_u_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); - // Set up Mat - // User Operator + // PETSc pointwise mult vectors + // -- Calculate multiplicity + { + ierr = VecSet(X_loc, 1.0); CHKERRQ(ierr); + + // Local-to-global + ierr = VecZeroEntries(X); CHKERRQ(ierr); + ierr = DMLocalToGlobal(dm, X_loc, ADD_VALUES, X); + CHKERRQ(ierr); + ierr = VecZeroEntries(X_loc); CHKERRQ(ierr); + + // Global-to-local + ierr = DMGlobalToLocal(dm, X, INSERT_VALUES, X_loc); + CHKERRQ(ierr); + ierr = VecZeroEntries(X); CHKERRQ(ierr); + + // CEED vector + PetscScalar *x; + ierr = VecGetArrayAndMemType(X_loc, &x, &mem_type); CHKERRQ(ierr); + CeedInt len; + CeedVectorGetLength(ceed_data->x_ceed, &len);) + CeedVectorCreate(ceed_data->ceed, len, &mult_ceed); + CeedVectorSetArray(mult_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, x); + + // Multiplicity + CeedVector e_vec; + CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, NULL, &e_vec); + CeedVectorSetValue(e_vec, 0.0); + CeedElemRestrictionApply(ceed_data->elem_restr_u, CEED_NOTRANSPOSE, mult_ceed, + e_vec, CEED_REQUEST_IMMEDIATE); + CeedVectorSetValue(mult_ceed, 0.0); + CeedElemRestrictionApply(ceed_data->elem_restr_u, CEED_TRANSPOSE, e_vec, + mult_ceed, CEED_REQUEST_IMMEDIATE); + CeedVectorSyncArray(mult_ceed, MemTypeP2C(mem_type)); + CeedVectorDestroy(&e_vec); + + // Restore vector + ierr = VecRestoreArrayAndMemType(X_loc, &x); CHKERRQ(ierr); + + // Multiplicity scaling + ierr = VecReciprocal(X_loc); + + // Copy to Ceed vector + ierr = VecGetArrayAndMemType(X_loc, &x, &mem_type); CHKERRQ(ierr); + CeedVectorSetArray(mult_ceed, MemTypeP2C(mem_type), CEED_COPY_VALUES, x); + ierr = VecRestoreArrayAndMemType(X_loc, &x); CHKERRQ(ierr); + ierr = VecZeroEntries(X_loc); CHKERRQ(ierr); + } + + // -- Masks for subdomains + { + CeedInt length_r; + CeedVectorGetLength(ceed_data_bddc->x_r_ceed, &length_r);) + ierr = VecDuplicate(X_loc, &mask_r); CHKERRQ(ierr); + ierr = VecSetSizes(mask_r, length_r, PETSC_DECIDE); CHKERRQ(ierr); + ierr = VecDuplicate(mask_r, &mask_Gamma); CHKERRQ(ierr); + ierr = VecDuplicate(mask_r, &mask_I); CHKERRQ(ierr); + + // Set mask contents + CeedScalar *mask_r_array, *mask_Gamma_array, *mask_I_array; + CeedInt num_elem, elem_size; + CeedElemRestrictionGetNumElements(ceed_data_bddc->elem_restr_r, &num_elem); + CeedElemRestrictionGetElementSize(ceed_data_bddc->elem_restr_r, &elem_size); + ierr = VecGetArray(mask_r_ceed, &mask_r_array); + ierr = VecGetArray(mask_Gamma_ceed, &mask_Gamma_array); + ierr = VecGetArray(mask_I_ceed, &mask_I_array); + for (CeedInt e=0; estrides[0]*n + ceed_data_bddc->strides[1]*c + + ceed_data_bddc->strides[2]*e; + mask_r_array[index] = r ? 1.0 : 0.0; + mask_Gamma_array[index] = Gamma ? 1.0 : 0.0; + mask_I_array[index] = I ? 1.0 : 0.0; + } + } + } + ierr = VecRestoreArray(mask_r_ceed, &mask_r_array); CHKERRQ(ierr); + ierr = VecRestoreArray(mask_Gamma_ceed, &mask_Gamma_array); CHKERRQ(ierr); + ierr = VecRestoreArray(mask_I_ceed, &mask_I_array); CHKERRQ(ierr); + } + + // Set up MatShell user data user_O->comm = comm; user_O->dm = dm; user_O->X_loc = X_loc; @@ -326,18 +409,9 @@ int main(int argc, char **argv) { user_O->op = ceed_data->op_apply; user_O->ceed = ceed; - // Injection/Restriction Operator - user_pr->comm = comm; - user_pr->dmf = dm; - user_pr->dmc = dm_vertex; - user_pr->loc_vec_c = X_loc; - user_pr->loc_vec_f = user_O->Y_loc; - user_pr->mult_vec = mult; - user_pr->ceed_vec_c = user_O->x_ceed; - user_pr->ceed_vec_f = user_O->y_ceed; - user_pr->op_prolong = ceed_data->op_prolong; - user_pr->op_restrict = ceed_data->op_restrict; - user_pr->ceed = ceed; +// TODO + // Set up PCShell user data +// TODO // Set up KSP ierr = KSPCreate(comm, &ksp); CHKERRQ(ierr); @@ -354,10 +428,10 @@ int main(int argc, char **argv) { ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr); { ierr = PCSetType(pc, PCSHELL); CHKERRQ(ierr); + ierr = PCShellSetContext(pc, user_bddc); CHKERRQ(ierr); + ierr = PCShellSetApply(pc, (void(*)(void))PCShellApply_BDDC); CHKERRQ(ierr); // TODO - ierr = PCShellSetContext(pc, ); CHKERRQ(ierr); - ierr = PCShellSetApply(pc, ); CHKERRQ(ierr); - ierr = PCShellSetSetup(pc, ); CHKERRQ(ierr); + //ierr = PCShellSetSetup(pc, ); CHKERRQ(ierr); // TODO } @@ -471,9 +545,9 @@ int main(int argc, char **argv) { ierr = MatDestroy(&mat_pr); CHKERRQ(ierr); ierr = PetscFree(user_pr); CHKERRQ(ierr); ierr = CeedDataDestroy(i, ceed_data); CHKERRQ(ierr); - ierr = CeedDataDestroy(i, ceed_data_vertex); CHKERRQ(ierr); + ierr = CeedDataBDDCDestroy(i, ceed_data_bddc); CHKERRQ(ierr); ierr = DMDestroy(&dm); CHKERRQ(ierr); - ierr = DMDestroy(&dm_vertex); CHKERRQ(ierr); + ierr = DMDestroy(&dm_Pi); CHKERRQ(ierr); ierr = VecDestroy(&rhs); CHKERRQ(ierr); ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr); ierr = KSPDestroy(&ksp); CHKERRQ(ierr); diff --git a/examples/petsc/include/matops.h b/examples/petsc/include/matops.h index 7582e44755..4215bfa74d 100644 --- a/examples/petsc/include/matops.h +++ b/examples/petsc/include/matops.h @@ -13,8 +13,7 @@ PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y); PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y); PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y); -PetscErrorCode MatMult_Inject(Mat A, Vec X, Vec Y); -PetscErrorCode MatMult_Inject_T(Mat A, Vec X, Vec Y); +PetscErrorCode PCShellApply_BDDC(Mat A, Vec X, Vec Y); PetscErrorCode ComputeErrorMax(UserO user, CeedOperator op_error, Vec X, CeedVector target, PetscReal *max_error); diff --git a/examples/petsc/include/structs.h b/examples/petsc/include/structs.h index 5144220b8d..e2dc010990 100644 --- a/examples/petsc/include/structs.h +++ b/examples/petsc/include/structs.h @@ -29,15 +29,14 @@ struct UserProlongRestr_ { CeedOperator op_prolong, op_restrict; Ceed ceed; }; - -// Data for PETSc Inject/Inject_T Matshells -typedef struct UserInject_ *UserInject; -struct UserInject_ { +// Data for PETSc PCshell +typedef struct UserBDDC_ *UserBDDC; +struct UserBDDC_ { MPI_Comm comm; - DM dmc, dmf; - Vec loc_vec_f, loc_vec_Pi, mult_vec; - CeedVector ceed_vec_f, ceed_vec_Pi, ceed_vec_r; - CeedOperator op_inject, op_inject_t; + DM dm, dm_Pi; + Vec X_loc, Y_loc, diag; + CeedVector x_ceed, y_ceed; + CeedOperator op; Ceed ceed; }; @@ -60,9 +59,10 @@ struct CeedData_ { typedef struct CeedDataBDDC_ *CeedDataBDDC; struct CeedDataBDDC_ { CeedBasis basis_Pi; + CeedInt strides[3]; CeedElemRestriction elem_restr_Pi, elem_restr_r; CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv; - CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mask_r, mask_Gamma, mask_I; + CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mult_ceed; }; // BP specific data diff --git a/examples/petsc/multigrid.c b/examples/petsc/multigrid.c index 5ec02e3266..af3c166155 100644 --- a/examples/petsc/multigrid.c +++ b/examples/petsc/multigrid.c @@ -40,6 +40,23 @@ const char help[] = "Solve CEED BPs using p-multigrid with PETSc and DMPlex\n"; #include "bps.h" +// The BDDC example uses vectors in three spaces +// +// Fine mesh: Broken mesh: Vertex mesh: +// x----x----x x----x x----x x x x +// | | | | | | | +// | | | | | | | +// | | | | | | | +// | | | | | | | +// x----x----x x----x x----x x x x +// +// Vectors are organized as follows +// - *_Pi : vector on the vertex mesh +// - *_r : vector on the broken mesh, all points but vertices +// - *_Gamma : vector on the broken mesh, face/vertex/edge points +// - *_I : vector on the broken mesh, interior points +// - * : all other vectors are on the fine mesh + int main(int argc, char **argv) { PetscInt ierr; MPI_Comm comm; diff --git a/examples/petsc/src/libceedsetup.c b/examples/petsc/src/libceedsetup.c index 0f0e4cf268..3f1c3e708e 100644 --- a/examples/petsc/src/libceedsetup.c +++ b/examples/petsc/src/libceedsetup.c @@ -263,7 +263,7 @@ BPData bp_data { CeedBasis basis_Pi, basis_u = data_fine->basis_u; CeedElemRestriction elem_restr_Pi, elem_restr_r; CeedOperator op_Pi_r, op_r_Pi, op_Pi_Pi, op_r_r, op_r_r_inv,; - CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mask_r, mask_Gamma, mask_I; + CeedVector x_Pi_ceed, y_Pi_ceed, x_r_ceed, y_r_ceed, mask_r_ceed, mask_Gamma_ceed, mask_I_ceed; CeedInt topo_dim, num_comp_u, P, Q, num_qpts, num_elem, elem_size, q_data_size = bp_data.q_data_size; @@ -303,38 +303,10 @@ BPData bp_data { ierr = DMPlexGetHeightStratum(dm_vertex, 0, &c_start, &c_end); CHKERRQ(ierr); num_elem = c_end - c_start; CeedInt strides = [num_comp_u, 1, num_comp_u*elem_size]; - CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, num_comp_u, + CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, num_comp_u, num_comp_u *num_elem*elem_size, strides, &elem_restr_r); - // -- Masks for subdomains - CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &mask_r); - CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &mask_Gamma); - CeedVectorCreate(ceed, num_comp_u *elem_size*num_elem, &mask_I); - CeedScalar *mask_r_array, *mask_Gamma_array, *mask_I_array; - CeedVectorGetArray(mask_r, CEED_MEM_HOST, &mask_r_array); - CeedVectorGetArray(mask_Gamma, CEED_MEM_HOST, &mask_Gamma_array); - CeedVectorGetArray(mask_I, CEED_MEM_HOST, &mask_I_array); - for (CeedInt e=0; ey_Pi_ceed = y_Pi_ceed; data_vertex->x_r_ceed = x_r_ceed; data_vertex->y_r_ceed = y_r_ceed; - data_vertex->mask_r = mask_r; - data_vertex->mask_Gamma = mask_Gamma; - data_vertex->mask_I = mask_I; PetscFunctionReturn(0); }; diff --git a/examples/petsc/src/matops.c b/examples/petsc/src/matops.c index d71366295a..4b0761a23b 100644 --- a/examples/petsc/src/matops.c +++ b/examples/petsc/src/matops.c @@ -205,6 +205,49 @@ PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) { PetscFunctionReturn(0); }; +// ----------------------------------------------------------------------------- +// This function uses libCEED to compute the action of the BDDC preconditioner +// ----------------------------------------------------------------------------- +PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y) { + PetscErrorCode ierr; + UserBDDC user; + + PetscFunctionBeginUser; + + // Inject to broken space + // -- Scaled injection, point multiply by 1/multiplicity + // -- Harmonic injection, scaled with jump map + // ---- A_I,I^-1 + // ---- A_Gamma,I + // ---- J^T (jump map) + // ---- X_r -= J^T A_Gamma,I A_I,I^-1 X_r + // ---- X_Pi copy nodal values from X_r + + // K_u^-1 - update nodal values from subdomain + // -- A_r,r^-1 + // -- A_Pi,r + // -- X_Pi -= A_Pi_r A_r,r^-1 X_Pi + + // P - subdomain and Schur compliment solve + // -- X_r = A_r,r^-1 X_r + // -- X_Pi = S_Pi^-1 + + // K_u^-T - update subdomain values from nodal + // -- A_r,Pi + // -- A_r,r^-1 + // -- X_r -= A_r,r^-1 A_r,Pi X_Pi + + // Restrict to fine space + // -- Scaled restriction, point multiply by 1/multiplicity + // -- Harmonic injection, scaled with jump map + // ---- J^T (jump map) + // ---- A_I,Gamma + // ---- A_I,I^-1 + // ---- X -= A_I,I^-1 A_Gamma,I J^T X_r + + PetscFunctionReturn(0); +}; + // ----------------------------------------------------------------------------- // This function calculates the error in the final solution // -----------------------------------------------------------------------------