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