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 = VecGetArrayRead(user->Xloc, (const PetscScalar **)&x); CHKERRQ(ierr); 38 ierr = VecGetArray(user->Yloc, &y); CHKERRQ(ierr); 39 CeedVectorSetArray(user->Xceed, CEED_MEM_HOST, CEED_USE_POINTER, x); 40 CeedVectorSetArray(user->Yceed, CEED_MEM_HOST, CEED_USE_POINTER, y); 41 42 // Apply CEED operator 43 CeedOperatorApply(user->op, user->Xceed, user->Yceed, CEED_REQUEST_IMMEDIATE); 44 CeedVectorSyncArray(user->Yceed, CEED_MEM_HOST); 45 46 // Restore PETSc vectors 47 ierr = VecRestoreArrayRead(user->Xloc, (const PetscScalar **)&x); 48 CHKERRQ(ierr); 49 ierr = VecRestoreArray(user->Yloc, &y); CHKERRQ(ierr); 50 51 // Local-to-global 52 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 53 ierr = DMLocalToGlobal(user->dm, user->Yloc, ADD_VALUES, Y); CHKERRQ(ierr); 54 55 PetscFunctionReturn(0); 56 }; 57 58 // This function uses libCEED to compute the non-linear residual 59 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 60 PetscErrorCode ierr; 61 UserMult user = (UserMult)ctx; 62 63 PetscFunctionBeginUser; 64 65 // Use computed BCs 66 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 67 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->Xloc, 68 user->loadIncrement, NULL, NULL, NULL); 69 CHKERRQ(ierr); 70 71 // libCEED for local action of residual evaluator 72 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 73 74 PetscFunctionReturn(0); 75 }; 76 77 // This function uses libCEED to apply the Jacobian for assembly via a SNES 78 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 79 PetscErrorCode ierr; 80 UserMult user = (UserMult)ctx; 81 82 PetscFunctionBeginUser; 83 84 // Use computed BCs 85 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 86 87 // libCEED for local action of residual evaluator 88 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 89 90 PetscFunctionReturn(0); 91 }; 92 93 // This function uses libCEED to compute the action of the Jacobian 94 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y) { 95 PetscErrorCode ierr; 96 UserMult user; 97 98 PetscFunctionBeginUser; 99 100 // Zero boundary values 101 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 102 ierr = VecZeroEntries(user->Xloc); CHKERRQ(ierr); 103 104 // libCEED for local action of Jacobian 105 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 106 107 PetscFunctionReturn(0); 108 }; 109 110 // This function uses libCEED to compute the action of the prolongation operator 111 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) { 112 PetscErrorCode ierr; 113 UserMultProlongRestr user; 114 PetscScalar *c, *f; 115 116 PetscFunctionBeginUser; 117 118 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 119 120 // Global-to-local 121 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 122 ierr = DMGlobalToLocal(user->dmC, X, INSERT_VALUES, user->locVecC); 123 CHKERRQ(ierr); 124 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr); 125 126 // Setup CEED vectors 127 ierr = VecGetArrayRead(user->locVecC, (const PetscScalar **)&c); 128 CHKERRQ(ierr); 129 ierr = VecGetArray(user->locVecF, &f); CHKERRQ(ierr); 130 CeedVectorSetArray(user->ceedVecC, CEED_MEM_HOST, CEED_USE_POINTER, c); 131 CeedVectorSetArray(user->ceedVecF, CEED_MEM_HOST, CEED_USE_POINTER, f); 132 133 // Apply CEED operator 134 CeedOperatorApply(user->opProlong, user->ceedVecC, user->ceedVecF, 135 CEED_REQUEST_IMMEDIATE); 136 CeedVectorSyncArray(user->ceedVecF, CEED_MEM_HOST); 137 138 // Restore PETSc vectors 139 ierr = VecRestoreArrayRead(user->locVecC, (const PetscScalar **)c); 140 CHKERRQ(ierr); 141 ierr = VecRestoreArray(user->locVecF, &f); CHKERRQ(ierr); 142 143 // Multiplicity 144 ierr = VecPointwiseMult(user->locVecF, user->locVecF, user->multVec); 145 146 // Local-to-global 147 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 148 ierr = DMLocalToGlobal(user->dmF, user->locVecF, ADD_VALUES, Y); 149 CHKERRQ(ierr); 150 151 PetscFunctionReturn(0); 152 } 153 154 // This function uses libCEED to compute the action of the restriction operator 155 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) { 156 PetscErrorCode ierr; 157 UserMultProlongRestr user; 158 PetscScalar *c, *f; 159 160 PetscFunctionBeginUser; 161 162 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 163 164 // Global-to-local 165 ierr = VecZeroEntries(user->locVecF); CHKERRQ(ierr); 166 ierr = DMGlobalToLocal(user->dmF, X, INSERT_VALUES, user->locVecF); 167 CHKERRQ(ierr); 168 ierr = VecZeroEntries(user->locVecC); CHKERRQ(ierr); 169 170 // Multiplicity 171 ierr = VecPointwiseMult(user->locVecF, user->locVecF, user->multVec); 172 CHKERRQ(ierr); 173 174 // Setup CEED vectors 175 ierr = VecGetArrayRead(user->locVecF, (const PetscScalar **)&f); CHKERRQ(ierr); 176 ierr = VecGetArray(user->locVecC, &c); CHKERRQ(ierr); 177 CeedVectorSetArray(user->ceedVecF, CEED_MEM_HOST, CEED_USE_POINTER, f); 178 CeedVectorSetArray(user->ceedVecC, CEED_MEM_HOST, CEED_USE_POINTER, c); 179 180 // Apply CEED operator 181 CeedOperatorApply(user->opRestrict, user->ceedVecF, user->ceedVecC, 182 CEED_REQUEST_IMMEDIATE); 183 CeedVectorSyncArray(user->ceedVecC, CEED_MEM_HOST); 184 185 // Restore PETSc vectors 186 ierr = VecRestoreArrayRead(user->locVecF, (const PetscScalar **)&f); 187 CHKERRQ(ierr); 188 ierr = VecRestoreArray(user->locVecC, &c); CHKERRQ(ierr); 189 190 // Local-to-global 191 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 192 ierr = DMLocalToGlobal(user->dmC, user->locVecC, ADD_VALUES, Y); 193 CHKERRQ(ierr); 194 195 PetscFunctionReturn(0); 196 }; 197 // This function returns the computed diagonal of the operator 198 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) { 199 PetscErrorCode ierr; 200 UserMult user; 201 202 PetscFunctionBeginUser; 203 204 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 205 206 // Compute Diagonal via libCEED 207 CeedVector ceedDiagVec; 208 const CeedScalar *diagArray; 209 210 // -- Compute Diagonal 211 CeedOperatorAssembleLinearDiagonal(user->op, &ceedDiagVec, 212 CEED_REQUEST_IMMEDIATE); 213 214 // -- Place in PETSc vector 215 CeedVectorGetArrayRead(ceedDiagVec, CEED_MEM_HOST, &diagArray); 216 ierr = VecPlaceArray(user->Xloc, diagArray); CHKERRQ(ierr); 217 218 // -- Local-to-Global 219 ierr = VecZeroEntries(D); CHKERRQ(ierr); 220 ierr = DMLocalToGlobal(user->dm, user->Xloc, ADD_VALUES, D); CHKERRQ(ierr); 221 222 // -- Cleanup 223 ierr = VecResetArray(user->Xloc); CHKERRQ(ierr); 224 CeedVectorRestoreArrayRead(ceedDiagVec, &diagArray); 225 CeedVectorDestroy(&ceedDiagVec); 226 227 PetscFunctionReturn(0); 228 }; 229