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 13 #include "../include/structs.h" 14 #include "../include/utils.h" 15 16 // ----------------------------------------------------------------------------- 17 // libCEED Operators for MatShell 18 // ----------------------------------------------------------------------------- 19 // This function uses libCEED to compute the local action of an operator 20 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user) { 21 PetscScalar *x, *y; 22 PetscMemType x_mem_type, y_mem_type; 23 24 PetscFunctionBeginUser; 25 26 // Global-to-local 27 PetscCall(DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc)); 28 PetscCall(VecZeroEntries(user->Y_loc)); 29 30 // Setup CEED vectors 31 PetscCall(VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, &x_mem_type)); 32 PetscCall(VecGetArrayAndMemType(user->Y_loc, &y, &y_mem_type)); 33 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 34 CeedVectorSetArray(user->y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 35 36 // Apply CEED operator 37 CeedOperatorApply(user->op, user->x_ceed, user->y_ceed, CEED_REQUEST_IMMEDIATE); 38 39 // Restore PETSc vectors 40 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 41 CeedVectorTakeArray(user->y_ceed, MemTypeP2C(y_mem_type), NULL); 42 PetscCall(VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x)); 43 PetscCall(VecRestoreArrayAndMemType(user->Y_loc, &y)); 44 45 // Local-to-global 46 PetscCall(VecZeroEntries(Y)); 47 PetscCall(DMLocalToGlobal(user->dm, user->Y_loc, ADD_VALUES, Y)); 48 49 PetscFunctionReturn(0); 50 }; 51 52 // This function uses libCEED to compute the non-linear residual 53 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 54 UserMult user = (UserMult)ctx; 55 56 PetscFunctionBeginUser; 57 58 // Use computed BCs 59 PetscCall(VecZeroEntries(user->X_loc)); 60 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, user->load_increment, NULL, NULL, NULL)); 61 62 // libCEED for local action of residual evaluator 63 PetscCall(ApplyLocalCeedOp(X, Y, user)); 64 65 // Neumann BCs 66 if (user->neumann_bcs) { 67 PetscCall(VecAXPY(Y, -user->load_increment, user->neumann_bcs)); 68 } 69 70 PetscFunctionReturn(0); 71 }; 72 73 // This function uses libCEED to apply the Jacobian for assembly via a SNES 74 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx) { 75 UserMult user = (UserMult)ctx; 76 77 PetscFunctionBeginUser; 78 79 // Zero boundary values 80 PetscCall(VecZeroEntries(user->X_loc)); 81 82 // libCEED for local action of residual evaluator 83 PetscCall(ApplyLocalCeedOp(X, Y, user)); 84 85 PetscFunctionReturn(0); 86 }; 87 88 // This function uses libCEED to compute the action of the Jacobian 89 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y) { 90 UserMult user; 91 92 PetscFunctionBeginUser; 93 94 // Zero boundary values 95 PetscCall(MatShellGetContext(A, &user)); 96 PetscCall(VecZeroEntries(user->X_loc)); 97 98 // libCEED for local action of Jacobian 99 PetscCall(ApplyLocalCeedOp(X, Y, user)); 100 101 PetscFunctionReturn(0); 102 }; 103 104 // This function uses libCEED to compute the action of the prolongation operator 105 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y) { 106 UserMultProlongRestr user; 107 PetscScalar *c, *f; 108 PetscMemType c_mem_type, f_mem_type; 109 110 PetscFunctionBeginUser; 111 112 PetscCall(MatShellGetContext(A, &user)); 113 114 // Global-to-local 115 PetscCall(VecZeroEntries(user->loc_vec_c)); 116 PetscCall(DMGlobalToLocal(user->dm_c, X, INSERT_VALUES, user->loc_vec_c)); 117 PetscCall(VecZeroEntries(user->loc_vec_f)); 118 119 // Setup CEED vectors 120 PetscCall(VecGetArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c, &c_mem_type)); 121 PetscCall(VecGetArrayAndMemType(user->loc_vec_f, &f, &f_mem_type)); 122 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, c); 123 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, f); 124 125 // Apply CEED operator 126 CeedOperatorApply(user->op_prolong, user->ceed_vec_c, user->ceed_vec_f, CEED_REQUEST_IMMEDIATE); 127 128 // Restore PETSc vectors 129 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 130 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 131 PetscCall(VecRestoreArrayReadAndMemType(user->loc_vec_c, (const PetscScalar **)&c)); 132 PetscCall(VecRestoreArrayAndMemType(user->loc_vec_f, &f)); 133 134 // Local-to-global 135 PetscCall(VecZeroEntries(Y)); 136 PetscCall(DMLocalToGlobal(user->dm_f, user->loc_vec_f, ADD_VALUES, Y)); 137 138 PetscFunctionReturn(0); 139 } 140 141 // This function uses libCEED to compute the action of the restriction operator 142 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y) { 143 UserMultProlongRestr user; 144 PetscScalar *c, *f; 145 PetscMemType c_mem_type, f_mem_type; 146 147 PetscFunctionBeginUser; 148 149 PetscCall(MatShellGetContext(A, &user)); 150 151 // Global-to-local 152 PetscCall(VecZeroEntries(user->loc_vec_f)); 153 PetscCall(DMGlobalToLocal(user->dm_f, X, INSERT_VALUES, user->loc_vec_f)); 154 PetscCall(VecZeroEntries(user->loc_vec_c)); 155 156 // Setup CEED vectors 157 PetscCall(VecGetArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f, &f_mem_type)); 158 PetscCall(VecGetArrayAndMemType(user->loc_vec_c, &c, &c_mem_type)); 159 CeedVectorSetArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), CEED_USE_POINTER, f); 160 CeedVectorSetArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), CEED_USE_POINTER, c); 161 162 // Apply CEED operator 163 CeedOperatorApply(user->op_restrict, user->ceed_vec_f, user->ceed_vec_c, CEED_REQUEST_IMMEDIATE); 164 165 // Restore PETSc vectors 166 CeedVectorTakeArray(user->ceed_vec_f, MemTypeP2C(f_mem_type), NULL); 167 CeedVectorTakeArray(user->ceed_vec_c, MemTypeP2C(c_mem_type), NULL); 168 PetscCall(VecRestoreArrayReadAndMemType(user->loc_vec_f, (const PetscScalar **)&f)); 169 PetscCall(VecRestoreArrayAndMemType(user->loc_vec_c, &c)); 170 171 // Local-to-global 172 PetscCall(VecZeroEntries(Y)); 173 PetscCall(DMLocalToGlobal(user->dm_c, user->loc_vec_c, ADD_VALUES, Y)); 174 175 PetscFunctionReturn(0); 176 }; 177 178 // This function returns the computed diagonal of the operator 179 PetscErrorCode GetDiag_Ceed(Mat A, Vec D) { 180 UserMult user; 181 182 PetscFunctionBeginUser; 183 184 PetscCall(MatShellGetContext(A, &user)); 185 186 // -- Set physics context 187 if (user->ctx_phys_smoother) CeedQFunctionSetContext(user->qf, user->ctx_phys_smoother); 188 189 // Compute Diagonal via libCEED 190 PetscScalar *x; 191 PetscMemType x_mem_type; 192 193 // -- Place PETSc vector in libCEED vector 194 PetscCall(VecGetArrayAndMemType(user->X_loc, &x, &x_mem_type)); 195 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 196 197 // -- Compute Diagonal 198 CeedOperatorLinearAssembleDiagonal(user->op, user->x_ceed, CEED_REQUEST_IMMEDIATE); 199 200 // -- Reset physics context 201 if (user->ctx_phys_smoother) CeedQFunctionSetContext(user->qf, user->ctx_phys); 202 203 // -- Local-to-Global 204 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 205 PetscCall(VecRestoreArrayAndMemType(user->X_loc, &x)); 206 PetscCall(VecZeroEntries(D)); 207 PetscCall(DMLocalToGlobal(user->dm, user->X_loc, ADD_VALUES, D)); 208 209 // Cleanup 210 PetscCall(VecZeroEntries(user->X_loc)); 211 212 PetscFunctionReturn(0); 213 }; 214 215 // This function calculates the strain energy in the final solution 216 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, CeedOperator op_energy, Vec X, PetscReal *energy) { 217 PetscScalar *x; 218 PetscMemType x_mem_type; 219 CeedInt length; 220 221 PetscFunctionBeginUser; 222 223 // Global-to-local 224 PetscCall(VecZeroEntries(user->X_loc)); 225 PetscCall(DMGlobalToLocal(user->dm, X, INSERT_VALUES, user->X_loc)); 226 PetscCall(DMPlexInsertBoundaryValues(user->dm, PETSC_TRUE, user->X_loc, user->load_increment, NULL, NULL, NULL)); 227 228 // Setup libCEED input vector 229 PetscCall(VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, &x_mem_type)); 230 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 231 232 // Setup libCEED output vector 233 Vec E_loc; 234 CeedVector e_loc; 235 PetscCall(DMCreateLocalVector(dmEnergy, &E_loc)); 236 PetscCall(VecGetSize(E_loc, &length)); 237 PetscCall(VecDestroy(&E_loc)); 238 CeedVectorCreate(user->ceed, length, &e_loc); 239 240 // Apply libCEED operator 241 CeedOperatorApply(op_energy, user->x_ceed, e_loc, CEED_REQUEST_IMMEDIATE); 242 243 // Restore PETSc vector 244 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 245 PetscCall(VecRestoreArrayRead(user->X_loc, (const PetscScalar **)&x)); 246 247 // Reduce max error 248 const CeedScalar *e; 249 CeedVectorGetArrayRead(e_loc, CEED_MEM_HOST, &e); 250 (*energy) = 0; 251 for (CeedInt i = 0; i < length; i++) (*energy) += e[i]; 252 CeedVectorRestoreArrayRead(e_loc, &e); 253 CeedVectorDestroy(&e_loc); 254 255 PetscCall(MPI_Allreduce(MPI_IN_PLACE, energy, 1, MPIU_REAL, MPIU_SUM, user->comm)); 256 257 PetscFunctionReturn(0); 258 }; 259