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