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 <stdbool.h> 21 #include <string.h> 22 23 #include <petsc.h> 24 #include <petscdmplex.h> 25 #include <petscksp.h> 26 #include <petscfe.h> 27 28 #include <ceed.h> 29 30 #if PETSC_VERSION_LT(3,14,0) 31 # 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) 32 #endif 33 34 #ifndef PHYSICS_STRUCT 35 #define PHYSICS_STRUCT 36 typedef struct Physics_private *Physics; 37 struct Physics_private { 38 CeedScalar nu; // Poisson's ratio 39 CeedScalar E; // Young's Modulus 40 }; 41 #endif 42 43 // ----------------------------------------------------------------------------- 44 // Command Line Options 45 // ----------------------------------------------------------------------------- 46 // Problem options 47 typedef enum { 48 ELAS_LIN = 0, ELAS_HYPER_SS = 1, ELAS_HYPER_FS = 2 49 } problemType; 50 static const char *const problemTypes[] = {"linElas", 51 "hyperSS", 52 "hyperFS", 53 "problemType","ELAS_",0 54 }; 55 static const char *const problemTypesForDisp[] = {"Linear elasticity", 56 "Hyper elasticity small strain", 57 "Hyper elasticity finite strain" 58 }; 59 60 // Forcing function options 61 typedef enum { 62 FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2 63 } forcingType; 64 static const char *const forcingTypes[] = {"none", 65 "constant", 66 "mms", 67 "forcingType","FORCE_",0 68 }; 69 static const char *const forcingTypesForDisp[] = {"None", 70 "Constant", 71 "Manufactured solution" 72 }; 73 74 // Multigrid options 75 typedef enum { 76 MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2 77 } multigridType; 78 static const char *const multigridTypes [] = {"logarithmic", 79 "uniform", 80 "none", 81 "multigridType","MULTIGRID",0 82 }; 83 static const char *const multigridTypesForDisp[] = {"P-multigrid, logarithmic coarsening", 84 "P-multigrind, uniform coarsening", 85 "No multigrid" 86 }; 87 88 typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt, 89 PetscScalar *, void *); 90 // Note: These variables should be updated if additional boundary conditions 91 // are added to boundary.c. 92 BCFunc BCMMS, BCZero, BCClamp; 93 94 // MemType Options 95 static const char *const memTypes[] = {"host","device","memType", 96 "CEED_MEM_",0 97 }; 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 PetscBool testMode; 119 PetscBool viewSoln; 120 PetscBool viewFinalSoln; 121 problemType problemChoice; 122 forcingType forcingChoice; 123 multigridType multigridChoice; 124 PetscScalar nuSmoother; 125 PetscInt degree; 126 PetscInt qextra; 127 PetscInt numLevels; 128 PetscInt *levelDegrees; 129 PetscInt numIncrements; // Number of steps 130 PetscInt bcClampCount; 131 PetscInt bcClampFaces[16]; 132 PetscScalar bcClampMax[16][7]; 133 PetscInt bcTractionCount; 134 PetscInt bcTractionFaces[16]; 135 PetscScalar bcTractionVector[16][3]; 136 PetscScalar forcingVector[3]; 137 PetscBool petscHaveCuda, setMemTypeRequest; 138 CeedMemType memTypeRequested; 139 }; 140 141 // Problem specific data 142 // *INDENT-OFF* 143 typedef struct { 144 CeedInt qdatasize; 145 CeedQFunctionUser setupgeo, apply, jacob, energy, diagnostic; 146 const char *setupgeofname, *applyfname, *jacobfname, *energyfname, 147 *diagnosticfname; 148 CeedQuadMode qmode; 149 } problemData; 150 // *INDENT-ON* 151 152 // Data specific to each problem option 153 extern problemData problemOptions[3]; 154 155 // Forcing function data 156 typedef struct { 157 CeedQFunctionUser setupforcing; 158 const char *setupforcingfname; 159 } forcingData; 160 161 extern forcingData forcingOptions[3]; 162 163 // Data for PETSc Matshell 164 typedef struct UserMult_private *UserMult; 165 struct UserMult_private { 166 MPI_Comm comm; 167 DM dm; 168 Vec Xloc, Yloc, NBCs; 169 CeedVector Xceed, Yceed; 170 CeedOperator op; 171 CeedQFunction qf; 172 Ceed ceed; 173 PetscScalar loadIncrement; 174 CeedQFunctionContext ctxPhys, ctxPhysSmoother; 175 CeedMemType memType; 176 int (*VecGetArray)(Vec, PetscScalar **); 177 int (*VecGetArrayRead)(Vec, const PetscScalar **); 178 int (*VecRestoreArray)(Vec, PetscScalar **); 179 int (*VecRestoreArrayRead)(Vec, const PetscScalar **); 180 }; 181 182 // Data for Jacobian setup routine 183 typedef struct FormJacobCtx_private *FormJacobCtx; 184 struct FormJacobCtx_private { 185 UserMult *jacobCtx; 186 PetscInt numLevels; 187 SNES snesCoarse; 188 Mat *jacobMat, jacobMatCoarse; 189 Vec Ucoarse; 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 dmC, dmF; 197 Vec locVecC, locVecF; 198 CeedVector ceedVecC, ceedVecF; 199 CeedOperator opProlong, opRestrict; 200 Ceed ceed; 201 CeedMemType memType; 202 int (*VecGetArray)(Vec, PetscScalar **); 203 int (*VecGetArrayRead)(Vec, const PetscScalar **); 204 int (*VecRestoreArray)(Vec, PetscScalar **); 205 int (*VecRestoreArrayRead)(Vec, const PetscScalar **); 206 }; 207 208 // libCEED data struct for level 209 typedef struct CeedData_private *CeedData; 210 struct CeedData_private { 211 Ceed ceed; 212 CeedBasis basisx, basisu, basisCtoF, basisEnergy, basisDiagnostic; 213 CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi, 214 ErestrictGradui, ErestrictEnergy, ErestrictDiagnostic, 215 ErestrictqdDiagnostici; 216 CeedQFunction qfApply, qfJacob, qfEnergy, qfDiagnostic; 217 CeedOperator opApply, opJacob, opRestrict, opProlong, opEnergy, 218 opDiagnostic; 219 CeedVector qdata, qdataDiagnostic, gradu, xceed, yceed, truesoln; 220 }; 221 222 // ----------------------------------------------------------------------------- 223 // Process command line options 224 // ----------------------------------------------------------------------------- 225 // Process general command line options 226 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx appCtx); 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 isSimplex, const char prefix[], 239 PetscInt order, PetscFE *fem); 240 241 // Read mesh and distribute DM in parallel 242 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx appCtx, DM *dm); 243 244 // Setup DM with FE space of appropriate degree 245 PetscErrorCode SetupDMByDegree(DM dm, AppCtx appCtx, PetscInt order, 246 PetscBool boundary, PetscInt ncompu); 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 domainLabel, CeedInt value, 260 CeedElemRestriction *Erestrict); 261 262 // Utility function to get Ceed Restriction for each domain 263 PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, 264 DMLabel domainLabel, PetscInt value, CeedInt P, CeedInt Q, CeedInt qdatasize, 265 CeedElemRestriction *restrictq, CeedElemRestriction *restrictx, 266 CeedElemRestriction *restrictqdi); 267 268 // Set up libCEED for a given degree 269 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dmEnergy, DM dmDiagnostic, 270 Ceed ceed, AppCtx appCtx, 271 CeedQFunctionContext physCtx, 272 CeedData *data, PetscInt fineLevel, 273 PetscInt ncompu, PetscInt Ugsz, 274 PetscInt Ulocsz, CeedVector forceCeed, 275 CeedVector neumannCeed); 276 277 // Set up libCEED multigrid level for a given degree 278 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx appCtx, 279 CeedData *data, PetscInt level, 280 PetscInt ncompu, PetscInt Ugsz, 281 PetscInt Ulocsz, CeedVector fineMult); 282 283 // Setup context data for Jacobian evaluation 284 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V, 285 Vec Vloc, CeedData ceedData, Ceed ceed, 286 CeedQFunctionContext ctxPhys, 287 CeedQFunctionContext ctxPhysSmoother, 288 UserMult jacobianCtx); 289 290 // Setup context data for prolongation and restriction operators 291 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx appCtx, DM dmC, 292 DM dmF, Vec VF, Vec VlocC, Vec VlocF, 293 CeedData ceedDataC, CeedData ceedDataF, 294 Ceed ceed, 295 UserMultProlongRestr prolongRestrCtx); 296 297 // ----------------------------------------------------------------------------- 298 // Jacobian setup 299 // ----------------------------------------------------------------------------- 300 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx); 301 302 // ----------------------------------------------------------------------------- 303 // Solution output 304 // ----------------------------------------------------------------------------- 305 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment, 306 PetscScalar loadIncrement); 307 308 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, 309 UserMult user, Vec U, 310 CeedElemRestriction ErestrictDiagnostic); 311 312 // ----------------------------------------------------------------------------- 313 // libCEED Operators for MatShell 314 // ----------------------------------------------------------------------------- 315 // This function uses libCEED to compute the local action of an operator 316 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user); 317 318 // This function uses libCEED to compute the non-linear residual 319 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 320 321 // This function uses libCEED to apply the Jacobian for assembly via a SNES 322 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 323 324 // This function uses libCEED to compute the action of the Jacobian 325 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y); 326 327 // This function uses libCEED to compute the action of the prolongation operator 328 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y); 329 330 // This function uses libCEED to compute the action of the restriction operator 331 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y); 332 333 // This function returns the computed diagonal of the operator 334 PetscErrorCode GetDiag_Ceed(Mat A, Vec D); 335 336 // This function calculates the strain energy in the final solution 337 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 338 CeedOperator opEnergy, Vec X, 339 PetscReal *energy); 340 341 // ----------------------------------------------------------------------------- 342 // Boundary Functions 343 // ----------------------------------------------------------------------------- 344 // Note: If additional boundary conditions are added, an update is needed in 345 // elasticity.h for the boundaryOptions variable. 346 347 // BCMMS - boundary function 348 // Values on all points of the mesh is set based on given solution below 349 // for u[0], u[1], u[2] 350 PetscErrorCode BCMMS(PetscInt dim, PetscReal loadIncrement, 351 const PetscReal coords[], PetscInt ncompu, 352 PetscScalar *u, void *ctx); 353 354 // BCClamp - fix boundary values with affine transformation at fraction of load 355 // increment 356 PetscErrorCode BCClamp(PetscInt dim, PetscReal loadIncrement, 357 const PetscReal coords[], PetscInt ncompu, 358 PetscScalar *u, void *ctx); 359 360 #endif //setup_h 361