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 /// Helper functions for solid mechanics example using PETSc 10 11 #include "../include/misc.h" 12 #include "../include/utils.h" 13 14 // ----------------------------------------------------------------------------- 15 // Create libCEED operator context 16 // ----------------------------------------------------------------------------- 17 // Setup context data for Jacobian evaluation 18 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx app_ctx, DM dm, Vec V, 19 Vec V_loc, CeedData ceed_data, Ceed ceed, 20 CeedQFunctionContext ctx_phys, 21 CeedQFunctionContext ctx_phys_smoother, 22 UserMult jacobian_ctx) { 23 PetscErrorCode ierr; 24 25 PetscFunctionBeginUser; 26 27 // PETSc objects 28 jacobian_ctx->comm = comm; 29 jacobian_ctx->dm = dm; 30 31 // Work vectors 32 jacobian_ctx->X_loc = V_loc; 33 ierr = VecDuplicate(V_loc, &jacobian_ctx->Y_loc); CHKERRQ(ierr); 34 jacobian_ctx->x_ceed = ceed_data->x_ceed; 35 jacobian_ctx->y_ceed = ceed_data->y_ceed; 36 37 // libCEED operator 38 jacobian_ctx->op = ceed_data->op_jacobian; 39 jacobian_ctx->qf = ceed_data->qf_jacobian; 40 41 // Ceed 42 jacobian_ctx->ceed = ceed; 43 44 // Physics 45 jacobian_ctx->ctx_phys = ctx_phys; 46 jacobian_ctx->ctx_phys_smoother = ctx_phys_smoother; 47 PetscFunctionReturn(0); 48 }; 49 50 // Setup context data for prolongation and restriction operators 51 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx app_ctx, DM dm_c, 52 DM dm_f, Vec V_f, Vec V_loc_c, Vec V_loc_f, 53 CeedData ceed_data_c, CeedData ceed_data_f, 54 Ceed ceed, UserMultProlongRestr prolong_restr_ctx) { 55 PetscFunctionBeginUser; 56 57 // PETSc objects 58 prolong_restr_ctx->comm = comm; 59 prolong_restr_ctx->dm_c = dm_c; 60 prolong_restr_ctx->dm_f = dm_f; 61 62 // Work vectors 63 prolong_restr_ctx->loc_vec_c = V_loc_c; 64 prolong_restr_ctx->loc_vec_f = V_loc_f; 65 prolong_restr_ctx->ceed_vec_c = ceed_data_c->x_ceed; 66 prolong_restr_ctx->ceed_vec_f = ceed_data_f->x_ceed; 67 68 // libCEED operators 69 prolong_restr_ctx->op_prolong = ceed_data_f->op_prolong; 70 prolong_restr_ctx->op_restrict = ceed_data_f->op_restrict; 71 72 // Ceed 73 prolong_restr_ctx->ceed = ceed; 74 PetscFunctionReturn(0); 75 }; 76 77 // ----------------------------------------------------------------------------- 78 // Jacobian setup 79 // ----------------------------------------------------------------------------- 80 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx) { 81 PetscErrorCode ierr; 82 83 PetscFunctionBeginUser; 84 85 // Context data 86 FormJacobCtx form_jacob_ctx = (FormJacobCtx)ctx; 87 PetscInt num_levels = form_jacob_ctx->num_levels; 88 Mat *jacob_mat = form_jacob_ctx->jacob_mat; 89 90 // Update Jacobian on each level 91 for (PetscInt level = 0; level < num_levels; level++) { 92 ierr = MatAssemblyBegin(jacob_mat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 93 ierr = MatAssemblyEnd(jacob_mat[level], MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 94 } 95 96 // Form coarse assembled matrix 97 CeedOperatorLinearAssemble(form_jacob_ctx->op_coarse, 98 form_jacob_ctx->coo_values); 99 const CeedScalar *values; 100 CeedVectorGetArrayRead(form_jacob_ctx->coo_values, CEED_MEM_HOST, &values); 101 ierr = MatSetValuesCOO(form_jacob_ctx->jacob_mat_coarse, values, ADD_VALUES); 102 CHKERRQ(ierr); 103 CeedVectorRestoreArrayRead(form_jacob_ctx->coo_values, &values); 104 105 // J_pre might be AIJ (e.g., when using coloring), so we need to assemble it 106 ierr = MatAssemblyBegin(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 107 ierr = MatAssemblyEnd(J_pre, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 108 if (J != J_pre) { 109 ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 110 ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 111 } 112 PetscFunctionReturn(0); 113 }; 114 115 // ----------------------------------------------------------------------------- 116 // Output solution for visualization 117 // ----------------------------------------------------------------------------- 118 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U, 119 PetscInt increment, PetscScalar load_increment) { 120 PetscErrorCode ierr; 121 DM dm; 122 PetscViewer viewer; 123 char output_filename[PETSC_MAX_PATH_LEN]; 124 PetscMPIInt rank; 125 126 PetscFunctionBeginUser; 127 128 // Create output directory 129 MPI_Comm_rank(comm, &rank); 130 if (!rank) {ierr = PetscMkdir(app_ctx->output_dir); CHKERRQ(ierr);} 131 132 // Build file name 133 ierr = PetscSNPrintf(output_filename, sizeof output_filename, 134 "%s/solution-%03" PetscInt_FMT ".vtu", app_ctx->output_dir, 135 increment); CHKERRQ(ierr); 136 137 // Increment sequence 138 ierr = VecGetDM(U, &dm); CHKERRQ(ierr); 139 ierr = DMSetOutputSequenceNumber(dm, increment, load_increment); CHKERRQ(ierr); 140 141 // Output solution vector 142 ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer); 143 CHKERRQ(ierr); 144 ierr = VecView(U, viewer); CHKERRQ(ierr); 145 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 146 147 PetscFunctionReturn(0); 148 }; 149 150 // ----------------------------------------------------------------------------- 151 // Output diagnostic quantities for visualization 152 // ----------------------------------------------------------------------------- 153 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, 154 UserMult user, AppCtx app_ctx, Vec U, 155 CeedElemRestriction elem_restr_diagnostic) { 156 PetscErrorCode ierr; 157 Vec Diagnostic, Y_loc, mult_vec; 158 CeedVector y_ceed; 159 CeedScalar *x, *y; 160 PetscMemType x_mem_type, y_mem_type; 161 PetscInt loc_size; 162 PetscViewer viewer; 163 char output_filename[PETSC_MAX_PATH_LEN]; 164 165 PetscFunctionBeginUser; 166 167 // --------------------------------------------------------------------------- 168 // PETSc and libCEED vectors 169 // --------------------------------------------------------------------------- 170 ierr = DMCreateGlobalVector(user->dm, &Diagnostic); CHKERRQ(ierr); 171 ierr = PetscObjectSetName((PetscObject)Diagnostic, ""); CHKERRQ(ierr); 172 ierr = DMCreateLocalVector(user->dm, &Y_loc); CHKERRQ(ierr); 173 ierr = VecGetSize(Y_loc, &loc_size); CHKERRQ(ierr); 174 CeedVectorCreate(user->ceed, loc_size, &y_ceed); 175 176 // --------------------------------------------------------------------------- 177 // Compute quantities 178 // --------------------------------------------------------------------------- 179 // -- Global-to-local 180 ierr = VecZeroEntries(user->X_loc); CHKERRQ(ierr); 181 ierr = DMPlexInsertBoundaryValues(dmU, PETSC_TRUE, user->X_loc, 182 user->load_increment, NULL, NULL, NULL); 183 CHKERRQ(ierr); 184 ierr = DMGlobalToLocal(dmU, U, INSERT_VALUES, user->X_loc); CHKERRQ(ierr); 185 ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr); 186 187 // -- Setup CEED vectors 188 ierr = VecGetArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x, 189 &x_mem_type); CHKERRQ(ierr); 190 ierr = VecGetArrayAndMemType(Y_loc, &y, &y_mem_type); CHKERRQ(ierr); 191 CeedVectorSetArray(user->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER, x); 192 CeedVectorSetArray(y_ceed, MemTypeP2C(y_mem_type), CEED_USE_POINTER, y); 193 194 // -- Apply CEED operator 195 CeedOperatorApply(user->op, user->x_ceed, y_ceed, CEED_REQUEST_IMMEDIATE); 196 197 // -- Restore PETSc vector; keep y_ceed viewing memory of Y_loc for use below 198 CeedVectorTakeArray(user->x_ceed, MemTypeP2C(x_mem_type), NULL); 199 ierr = VecRestoreArrayReadAndMemType(user->X_loc, (const PetscScalar **)&x); 200 CHKERRQ(ierr); 201 202 // -- Local-to-global 203 ierr = VecZeroEntries(Diagnostic); CHKERRQ(ierr); 204 ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, Diagnostic); 205 CHKERRQ(ierr); 206 207 // --------------------------------------------------------------------------- 208 // Scale for multiplicity 209 // --------------------------------------------------------------------------- 210 // -- Setup vectors 211 ierr = VecDuplicate(Diagnostic, &mult_vec); CHKERRQ(ierr); 212 ierr = VecZeroEntries(Y_loc); CHKERRQ(ierr); 213 214 // -- Compute multiplicity 215 CeedElemRestrictionGetMultiplicity(elem_restr_diagnostic, y_ceed); 216 217 // -- Restore vectors 218 CeedVectorTakeArray(y_ceed, MemTypeP2C(y_mem_type), NULL); 219 ierr = VecRestoreArrayAndMemType(Y_loc, &y); CHKERRQ(ierr); 220 221 // -- Local-to-global 222 ierr = VecZeroEntries(mult_vec); CHKERRQ(ierr); 223 ierr = DMLocalToGlobal(user->dm, Y_loc, ADD_VALUES, mult_vec); 224 CHKERRQ(ierr); 225 226 // -- Scale 227 ierr = VecReciprocal(mult_vec); CHKERRQ(ierr); 228 ierr = VecPointwiseMult(Diagnostic, Diagnostic, mult_vec); 229 230 // --------------------------------------------------------------------------- 231 // Output solution vector 232 // --------------------------------------------------------------------------- 233 ierr = PetscSNPrintf(output_filename, sizeof output_filename, 234 "%s/diagnostic_quantities.vtu", 235 app_ctx->output_dir); CHKERRQ(ierr); 236 ierr = PetscViewerVTKOpen(comm, output_filename, FILE_MODE_WRITE, &viewer); 237 CHKERRQ(ierr); 238 ierr = VecView(Diagnostic, viewer); CHKERRQ(ierr); 239 ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr); 240 241 // --------------------------------------------------------------------------- 242 // Cleanup 243 // --------------------------------------------------------------------------- 244 ierr = VecDestroy(&Diagnostic); CHKERRQ(ierr); 245 ierr = VecDestroy(&mult_vec); CHKERRQ(ierr); 246 ierr = VecDestroy(&Y_loc); CHKERRQ(ierr); 247 CeedVectorDestroy(&y_ceed); 248 249 PetscFunctionReturn(0); 250 }; 251 252 // ----------------------------------------------------------------------------- 253 // Regression testing 254 // ----------------------------------------------------------------------------- 255 // test option change. could remove the loading step. Run only with one loading step and compare relatively to ref file 256 // option: expect_final_strain_energy and check against the relative error to ref is within tolerance (10^-5) I.e. one Newton solve then check final energy 257 PetscErrorCode RegressionTests_solids(AppCtx app_ctx, PetscReal energy) { 258 PetscFunctionBegin; 259 260 if (app_ctx->expect_final_strain >= 0.) { 261 PetscReal energy_ref = app_ctx->expect_final_strain; 262 PetscReal error = PetscAbsReal(energy - energy_ref) / energy_ref; 263 264 if (error > app_ctx->test_tol) { 265 PetscErrorCode ierr; 266 ierr = PetscPrintf(PETSC_COMM_WORLD, 267 "Energy %e does not match expected energy %e: relative tolerance %e > %e\n", 268 (double)energy, (double)energy_ref, (double)error, app_ctx->test_tol); 269 CHKERRQ(ierr); 270 } 271 } 272 PetscFunctionReturn(0); 273 }; 274