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 memtype; 18 19 // -- Place PETSc vector in libCEED vector 20 ierr = VecGetArrayAndMemType(user->Xloc, &x, &memtype); CHKERRQ(ierr); 21 CeedVectorSetArray(user->Xceed, MemTypeP2C(memtype), CEED_USE_POINTER, x); 22 23 // -- Compute Diagonal 24 CeedOperatorLinearAssembleDiagonal(user->op, user->Xceed, 25 CEED_REQUEST_IMMEDIATE); 26 27 // -- Local-to-Global 28 CeedVectorTakeArray(user->Xceed, MemTypeP2C(memtype), NULL); 29 ierr = VecRestoreArrayAndMemType(user->Xloc, &x); CHKERRQ(ierr); 30 ierr = VecZeroEntries(D); CHKERRQ(ierr); 31 ierr = DMLocalToGlobal(user->dm, user->Xloc, ADD_VALUES, D); CHKERRQ(ierr); 32 33 // Cleanup 34 ierr = VecZeroEntries(user->Xloc); 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 xmemtype, ymemtype; 47 48 PetscFunctionBeginUser; 49 50 // Global-to-local 51 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 52 53 // Setup libCEED vectors 54 ierr = VecGetArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x, 55 &xmemtype); CHKERRQ(ierr); 56 ierr = VecGetArrayAndMemType(user->Yloc, &y, &ymemtype); CHKERRQ(ierr); 57 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x); 58 CeedVectorSetArray(user->Yceed, MemTypeP2C(ymemtype), CEED_USE_POINTER, y); 59 60 // Apply libCEED operator 61 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE); 62 63 // Restore PETSc vectors 64 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL); 65 CeedVectorTakeArray(user->Yceed, MemTypeP2C(ymemtype), NULL); 66 ierr = VecRestoreArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x); 67 CHKERRQ(ierr); 68 ierr = VecRestoreArrayAndMemType(user->Yloc, &y); CHKERRQ(ierr); 69 70 // Local-to-global 71 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 72 ierr = DMLocalToGlobal(user->dm, user->Yloc, 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 cmemtype, fmemtype; 117 118 PetscFunctionBeginUser; 119 120 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 121 122 // Global-to-local 123 ierr = VecZeroEntries(user->locvecc); CHKERRQ(ierr); 124 ierr = DMGlobalToLocal(user->dmc, X, INSERT_VALUES, user->locvecc); 125 CHKERRQ(ierr); 126 127 // Setup libCEED vectors 128 ierr = VecGetArrayReadAndMemType(user->locvecc, (const PetscScalar **)&c, 129 &cmemtype); CHKERRQ(ierr); 130 ierr = VecGetArrayAndMemType(user->locvecf, &f, &fmemtype); CHKERRQ(ierr); 131 CeedVectorSetArray(user->ceedvecc, MemTypeP2C(cmemtype), CEED_USE_POINTER, c); 132 CeedVectorSetArray(user->ceedvecf, MemTypeP2C(fmemtype), CEED_USE_POINTER, f); 133 134 // Apply libCEED operator 135 CeedOperatorApply(user->opProlong, user->ceedvecc, user->ceedvecf, 136 CEED_REQUEST_IMMEDIATE); 137 138 // Restore PETSc vectors 139 CeedVectorTakeArray(user->ceedvecc, MemTypeP2C(cmemtype), NULL); 140 CeedVectorTakeArray(user->ceedvecf, MemTypeP2C(fmemtype), NULL); 141 ierr = VecRestoreArrayReadAndMemType(user->locvecc, (const PetscScalar **)&c); 142 CHKERRQ(ierr); 143 ierr = VecRestoreArrayAndMemType(user->locvecf, &f); CHKERRQ(ierr); 144 145 // Multiplicity 146 ierr = VecPointwiseMult(user->locvecf, user->locvecf, user->multvec); 147 148 // Local-to-global 149 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 150 ierr = DMLocalToGlobal(user->dmf, user->locvecf, ADD_VALUES, Y); 151 CHKERRQ(ierr); 152 153 PetscFunctionReturn(0); 154 }; 155 156 // ----------------------------------------------------------------------------- 157 // This function uses libCEED to compute the action of the restriction operator 158 // ----------------------------------------------------------------------------- 159 PetscErrorCode MatMult_Restrict(Mat A, Vec X, Vec Y) { 160 PetscErrorCode ierr; 161 UserProlongRestr user; 162 PetscScalar *c, *f; 163 PetscMemType cmemtype, fmemtype; 164 165 PetscFunctionBeginUser; 166 167 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 168 169 // Global-to-local 170 ierr = VecZeroEntries(user->locvecf); CHKERRQ(ierr); 171 ierr = DMGlobalToLocal(user->dmf, X, INSERT_VALUES, user->locvecf); 172 CHKERRQ(ierr); 173 174 // Multiplicity 175 ierr = VecPointwiseMult(user->locvecf, user->locvecf, user->multvec); 176 CHKERRQ(ierr); 177 178 // Setup libCEED vectors 179 ierr = VecGetArrayReadAndMemType(user->locvecf, (const PetscScalar **)&f, 180 &fmemtype); CHKERRQ(ierr); 181 ierr = VecGetArrayAndMemType(user->locvecc, &c, &cmemtype); CHKERRQ(ierr); 182 CeedVectorSetArray(user->ceedvecf, MemTypeP2C(fmemtype), CEED_USE_POINTER, f); 183 CeedVectorSetArray(user->ceedvecc, MemTypeP2C(cmemtype), CEED_USE_POINTER, c); 184 185 // Apply CEED operator 186 CeedOperatorApply(user->opRestrict, user->ceedvecf, user->ceedvecc, 187 CEED_REQUEST_IMMEDIATE); 188 189 // Restore PETSc vectors 190 CeedVectorTakeArray(user->ceedvecc, MemTypeP2C(cmemtype), NULL); 191 CeedVectorTakeArray(user->ceedvecf, MemTypeP2C(fmemtype), NULL); 192 ierr = VecRestoreArrayReadAndMemType(user->locvecf, (const PetscScalar **)&f); 193 CHKERRQ(ierr); 194 ierr = VecRestoreArrayAndMemType(user->locvecc, &c); CHKERRQ(ierr); 195 196 // Local-to-global 197 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 198 ierr = DMLocalToGlobal(user->dmc, user->locvecc, ADD_VALUES, Y); 199 CHKERRQ(ierr); 200 201 PetscFunctionReturn(0); 202 }; 203 204 // ----------------------------------------------------------------------------- 205 // This function calculates the error in the final solution 206 // ----------------------------------------------------------------------------- 207 PetscErrorCode ComputeErrorMax(UserO user, CeedOperator opError, 208 Vec X, CeedVector target, 209 PetscScalar *maxerror) { 210 PetscErrorCode ierr; 211 PetscScalar *x; 212 PetscMemType memtype; 213 CeedVector collocated_error; 214 CeedInt length; 215 216 PetscFunctionBeginUser; 217 CeedVectorGetLength(target, &length); 218 CeedVectorCreate(user->ceed, length, &collocated_error); 219 220 // Global-to-local 221 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 222 223 // Setup CEED vector 224 ierr = VecGetArrayAndMemType(user->Xloc, &x, &memtype); CHKERRQ(ierr); 225 CeedVectorSetArray(user->Xceed, MemTypeP2C(memtype), CEED_USE_POINTER, x); 226 227 // Apply CEED operator 228 CeedOperatorApply(opError, user->Xceed, collocated_error, 229 CEED_REQUEST_IMMEDIATE); 230 231 // Restore PETSc vector 232 CeedVectorTakeArray(user->Xceed, MemTypeP2C(memtype), NULL); 233 ierr = VecRestoreArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x); 234 CHKERRQ(ierr); 235 236 // Reduce max error 237 *maxerror = 0; 238 const CeedScalar *e; 239 CeedVectorGetArrayRead(collocated_error, CEED_MEM_HOST, &e); 240 for (CeedInt i=0; i<length; i++) { 241 *maxerror = PetscMax(*maxerror, PetscAbsScalar(e[i])); 242 } 243 CeedVectorRestoreArrayRead(collocated_error, &e); 244 ierr = MPI_Allreduce(MPI_IN_PLACE, maxerror, 1, MPIU_REAL, MPIU_MAX, 245 user->comm); 246 CHKERRQ(ierr); 247 248 // Cleanup 249 CeedVectorDestroy(&collocated_error); 250 251 PetscFunctionReturn(0); 252 }; 253 254 // ----------------------------------------------------------------------------- 255