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