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_LIN = 0, ELAS_HYPER_SS = 1, ELAS_HYPER_FS = 2 47 } problemType; 48 static const char *const problemTypes[] = {"linElas", 49 "hyperSS", 50 "hyperFS", 51 "problemType","ELAS_",0 52 }; 53 static const char *const problemTypesForDisp[] = {"Linear elasticity", 54 "Hyper elasticity small strain", 55 "Hyper elasticity finite strain" 56 }; 57 58 // Forcing function options 59 typedef enum { 60 FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2 61 } forcingType; 62 static const char *const forcingTypes[] = {"none", 63 "constant", 64 "mms", 65 "forcingType","FORCE_",0 66 }; 67 static const char *const forcingTypesForDisp[] = {"None", 68 "Constant", 69 "Manufactured solution" 70 }; 71 72 // Multigrid options 73 typedef enum { 74 MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2 75 } multigridType; 76 static const char *const multigridTypes [] = {"logarithmic", 77 "uniform", 78 "none", 79 "multigridType","MULTIGRID",0 80 }; 81 static const char *const multigridTypesForDisp[] = {"P-multigrid, logarithmic coarsening", 82 "P-multigrind, uniform coarsening", 83 "No multigrid" 84 }; 85 86 typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt, 87 PetscScalar *, void *); 88 // Note: These variables should be updated if additional boundary conditions 89 // are added to boundary.c. 90 BCFunc BCMMS, BCZero, BCClamp; 91 92 // ----------------------------------------------------------------------------- 93 // Structs 94 // ----------------------------------------------------------------------------- 95 // Units 96 typedef struct Units_private *Units; 97 struct Units_private { 98 // Fundamental units 99 PetscScalar meter; 100 PetscScalar kilogram; 101 PetscScalar second; 102 // Derived unit 103 PetscScalar Pascal; 104 }; 105 106 // Application context from user command line options 107 typedef struct AppCtx_private *AppCtx; 108 struct AppCtx_private { 109 char ceedResource[PETSC_MAX_PATH_LEN]; // libCEED backend 110 char meshFile[PETSC_MAX_PATH_LEN]; // exodusII mesh file 111 char outputdir[PETSC_MAX_PATH_LEN]; 112 PetscBool testMode; 113 PetscBool viewSoln; 114 PetscBool viewFinalSoln; 115 problemType problemChoice; 116 forcingType forcingChoice; 117 multigridType multigridChoice; 118 PetscScalar nuSmoother; 119 PetscInt degree; 120 PetscInt qextra; 121 PetscInt numLevels; 122 PetscInt *levelDegrees; 123 PetscInt numIncrements; // Number of steps 124 PetscInt bcClampCount; 125 PetscInt bcClampFaces[16]; 126 // [translation; 3] [rotation axis; 3] [rotation magnitude c_0, c_1] 127 // The rotations are (c_0 + c_1 s) \pi, where s = x · axis 128 PetscScalar bcClampMax[16][8]; 129 PetscInt bcTractionCount; 130 PetscInt bcTractionFaces[16]; 131 PetscScalar bcTractionVector[16][3]; 132 PetscScalar forcingVector[3]; 133 }; 134 135 // Problem specific data 136 // *INDENT-OFF* 137 typedef struct { 138 CeedInt qdatasize; 139 CeedQFunctionUser setupgeo, apply, jacob, energy, diagnostic; 140 const char *setupgeofname, *applyfname, *jacobfname, *energyfname, 141 *diagnosticfname; 142 CeedQuadMode qmode; 143 } problemData; 144 // *INDENT-ON* 145 146 // Data specific to each problem option 147 extern problemData problemOptions[3]; 148 149 // Forcing function data 150 typedef struct { 151 CeedQFunctionUser setupforcing; 152 const char *setupforcingfname; 153 } forcingData; 154 155 extern forcingData forcingOptions[3]; 156 157 // Data for PETSc Matshell 158 typedef struct UserMult_private *UserMult; 159 struct UserMult_private { 160 MPI_Comm comm; 161 DM dm; 162 Vec Xloc, Yloc, NBCs; 163 CeedVector Xceed, Yceed; 164 CeedOperator op; 165 CeedQFunction qf; 166 Ceed ceed; 167 PetscScalar loadIncrement; 168 CeedQFunctionContext ctxPhys, ctxPhysSmoother; 169 }; 170 171 // Data for Jacobian setup routine 172 typedef struct FormJacobCtx_private *FormJacobCtx; 173 struct FormJacobCtx_private { 174 UserMult *jacobCtx; 175 PetscInt numLevels; 176 SNES snesCoarse; 177 Mat *jacobMat, jacobMatCoarse; 178 Vec Ucoarse; 179 }; 180 181 // Data for PETSc Prolongation/Restriction Matshell 182 typedef struct UserMultProlongRestr_private *UserMultProlongRestr; 183 struct UserMultProlongRestr_private { 184 MPI_Comm comm; 185 DM dmC, dmF; 186 Vec locVecC, locVecF; 187 CeedVector ceedVecC, ceedVecF; 188 CeedOperator opProlong, opRestrict; 189 Ceed ceed; 190 }; 191 192 // libCEED data struct for level 193 typedef struct CeedData_private *CeedData; 194 struct CeedData_private { 195 Ceed ceed; 196 CeedBasis basisx, basisu, basisCtoF, basisEnergy, basisDiagnostic; 197 CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi, 198 ErestrictGradui, ErestrictEnergy, ErestrictDiagnostic, 199 ErestrictqdDiagnostici; 200 CeedQFunction qfApply, qfJacob, qfEnergy, qfDiagnostic; 201 CeedOperator opApply, opJacob, opRestrict, opProlong, opEnergy, 202 opDiagnostic; 203 CeedVector qdata, qdataDiagnostic, gradu, xceed, yceed, truesoln; 204 }; 205 206 // Translate PetscMemType to CeedMemType 207 static inline CeedMemType MemTypeP2C(PetscMemType mtype) { 208 return PetscMemTypeDevice(mtype) ? CEED_MEM_DEVICE : CEED_MEM_HOST; 209 } 210 211 // ----------------------------------------------------------------------------- 212 // Process command line options 213 // ----------------------------------------------------------------------------- 214 // Process general command line options 215 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx appCtx); 216 217 // Process physics options 218 PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units); 219 220 // ----------------------------------------------------------------------------- 221 // Setup DM 222 // ----------------------------------------------------------------------------- 223 PetscErrorCode CreateBCLabel(DM dm, const char name[]); 224 225 // Create FE by degree 226 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 227 PetscBool isSimplex, const char prefix[], 228 PetscInt order, PetscFE *fem); 229 230 // Read mesh and distribute DM in parallel 231 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx appCtx, DM *dm); 232 233 // Setup DM with FE space of appropriate degree 234 PetscErrorCode SetupDMByDegree(DM dm, AppCtx appCtx, PetscInt order, 235 PetscBool boundary, PetscInt ncompu); 236 237 // ----------------------------------------------------------------------------- 238 // libCEED Functions 239 // ----------------------------------------------------------------------------- 240 // Destroy libCEED objects 241 PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data); 242 243 // Utility function - essential BC dofs are encoded in closure indices as -(i+1) 244 PetscInt Involute(PetscInt i); 245 246 // Utility function to create local CEED restriction from DMPlex 247 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt P, 248 CeedInt height, DMLabel domainLabel, CeedInt value, 249 CeedElemRestriction *Erestrict); 250 251 // Utility function to get Ceed Restriction for each domain 252 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 253 DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize, 254 CeedElemRestriction *restrictq, CeedElemRestriction *restrictx, 255 CeedElemRestriction *restrictqdi); 256 257 // Set up libCEED for a given degree 258 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dmEnergy, DM dmDiagnostic, 259 Ceed ceed, AppCtx appCtx, 260 CeedQFunctionContext physCtx, 261 CeedData *data, PetscInt fineLevel, 262 PetscInt ncompu, PetscInt Ugsz, 263 PetscInt Ulocsz, CeedVector forceCeed, 264 CeedVector neumannCeed); 265 266 // Set up libCEED multigrid level for a given degree 267 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx appCtx, 268 CeedData *data, PetscInt level, 269 PetscInt ncompu, PetscInt Ugsz, 270 PetscInt Ulocsz, CeedVector fineMult); 271 272 // Setup context data for Jacobian evaluation 273 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V, 274 Vec Vloc, CeedData ceedData, Ceed ceed, 275 CeedQFunctionContext ctxPhys, 276 CeedQFunctionContext ctxPhysSmoother, 277 UserMult jacobianCtx); 278 279 // Setup context data for prolongation and restriction operators 280 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx appCtx, DM dmC, 281 DM dmF, Vec VF, Vec VlocC, Vec VlocF, 282 CeedData ceedDataC, CeedData ceedDataF, 283 Ceed ceed, 284 UserMultProlongRestr prolongRestrCtx); 285 286 // ----------------------------------------------------------------------------- 287 // Jacobian setup 288 // ----------------------------------------------------------------------------- 289 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx); 290 291 // ----------------------------------------------------------------------------- 292 // Solution output 293 // ----------------------------------------------------------------------------- 294 PetscErrorCode ViewSolution(MPI_Comm comm, AppCtx appCtx, Vec U, 295 PetscInt increment, 296 PetscScalar loadIncrement); 297 298 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, 299 UserMult user, AppCtx appCtx, Vec U, 300 CeedElemRestriction ErestrictDiagnostic); 301 302 // ----------------------------------------------------------------------------- 303 // libCEED Operators for MatShell 304 // ----------------------------------------------------------------------------- 305 // This function uses libCEED to compute the local action of an operator 306 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user); 307 308 // This function uses libCEED to compute the non-linear residual 309 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 310 311 // This function uses libCEED to apply the Jacobian for assembly via a SNES 312 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 313 314 // This function uses libCEED to compute the action of the Jacobian 315 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y); 316 317 // This function uses libCEED to compute the action of the prolongation operator 318 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y); 319 320 // This function uses libCEED to compute the action of the restriction operator 321 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y); 322 323 // This function returns the computed diagonal of the operator 324 PetscErrorCode GetDiag_Ceed(Mat A, Vec D); 325 326 // This function calculates the strain energy in the final solution 327 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 328 CeedOperator opEnergy, Vec X, 329 PetscReal *energy); 330 331 // ----------------------------------------------------------------------------- 332 // Boundary Functions 333 // ----------------------------------------------------------------------------- 334 // Note: If additional boundary conditions are added, an update is needed in 335 // elasticity.h for the boundaryOptions variable. 336 337 // BCMMS - boundary function 338 // Values on all points of the mesh is set based on given solution below 339 // for u[0], u[1], u[2] 340 PetscErrorCode BCMMS(PetscInt dim, PetscReal loadIncrement, 341 const PetscReal coords[], PetscInt ncompu, 342 PetscScalar *u, void *ctx); 343 344 // BCClamp - fix boundary values with affine transformation at fraction of load 345 // increment 346 PetscErrorCode BCClamp(PetscInt dim, PetscReal loadIncrement, 347 const PetscReal coords[], PetscInt ncompu, 348 PetscScalar *u, void *ctx); 349 350 #endif //setup_h 351