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