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