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 /// Command line option processing for solid mechanics example using PETSc 10 11 #include "../include/cl-options.h" 12 13 #include <ceed.h> 14 #include <petscsys.h> 15 16 // ----------------------------------------------------------------------------- 17 // Process command line options 18 // ----------------------------------------------------------------------------- 19 // Process general command line options 20 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx) { 21 PetscBool ceed_flag = PETSC_FALSE; 22 23 PetscFunctionBeginUser; 24 25 PetscOptionsBegin(comm, NULL, "Elasticity / Hyperelasticity in PETSc with libCEED", NULL); 26 27 PetscCall(PetscOptionsString("-ceed", "CEED resource specifier", NULL, app_ctx->ceed_resource, app_ctx->ceed_resource, 28 sizeof(app_ctx->ceed_resource), &ceed_flag)); 29 30 PetscCall(PetscStrncpy(app_ctx->output_dir, ".", 2)); // Default - current directory 31 PetscCall(PetscOptionsString("-output_dir", "Output directory", NULL, app_ctx->output_dir, app_ctx->output_dir, sizeof(app_ctx->output_dir), NULL)); 32 33 app_ctx->degree = 3; 34 PetscCall(PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", NULL, app_ctx->degree, &app_ctx->degree, NULL)); 35 36 app_ctx->q_extra = 0; 37 PetscCall(PetscOptionsInt("-q_extra", "Number of extra quadrature points", NULL, app_ctx->q_extra, &app_ctx->q_extra, NULL)); 38 39 PetscCall(PetscOptionsString("-mesh", "Read mesh from file", NULL, app_ctx->mesh_file, app_ctx->mesh_file, sizeof(app_ctx->mesh_file), NULL)); 40 41 app_ctx->problem_choice = ELAS_LINEAR; // Default - Linear Elasticity 42 PetscCall(PetscOptionsEnum("-problem", "Solves Elasticity & Hyperelasticity Problems", NULL, problemTypes, (PetscEnum)app_ctx->problem_choice, 43 (PetscEnum *)&app_ctx->problem_choice, NULL)); 44 app_ctx->name = problemTypes[app_ctx->problem_choice]; 45 app_ctx->name_for_disp = problemTypesForDisp[app_ctx->problem_choice]; 46 47 app_ctx->num_increments = app_ctx->problem_choice == ELAS_LINEAR ? 1 : 10; 48 PetscCall(PetscOptionsInt("-num_steps", "Number of pseudo-time steps", NULL, app_ctx->num_increments, &app_ctx->num_increments, NULL)); 49 50 app_ctx->forcing_choice = FORCE_NONE; // Default - no forcing term 51 PetscCall(PetscOptionsEnum("-forcing", "Set forcing function option", NULL, forcing_types, (PetscEnum)app_ctx->forcing_choice, 52 (PetscEnum *)&app_ctx->forcing_choice, NULL)); 53 54 PetscInt max_n = 3; 55 app_ctx->forcing_vector[0] = 0; 56 app_ctx->forcing_vector[1] = -1; 57 app_ctx->forcing_vector[2] = 0; 58 PetscCall(PetscOptionsScalarArray("-forcing_vec", "Direction to apply constant force", NULL, app_ctx->forcing_vector, &max_n, NULL)); 59 60 if ((app_ctx->problem_choice == ELAS_FSInitial_NH1 || app_ctx->problem_choice == ELAS_FSInitial_NH2 || 61 app_ctx->problem_choice == ELAS_FSCurrent_NH1 || app_ctx->problem_choice == ELAS_FSCurrent_NH2 || 62 app_ctx->problem_choice == ELAS_FSInitial_MR1) && 63 app_ctx->forcing_choice == FORCE_CONST) { 64 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 65 "Cannot use constant forcing and finite strain formulation. " 66 "Constant forcing in reference frame currently unavailable."); 67 } 68 69 // Dirichlet boundary conditions 70 app_ctx->bc_clamp_count = 16; 71 PetscCall( 72 PetscOptionsIntArray("-bc_clamp", "Face IDs to apply incremental Dirichlet BC", NULL, app_ctx->bc_clamp_faces, &app_ctx->bc_clamp_count, NULL)); 73 // Set vector for each clamped BC 74 for (PetscInt i = 0; i < app_ctx->bc_clamp_count; i++) { 75 // Translation vector 76 char option_name[25]; 77 const size_t nclamp_params = sizeof(app_ctx->bc_clamp_max[0]) / sizeof(app_ctx->bc_clamp_max[0][0]); 78 for (PetscInt j = 0; j < nclamp_params; j++) app_ctx->bc_clamp_max[i][j] = 0.; 79 80 snprintf(option_name, sizeof option_name, "-bc_clamp_%" PetscInt_FMT "_translate", app_ctx->bc_clamp_faces[i]); 81 max_n = 3; 82 PetscCall(PetscOptionsScalarArray(option_name, "Vector to translate clamped end by", NULL, app_ctx->bc_clamp_max[i], &max_n, NULL)); 83 84 // Rotation vector 85 max_n = 5; 86 snprintf(option_name, sizeof option_name, "-bc_clamp_%" PetscInt_FMT "_rotate", app_ctx->bc_clamp_faces[i]); 87 PetscCall(PetscOptionsScalarArray(option_name, "Vector with axis of rotation and rotation, in radians", NULL, &app_ctx->bc_clamp_max[i][3], 88 &max_n, NULL)); 89 90 // Normalize 91 PetscScalar norm = sqrt(app_ctx->bc_clamp_max[i][3] * app_ctx->bc_clamp_max[i][3] + app_ctx->bc_clamp_max[i][4] * app_ctx->bc_clamp_max[i][4] + 92 app_ctx->bc_clamp_max[i][5] * app_ctx->bc_clamp_max[i][5]); 93 if (fabs(norm) < 1e-16) norm = 1; 94 for (PetscInt j = 0; j < 3; j++) app_ctx->bc_clamp_max[i][3 + j] /= norm; 95 } 96 97 // Neumann boundary conditions 98 app_ctx->bc_traction_count = 16; 99 PetscCall(PetscOptionsIntArray("-bc_traction", "Face IDs to apply traction (Neumann) BC", NULL, app_ctx->bc_traction_faces, 100 &app_ctx->bc_traction_count, NULL)); 101 // Set vector for each traction BC 102 for (PetscInt i = 0; i < app_ctx->bc_traction_count; i++) { 103 // Translation vector 104 char option_name[25]; 105 for (PetscInt j = 0; j < 3; j++) app_ctx->bc_traction_vector[i][j] = 0.; 106 107 snprintf(option_name, sizeof option_name, "-bc_traction_%" PetscInt_FMT, app_ctx->bc_traction_faces[i]); 108 max_n = 3; 109 PetscBool set = false; 110 PetscCall(PetscOptionsScalarArray(option_name, "Traction vector for constrained face", NULL, app_ctx->bc_traction_vector[i], &max_n, &set)); 111 112 if (!set) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Traction vector must be set for all traction boundary conditions."); 113 } 114 115 app_ctx->multigrid_choice = MULTIGRID_LOGARITHMIC; 116 PetscCall(PetscOptionsEnum("-multigrid", "Set multigrid type option", NULL, multigrid_types, (PetscEnum)app_ctx->multigrid_choice, 117 (PetscEnum *)&app_ctx->multigrid_choice, NULL)); 118 119 app_ctx->test_mode = PETSC_FALSE; 120 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, app_ctx->test_mode, &(app_ctx->test_mode), NULL)); 121 122 app_ctx->expect_final_strain = -1.; 123 PetscCall(PetscOptionsReal("-expect_final_strain_energy", "Expect final strain energy close to this value.", NULL, app_ctx->expect_final_strain, 124 &app_ctx->expect_final_strain, NULL)); 125 126 app_ctx->test_tol = 1e-8; 127 PetscCall(PetscOptionsReal("-expect_final_state_rtol", "Relative tolerance for final strain energy test", NULL, app_ctx->test_tol, 128 &app_ctx->test_tol, NULL)); 129 130 app_ctx->view_soln = PETSC_FALSE; 131 PetscCall(PetscOptionsBool("-view_soln", "Write out solution vector for viewing", NULL, app_ctx->view_soln, &(app_ctx->view_soln), NULL)); 132 133 app_ctx->view_final_soln = PETSC_FALSE; 134 PetscCall(PetscOptionsBool("-view_final_soln", "Write out final solution vector for viewing", NULL, app_ctx->view_final_soln, 135 &(app_ctx->view_final_soln), NULL)); 136 137 PetscBool set; 138 char energy_viewer_filename[PETSC_MAX_PATH_LEN] = ""; 139 PetscCall(PetscOptionsString("-strain_energy_monitor", "Print out current strain energy at every load increment", NULL, energy_viewer_filename, 140 energy_viewer_filename, sizeof(energy_viewer_filename), &set)); 141 if (set) { 142 PetscCall(PetscViewerASCIIOpen(comm, energy_viewer_filename, &app_ctx->energy_viewer)); 143 PetscCall(PetscViewerASCIIPrintf(app_ctx->energy_viewer, "increment,energy\n")); 144 // Initial configuration is base energy state; this may not be true if we extend in the future to initially loaded configurations (because a truly 145 // at-rest initial state may not be realizable). 146 PetscCall(PetscViewerASCIIPrintf(app_ctx->energy_viewer, "%f,%e\n", 0., 0.)); 147 } 148 PetscOptionsEnd(); // End of setting AppCtx 149 150 // Check for all required values set 151 if (app_ctx->test_mode) { 152 if (app_ctx->forcing_choice == FORCE_NONE && !app_ctx->bc_clamp_count) app_ctx->forcing_choice = FORCE_MMS; 153 } 154 if (!app_ctx->bc_clamp_count && app_ctx->forcing_choice != FORCE_MMS) { 155 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "-boundary options needed"); 156 } 157 158 // Provide default ceed resource if not specified 159 if (!ceed_flag) { 160 const char *ceed_resource = "/cpu/self"; 161 strncpy(app_ctx->ceed_resource, ceed_resource, 10); 162 } 163 164 // Determine number of levels 165 switch (app_ctx->multigrid_choice) { 166 case MULTIGRID_LOGARITHMIC: 167 app_ctx->num_levels = ceil(log(app_ctx->degree) / log(2)) + 1; 168 break; 169 case MULTIGRID_UNIFORM: 170 app_ctx->num_levels = app_ctx->degree; 171 break; 172 case MULTIGRID_NONE: 173 app_ctx->num_levels = 1; 174 break; 175 } 176 177 // Populate array of degrees for each level for multigrid 178 PetscCall(PetscMalloc1(app_ctx->num_levels, &(app_ctx->level_degrees))); 179 180 switch (app_ctx->multigrid_choice) { 181 case MULTIGRID_LOGARITHMIC: 182 for (int i = 0; i < app_ctx->num_levels - 1; i++) app_ctx->level_degrees[i] = pow(2, i); 183 app_ctx->level_degrees[app_ctx->num_levels - 1] = app_ctx->degree; 184 break; 185 case MULTIGRID_UNIFORM: 186 for (int i = 0; i < app_ctx->num_levels; i++) app_ctx->level_degrees[i] = i + 1; 187 break; 188 case MULTIGRID_NONE: 189 app_ctx->level_degrees[0] = app_ctx->degree; 190 break; 191 } 192 193 PetscFunctionReturn(PETSC_SUCCESS); 194 }; 195