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