1 #include "../include/matops.h" 2 3 #include "../include/petscutils.h" 4 5 // ----------------------------------------------------------------------------- 6 // Setup apply operator context data 7 // ----------------------------------------------------------------------------- 8 PetscErrorCode SetupApplyOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, Vec X_loc, OperatorApplyContext op_apply_ctx) { 9 PetscFunctionBeginUser; 10 op_apply_ctx->comm = comm; 11 op_apply_ctx->dm = dm; 12 op_apply_ctx->X_loc = X_loc; 13 PetscCall(VecDuplicate(X_loc, &op_apply_ctx->Y_loc)); 14 op_apply_ctx->x_ceed = ceed_data->x_ceed; 15 op_apply_ctx->y_ceed = ceed_data->y_ceed; 16 op_apply_ctx->op = ceed_data->op_apply; 17 op_apply_ctx->ceed = ceed; 18 PetscFunctionReturn(PETSC_SUCCESS); 19 } 20 21 // ----------------------------------------------------------------------------- 22 // Setup error operator context data 23 // ----------------------------------------------------------------------------- 24 PetscErrorCode SetupErrorOperatorCtx(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, Vec X_loc, CeedOperator op_error, 25 OperatorApplyContext op_error_ctx) { 26 PetscFunctionBeginUser; 27 op_error_ctx->comm = comm; 28 op_error_ctx->dm = dm; 29 op_error_ctx->X_loc = X_loc; 30 PetscCall(VecDuplicate(X_loc, &op_error_ctx->Y_loc)); 31 op_error_ctx->x_ceed = ceed_data->x_ceed; 32 op_error_ctx->y_ceed = ceed_data->y_ceed; 33 op_error_ctx->op = op_error; 34 op_error_ctx->ceed = ceed; 35 PetscFunctionReturn(PETSC_SUCCESS); 36 } 37 38 // ----------------------------------------------------------------------------- 39 // This function returns the computed diagonal of the operator 40 // ----------------------------------------------------------------------------- 41 PetscErrorCode MatGetDiag(Mat A, Vec D) { 42 OperatorApplyContext op_apply_ctx; 43 44 PetscFunctionBeginUser; 45 PetscCall(MatShellGetContext(A, &op_apply_ctx)); 46 47 // Compute Diagonal via libCEED 48 PetscMemType mem_type; 49 50 // Place PETSc vector in libCEED vector 51 PetscCall(VecP2C(op_apply_ctx->Y_loc, &mem_type, op_apply_ctx->y_ceed)); 52 53 // Compute Diagonal 54 CeedOperatorLinearAssembleDiagonal(op_apply_ctx->op, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE); 55 56 // Local-to-Global 57 PetscCall(VecC2P(op_apply_ctx->y_ceed, mem_type, op_apply_ctx->Y_loc)); 58 PetscCall(VecZeroEntries(D)); 59 PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, D)); 60 PetscFunctionReturn(PETSC_SUCCESS); 61 } 62 63 // ----------------------------------------------------------------------------- 64 // This function uses libCEED to compute the action of the Laplacian with Dirichlet boundary conditions 65 // ----------------------------------------------------------------------------- 66 PetscErrorCode ApplyLocal_Ceed(Vec X, Vec Y, OperatorApplyContext op_apply_ctx) { 67 PetscMemType x_mem_type, y_mem_type; 68 69 PetscFunctionBeginUser; 70 // Global-to-local 71 PetscCall(DMGlobalToLocal(op_apply_ctx->dm, X, INSERT_VALUES, op_apply_ctx->X_loc)); 72 73 // Setup libCEED vectors 74 PetscCall(VecReadP2C(op_apply_ctx->X_loc, &x_mem_type, op_apply_ctx->x_ceed)); 75 PetscCall(VecP2C(op_apply_ctx->Y_loc, &y_mem_type, op_apply_ctx->y_ceed)); 76 77 // Apply libCEED operator 78 CeedOperatorApply(op_apply_ctx->op, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE); 79 80 // Restore PETSc vectors 81 PetscCall(VecReadC2P(op_apply_ctx->x_ceed, x_mem_type, op_apply_ctx->X_loc)); 82 PetscCall(VecC2P(op_apply_ctx->y_ceed, y_mem_type, op_apply_ctx->Y_loc)); 83 84 // Local-to-global 85 PetscCall(VecZeroEntries(Y)); 86 PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, Y)); 87 PetscFunctionReturn(PETSC_SUCCESS); 88 } 89 90 // ----------------------------------------------------------------------------- 91 // This function wraps the libCEED operator for a MatShell 92 // ----------------------------------------------------------------------------- 93 PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y) { 94 OperatorApplyContext op_apply_ctx; 95 96 PetscFunctionBeginUser; 97 PetscCall(MatShellGetContext(A, &op_apply_ctx)); 98 99 // libCEED for local action of residual evaluator 100 PetscCall(ApplyLocal_Ceed(X, Y, op_apply_ctx)); 101 PetscFunctionReturn(PETSC_SUCCESS); 102 } 103 104 // ----------------------------------------------------------------------------- 105 // This function uses libCEED to compute the action of the prolongation operator 106 // ----------------------------------------------------------------------------- 107 PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y) { 108 ProlongRestrContext pr_restr_ctx; 109 PetscMemType c_mem_type, f_mem_type; 110 111 PetscFunctionBeginUser; 112 PetscCall(MatShellGetContext(A, &pr_restr_ctx)); 113 114 // Global-to-local 115 PetscCall(VecZeroEntries(pr_restr_ctx->loc_vec_c)); 116 PetscCall(DMGlobalToLocal(pr_restr_ctx->dmc, X, INSERT_VALUES, pr_restr_ctx->loc_vec_c)); 117 118 // Setup libCEED vectors 119 PetscCall(VecReadP2C(pr_restr_ctx->loc_vec_c, &c_mem_type, pr_restr_ctx->ceed_vec_c)); 120 PetscCall(VecP2C(pr_restr_ctx->loc_vec_f, &f_mem_type, pr_restr_ctx->ceed_vec_f)); 121 122 // Apply libCEED operator 123 CeedOperatorApply(pr_restr_ctx->op_prolong, pr_restr_ctx->ceed_vec_c, pr_restr_ctx->ceed_vec_f, CEED_REQUEST_IMMEDIATE); 124 125 // Restore PETSc vectors 126 PetscCall(VecReadC2P(pr_restr_ctx->ceed_vec_c, c_mem_type, pr_restr_ctx->loc_vec_c)); 127 PetscCall(VecC2P(pr_restr_ctx->ceed_vec_f, f_mem_type, pr_restr_ctx->loc_vec_f)); 128 129 // Multiplicity 130 PetscCall(VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, pr_restr_ctx->mult_vec)); 131 132 // Local-to-global 133 PetscCall(VecZeroEntries(Y)); 134 PetscCall(DMLocalToGlobal(pr_restr_ctx->dmf, pr_restr_ctx->loc_vec_f, ADD_VALUES, Y)); 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 // ----------------------------------------------------------------------------- 139 // This function uses libCEED to compute the action of the restriction operator 140 // ----------------------------------------------------------------------------- 141 PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) { 142 ProlongRestrContext pr_restr_ctx; 143 PetscMemType c_mem_type, f_mem_type; 144 145 PetscFunctionBeginUser; 146 PetscCall(MatShellGetContext(A, &pr_restr_ctx)); 147 148 // Global-to-local 149 PetscCall(VecZeroEntries(pr_restr_ctx->loc_vec_f)); 150 PetscCall(DMGlobalToLocal(pr_restr_ctx->dmf, X, INSERT_VALUES, pr_restr_ctx->loc_vec_f)); 151 152 // Multiplicity 153 PetscCall(VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, pr_restr_ctx->mult_vec)); 154 155 // Setup libCEED vectors 156 PetscCall(VecReadP2C(pr_restr_ctx->loc_vec_f, &f_mem_type, pr_restr_ctx->ceed_vec_f)); 157 PetscCall(VecP2C(pr_restr_ctx->loc_vec_c, &c_mem_type, pr_restr_ctx->ceed_vec_c)); 158 159 // Apply CEED operator 160 CeedOperatorApply(pr_restr_ctx->op_restrict, pr_restr_ctx->ceed_vec_f, pr_restr_ctx->ceed_vec_c, CEED_REQUEST_IMMEDIATE); 161 162 // Restore PETSc vectors 163 PetscCall(VecReadC2P(pr_restr_ctx->ceed_vec_f, f_mem_type, pr_restr_ctx->loc_vec_f)); 164 PetscCall(VecC2P(pr_restr_ctx->ceed_vec_c, c_mem_type, pr_restr_ctx->loc_vec_c)); 165 166 // Local-to-global 167 PetscCall(VecZeroEntries(Y)); 168 PetscCall(DMLocalToGlobal(pr_restr_ctx->dmc, pr_restr_ctx->loc_vec_c, ADD_VALUES, Y)); 169 PetscFunctionReturn(PETSC_SUCCESS); 170 } 171 172 // ----------------------------------------------------------------------------- 173 // This function calculates the error in the final solution 174 // ----------------------------------------------------------------------------- 175 PetscErrorCode ComputeL2Error(Vec X, PetscScalar *l2_error, OperatorApplyContext op_error_ctx) { 176 Vec E; 177 178 PetscFunctionBeginUser; 179 PetscCall(VecDuplicate(X, &E)); 180 PetscCall(ApplyLocal_Ceed(X, E, op_error_ctx)); 181 PetscScalar error_sq = 1.0; 182 PetscCall(VecSum(E, &error_sq)); 183 *l2_error = sqrt(error_sq); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 // ----------------------------------------------------------------------------- 188