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), Cinv, 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 forcingTypes[] = {"none", 70 "constant", 71 "mms", 72 "forcingType","FORCE_",0 73 }; 74 static const char *const forcingTypesForDisp[] = {"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 multigridTypes [] = {"logarithmic", 84 "uniform", 85 "none", 86 "multigridType","MULTIGRID",0 87 }; 88 static const char *const multigridTypesForDisp[] = {"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 ceedResource[PETSC_MAX_PATH_LEN]; // libCEED backend 117 char meshFile[PETSC_MAX_PATH_LEN]; // exodusII mesh file 118 char outputdir[PETSC_MAX_PATH_LEN]; 119 PetscBool testMode; 120 PetscBool viewSoln; 121 PetscBool viewFinalSoln; 122 problemType problemChoice; 123 forcingType forcingChoice; 124 multigridType multigridChoice; 125 PetscScalar nuSmoother; 126 PetscInt degree; 127 PetscInt qextra; 128 PetscInt numLevels; 129 PetscInt *levelDegrees; 130 PetscInt numIncrements; // Number of steps 131 PetscInt bcClampCount; 132 PetscInt bcClampFaces[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 bcClampMax[16][8]; 136 PetscInt bcTractionCount; 137 PetscInt bcTractionFaces[16]; 138 PetscScalar bcTractionVector[16][3]; 139 PetscScalar forcingVector[3]; 140 }; 141 142 // Problem specific data 143 // *INDENT-OFF* 144 typedef struct { 145 CeedInt qdatasize; 146 CeedQFunctionUser setupgeo, apply, jacob, energy, diagnostic; 147 const char *setupgeofname, *applyfname, *jacobfname, *energyfname, 148 *diagnosticfname; 149 CeedQuadMode qmode; 150 } problemData; 151 // *INDENT-ON* 152 153 // Data specific to each problem option 154 extern problemData problemOptions[6]; 155 156 // Forcing function data 157 typedef struct { 158 CeedQFunctionUser setupforcing; 159 const char *setupforcingfname; 160 } forcingData; 161 162 extern forcingData forcingOptions[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 Xloc, Yloc, NBCs; 170 CeedVector Xceed, Yceed; 171 CeedOperator op; 172 CeedQFunction qf; 173 Ceed ceed; 174 PetscScalar loadIncrement; 175 CeedQFunctionContext ctxPhys, ctxPhysSmoother; 176 }; 177 178 // Data for Jacobian setup routine 179 typedef struct FormJacobCtx_private *FormJacobCtx; 180 struct FormJacobCtx_private { 181 UserMult *jacobCtx; 182 PetscInt numLevels; 183 SNES snesCoarse; 184 Mat *jacobMat, jacobMatCoarse; 185 Vec Ucoarse; 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 dmC, dmF; 193 Vec locVecC, locVecF; 194 CeedVector ceedVecC, ceedVecF; 195 CeedOperator opProlong, opRestrict; 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 basisx, basisu, basisCtoF, basisEnergy, basisDiagnostic; 204 CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi,ErestrictGradui, 205 ErestrictEnergy, ErestrictDiagnostic, 206 ErestrictdXdx, Erestricttau, ErestrictCc1, 207 ErestrictCinv, ErestrictlamlogJ, ErestrictqdDiagnostici; 208 CeedQFunction qfApply, qfJacob, qfEnergy, qfDiagnostic; 209 CeedOperator opApply, opJacob, opRestrict, opProlong, opEnergy, 210 opDiagnostic; 211 CeedVector qdata, qdataDiagnostic, gradu, xceed, 212 yceed, truesoln, dXdx, tau, Cc1, Cinv, lamlogJ; 213 }; 214 215 // Translate PetscMemType to CeedMemType 216 static inline CeedMemType MemTypeP2C(PetscMemType mtype) { 217 return PetscMemTypeDevice(mtype) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 218 } 219 220 // ----------------------------------------------------------------------------- 221 // Process command line options 222 // ----------------------------------------------------------------------------- 223 // Process general command line options 224 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx appCtx); 225 226 // Process physics options 227 PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units); 228 229 // ----------------------------------------------------------------------------- 230 // Setup DM 231 // ----------------------------------------------------------------------------- 232 PetscErrorCode CreateBCLabel(DM dm, const char name[]); 233 234 // Create FE by degree 235 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 236 PetscBool isSimplex, const char prefix[], 237 PetscInt order, PetscFE *fem); 238 239 // Read mesh and distribute DM in parallel 240 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx appCtx, DM *dm); 241 242 // Setup DM with FE space of appropriate degree 243 PetscErrorCode SetupDMByDegree(DM dm, AppCtx appCtx, PetscInt order, 244 PetscBool boundary, PetscInt ncompu); 245 246 // ----------------------------------------------------------------------------- 247 // libCEED Functions 248 // ----------------------------------------------------------------------------- 249 // Destroy libCEED objects 250 PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data); 251 252 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 253 PetscInt Involute(PetscInt i); 254 255 // Utility function to create local CEED restriction from DMPlex 256 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 257 CeedInt height, DMLabel domainLabel, CeedInt value, 258 CeedElemRestriction *Erestrict); 259 260 // Utility function to get Ceed Restriction for each domain 261 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 262 DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize, 263 CeedElemRestriction *restrictq, CeedElemRestriction *restrictx, 264 CeedElemRestriction *restrictqdi); 265 266 // Set up libCEED for a given degree 267 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dmEnergy, DM dmDiagnostic, 268 Ceed ceed, AppCtx appCtx, 269 CeedQFunctionContext physCtx, 270 CeedData *data, PetscInt fineLevel, 271 PetscInt ncompu, PetscInt Ugsz, 272 PetscInt Ulocsz, CeedVector forceCeed, 273 CeedVector neumannCeed); 274 275 // Set up libCEED multigrid level for a given degree 276 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx appCtx, 277 CeedData *data, PetscInt level, 278 PetscInt ncompu, PetscInt Ugsz, 279 PetscInt Ulocsz, CeedVector fineMult); 280 281 // Setup context data for Jacobian evaluation 282 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V, 283 Vec Vloc, CeedData ceedData, Ceed ceed, 284 CeedQFunctionContext ctxPhys, 285 CeedQFunctionContext ctxPhysSmoother, 286 UserMult jacobianCtx); 287 288 // Setup context data for prolongation and restriction operators 289 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx appCtx, DM dmC, 290 DM dmF, Vec VF, Vec VlocC, Vec VlocF, 291 CeedData ceedDataC, CeedData ceedDataF, 292 Ceed ceed, 293 UserMultProlongRestr prolongRestrCtx); 294 295 // ----------------------------------------------------------------------------- 296 // Jacobian setup 297 // ----------------------------------------------------------------------------- 298 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx); 299 300 // ----------------------------------------------------------------------------- 301 // Solution output 302 // ----------------------------------------------------------------------------- 303 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx appCtx, Vec U, 304 PetscInt increment, 305 PetscScalar loadIncrement); 306 307 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, 308 UserMult user, AppCtx appCtx, Vec U, 309 CeedElemRestriction ErestrictDiagnostic); 310 311 // ----------------------------------------------------------------------------- 312 // libCEED Operators for MatShell 313 // ----------------------------------------------------------------------------- 314 // This function uses libCEED to compute the local action of an operator 315 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user); 316 317 // This function uses libCEED to compute the non-linear residual 318 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 319 320 // This function uses libCEED to apply the Jacobian for assembly via a SNES 321 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 322 323 // This function uses libCEED to compute the action of the Jacobian 324 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y); 325 326 // This function uses libCEED to compute the action of the prolongation operator 327 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y); 328 329 // This function uses libCEED to compute the action of the restriction operator 330 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y); 331 332 // This function returns the computed diagonal of the operator 333 PetscErrorCode GetDiag_Ceed(Mat A, Vec D); 334 335 // This function calculates the strain energy in the final solution 336 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 337 CeedOperator opEnergy, Vec X, 338 PetscReal *energy); 339 340 // ----------------------------------------------------------------------------- 341 // Boundary Functions 342 // ----------------------------------------------------------------------------- 343 // Note: If additional boundary conditions are added, an update is needed in 344 // elasticity.h for the boundaryOptions variable. 345 346 // BCMMS - boundary function 347 // Values on all points of the mesh is set based on given solution below 348 // for u[0], u[1], u[2] 349 PetscErrorCode BCMMS(PetscInt dim, PetscReal loadIncrement, 350 const PetscReal coords[], PetscInt ncompu, 351 PetscScalar *u, void *ctx); 352 353 // BCClamp - fix boundary values with affine transformation at fraction of load 354 // increment 355 PetscErrorCode BCClamp(PetscInt dim, PetscReal loadIncrement, 356 const PetscReal coords[], PetscInt ncompu, 357 PetscScalar *u, void *ctx); 358 359 #endif //setup_h 360