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 30 PetscFunctionBeginUser; 31 32 // Global-to-local 33 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 34 ierr = VecZeroEntries(user->Yloc); CHKERRQ(ierr); 35 36 // Setup CEED vectors 37 ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); 38 CHKERRQ(ierr); 39 ierr = user->VecGetArray(user->Yloc, &y); CHKERRQ(ierr); 40 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x); 41 CeedVectorSetArray(user->Yceed, user->memType, CEED_USE_POINTER, y); 42 43 // Apply CEED operator 44 // Note: We could use VecGetArrayInPlace. Instead, we use SetArray/TakeArray 45 // so we can request host memory for easier debugging. 46 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE); 47 48 // Restore PETSc vectors 49 CeedVectorTakeArray(user->Xceed, user->memType, NULL); 50 CeedVectorTakeArray(user->Yceed, user->memType, NULL); 51 ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); 52 CHKERRQ(ierr); 53 ierr = user->VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr); 54 55 // Local-to-global 56 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 57 ierr = DMLocalToGlobal(user->dm, user->Yloc, ADD_VALUES, Y); CHKERRQ(ierr); 58 59 PetscFunctionReturn(0); 60 }; 61 62 // This function uses libCEED to compute the non-linear residual 63 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 64 PetscErrorCode ierr; 65 UserMult user = (UserMult)ctx; 66 67 PetscFunctionBeginUser; 68 69 // Use computed BCs 70 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 71 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 72 user->loadIncrement, NULL, NULL, NULL); 73 CHKERRQ(ierr); 74 75 // libCEED for local action of residual evaluator 76 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 77 78 // Neumann BCs 79 if (user->NBCs) { 80 ierr = VecAXPY(Y, -user->loadIncrement, user->NBCs); CHKERRQ(ierr); 81 } 82 83 PetscFunctionReturn(0); 84 }; 85 86 // This function uses libCEED to apply the Jacobian for assembly via a SNES 87 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 88 PetscErrorCode ierr; 89 UserMult user = (UserMult)ctx; 90 91 PetscFunctionBeginUser; 92 93 // Zero boundary values 94 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 95 96 // libCEED for local action of residual evaluator 97 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 98 99 PetscFunctionReturn(0); 100 }; 101 102 // This function uses libCEED to compute the action of the Jacobian 103 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y) { 104 PetscErrorCode ierr; 105 UserMult user; 106 107 PetscFunctionBeginUser; 108 109 // Zero boundary values 110 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 111 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 112 113 // libCEED for local action of Jacobian 114 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 115 116 PetscFunctionReturn(0); 117 }; 118 119 // This function uses libCEED to compute the action of the prolongation operator 120 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) { 121 PetscErrorCode ierr; 122 UserMultProlongRestr user; 123 PetscScalar *c, *f; 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 = user->VecGetArrayRead(user->locVecC, (const PetscScalar **)&c); 137 CHKERRQ(ierr); 138 ierr = user->VecGetArray(user->locVecF, &f); CHKERRQ(ierr); 139 CeedVectorSetArray(user->ceedVecC, user->memType, CEED_USE_POINTER, c); 140 CeedVectorSetArray(user->ceedVecF, user->memType, 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, user->memType, NULL); 148 CeedVectorTakeArray(user->ceedVecF, user->memType, NULL); 149 ierr = user->VecRestoreArrayRead(user->locVecC, (const PetscScalar **)&c); 150 CHKERRQ(ierr); 151 ierr = user->VecRestoreArray(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 167 PetscFunctionBeginUser; 168 169 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 170 171 // Global-to-local 172 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr); 173 ierr = DMGlobalToLocal(user->dmF, X, INSERT_VALUES, user->locVecF); 174 CHKERRQ(ierr); 175 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 176 177 // Setup CEED vectors 178 ierr = user->VecGetArrayRead(user->locVecF, (const PetscScalar **)&f); 179 CHKERRQ(ierr); 180 ierr = user->VecGetArray(user->locVecC, &c); CHKERRQ(ierr); 181 CeedVectorSetArray(user->ceedVecF, user->memType, CEED_USE_POINTER, f); 182 CeedVectorSetArray(user->ceedVecC, user->memType, CEED_USE_POINTER, c); 183 184 // Apply CEED operator 185 CeedOperatorApply(user->opRestrict, user->ceedVecF, user->ceedVecC, 186 CEED_REQUEST_IMMEDIATE); 187 188 // Restore PETSc vectors 189 CeedVectorTakeArray(user->ceedVecF, user->memType, NULL); 190 CeedVectorTakeArray(user->ceedVecC, user->memType, NULL); 191 ierr = user->VecRestoreArrayRead(user->locVecF, (const PetscScalar **)&f); 192 CHKERRQ(ierr); 193 ierr = user->VecRestoreArray(user->locVecC, &c); CHKERRQ(ierr); 194 195 // Local-to-global 196 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 197 ierr = DMLocalToGlobal(user->dmC, user->locVecC, ADD_VALUES, Y); 198 CHKERRQ(ierr); 199 200 PetscFunctionReturn(0); 201 }; 202 203 // This function returns the computed diagonal of the operator 204 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) { 205 PetscErrorCode ierr; 206 UserMult user; 207 208 PetscFunctionBeginUser; 209 210 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 211 212 // -- Set physics context 213 if (user->ctxPhysSmoother) 214 CeedQFunctionSetContext(user->qf, user->ctxPhysSmoother); 215 216 // Compute Diagonal via libCEED 217 PetscScalar *x; 218 219 // -- Place PETSc vector in libCEED vector 220 ierr = user->VecGetArray(user->Xloc, &x); CHKERRQ(ierr); 221 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x); 222 223 // -- Compute Diagonal 224 CeedOperatorLinearAssembleDiagonal(user->op, user->Xceed, 225 CEED_REQUEST_IMMEDIATE); 226 227 // -- Reset physics context 228 if (user->ctxPhysSmoother) 229 CeedQFunctionSetContext(user->qf, user->ctxPhys); 230 231 // -- Local-to-Global 232 CeedVectorTakeArray(user->Xceed, user->memType, NULL); 233 ierr = user->VecRestoreArray(user->Xloc, &x); CHKERRQ(ierr); 234 ierr = VecZeroEntries(D); CHKERRQ(ierr); 235 ierr = DMLocalToGlobal(user->dm, user->Xloc, ADD_VALUES, D); CHKERRQ(ierr); 236 237 // Cleanup 238 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 239 240 PetscFunctionReturn(0); 241 }; 242 243 // This function calculates the strain energy in the final solution 244 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 245 CeedOperator opEnergy, Vec X, 246 PetscReal *energy) { 247 PetscErrorCode ierr; 248 PetscScalar *x; 249 CeedInt length; 250 251 PetscFunctionBeginUser; 252 253 // Global-to-local 254 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 255 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 256 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 257 user->loadIncrement, NULL, NULL, NULL); 258 CHKERRQ(ierr); 259 260 // Setup libCEED input vector 261 ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); 262 CHKERRQ(ierr); 263 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x); 264 265 // Setup libCEED output vector 266 Vec Eloc; 267 CeedVector eloc; 268 ierr = DMCreateLocalVector(dmEnergy, &Eloc); CHKERRQ(ierr); 269 ierr = VecGetSize(Eloc, &length); CHKERRQ(ierr); 270 ierr = VecDestroy(&Eloc); CHKERRQ(ierr); 271 CeedVectorCreate(user->ceed, length, &eloc); 272 273 // Apply libCEED operator 274 CeedOperatorApply(opEnergy, user->Xceed, eloc, CEED_REQUEST_IMMEDIATE); 275 276 // Restore PETSc vector 277 ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); 278 CHKERRQ(ierr); 279 280 // Reduce max error 281 const CeedScalar *e; 282 CeedVectorGetArrayRead(eloc, CEED_MEM_HOST, &e); 283 (*energy) = 0; 284 for (CeedInt i=0; i<length; i++) 285 (*energy) += e[i]; 286 CeedVectorRestoreArrayRead(eloc, &e); 287 CeedVectorDestroy(&eloc); 288 289 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, 290 user->comm); CHKERRQ(ierr); 291 292 PetscFunctionReturn(0); 293 }; 294