1 #include "../include/matops.h" 2 #include "../include/petscutils.h" 3 4 // ----------------------------------------------------------------------------- 5 // This function returns the computed diagonal of the operator 6 // ----------------------------------------------------------------------------- 7 PetscErrorCode MatGetDiag(Mat A, Vec D) { 8 PetscErrorCode ierr; 9 OperatorApplyContext op_apply_ctx; 10 11 PetscFunctionBeginUser; 12 13 ierr = MatShellGetContext(A, &op_apply_ctx); CHKERRQ(ierr); 14 15 // Compute Diagonal via libCEED 16 PetscScalar *y; 17 PetscMemType mem_type; 18 19 // -- Place PETSc vector in libCEED vector 20 ierr = VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &mem_type); CHKERRQ(ierr); 21 CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(mem_type), 22 CEED_USE_POINTER, y); 23 24 // -- Compute Diagonal 25 CeedOperatorLinearAssembleDiagonal(op_apply_ctx->op, op_apply_ctx->y_ceed, 26 CEED_REQUEST_IMMEDIATE); 27 28 // -- Local-to-Global 29 CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(mem_type), NULL); 30 ierr = VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y); CHKERRQ(ierr); 31 ierr = VecZeroEntries(D); CHKERRQ(ierr); 32 ierr = DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, D); 33 CHKERRQ(ierr); 34 35 PetscFunctionReturn(0); 36 }; 37 38 // ----------------------------------------------------------------------------- 39 // This function uses libCEED to compute the action of the Laplacian with 40 // Dirichlet boundary conditions 41 // ----------------------------------------------------------------------------- 42 PetscErrorCode ApplyLocal_Ceed(Vec X, Vec Y, 43 OperatorApplyContext op_apply_ctx) { 44 PetscErrorCode ierr; 45 PetscScalar *x, *y; 46 PetscMemType x_mem_type, y_mem_type; 47 48 PetscFunctionBeginUser; 49 50 // Global-to-local 51 ierr = DMGlobalToLocal(op_apply_ctx->dm, X, INSERT_VALUES, op_apply_ctx->X_loc); 52 CHKERRQ(ierr); 53 54 // Setup libCEED vectors 55 ierr = VecGetArrayReadAndMemType(op_apply_ctx->X_loc, (const PetscScalar **)&x, 56 &x_mem_type); CHKERRQ(ierr); 57 ierr = VecGetArrayAndMemType(op_apply_ctx->Y_loc, &y, &y_mem_type); 58 CHKERRQ(ierr); 59 CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), 60 CEED_USE_POINTER, x); 61 CeedVectorSetArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), 62 CEED_USE_POINTER, y); 63 64 // Apply libCEED operator 65 CeedOperatorApply(op_apply_ctx->op, op_apply_ctx->x_ceed, op_apply_ctx->y_ceed, 66 CEED_REQUEST_IMMEDIATE); 67 68 // Restore PETSc vectors 69 CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(x_mem_type), NULL); 70 CeedVectorTakeArray(op_apply_ctx->y_ceed, MemTypeP2C(y_mem_type), NULL); 71 ierr = VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, 72 (const PetscScalar **)&x); 73 CHKERRQ(ierr); 74 ierr = VecRestoreArrayAndMemType(op_apply_ctx->Y_loc, &y); CHKERRQ(ierr); 75 76 // Local-to-global 77 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 78 ierr = DMLocalToGlobal(op_apply_ctx->dm, op_apply_ctx->Y_loc, ADD_VALUES, Y); 79 CHKERRQ(ierr); 80 81 PetscFunctionReturn(0); 82 }; 83 84 // ----------------------------------------------------------------------------- 85 // This function wraps the libCEED operator for a MatShell 86 // ----------------------------------------------------------------------------- 87 PetscErrorCode MatMult_Ceed(Mat A, Vec X, Vec Y) { 88 PetscErrorCode ierr; 89 OperatorApplyContext op_apply_ctx; 90 91 PetscFunctionBeginUser; 92 93 ierr = MatShellGetContext(A, &op_apply_ctx); CHKERRQ(ierr); 94 95 // libCEED for local action of residual evaluator 96 ierr = ApplyLocal_Ceed(X, Y, op_apply_ctx); CHKERRQ(ierr); 97 98 PetscFunctionReturn(0); 99 }; 100 101 // ----------------------------------------------------------------------------- 102 // This function wraps the libCEED operator for a SNES residual evaluation 103 // ----------------------------------------------------------------------------- 104 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 105 PetscErrorCode ierr; 106 OperatorApplyContext op_apply_ctx = (OperatorApplyContext)ctx; 107 108 PetscFunctionBeginUser; 109 110 // libCEED for local action of residual evaluator 111 ierr = ApplyLocal_Ceed(X, Y, op_apply_ctx); CHKERRQ(ierr); 112 113 PetscFunctionReturn(0); 114 }; 115 116 // ----------------------------------------------------------------------------- 117 // This function uses libCEED to compute the action of the prolongation operator 118 // ----------------------------------------------------------------------------- 119 PetscErrorCode MatMult_Prolong(Mat A, Vec X, Vec Y) { 120 PetscErrorCode ierr; 121 ProlongRestrContext pr_restr_ctx; 122 PetscScalar *c, *f; 123 PetscMemType c_mem_type, f_mem_type; 124 125 PetscFunctionBeginUser; 126 127 ierr = MatShellGetContext(A, &pr_restr_ctx); CHKERRQ(ierr); 128 129 // Global-to-local 130 ierr = VecZeroEntries(pr_restr_ctx->loc_vec_c); CHKERRQ(ierr); 131 ierr = DMGlobalToLocal(pr_restr_ctx->dmc, X, INSERT_VALUES, 132 pr_restr_ctx->loc_vec_c); 133 CHKERRQ(ierr); 134 135 // Setup libCEED vectors 136 ierr = VecGetArrayReadAndMemType(pr_restr_ctx->loc_vec_c, 137 (const PetscScalar **)&c, 138 &c_mem_type); CHKERRQ(ierr); 139 ierr = VecGetArrayAndMemType(pr_restr_ctx->loc_vec_f, &f, &f_mem_type); 140 CHKERRQ(ierr); 141 CeedVectorSetArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), 142 CEED_USE_POINTER, c); 143 CeedVectorSetArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), 144 CEED_USE_POINTER, f); 145 146 // Apply libCEED operator 147 CeedOperatorApply(pr_restr_ctx->op_prolong, pr_restr_ctx->ceed_vec_c, 148 pr_restr_ctx->ceed_vec_f, CEED_REQUEST_IMMEDIATE); 149 150 // Restore PETSc vectors 151 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 152 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 153 ierr = VecRestoreArrayReadAndMemType(pr_restr_ctx->loc_vec_c, 154 (const PetscScalar **)&c); 155 CHKERRQ(ierr); 156 ierr = VecRestoreArrayAndMemType(pr_restr_ctx->loc_vec_f, &f); CHKERRQ(ierr); 157 158 // Multiplicity 159 ierr = VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, 160 pr_restr_ctx->mult_vec); 161 162 // Local-to-global 163 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 164 ierr = DMLocalToGlobal(pr_restr_ctx->dmf, pr_restr_ctx->loc_vec_f, ADD_VALUES, 165 Y); CHKERRQ(ierr); 166 167 PetscFunctionReturn(0); 168 }; 169 170 // ----------------------------------------------------------------------------- 171 // This function uses libCEED to compute the action of the restriction operator 172 // ----------------------------------------------------------------------------- 173 PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) { 174 PetscErrorCode ierr; 175 ProlongRestrContext pr_restr_ctx; 176 PetscScalar *c, *f; 177 PetscMemType c_mem_type, f_mem_type; 178 179 PetscFunctionBeginUser; 180 181 ierr = MatShellGetContext(A, &pr_restr_ctx); CHKERRQ(ierr); 182 183 // Global-to-local 184 ierr = VecZeroEntries(pr_restr_ctx->loc_vec_f); CHKERRQ(ierr); 185 ierr = DMGlobalToLocal(pr_restr_ctx->dmf, X, INSERT_VALUES, 186 pr_restr_ctx->loc_vec_f); 187 CHKERRQ(ierr); 188 189 // Multiplicity 190 ierr = VecPointwiseMult(pr_restr_ctx->loc_vec_f, pr_restr_ctx->loc_vec_f, 191 pr_restr_ctx->mult_vec); 192 CHKERRQ(ierr); 193 194 // Setup libCEED vectors 195 ierr = VecGetArrayReadAndMemType(pr_restr_ctx->loc_vec_f, 196 (const PetscScalar **)&f, 197 &f_mem_type); CHKERRQ(ierr); 198 ierr = VecGetArrayAndMemType(pr_restr_ctx->loc_vec_c, &c, &c_mem_type); 199 CHKERRQ(ierr); 200 CeedVectorSetArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), 201 CEED_USE_POINTER, f); 202 CeedVectorSetArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), 203 CEED_USE_POINTER, c); 204 205 // Apply CEED operator 206 CeedOperatorApply(pr_restr_ctx->op_restrict, pr_restr_ctx->ceed_vec_f, 207 pr_restr_ctx->ceed_vec_c, CEED_REQUEST_IMMEDIATE); 208 209 // Restore PETSc vectors 210 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 211 CeedVectorTakeArray(pr_restr_ctx->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 212 ierr = VecRestoreArrayReadAndMemType(pr_restr_ctx->loc_vec_f, 213 (const PetscScalar **)&f); CHKERRQ(ierr); 214 ierr = VecRestoreArrayAndMemType(pr_restr_ctx->loc_vec_c, &c); CHKERRQ(ierr); 215 216 // Local-to-global 217 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 218 ierr = DMLocalToGlobal(pr_restr_ctx->dmc, pr_restr_ctx->loc_vec_c, ADD_VALUES, 219 Y); 220 CHKERRQ(ierr); 221 222 PetscFunctionReturn(0); 223 }; 224 225 // ----------------------------------------------------------------------------- 226 // This function calculates the error in the final solution 227 // ----------------------------------------------------------------------------- 228 PetscErrorCode ComputeErrorMax(OperatorApplyContext op_apply_ctx, 229 CeedOperator op_error, 230 Vec X, CeedVector target, 231 PetscScalar *max_error) { 232 PetscErrorCode ierr; 233 PetscScalar *x; 234 PetscMemType mem_type; 235 CeedVector collocated_error; 236 CeedSize length; 237 238 PetscFunctionBeginUser; 239 CeedVectorGetLength(target, &length); 240 CeedVectorCreate(op_apply_ctx->ceed, length, &collocated_error); 241 242 // Global-to-local 243 ierr = DMGlobalToLocal(op_apply_ctx->dm, X, INSERT_VALUES, op_apply_ctx->X_loc); 244 CHKERRQ(ierr); 245 246 // Setup CEED vector 247 ierr = VecGetArrayAndMemType(op_apply_ctx->X_loc, &x, &mem_type); CHKERRQ(ierr); 248 CeedVectorSetArray(op_apply_ctx->x_ceed, MemTypeP2C(mem_type), 249 CEED_USE_POINTER, x); 250 251 // Apply CEED operator 252 CeedOperatorApply(op_error, op_apply_ctx->x_ceed, collocated_error, 253 CEED_REQUEST_IMMEDIATE); 254 255 // Restore PETSc vector 256 CeedVectorTakeArray(op_apply_ctx->x_ceed, MemTypeP2C(mem_type), NULL); 257 ierr = VecRestoreArrayReadAndMemType(op_apply_ctx->X_loc, 258 (const PetscScalar **)&x); 259 CHKERRQ(ierr); 260 261 // Reduce max error 262 *max_error = 0; 263 const CeedScalar *e; 264 CeedVectorGetArrayRead(collocated_error, CEED_MEM_HOST, &e); 265 for (CeedInt i=0; i<length; i++) { 266 *max_error = PetscMax(*max_error, PetscAbsScalar(e[i])); 267 } 268 CeedVectorRestoreArrayRead(collocated_error, &e); 269 ierr = MPI_Allreduce(MPI_IN_PLACE, max_error, 1, MPIU_REAL, MPIU_MAX, 270 op_apply_ctx->comm); 271 CHKERRQ(ierr); 272 273 // Cleanup 274 CeedVectorDestroy(&collocated_error); 275 276 PetscFunctionReturn(0); 277 }; 278 279 // ----------------------------------------------------------------------------- 280