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