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 #ifndef PHYSICS_STRUCT 31 #define PHYSICS_STRUCT 32 typedef struct Physics_private *Physics; 33 struct Physics_private { 34 CeedScalar nu; // Poisson's ratio 35 CeedScalar E; // Young's Modulus 36 }; 37 #endif 38 39 // ----------------------------------------------------------------------------- 40 // Command Line Options 41 // ----------------------------------------------------------------------------- 42 // Problem options 43 typedef enum { 44 ELAS_LIN = 0, ELAS_HYPER_SS = 1, ELAS_HYPER_FS = 2 45 } problemType; 46 static const char *const problemTypes[] = {"linElas", 47 "hyperSS", 48 "hyperFS", 49 "problemType","ELAS_",0 50 }; 51 static const char *const problemTypesForDisp[] = {"Linear elasticity", 52 "Hyper elasticity small strain", 53 "Hyper elasticity finite strain" 54 }; 55 56 // Forcing function options 57 typedef enum { 58 FORCE_NONE = 0, FORCE_CONST = 1, FORCE_MMS = 2 59 } forcingType; 60 static const char *const forcingTypes[] = {"none", 61 "constant", 62 "mms", 63 "forcingType","FORCE_",0 64 }; 65 static const char *const forcingTypesForDisp[] = {"None", 66 "Constant", 67 "Manufactured solution" 68 }; 69 70 // Multigrid options 71 typedef enum { 72 MULTIGRID_LOGARITHMIC = 0, MULTIGRID_UNIFORM = 1, MULTIGRID_NONE = 2 73 } multigridType; 74 static const char *const multigridTypes [] = {"logarithmic", 75 "uniform", 76 "none", 77 "multigridType","MULTIGRID",0 78 }; 79 static const char *const multigridTypesForDisp[] = {"P-multigrid, logarithmic coarsening", 80 "P-multigrind, uniform coarsening", 81 "No multigrid" 82 }; 83 84 typedef PetscErrorCode BCFunc(PetscInt, PetscReal, const PetscReal *, PetscInt, 85 PetscScalar *, void *); 86 // Note: These variables should be updated if additional boundary conditions 87 // are added to boundary.c. 88 BCFunc BCMMS, BCZero, BCClamp; 89 90 // MemType Options 91 static const char *const memTypes[] = {"host","device","memType", 92 "CEED_MEM_",0 93 }; 94 95 // ----------------------------------------------------------------------------- 96 // Structs 97 // ----------------------------------------------------------------------------- 98 // Units 99 typedef struct Units_private *Units; 100 struct Units_private { 101 // Fundamental units 102 PetscScalar meter; 103 PetscScalar kilogram; 104 PetscScalar second; 105 // Derived unit 106 PetscScalar Pascal; 107 }; 108 109 // Application context from user command line options 110 typedef struct AppCtx_private *AppCtx; 111 struct AppCtx_private { 112 char ceedResource[PETSC_MAX_PATH_LEN]; // libCEED backend 113 char ceedResourceFine[PETSC_MAX_PATH_LEN]; // libCEED for fine grid 114 char meshFile[PETSC_MAX_PATH_LEN]; // exodusII mesh file 115 PetscBool testMode; 116 PetscBool viewSoln; 117 PetscBool viewFinalSoln; 118 problemType problemChoice; 119 forcingType forcingChoice; 120 multigridType multigridChoice; 121 PetscScalar nuSmoother; 122 PetscInt degree; 123 PetscInt qextra; 124 PetscInt numLevels; 125 PetscInt *levelDegrees; 126 PetscInt numIncrements; // Number of steps 127 PetscInt bcClampFaces[16]; 128 PetscInt bcClampCount; 129 PetscScalar bcClampMax[16][7]; 130 PetscScalar forcingVector[3]; 131 PetscBool petscHaveCuda, setMemTypeRequest; 132 CeedMemType memTypeRequested; 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; 163 CeedVector Xceed, Yceed; 164 CeedOperator op; 165 CeedQFunction qf; 166 Ceed ceed; 167 PetscScalar loadIncrement; 168 Physics phys, physSmoother; 169 CeedMemType memType; 170 int (*VecGetArray)(Vec, PetscScalar **); 171 int (*VecGetArrayRead)(Vec, const PetscScalar **); 172 int (*VecRestoreArray)(Vec, PetscScalar **); 173 int (*VecRestoreArrayRead)(Vec, const PetscScalar **); 174 }; 175 176 // Data for Jacobian setup routine 177 typedef struct FormJacobCtx_private *FormJacobCtx; 178 struct FormJacobCtx_private { 179 UserMult *jacobCtx; 180 PetscInt numLevels; 181 SNES snesCoarse; 182 Mat *jacobMat, jacobMatCoarse; 183 Vec Ucoarse; 184 }; 185 186 // Data for PETSc Prolongation/Restriction Matshell 187 typedef struct UserMultProlongRestr_private *UserMultProlongRestr; 188 struct UserMultProlongRestr_private { 189 MPI_Comm comm; 190 DM dmC, dmF; 191 Vec locVecC, locVecF, multVec; 192 CeedVector ceedVecC, ceedVecF; 193 CeedOperator opProlong, opRestrict; 194 Ceed ceed; 195 CeedMemType memType; 196 int (*VecGetArray)(Vec, PetscScalar **); 197 int (*VecGetArrayRead)(Vec, const PetscScalar **); 198 int (*VecRestoreArray)(Vec, PetscScalar **); 199 int (*VecRestoreArrayRead)(Vec, const PetscScalar **); 200 }; 201 202 // libCEED data struct for level 203 typedef struct CeedData_private *CeedData; 204 struct CeedData_private { 205 Ceed ceed; 206 CeedBasis basisx, basisu, basisCtoF, basisEnergy, basisDiagnostic; 207 CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi, 208 ErestrictGradui, ErestrictEnergy, ErestrictDiagnostic, 209 ErestrictqdDiagnostici; 210 CeedQFunction qfApply, qfJacob, qfEnergy, qfDiagnostic; 211 CeedOperator opApply, opJacob, opRestrict, opProlong, opEnergy, 212 opDiagnostic; 213 CeedVector qdata, qdataDiagnostic, gradu, xceed, yceed, truesoln; 214 }; 215 216 // ----------------------------------------------------------------------------- 217 // Process command line options 218 // ----------------------------------------------------------------------------- 219 // Process general command line options 220 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx appCtx); 221 222 // Process physics options 223 PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units); 224 225 // ----------------------------------------------------------------------------- 226 // Setup DM 227 // ----------------------------------------------------------------------------- 228 PetscErrorCode CreateBCLabel(DM dm, const char name[]); 229 230 // Create FE by degree 231 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 232 PetscBool isSimplex, const char prefix[], 233 PetscInt order, PetscFE *fem); 234 235 // Read mesh and distribute DM in parallel 236 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx appCtx, DM *dm); 237 238 // Setup DM with FE space of appropriate degree 239 PetscErrorCode SetupDMByDegree(DM dm, AppCtx appCtx, PetscInt order, 240 PetscBool boundary, PetscInt ncompu); 241 242 // ----------------------------------------------------------------------------- 243 // libCEED Functions 244 // ----------------------------------------------------------------------------- 245 // Destroy libCEED objects 246 PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data); 247 248 // Get libCEED restriction data from DMPlex 249 PetscErrorCode CreateRestrictionPlex(Ceed ceed, CeedInt P, CeedInt ncomp, 250 CeedElemRestriction *Erestrict, DM dm); 251 252 // Set up libCEED for a given degree 253 PetscErrorCode SetupLibceedFineLevel(DM dm, DM dmEnergy, DM dmDiagnostic, 254 Ceed ceed, AppCtx appCtx, Physics phys, 255 CeedData *data, PetscInt fineLevel, 256 PetscInt ncompu, PetscInt Ugsz, 257 PetscInt Ulocsz, CeedVector forceCeed, 258 CeedQFunction qfRestrict, 259 CeedQFunction qfProlong); 260 261 // Set up libCEED for a given degree 262 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx appCtx, Physics phys, 263 CeedData *data, PetscInt level, 264 PetscInt ncompu, PetscInt Ugsz, 265 PetscInt Ulocsz, CeedVector forceCeed, 266 CeedQFunction qfRestrict, 267 CeedQFunction qfProlong); 268 269 // Setup context data for Jacobian evaluation 270 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V, 271 Vec Vloc, CeedData ceedData, Ceed ceed, 272 Physics phys, Physics physSmoother, 273 UserMult jacobianCtx); 274 275 // Setup context data for prolongation and restriction operators 276 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, AppCtx appCtx, DM dmC, 277 DM dmF, Vec VF, Vec VlocC, Vec VlocF, 278 CeedData ceedDataC, CeedData ceedDataF, 279 Ceed ceed, 280 UserMultProlongRestr prolongRestrCtx); 281 282 // ----------------------------------------------------------------------------- 283 // Jacobian setup 284 // ----------------------------------------------------------------------------- 285 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx); 286 287 // ----------------------------------------------------------------------------- 288 // Solution output 289 // ----------------------------------------------------------------------------- 290 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment, 291 PetscScalar loadIncrement); 292 293 PetscErrorCode ViewDiagnosticQuantities(MPI_Comm comm, DM dmU, 294 UserMult user, Vec U, 295 CeedElemRestriction ErestrictDiagnostic); 296 297 // ----------------------------------------------------------------------------- 298 // libCEED Operators for MatShell 299 // ----------------------------------------------------------------------------- 300 // This function uses libCEED to compute the local action of an operator 301 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user); 302 303 // This function uses libCEED to compute the non-linear residual 304 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 305 306 // This function uses libCEED to apply the Jacobian for assembly via a SNES 307 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 308 309 // This function uses libCEED to compute the action of the Jacobian 310 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y); 311 312 // This function uses libCEED to compute the action of the prolongation operator 313 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y); 314 315 // This function uses libCEED to compute the action of the restriction operator 316 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y); 317 318 // This function returns the computed diagonal of the operator 319 PetscErrorCode GetDiag_Ceed(Mat A, Vec D); 320 321 // This function calculates the strain energy in the final solution 322 PetscErrorCode ComputeStrainEnergy(DM dmEnergy, UserMult user, 323 CeedOperator opEnergy, Vec X, 324 PetscReal *energy); 325 326 // ----------------------------------------------------------------------------- 327 // Boundary Functions 328 // ----------------------------------------------------------------------------- 329 // Note: If additional boundary conditions are added, an update is needed in 330 // elasticity.h for the boundaryOptions variable. 331 332 // BCMMS - boundary function 333 // Values on all points of the mesh is set based on given solution below 334 // for u[0], u[1], u[2] 335 PetscErrorCode BCMMS(PetscInt dim, PetscReal loadIncrement, 336 const PetscReal coords[], PetscInt ncompu, 337 PetscScalar *u, void *ctx); 338 339 // BCClamp - fix boundary values with affine transformation at fraction of load 340 // increment 341 PetscErrorCode BCClamp(PetscInt dim, PetscReal loadIncrement, 342 const PetscReal coords[], PetscInt ncompu, 343 PetscScalar *u, void *ctx); 344 345 #endif //setup_h 346