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(0); 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(0); 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 PetscScalar *y; 54 PetscMemType mem_type; 55 56 // -- Place PETSc vector in libCEED vector 57 PetscCall(VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &mem_type)); 58 CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, y); 59 60 // -- Compute Diagonal 61 CeedOperatorLinearAssembleDiagonal(op_apply_ctx->op, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE); 62 63 // -- Local-to-Global 64 CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(mem_type), NULL); 65 PetscCall(VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y)); 66 PetscCall(VecZeroEntries(D)); 67 PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, D)); 68 69 PetscFunctionReturn(0); 70 }; 71 72 // ----------------------------------------------------------------------------- 73 // This function uses libCEED to compute the action of the Laplacian with 74 // Dirichlet boundary conditions 75 // ----------------------------------------------------------------------------- 76 PetscErrorCode ApplyLocal_Ceed(Vec X, Vec Y, OperatorApplyContext op_apply_ctx) { 77 PetscScalar *x, *y; 78 PetscMemType x_mem_type, y_mem_type; 79 80 PetscFunctionBeginUser; 81 82 // Global-to-local 83 PetscCall(DMGlobalToLocal(op_apply_ctx->dm, X, INSERT_VALUES, op_apply_ctx->X_loc)); 84 85 // Setup libCEED vectors 86 PetscCall(VecGetArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x, &x_mem_type)); 87 PetscCall(VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &y_mem_type)); 88 CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 89 CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 90 91 // Apply libCEED operator 92 CeedOperatorApply(op_apply_ctx->op, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, CEED_REQUEST_IMMEDIATE); 93 94 // Restore PETSc vectors 95 CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), NULL); 96 CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), NULL); 97 PetscCall(VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x)); 98 PetscCall(VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y)); 99 100 // Local-to-global 101 PetscCall(VecZeroEntries(Y)); 102 PetscCall(DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, Y)); 103 104 PetscFunctionReturn(0); 105 }; 106 107 // ----------------------------------------------------------------------------- 108 // This function wraps the libCEED operator for a MatShell 109 // ----------------------------------------------------------------------------- 110 PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y) { 111 OperatorApplyContext op_apply_ctx; 112 113 PetscFunctionBeginUser; 114 115 PetscCall(MatShellGetContext(A, &op_apply_ctx)); 116 117 // libCEED for local action of residual evaluator 118 PetscCall(ApplyLocal_Ceed(X, Y, op_apply_ctx)); 119 120 PetscFunctionReturn(0); 121 }; 122 123 // ----------------------------------------------------------------------------- 124 // This function uses libCEED to compute the action of the prolongation operator 125 // ----------------------------------------------------------------------------- 126 PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y) { 127 ProlongRestrContext pr_restr_ctx; 128 PetscScalar *c, *f; 129 PetscMemType c_mem_type, f_mem_type; 130 131 PetscFunctionBeginUser; 132 133 PetscCall(MatShellGetContext(A, &pr_restr_ctx)); 134 135 // Global-to-local 136 PetscCall(VecZeroEntries(pr_restr_ctx->loc_vec_c)); 137 PetscCall(DMGlobalToLocal(pr_restr_ctx->dmc, X, INSERT_VALUES, pr_restr_ctx->loc_vec_c)); 138 139 // Setup libCEED vectors 140 PetscCall(VecGetArrayReadAndMemType(pr_restr_ctx->loc_vec_c, (const PetscScalar **)&c, &c_mem_type)); 141 PetscCall(VecGetArrayAndMemType(pr_restr_ctx->loc_vec_f, &f, &f_mem_type)); 142 CeedVectorSetArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, c); 143 CeedVectorSetArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, f); 144 145 // Apply libCEED operator 146 CeedOperatorApply(pr_restr_ctx->op_prolong, pr_restr_ctx->ceed_vec_c, pr_restr_ctx->ceed_vec_f, CEED_REQUEST_IMMEDIATE); 147 148 // Restore PETSc vectors 149 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 150 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 151 PetscCall(VecRestoreArrayReadAndMemType(pr_restr_ctx->loc_vec_c, (const PetscScalar **)&c)); 152 PetscCall(VecRestoreArrayAndMemType(pr_restr_ctx->loc_vec_f, &f)); 153 154 // Multiplicity 155 PetscCall(VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, pr_restr_ctx->mult_vec)); 156 157 // Local-to-global 158 PetscCall(VecZeroEntries(Y)); 159 PetscCall(DMLocalToGlobal(pr_restr_ctx->dmf, pr_restr_ctx->loc_vec_f, ADD_VALUES, Y)); 160 161 PetscFunctionReturn(0); 162 }; 163 164 // ----------------------------------------------------------------------------- 165 // This function uses libCEED to compute the action of the restriction operator 166 // ----------------------------------------------------------------------------- 167 PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) { 168 ProlongRestrContext pr_restr_ctx; 169 PetscScalar *c, *f; 170 PetscMemType c_mem_type, f_mem_type; 171 172 PetscFunctionBeginUser; 173 174 PetscCall(MatShellGetContext(A, &pr_restr_ctx)); 175 176 // Global-to-local 177 PetscCall(VecZeroEntries(pr_restr_ctx->loc_vec_f)); 178 PetscCall(DMGlobalToLocal(pr_restr_ctx->dmf, X, INSERT_VALUES, pr_restr_ctx->loc_vec_f)); 179 180 // Multiplicity 181 PetscCall(VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, pr_restr_ctx->mult_vec)); 182 183 // Setup libCEED vectors 184 PetscCall(VecGetArrayReadAndMemType(pr_restr_ctx->loc_vec_f, (const PetscScalar **)&f, &f_mem_type)); 185 PetscCall(VecGetArrayAndMemType(pr_restr_ctx->loc_vec_c, &c, &c_mem_type)); 186 CeedVectorSetArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, f); 187 CeedVectorSetArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, c); 188 189 // Apply CEED operator 190 CeedOperatorApply(pr_restr_ctx->op_restrict, pr_restr_ctx->ceed_vec_f, pr_restr_ctx->ceed_vec_c, CEED_REQUEST_IMMEDIATE); 191 192 // Restore PETSc vectors 193 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 194 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 195 PetscCall(VecRestoreArrayReadAndMemType(pr_restr_ctx->loc_vec_f, (const PetscScalar **)&f)); 196 PetscCall(VecRestoreArrayAndMemType(pr_restr_ctx->loc_vec_c, &c)); 197 198 // Local-to-global 199 PetscCall(VecZeroEntries(Y)); 200 PetscCall(DMLocalToGlobal(pr_restr_ctx->dmc, pr_restr_ctx->loc_vec_c, ADD_VALUES, Y)); 201 202 PetscFunctionReturn(0); 203 }; 204 205 // ----------------------------------------------------------------------------- 206 // This function calculates the error in the final solution 207 // ----------------------------------------------------------------------------- 208 PetscErrorCode ComputeL2Error(Vec X, PetscScalar *l2_error, OperatorApplyContext op_error_ctx) { 209 Vec E; 210 PetscFunctionBeginUser; 211 PetscCall(VecDuplicate(X, &E)); 212 PetscCall(ApplyLocal_Ceed(X, E, op_error_ctx)); 213 PetscScalar error_sq = 1.0; 214 PetscCall(VecSum(E, &error_sq)); 215 *l2_error = sqrt(error_sq); 216 217 PetscFunctionReturn(0); 218 }; 219 220 // ----------------------------------------------------------------------------- 221