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: CeedVectorSyncArray is not required when the backend memtype and 45 // the PETSc vector agree. We could in that case use 46 // VecGetArrayInPlace. Instead, we use SetArray/SyncArray so we can 47 // request host memory for easier debugging. 48 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE); 49 CeedVectorSyncArray(user->Yceed, user->memType); 50 51 // Restore PETSc vectors 52 ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); 53 CHKERRQ(ierr); 54 ierr = user->VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr); 55 56 // Local-to-global 57 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 58 ierr = DMLocalToGlobal(user->dm, user->Yloc, ADD_VALUES, Y); CHKERRQ(ierr); 59 60 PetscFunctionReturn(0); 61 }; 62 63 // This function uses libCEED to compute the non-linear residual 64 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 65 PetscErrorCode ierr; 66 UserMult user = (UserMult)ctx; 67 68 PetscFunctionBeginUser; 69 70 // Use computed BCs 71 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 72 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 73 user->loadIncrement, NULL, NULL, NULL); 74 CHKERRQ(ierr); 75 76 // libCEED for local action of residual evaluator 77 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 78 79 PetscFunctionReturn(0); 80 }; 81 82 // This function uses libCEED to apply the Jacobian for assembly via a SNES 83 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 84 PetscErrorCode ierr; 85 UserMult user = (UserMult)ctx; 86 87 PetscFunctionBeginUser; 88 89 // Zero boundary values 90 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 91 92 // libCEED for local action of residual evaluator 93 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 94 95 PetscFunctionReturn(0); 96 }; 97 98 // This function uses libCEED to compute the action of the Jacobian 99 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y) { 100 PetscErrorCode ierr; 101 UserMult user; 102 103 PetscFunctionBeginUser; 104 105 // Zero boundary values 106 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 107 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 108 109 // libCEED for local action of Jacobian 110 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 111 112 PetscFunctionReturn(0); 113 }; 114 115 // This function uses libCEED to compute the action of the prolongation operator 116 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) { 117 PetscErrorCode ierr; 118 UserMultProlongRestr user; 119 PetscScalar *c, *f; 120 121 PetscFunctionBeginUser; 122 123 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 124 125 // Global-to-local 126 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 127 ierr = DMGlobalToLocal(user->dmC, X, INSERT_VALUES, user->locVecC); 128 CHKERRQ(ierr); 129 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr); 130 131 // Setup CEED vectors 132 ierr = user->VecGetArrayRead(user->locVecC, (const PetscScalar **)&c); 133 CHKERRQ(ierr); 134 ierr = user->VecGetArray(user->locVecF, &f); CHKERRQ(ierr); 135 CeedVectorSetArray(user->ceedVecC, user->memType, CEED_USE_POINTER, c); 136 CeedVectorSetArray(user->ceedVecF, user->memType, CEED_USE_POINTER, f); 137 138 // Apply CEED operator 139 CeedOperatorApply(user->opProlong, user->ceedVecC, user->ceedVecF, 140 CEED_REQUEST_IMMEDIATE); 141 CeedVectorSyncArray(user->ceedVecF, user->memType); 142 143 // Restore PETSc vectors 144 ierr = user->VecRestoreArrayRead(user->locVecC, (const PetscScalar **)&c); 145 CHKERRQ(ierr); 146 ierr = user->VecRestoreArray(user->locVecF, &f); CHKERRQ(ierr); 147 148 // Multiplicity 149 ierr = VecPointwiseMult(user->locVecF, user->locVecF, user->multVec); 150 151 // Local-to-global 152 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 153 ierr = DMLocalToGlobal(user->dmF, user->locVecF, ADD_VALUES, Y); 154 CHKERRQ(ierr); 155 156 PetscFunctionReturn(0); 157 } 158 159 // This function uses libCEED to compute the action of the restriction operator 160 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) { 161 PetscErrorCode ierr; 162 UserMultProlongRestr user; 163 PetscScalar *c, *f; 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 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 174 175 // Multiplicity 176 ierr = VecPointwiseMult(user->locVecF, user->locVecF, user->multVec); 177 CHKERRQ(ierr); 178 179 // Setup CEED vectors 180 ierr = user->VecGetArrayRead(user->locVecF, (const PetscScalar **)&f); 181 CHKERRQ(ierr); 182 ierr = user->VecGetArray(user->locVecC, &c); CHKERRQ(ierr); 183 CeedVectorSetArray(user->ceedVecF, user->memType, CEED_USE_POINTER, f); 184 CeedVectorSetArray(user->ceedVecC, user->memType, CEED_USE_POINTER, c); 185 186 // Apply CEED operator 187 CeedOperatorApply(user->opRestrict, user->ceedVecF, user->ceedVecC, 188 CEED_REQUEST_IMMEDIATE); 189 CeedVectorSyncArray(user->ceedVecC, user->memType); 190 191 // Restore PETSc vectors 192 ierr = user->VecRestoreArrayRead(user->locVecF, (const PetscScalar **)&f); 193 CHKERRQ(ierr); 194 ierr = user->VecRestoreArray(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->physSmoother) 215 CeedQFunctionSetContext(user->qf, user->physSmoother, 216 sizeof(*user->physSmoother)); 217 218 // Compute Diagonal via libCEED 219 PetscScalar *x; 220 221 // -- Place PETSc vector in libCEED vector 222 ierr = user->VecGetArray(user->Xloc, &x); CHKERRQ(ierr); 223 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x); 224 225 // -- Compute Diagonal 226 CeedOperatorLinearAssembleDiagonal(user->op, user->Xceed, 227 CEED_REQUEST_IMMEDIATE); 228 CeedVectorSyncArray(user->Xceed, user->memType); 229 230 // -- Reset physics context 231 if (user->physSmoother) 232 CeedQFunctionSetContext(user->qf, user->phys, sizeof(*user->phys)); 233 234 // -- Local-to-Global 235 ierr = user->VecRestoreArray(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 CeedInt length; 252 253 PetscFunctionBeginUser; 254 255 // Global-to-local 256 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 257 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->Xloc); CHKERRQ(ierr); 258 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 259 user->loadIncrement, NULL, NULL, NULL); 260 CHKERRQ(ierr); 261 262 // Setup libCEED input vector 263 ierr = user->VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); 264 CHKERRQ(ierr); 265 CeedVectorSetArray(user->Xceed, user->memType, CEED_USE_POINTER, x); 266 267 // Setup libCEED output vector 268 Vec Eloc; 269 CeedVector eloc; 270 ierr = DMCreateLocalVector(dmEnergy, &Eloc); CHKERRQ(ierr); 271 ierr = VecGetSize(Eloc, &length); CHKERRQ(ierr); 272 ierr = VecDestroy(&Eloc); CHKERRQ(ierr); 273 CeedVectorCreate(user->ceed, length, &eloc); 274 275 // Apply libCEED operator 276 CeedOperatorApply(opEnergy, user->Xceed, eloc, CEED_REQUEST_IMMEDIATE); 277 278 // Restore PETSc vector 279 ierr = user->VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); 280 CHKERRQ(ierr); 281 282 // Reduce max error 283 const CeedScalar *e; 284 CeedVectorGetArrayRead(eloc, CEED_MEM_HOST, &e); 285 (*energy) = 0; 286 for (CeedInt i=0; i<length; i++) 287 (*energy) += e[i]; 288 CeedVectorRestoreArrayRead(eloc, &e); 289 CeedVectorDestroy(&eloc); 290 291 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, 292 user->comm); CHKERRQ(ierr); 293 294 PetscFunctionReturn(0); 295 }; 296