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