1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #ifndef setup_h 18 #define setup_h 19 20 #include <ceed.h> 21 #include <petsc.h> 22 #include <petscdmplex.h> 23 #include <petscfe.h> 24 #include <petscksp.h> 25 #include <stdbool.h> 26 #include <string.h> 27 28 #if PETSC_VERSION_LT(3,14,0) 29 # define DMAddBoundary(a,b,c,d,e,f,g,h,i,j,k,l) DMAddBoundary(a,b,c,d,e,f,g,h,j,k,l) 30 #endif 31 32 #ifndef PHYSICS_STRUCT 33 #define PHYSICS_STRUCT 34 typedef struct Physics_private *Physics; 35 struct Physics_private { 36 CeedScalar nu; // Poisson's ratio 37 CeedScalar E; // Young's Modulus 38 }; 39 #endif 40 41 // ----------------------------------------------------------------------------- 42 // Command Line Options 43 // ----------------------------------------------------------------------------- 44 // Problem options 45 typedef enum { 46 ELAS_LINEAR = 0, ELAS_SS_NH = 1, ELAS_FSInitial_NH1 = 2, ELAS_FSInitial_NH2 = 3, 47 ELAS_FSCurrent_NH1 = 4, ELAS_FSCurrent_NH2 = 5 48 } problemType; 49 static const char *const problemTypes[] = {"Linear", 50 "SS-NH", 51 "FSInitial-NH1", 52 "FSInitial-NH2", 53 "FSCurrent-NH1", 54 "FSCurrent-NH2", 55 "problemType","ELAS_",0 56 }; 57 static const char *const problemTypesForDisp[] = {"Linear elasticity", 58 "Hyperelasticity small strain, Neo-Hookean", 59 "Hyperelasticity finite strain Initial config Neo-Hookean w/ dXref_dxinit, Grad(u) storage", 60 "Hyperelasticity finite strain Initial config Neo-Hookean w/ dXref_dxinit, Grad(u), C_inv, constant storage", 61 "Hyperelasticity finite strain Current config Neo-Hookean w/ dXref_dxinit, Grad(u) storage", 62 "Hyperelasticity finite strain Current config Neo-Hookean w/ dXref_dxcurr, tau, constant storage", 63 }; 64 65 // Forcing function options 66 typedef enum { 67 FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2 68 } forcingType; 69 static const char *const forcing_types[] = {"none", 70 "constant", 71 "mms", 72 "forcingType","FORCE_",0 73 }; 74 static const char *const forcing_types_for_disp[] = {"None", 75 "Constant", 76 "Manufactured solution" 77 }; 78 79 // Multigrid options 80 typedef enum { 81 MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2 82 } multigridType; 83 static const char *const multigrid_types [] = {"logarithmic", 84 "uniform", 85 "none", 86 "multigridType","MULTIGRID",0 87 }; 88 static const char *const multigrid_types_for_disp[] = {"P-multigrid, logarithmic coarsening", 89 "P-multigrind, uniform coarsening", 90 "No multigrid" 91 }; 92 93 typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt, 94 PetscScalar *, void *); 95 // Note: These variables should be updated if additional boundary conditions 96 // are added to boundary.c. 97 BCFunc BCMMS, BCZero, BCClamp; 98 99 // ----------------------------------------------------------------------------- 100 // Structs 101 // ----------------------------------------------------------------------------- 102 // Units 103 typedef struct Units_private *Units; 104 struct Units_private { 105 // Fundamental units 106 PetscScalar meter; 107 PetscScalar kilogram; 108 PetscScalar second; 109 // Derived unit 110 PetscScalar Pascal; 111 }; 112 113 // Application context from user command line options 114 typedef struct AppCtx_private *AppCtx; 115 struct AppCtx_private { 116 char ceed_resource[PETSC_MAX_PATH_LEN]; // libCEED backend 117 char mesh_file[PETSC_MAX_PATH_LEN]; // exodusII mesh file 118 char output_dir[PETSC_MAX_PATH_LEN]; 119 PetscBool test_mode; 120 PetscBool view_soln; 121 PetscBool view_final_soln; 122 problemType problem_choice; 123 forcingType forcing_choice; 124 multigridType multigrid_choice; 125 PetscScalar nu_smoother; 126 PetscInt degree; 127 PetscInt q_extra; 128 PetscInt num_levels; 129 PetscInt *level_degrees; 130 PetscInt num_increments; // Number of steps 131 PetscInt bc_clamp_count; 132 PetscInt bc_clamp_faces[16]; 133 // [translation; 3] [rotation axis; 3] [rotation magnitude c_0, c_1] 134 // The rotations are (c_0 + c_1 s) \pi, where s = x · axis 135 PetscScalar bc_clamp_max[16][8]; 136 PetscInt bc_traction_count; 137 PetscInt bc_traction_faces[16]; 138 PetscScalar bc_traction_vector[16][3]; 139 PetscScalar forcing_vector[3]; 140 }; 141 142 // Problem specific data 143 // *INDENT-OFF* 144 typedef struct { 145 CeedInt q_data_size; 146 CeedQFunctionUser setup_geo, apply, jacob, energy, diagnostic; 147 const char *setup_geo_loc, *apply_loc, *jacob_loc, *energy_loc, 148 *diagnostic_loc; 149 CeedQuadMode quad_mode; 150 } problemData; 151 // *INDENT-ON* 152 153 // Data specific to each problem option 154 extern problemData problem_options[6]; 155 156 // Forcing function data 157 typedef struct { 158 CeedQFunctionUser setup_forcing; 159 const char *setup_forcing_loc; 160 } forcingData; 161 162 extern forcingData forcing_options[3]; 163 164 // Data for PETSc Matshell 165 typedef struct UserMult_private *UserMult; 166 struct UserMult_private { 167 MPI_Comm comm; 168 DM dm; 169 Vec X_loc, Y_loc, neumann_bcs; 170 CeedVector x_ceed, y_ceed; 171 CeedOperator op; 172 CeedQFunction qf; 173 Ceed ceed; 174 PetscScalar load_increment; 175 CeedQFunctionContext ctx_phys, ctx_phys_smoother; 176 }; 177 178 // Data for Jacobian setup routine 179 typedef struct FormJacobCtx_private *FormJacobCtx; 180 struct FormJacobCtx_private { 181 UserMult *jacob_ctx; 182 PetscInt num_levels; 183 SNES snes_coarse; 184 Mat *jacob_mat, jacob_mat_coarse; 185 Vec u_coarse; 186 }; 187 188 // Data for PETSc Prolongation/Restriction Matshell 189 typedef struct UserMultProlongRestr_private *UserMultProlongRestr; 190 struct UserMultProlongRestr_private { 191 MPI_Comm comm; 192 DM dm_c, dm_f; 193 Vec loc_vec_c, loc_vec_f; 194 CeedVector ceed_vec_c, ceed_vec_f; 195 CeedOperator op_prolong, op_restrict; 196 Ceed ceed; 197 }; 198 199 // libCEED data struct for level 200 typedef struct CeedData_private *CeedData; 201 struct CeedData_private { 202 Ceed ceed; 203 CeedBasis basis_x, basis_u, basis_c_to_f, basis_energy, 204 basis_diagnostic; 205 CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_qd_i, 206 elem_restr_gradu_i, 207 elem_restr_energy, elem_restr_diagnostic, 208 elem_restr_dXdx, elem_restr_tau, elem_restr_Cc1, 209 elem_restr_C_inv, elem_restr_lam_log_J, elem_restr_qd_diagnostic_i; 210 CeedQFunction qf_apply, qf_jacob, qf_energy, qf_diagnostic; 211 CeedOperator op_apply, op_jacob, op_restrict, op_prolong, op_energy, 212 op_diagnostic; 213 CeedVector q_data, q_data_diagnostic, grad_u, x_ceed, 214 y_ceed, true_soln, dXdx, tau, Cc1, C_inv, lam_log_J; 215 }; 216 217 // Translate PetscMemType to CeedMemType 218 static inline CeedMemType MemTypeP2C(PetscMemType mem_type) { 219 return PetscMemTypeDevice(mem_type) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 220 } 221 222 // ----------------------------------------------------------------------------- 223 // Process command line options 224 // ----------------------------------------------------------------------------- 225 // Process general command line options 226 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx); 227 228 // Process physics options 229 PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units); 230 231 // ----------------------------------------------------------------------------- 232 // Setup DM 233 // ----------------------------------------------------------------------------- 234 PetscErrorCode CreateBCLabel(DM dm, const char name[]); 235 236 // Create FE by degree 237 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 238 PetscBool is_simplex, const char prefix[], 239 PetscInt order, PetscFE *fem); 240 241 // Read mesh and distribute DM in parallel 242 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx app_ctx, DM *dm); 243 244 // Setup DM with FE space of appropriate degree 245 PetscErrorCode SetupDMByDegree(DM dm, AppCtx app_ctx, PetscInt order, 246 PetscBool boundary, PetscInt num_comp_u); 247 248 // ----------------------------------------------------------------------------- 249 // libCEED Functions 250 // ----------------------------------------------------------------------------- 251 // Destroy libCEED objects 252 PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data); 253 254 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 255 PetscInt Involute(PetscInt i); 256 257 // Utility function to create local CEED restriction from DMPlex 258 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 259 CeedInt height, DMLabel domain_label, CeedInt value, 260 CeedElemRestriction *elem_restr); 261 262 // Utility function to get Ceed Restriction for each domain 263 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 264 DMLabel domain_label, PetscInt value, CeedInt P, CeedInt Q, CeedInt q_data_size, 265 CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x, 266 CeedElemRestriction *elem_restr_qd_i); 267 268 // Set up libCEED for a given degree 269 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dm_energy, DM dm_diagnostic, 270 Ceed ceed, AppCtx app_ctx, 271 CeedQFunctionContext phys_ctx, 272 CeedData *data, PetscInt fine_level, 273 PetscInt num_comp_u, PetscInt U_g_size, 274 PetscInt u_loc_size, CeedVector force_ceed, 275 CeedVector neumann_ceed); 276 277 // Set up libCEED multigrid level for a given degree 278 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx app_ctx, 279 CeedData *data, PetscInt level, 280 PetscInt num_comp_u, PetscInt U_g_size, 281 PetscInt u_loc_size, CeedVector fine_mult); 282 283 // Setup context data for Jacobian evaluation 284 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx app_ctx, DM dm, Vec V, 285 Vec V_loc, CeedData ceed_data, Ceed ceed, 286 CeedQFunctionContext ctx_phys, 287 CeedQFunctionContext ctx_phys_smoother, 288 UserMult jacobian_ctx); 289 290 // Setup context data for prolongation and restriction operators 291 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx app_ctx, DM dm_c, 292 DM dm_f, Vec V_f, Vec V_loc_c, Vec V_loc_f, 293 CeedData ceed_data_c, CeedData ceed_data_f, 294 Ceed ceed, 295 UserMultProlongRestr prolong_restr_ctx); 296 297 // ----------------------------------------------------------------------------- 298 // Jacobian setup 299 // ----------------------------------------------------------------------------- 300 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx); 301 302 // ----------------------------------------------------------------------------- 303 // Solution output 304 // ----------------------------------------------------------------------------- 305 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx app_ctx, Vec U, 306 PetscInt increment, 307 PetscScalar load_increment); 308 309 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dm_U, 310 UserMult user, AppCtx app_ctx, Vec U, 311 CeedElemRestriction elem_restr_diagnostic); 312 313 // ----------------------------------------------------------------------------- 314 // libCEED Operators for MatShell 315 // ----------------------------------------------------------------------------- 316 // This function uses libCEED to compute the local action of an operator 317 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user); 318 319 // This function uses libCEED to compute the non-linear residual 320 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 321 322 // This function uses libCEED to apply the Jacobian for assembly via a SNES 323 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 324 325 // This function uses libCEED to compute the action of the Jacobian 326 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y); 327 328 // This function uses libCEED to compute the action of the prolongation operator 329 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y); 330 331 // This function uses libCEED to compute the action of the restriction operator 332 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y); 333 334 // This function returns the computed diagonal of the operator 335 PetscErrorCode GetDiag_Ceed(Mat A, Vec D); 336 337 // This function calculates the strain energy in the final solution 338 PetscErrorCode ComputeStrainEnergy(DM dm_energy, UserMult user, 339 CeedOperator op_energy, Vec X, 340 PetscReal *energy); 341 342 // ----------------------------------------------------------------------------- 343 // Boundary Functions 344 // ----------------------------------------------------------------------------- 345 // Note: If additional boundary conditions are added, an update is needed in 346 // elasticity.h for the boundaryOptions variable. 347 348 // BCMMS - boundary function 349 // Values on all points of the mesh is set based on given solution below 350 // for u[0], u[1], u[2] 351 PetscErrorCode BCMMS(PetscInt dim, PetscReal load_increment, 352 const PetscReal coords[], PetscInt num_comp_u, 353 PetscScalar *u, void *ctx); 354 355 // BCClamp - fix boundary values with affine transformation at fraction of load 356 // increment 357 PetscErrorCode BCClamp(PetscInt dim, PetscReal load_increment, 358 const PetscReal coords[], PetscInt num_comp_u, 359 PetscScalar *u, void *ctx); 360 361 #endif //setup_h 362