1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 /// @file 18 /// Matrix shell operations for solid mechanics example using PETSc 19 20 #include "../elasticity.h" 21 22 // ----------------------------------------------------------------------------- 23 // libCEED Operators for MatShell 24 // ----------------------------------------------------------------------------- 25 // This function uses libCEED to compute the local action of an operator 26 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user) { 27 PetscErrorCode ierr; 28 PetscScalar *x, *y; 29 PetscMemType xmemtype, ymemtype; 30 31 PetscFunctionBeginUser; 32 33 // Global-to-local 34 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 35 ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr); 36 37 // Setup CEED vectors 38 ierr = VecGetArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x, 39 &xmemtype); CHKERRQ(ierr); 40 ierr = VecGetArrayAndMemType(user->Yloc, &y, &ymemtype); CHKERRQ(ierr); 41 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x); 42 CeedVectorSetArray(user->Yceed, MemTypeP2C(ymemtype), CEED_USE_POINTER, y); 43 44 // Apply CEED operator 45 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE); 46 47 // Restore PETSc vectors 48 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL); 49 CeedVectorTakeArray(user->Yceed, MemTypeP2C(ymemtype), NULL); 50 ierr = VecRestoreArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x); 51 CHKERRQ(ierr); 52 ierr = VecRestoreArrayAndMemType(user->Yloc, &y); CHKERRQ(ierr); 53 54 // Local-to-global 55 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 56 ierr = DMLocalToGlobal(user->dm, user->Yloc, ADD_VALUES, Y); CHKERRQ(ierr); 57 58 PetscFunctionReturn(0); 59 }; 60 61 // This function uses libCEED to compute the non-linear residual 62 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 63 PetscErrorCode ierr; 64 UserMult user = (UserMult)ctx; 65 66 PetscFunctionBeginUser; 67 68 // Use computed BCs 69 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 70 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 71 user->loadIncrement, NULL, NULL, NULL); 72 CHKERRQ(ierr); 73 74 // libCEED for local action of residual evaluator 75 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 76 77 // Neumann BCs 78 if (user->NBCs) { 79 ierr = VecAXPY(Y, -user->loadIncrement, user->NBCs); CHKERRQ(ierr); 80 } 81 82 PetscFunctionReturn(0); 83 }; 84 85 // This function uses libCEED to apply the Jacobian for assembly via a SNES 86 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 87 PetscErrorCode ierr; 88 UserMult user = (UserMult)ctx; 89 90 PetscFunctionBeginUser; 91 92 // Zero boundary values 93 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 94 95 // libCEED for local action of residual evaluator 96 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 97 98 PetscFunctionReturn(0); 99 }; 100 101 // This function uses libCEED to compute the action of the Jacobian 102 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y) { 103 PetscErrorCode ierr; 104 UserMult user; 105 106 PetscFunctionBeginUser; 107 108 // Zero boundary values 109 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 110 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 111 112 // libCEED for local action of Jacobian 113 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 114 115 PetscFunctionReturn(0); 116 }; 117 118 // This function uses libCEED to compute the action of the prolongation operator 119 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) { 120 PetscErrorCode ierr; 121 UserMultProlongRestr user; 122 PetscScalar *c, *f; 123 PetscMemType cmemtype, fmemtype; 124 125 PetscFunctionBeginUser; 126 127 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 128 129 // Global-to-local 130 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 131 ierr = DMGlobalToLocal(user->dmC, X, INSERT_VALUES, user->locVecC); 132 CHKERRQ(ierr); 133 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr); 134 135 // Setup CEED vectors 136 ierr = VecGetArrayReadAndMemType(user->locVecC, (const PetscScalar **)&c, 137 &cmemtype); CHKERRQ(ierr); 138 ierr = VecGetArrayAndMemType(user->locVecF, &f, &fmemtype); CHKERRQ(ierr); 139 CeedVectorSetArray(user->ceedVecC, MemTypeP2C(cmemtype), CEED_USE_POINTER, c); 140 CeedVectorSetArray(user->ceedVecF, MemTypeP2C(fmemtype), CEED_USE_POINTER, f); 141 142 // Apply CEED operator 143 CeedOperatorApply(user->opProlong, user->ceedVecC, user->ceedVecF, 144 CEED_REQUEST_IMMEDIATE); 145 146 // Restore PETSc vectors 147 CeedVectorTakeArray(user->ceedVecC, MemTypeP2C(cmemtype), NULL); 148 CeedVectorTakeArray(user->ceedVecF, MemTypeP2C(fmemtype), NULL); 149 ierr = VecRestoreArrayReadAndMemType(user->locVecC, (const PetscScalar **)&c); 150 CHKERRQ(ierr); 151 ierr = VecRestoreArrayAndMemType(user->locVecF, &f); CHKERRQ(ierr); 152 153 // Local-to-global 154 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 155 ierr = DMLocalToGlobal(user->dmF, user->locVecF, ADD_VALUES, Y); 156 CHKERRQ(ierr); 157 158 PetscFunctionReturn(0); 159 } 160 161 // This function uses libCEED to compute the action of the restriction operator 162 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) { 163 PetscErrorCode ierr; 164 UserMultProlongRestr user; 165 PetscScalar *c, *f; 166 PetscMemType cmemtype, fmemtype; 167 168 PetscFunctionBeginUser; 169 170 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 171 172 // Global-to-local 173 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr); 174 ierr = DMGlobalToLocal(user->dmF, X, INSERT_VALUES, user->locVecF); 175 CHKERRQ(ierr); 176 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 177 178 // Setup CEED 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->ceedVecF, MemTypeP2C(fmemtype), NULL); 191 CeedVectorTakeArray(user->ceedVecC, MemTypeP2C(cmemtype), 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 // This function returns the computed diagonal of the operator 205 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) { 206 PetscErrorCode ierr; 207 UserMult user; 208 209 PetscFunctionBeginUser; 210 211 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 212 213 // -- Set physics context 214 if (user->ctxPhysSmoother) 215 CeedQFunctionSetContext(user->qf, user->ctxPhysSmoother); 216 217 // Compute Diagonal via libCEED 218 PetscScalar *x; 219 PetscMemType xmemtype; 220 221 // -- Place PETSc vector in libCEED vector 222 ierr = VecGetArrayAndMemType(user->Xloc, &x, &xmemtype); CHKERRQ(ierr); 223 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x); 224 225 // -- Compute Diagonal 226 CeedOperatorLinearAssembleDiagonal(user->op, user->Xceed, 227 CEED_REQUEST_IMMEDIATE); 228 229 // -- Reset physics context 230 if (user->ctxPhysSmoother) 231 CeedQFunctionSetContext(user->qf, user->ctxPhys); 232 233 // -- Local-to-Global 234 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL); 235 ierr = VecRestoreArrayAndMemType(user->Xloc, &x); CHKERRQ(ierr); 236 ierr = VecZeroEntries(D); CHKERRQ(ierr); 237 ierr = DMLocalToGlobal(user->dm, user->Xloc, ADD_VALUES, D); CHKERRQ(ierr); 238 239 // Cleanup 240 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 241 242 PetscFunctionReturn(0); 243 }; 244 245 // This function calculates the strain energy in the final solution 246 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 247 CeedOperator opEnergy, Vec X, 248 PetscReal *energy) { 249 PetscErrorCode ierr; 250 PetscScalar *x; 251 PetscMemType xmemtype; 252 CeedInt length; 253 254 PetscFunctionBeginUser; 255 256 // Global-to-local 257 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 258 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 259 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 260 user->loadIncrement, NULL, NULL, NULL); 261 CHKERRQ(ierr); 262 263 // Setup libCEED input vector 264 ierr = VecGetArrayReadAndMemType(user->Xloc, (const PetscScalar **)&x, 265 &xmemtype); CHKERRQ(ierr); 266 CeedVectorSetArray(user->Xceed, MemTypeP2C(xmemtype), CEED_USE_POINTER, x); 267 268 // Setup libCEED output vector 269 Vec Eloc; 270 CeedVector eloc; 271 ierr = DMCreateLocalVector(dmEnergy, &Eloc); CHKERRQ(ierr); 272 ierr = VecGetSize(Eloc, &length); CHKERRQ(ierr); 273 ierr = VecDestroy(&Eloc); CHKERRQ(ierr); 274 CeedVectorCreate(user->ceed, length, &eloc); 275 276 // Apply libCEED operator 277 CeedOperatorApply(opEnergy, user->Xceed, eloc, CEED_REQUEST_IMMEDIATE); 278 279 // Restore PETSc vector 280 CeedVectorTakeArray(user->Xceed, MemTypeP2C(xmemtype), NULL); 281 ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); 282 CHKERRQ(ierr); 283 284 // Reduce max error 285 const CeedScalar *e; 286 CeedVectorGetArrayRead(eloc, CEED_MEM_HOST, &e); 287 (*energy) = 0; 288 for (CeedInt i=0; i<length; i++) 289 (*energy) += e[i]; 290 CeedVectorRestoreArrayRead(eloc, &e); 291 CeedVectorDestroy(&eloc); 292 293 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, 294 user->comm); CHKERRQ(ierr); 295 296 PetscFunctionReturn(0); 297 }; 298