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