1 // Copyright (c) 2017-2026, 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
ProcessCommandLineOptions(MPI_Comm comm,AppCtx app_ctx)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_FS_NH || app_ctx->problem_choice == ELAS_FS_MR) && app_ctx->forcing_choice == FORCE_CONST) {
61 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
62 "Cannot use constant forcing and finite strain formulation. "
63 "Constant forcing in reference frame currently unavailable.");
64 }
65
66 // Dirichlet boundary conditions
67 app_ctx->bc_clamp_count = 16;
68 PetscCall(PetscOptionsIntArray("-bc_clamp", "Face IDs to apply incremental Dirichlet BC", NULL, app_ctx->bc_clamp_faces, &app_ctx->bc_clamp_count,
69 NULL));
70 // Set vector for each clamped BC
71 for (PetscInt i = 0; i < app_ctx->bc_clamp_count; i++) {
72 // Translation vector
73 char option_name[25];
74 const size_t nclamp_params = sizeof(app_ctx->bc_clamp_max[0]) / sizeof(app_ctx->bc_clamp_max[0][0]);
75 for (PetscInt j = 0; j < nclamp_params; j++) app_ctx->bc_clamp_max[i][j] = 0.;
76
77 snprintf(option_name, sizeof option_name, "-bc_clamp_%" PetscInt_FMT "_translate", app_ctx->bc_clamp_faces[i]);
78 max_n = 3;
79 PetscCall(PetscOptionsScalarArray(option_name, "Vector to translate clamped end by", NULL, app_ctx->bc_clamp_max[i], &max_n, NULL));
80
81 // Rotation vector
82 max_n = 5;
83 snprintf(option_name, sizeof option_name, "-bc_clamp_%" PetscInt_FMT "_rotate", app_ctx->bc_clamp_faces[i]);
84 PetscCall(PetscOptionsScalarArray(option_name, "Vector with axis of rotation and rotation, in radians", NULL, &app_ctx->bc_clamp_max[i][3],
85 &max_n, NULL));
86
87 // Normalize
88 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] +
89 app_ctx->bc_clamp_max[i][5] * app_ctx->bc_clamp_max[i][5]);
90 if (fabs(norm) < 1e-16) norm = 1;
91 for (PetscInt j = 0; j < 3; j++) app_ctx->bc_clamp_max[i][3 + j] /= norm;
92 }
93
94 // Neumann boundary conditions
95 app_ctx->bc_traction_count = 16;
96 PetscCall(PetscOptionsIntArray("-bc_traction", "Face IDs to apply traction (Neumann) BC", NULL, app_ctx->bc_traction_faces,
97 &app_ctx->bc_traction_count, NULL));
98 // Set vector for each traction BC
99 for (PetscInt i = 0; i < app_ctx->bc_traction_count; i++) {
100 // Translation vector
101 char option_name[25];
102 for (PetscInt j = 0; j < 3; j++) app_ctx->bc_traction_vector[i][j] = 0.;
103
104 snprintf(option_name, sizeof option_name, "-bc_traction_%" PetscInt_FMT, app_ctx->bc_traction_faces[i]);
105 max_n = 3;
106 PetscBool set = false;
107 PetscCall(PetscOptionsScalarArray(option_name, "Traction vector for constrained face", NULL, app_ctx->bc_traction_vector[i], &max_n, &set));
108
109 if (!set) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Traction vector must be set for all traction boundary conditions.");
110 }
111
112 app_ctx->multigrid_choice = MULTIGRID_LOGARITHMIC;
113 PetscCall(PetscOptionsEnum("-multigrid", "Set multigrid type option", NULL, multigrid_types, (PetscEnum)app_ctx->multigrid_choice,
114 (PetscEnum *)&app_ctx->multigrid_choice, NULL));
115
116 app_ctx->test_mode = PETSC_FALSE;
117 PetscCall(PetscOptionsBool("-test", "Testing mode (do not print unless error is large)", NULL, app_ctx->test_mode, &(app_ctx->test_mode), NULL));
118
119 app_ctx->expect_final_strain = -1.;
120 PetscCall(PetscOptionsReal("-expect_final_strain_energy", "Expect final strain energy close to this value.", NULL, app_ctx->expect_final_strain,
121 &app_ctx->expect_final_strain, NULL));
122
123 app_ctx->test_tol = 1e-8;
124 PetscCall(PetscOptionsReal("-expect_final_state_rtol", "Relative tolerance for final strain energy test", NULL, app_ctx->test_tol,
125 &app_ctx->test_tol, NULL));
126
127 app_ctx->view_soln = PETSC_FALSE;
128 PetscCall(PetscOptionsBool("-view_soln", "Write out solution vector for viewing", NULL, app_ctx->view_soln, &(app_ctx->view_soln), NULL));
129
130 app_ctx->view_final_soln = PETSC_FALSE;
131 PetscCall(PetscOptionsBool("-view_final_soln", "Write out final solution vector for viewing", NULL, app_ctx->view_final_soln,
132 &(app_ctx->view_final_soln), NULL));
133
134 PetscBool set;
135 char energy_viewer_filename[PETSC_MAX_PATH_LEN] = "";
136 PetscCall(PetscOptionsString("-strain_energy_monitor", "Print out current strain energy at every load increment", NULL, energy_viewer_filename,
137 energy_viewer_filename, sizeof(energy_viewer_filename), &set));
138 if (set) {
139 PetscCall(PetscViewerASCIIOpen(comm, energy_viewer_filename, &app_ctx->energy_viewer));
140 PetscCall(PetscViewerASCIIPrintf(app_ctx->energy_viewer, "increment,energy\n"));
141 // Initial configuration is base energy state; this may not be true if we extend in the future to initially loaded configurations (because a truly
142 // at-rest initial state may not be realizable).
143 PetscCall(PetscViewerASCIIPrintf(app_ctx->energy_viewer, "%f,%e\n", 0., 0.));
144 }
145 PetscOptionsEnd(); // End of setting AppCtx
146
147 // Check for all required values set
148 if (app_ctx->test_mode) {
149 if (app_ctx->forcing_choice == FORCE_NONE && !app_ctx->bc_clamp_count) app_ctx->forcing_choice = FORCE_MMS;
150 }
151 if (!app_ctx->bc_clamp_count && app_ctx->forcing_choice != FORCE_MMS) {
152 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "-boundary options needed");
153 }
154
155 // Provide default ceed resource if not specified
156 if (!ceed_flag) {
157 const char *ceed_resource = "/cpu/self";
158 strncpy(app_ctx->ceed_resource, ceed_resource, 10);
159 }
160
161 // Determine number of levels
162 switch (app_ctx->multigrid_choice) {
163 case MULTIGRID_LOGARITHMIC:
164 app_ctx->num_levels = ceil(log(app_ctx->degree) / log(2)) + 1;
165 break;
166 case MULTIGRID_UNIFORM:
167 app_ctx->num_levels = app_ctx->degree;
168 break;
169 case MULTIGRID_NONE:
170 app_ctx->num_levels = 1;
171 break;
172 }
173
174 // Populate array of degrees for each level for multigrid
175 PetscCall(PetscMalloc1(app_ctx->num_levels, &(app_ctx->level_degrees)));
176
177 switch (app_ctx->multigrid_choice) {
178 case MULTIGRID_LOGARITHMIC:
179 for (int i = 0; i < app_ctx->num_levels - 1; i++) app_ctx->level_degrees[i] = pow(2, i);
180 app_ctx->level_degrees[app_ctx->num_levels - 1] = app_ctx->degree;
181 break;
182 case MULTIGRID_UNIFORM:
183 for (int i = 0; i < app_ctx->num_levels; i++) app_ctx->level_degrees[i] = i + 1;
184 break;
185 case MULTIGRID_NONE:
186 app_ctx->level_degrees[0] = app_ctx->degree;
187 break;
188 }
189
190 PetscFunctionReturn(PETSC_SUCCESS);
191 };
192