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