1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Matrix shell operations for solid mechanics example using PETSc 10 11 #include "../include/matops.h" 12 #include "../include/structs.h" 13 #include "../include/utils.h" 14 15 // ----------------------------------------------------------------------------- 16 // libCEED Operators for MatShell 17 // ----------------------------------------------------------------------------- 18 // This function uses libCEED to compute the local action of an operator 19 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user) { 20 PetscErrorCode ierr; 21 PetscScalar *x, *y; 22 PetscMemType x_mem_type, y_mem_type; 23 24 PetscFunctionBeginUser; 25 26 // Global-to-local 27 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc); CHKERRQ(ierr); 28 ierr = VecZeroEntries(user->Y_loc); CHKERRQ(ierr); 29 30 // Setup CEED vectors 31 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 32 &x_mem_type); CHKERRQ(ierr); 33 ierr = VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type); CHKERRQ(ierr); 34 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 35 CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 36 37 // Apply CEED operator 38 CeedOperatorApply(user->op, user->x_ceed, user->y_ceed, CEED_REQUEST_IMMEDIATE); 39 40 // Restore PETSc vectors 41 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 42 CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL); 43 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 44 CHKERRQ(ierr); 45 ierr = VecRestoreArrayAndMemType(user->Y_loc, &y); CHKERRQ(ierr); 46 47 // Local-to-global 48 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 49 ierr = DMLocalToGlobal(user->dm, user->Y_loc, ADD_VALUES, Y); CHKERRQ(ierr); 50 51 PetscFunctionReturn(0); 52 }; 53 54 // This function uses libCEED to compute the non-linear residual 55 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 56 PetscErrorCode ierr; 57 UserMult user = (UserMult)ctx; 58 59 PetscFunctionBeginUser; 60 61 // Use computed BCs 62 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 63 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, 64 user->load_increment, NULL, NULL, NULL); 65 CHKERRQ(ierr); 66 67 // libCEED for local action of residual evaluator 68 ierr = ApplyLocalCeedOp(X, Y, user); CHKERRQ(ierr); 69 70 // Neumann BCs 71 if (user->neumann_bcs) { 72 ierr = VecAXPY(Y, -user->load_increment, user->neumann_bcs); CHKERRQ(ierr); 73 } 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->X_loc); 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->X_loc); 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 PetscMemType c_mem_type, f_mem_type; 117 118 PetscFunctionBeginUser; 119 120 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 121 122 // Global-to-local 123 ierr = VecZeroEntries(user->loc_vec_c); CHKERRQ(ierr); 124 ierr = DMGlobalToLocal(user->dm_c, X, INSERT_VALUES, user->loc_vec_c); 125 CHKERRQ(ierr); 126 ierr = VecZeroEntries(user->loc_vec_f); CHKERRQ(ierr); 127 128 // Setup CEED vectors 129 ierr = VecGetArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c, 130 &c_mem_type); CHKERRQ(ierr); 131 ierr = VecGetArrayAndMemType(user->loc_vec_f, &f, &f_mem_type); CHKERRQ(ierr); 132 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, 133 c); 134 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, 135 f); 136 137 // Apply CEED operator 138 CeedOperatorApply(user->op_prolong, user->ceed_vec_c, user->ceed_vec_f, 139 CEED_REQUEST_IMMEDIATE); 140 141 // Restore PETSc vectors 142 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 143 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 144 ierr = VecRestoreArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c); 145 CHKERRQ(ierr); 146 ierr = VecRestoreArrayAndMemType(user->loc_vec_f, &f); CHKERRQ(ierr); 147 148 // Local-to-global 149 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 150 ierr = DMLocalToGlobal(user->dm_f, user->loc_vec_f, ADD_VALUES, Y); 151 CHKERRQ(ierr); 152 153 PetscFunctionReturn(0); 154 } 155 156 // This function uses libCEED to compute the action of the restriction operator 157 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) { 158 PetscErrorCode ierr; 159 UserMultProlongRestr user; 160 PetscScalar *c, *f; 161 PetscMemType c_mem_type, f_mem_type; 162 163 PetscFunctionBeginUser; 164 165 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 166 167 // Global-to-local 168 ierr = VecZeroEntries(user->loc_vec_f); CHKERRQ(ierr); 169 ierr = DMGlobalToLocal(user->dm_f, X, INSERT_VALUES, user->loc_vec_f); 170 CHKERRQ(ierr); 171 ierr = VecZeroEntries(user->loc_vec_c); CHKERRQ(ierr); 172 173 // Setup CEED vectors 174 ierr = VecGetArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f, 175 &f_mem_type); CHKERRQ(ierr); 176 ierr = VecGetArrayAndMemType(user->loc_vec_c, &c, &c_mem_type); CHKERRQ(ierr); 177 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, 178 f); 179 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, 180 c); 181 182 // Apply CEED operator 183 CeedOperatorApply(user->op_restrict, user->ceed_vec_f, user->ceed_vec_c, 184 CEED_REQUEST_IMMEDIATE); 185 186 // Restore PETSc vectors 187 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 188 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 189 ierr = VecRestoreArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f); 190 CHKERRQ(ierr); 191 ierr = VecRestoreArrayAndMemType(user->loc_vec_c, &c); CHKERRQ(ierr); 192 193 // Local-to-global 194 ierr = VecZeroEntries(Y); CHKERRQ(ierr); 195 ierr = DMLocalToGlobal(user->dm_c, user->loc_vec_c, ADD_VALUES, Y); 196 CHKERRQ(ierr); 197 198 PetscFunctionReturn(0); 199 }; 200 201 // This function returns the computed diagonal of the operator 202 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) { 203 PetscErrorCode ierr; 204 UserMult user; 205 206 PetscFunctionBeginUser; 207 208 ierr = MatShellGetContext(A, &user); CHKERRQ(ierr); 209 210 // -- Set physics context 211 if (user->ctx_phys_smoother) 212 CeedQFunctionSetContext(user->qf, user->ctx_phys_smoother); 213 214 // Compute Diagonal via libCEED 215 PetscScalar *x; 216 PetscMemType x_mem_type; 217 218 // -- Place PETSc vector in libCEED vector 219 ierr = VecGetArrayAndMemType(user->X_loc, &x, &x_mem_type); CHKERRQ(ierr); 220 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 221 222 // -- Compute Diagonal 223 CeedOperatorLinearAssembleDiagonal(user->op, user->x_ceed, 224 CEED_REQUEST_IMMEDIATE); 225 226 // -- Reset physics context 227 if (user->ctx_phys_smoother) 228 CeedQFunctionSetContext(user->qf, user->ctx_phys); 229 230 // -- Local-to-Global 231 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 232 ierr = VecRestoreArrayAndMemType(user->X_loc, &x); CHKERRQ(ierr); 233 ierr = VecZeroEntries(D); CHKERRQ(ierr); 234 ierr = DMLocalToGlobal(user->dm, user->X_loc, ADD_VALUES, D); CHKERRQ(ierr); 235 236 // Cleanup 237 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 238 239 PetscFunctionReturn(0); 240 }; 241 242 // This function calculates the strain energy in the final solution 243 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 244 CeedOperator op_energy, Vec X, 245 PetscReal *energy) { 246 PetscErrorCode ierr; 247 PetscScalar *x; 248 PetscMemType x_mem_type; 249 CeedInt length; 250 251 PetscFunctionBeginUser; 252 253 // Global-to-local 254 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 255 ierr = DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc); CHKERRQ(ierr); 256 ierr = DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, 257 user->load_increment, NULL, NULL, NULL); 258 CHKERRQ(ierr); 259 260 // Setup libCEED input vector 261 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 262 &x_mem_type); CHKERRQ(ierr); 263 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 264 265 // Setup libCEED output vector 266 Vec E_loc; 267 CeedVector e_loc; 268 ierr = DMCreateLocalVector(dmEnergy, &E_loc); CHKERRQ(ierr); 269 ierr = VecGetSize(E_loc, &length); CHKERRQ(ierr); 270 ierr = VecDestroy(&E_loc); CHKERRQ(ierr); 271 CeedVectorCreate(user->ceed, length, &e_loc); 272 273 // Apply libCEED operator 274 CeedOperatorApply(op_energy, user->x_ceed, e_loc, CEED_REQUEST_IMMEDIATE); 275 276 // Restore PETSc vector 277 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 278 ierr = VecRestoreArrayRead(user->X_loc, (const PetscScalar **)&x); 279 CHKERRQ(ierr); 280 281 // Reduce max error 282 const CeedScalar *e; 283 CeedVectorGetArrayRead(e_loc, CEED_MEM_HOST, &e); 284 (*energy) = 0; 285 for (CeedInt i=0; i<length; i++) 286 (*energy) += e[i]; 287 CeedVectorRestoreArrayRead(e_loc, &e); 288 CeedVectorDestroy(&e_loc); 289 290 ierr = MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, 291 user->comm); CHKERRQ(ierr); 292 293 PetscFunctionReturn(0); 294 }; 295