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