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