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 x_mem_type, y_mem_type; 30 31 PetscFunctionBeginUser; 32 33 // Global-to-local 34 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc); CHKERRQ(ierr); 35 ierr = VecZeroEntries(user->Y_loc); CHKERRQ(ierr); 36 37 // Setup CEED vectors 38 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 39 &x_mem_type); CHKERRQ(ierr); 40 ierr = VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type); CHKERRQ(ierr); 41 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 42 CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 43 44 // Apply CEED operator 45 CeedOperatorApply(user->op, user->x_ceed, user->y_ceed, CEED_REQUEST_IMMEDIATE); 46 47 // Restore PETSc vectors 48 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 49 CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL); 50 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 51 CHKERRQ(ierr); 52 ierr = VecRestoreArrayAndMemType(user->Y_loc, &y); CHKERRQ(ierr); 53 54 // Local-to-global 55 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 56 ierr = DMLocalToGlobal(user->dm, user->Y_loc, 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->X_loc); CHKERRQ(ierr); 70 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, 71 user->load_increment, 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->neumann_bcs) { 79 ierr = VecAXPY(Y, -user->load_increment, user->neumann_bcs); 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->X_loc); 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->X_loc); 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 c_mem_type, f_mem_type; 124 125 PetscFunctionBeginUser; 126 127 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 128 129 // Global-to-local 130 ierr = VecZeroEntries(user->loc_vec_c); CHKERRQ(ierr); 131 ierr = DMGlobalToLocal(user->dm_c, X, INSERT_VALUES, user->loc_vec_c); 132 CHKERRQ(ierr); 133 ierr = VecZeroEntries(user->loc_vec_f); CHKERRQ(ierr); 134 135 // Setup CEED vectors 136 ierr = VecGetArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c, 137 &c_mem_type); CHKERRQ(ierr); 138 ierr = VecGetArrayAndMemType(user->loc_vec_f, &f, &f_mem_type); CHKERRQ(ierr); 139 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, 140 c); 141 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, 142 f); 143 144 // Apply CEED operator 145 CeedOperatorApply(user->op_prolong, user->ceed_vec_c, user->ceed_vec_f, 146 CEED_REQUEST_IMMEDIATE); 147 148 // Restore PETSc vectors 149 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 150 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 151 ierr = VecRestoreArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c); 152 CHKERRQ(ierr); 153 ierr = VecRestoreArrayAndMemType(user->loc_vec_f, &f); CHKERRQ(ierr); 154 155 // Local-to-global 156 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 157 ierr = DMLocalToGlobal(user->dm_f, user->loc_vec_f, ADD_VALUES, Y); 158 CHKERRQ(ierr); 159 160 PetscFunctionReturn(0); 161 } 162 163 // This function uses libCEED to compute the action of the restriction operator 164 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) { 165 PetscErrorCode ierr; 166 UserMultProlongRestr user; 167 PetscScalar *c, *f; 168 PetscMemType c_mem_type, f_mem_type; 169 170 PetscFunctionBeginUser; 171 172 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 173 174 // Global-to-local 175 ierr = VecZeroEntries(user->loc_vec_f); CHKERRQ(ierr); 176 ierr = DMGlobalToLocal(user->dm_f, X, INSERT_VALUES, user->loc_vec_f); 177 CHKERRQ(ierr); 178 ierr = VecZeroEntries(user->loc_vec_c); CHKERRQ(ierr); 179 180 // Setup CEED vectors 181 ierr = VecGetArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f, 182 &f_mem_type); CHKERRQ(ierr); 183 ierr = VecGetArrayAndMemType(user->loc_vec_c, &c, &c_mem_type); CHKERRQ(ierr); 184 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, 185 f); 186 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, 187 c); 188 189 // Apply CEED operator 190 CeedOperatorApply(user->op_restrict, user->ceed_vec_f, user->ceed_vec_c, 191 CEED_REQUEST_IMMEDIATE); 192 193 // Restore PETSc vectors 194 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 195 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 196 ierr = VecRestoreArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f); 197 CHKERRQ(ierr); 198 ierr = VecRestoreArrayAndMemType(user->loc_vec_c, &c); CHKERRQ(ierr); 199 200 // Local-to-global 201 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 202 ierr = DMLocalToGlobal(user->dm_c, user->loc_vec_c, ADD_VALUES, Y); 203 CHKERRQ(ierr); 204 205 PetscFunctionReturn(0); 206 }; 207 208 // This function returns the computed diagonal of the operator 209 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) { 210 PetscErrorCode ierr; 211 UserMult user; 212 213 PetscFunctionBeginUser; 214 215 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 216 217 // -- Set physics context 218 if (user->ctx_phys_smoother) 219 CeedQFunctionSetContext(user->qf, user->ctx_phys_smoother); 220 221 // Compute Diagonal via libCEED 222 PetscScalar *x; 223 PetscMemType x_mem_type; 224 225 // -- Place PETSc vector in libCEED vector 226 ierr = VecGetArrayAndMemType(user->X_loc, &x, &x_mem_type); CHKERRQ(ierr); 227 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 228 229 // -- Compute Diagonal 230 CeedOperatorLinearAssembleDiagonal(user->op, user->x_ceed, 231 CEED_REQUEST_IMMEDIATE); 232 233 // -- Reset physics context 234 if (user->ctx_phys_smoother) 235 CeedQFunctionSetContext(user->qf, user->ctx_phys); 236 237 // -- Local-to-Global 238 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 239 ierr = VecRestoreArrayAndMemType(user->X_loc, &x); CHKERRQ(ierr); 240 ierr = VecZeroEntries(D); CHKERRQ(ierr); 241 ierr = DMLocalToGlobal(user->dm, user->X_loc, ADD_VALUES, D); CHKERRQ(ierr); 242 243 // Cleanup 244 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 245 246 PetscFunctionReturn(0); 247 }; 248 249 // This function calculates the strain energy in the final solution 250 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 251 CeedOperator op_energy, Vec X, 252 PetscReal *energy) { 253 PetscErrorCode ierr; 254 PetscScalar *x; 255 PetscMemType x_mem_type; 256 CeedInt length; 257 258 PetscFunctionBeginUser; 259 260 // Global-to-local 261 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 262 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc); CHKERRQ(ierr); 263 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, 264 user->load_increment, NULL, NULL, NULL); 265 CHKERRQ(ierr); 266 267 // Setup libCEED input vector 268 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 269 &x_mem_type); CHKERRQ(ierr); 270 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 271 272 // Setup libCEED output vector 273 Vec E_loc; 274 CeedVector e_loc; 275 ierr = DMCreateLocalVector(dmEnergy, &E_loc); CHKERRQ(ierr); 276 ierr = VecGetSize(E_loc, &length); CHKERRQ(ierr); 277 ierr = VecDestroy(&E_loc); CHKERRQ(ierr); 278 CeedVectorCreate(user->ceed, length, &e_loc); 279 280 // Apply libCEED operator 281 CeedOperatorApply(op_energy, user->x_ceed, e_loc, CEED_REQUEST_IMMEDIATE); 282 283 // Restore PETSc vector 284 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 285 ierr = VecRestoreArrayRead(user->X_loc, (const PetscScalar **)&x); 286 CHKERRQ(ierr); 287 288 // Reduce max error 289 const CeedScalar *e; 290 CeedVectorGetArrayRead(e_loc, CEED_MEM_HOST, &e); 291 (*energy) = 0; 292 for (CeedInt i=0; i<length; i++) 293 (*energy) += e[i]; 294 CeedVectorRestoreArrayRead(e_loc, &e); 295 CeedVectorDestroy(&e_loc); 296 297 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, 298 user->comm); CHKERRQ(ierr); 299 300 PetscFunctionReturn(0); 301 }; 302