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 // ----------------------------------------------------------------------------- 91 // Structs 92 // ----------------------------------------------------------------------------- 93 // Units 94 typedef struct Units_private *Units; 95 struct Units_private { 96 // Fundamental units 97 PetscScalar meter; 98 PetscScalar kilogram; 99 PetscScalar second; 100 // Derived unit 101 PetscScalar Pascal; 102 }; 103 104 // Application context from user command line options 105 typedef struct AppCtx_private *AppCtx; 106 struct AppCtx_private { 107 char ceedResource[PETSC_MAX_PATH_LEN]; // libCEED backend 108 char ceedResourceFine[PETSC_MAX_PATH_LEN]; // libCEED for fine grid 109 char meshFile[PETSC_MAX_PATH_LEN]; // exodusII mesh file 110 PetscBool testMode; 111 PetscBool viewSoln; 112 problemType problemChoice; 113 forcingType forcingChoice; 114 multigridType multigridChoice; 115 PetscInt degree; 116 PetscInt numLevels; 117 PetscInt *levelDegrees; 118 PetscInt numIncrements; // Number of steps 119 PetscInt bcZeroFaces[16], bcClampFaces[16]; 120 PetscInt bcZeroCount, bcClampCount; 121 PetscBool bcClampTranslate[16]; 122 PetscScalar bcClampMax[16][4]; 123 PetscScalar forcingVector[3]; 124 }; 125 126 // Problem specific data 127 typedef struct { 128 CeedInt qdatasize; 129 CeedQFunctionUser setupgeo, apply, jacob, energy; 130 const char *setupgeofname, *applyfname, *jacobfname, *energyfname; 131 CeedQuadMode qmode; 132 } problemData; 133 134 // Data specific to each problem option 135 problemData problemOptions[3]; 136 137 // Forcing function data 138 typedef struct { 139 CeedQFunctionUser setupforcing; 140 const char *setupforcingfname; 141 } forcingData; 142 143 forcingData forcingOptions[3]; 144 145 // Data for PETSc Matshell 146 typedef struct UserMult_private *UserMult; 147 struct UserMult_private { 148 MPI_Comm comm; 149 DM dm; 150 Vec Xloc, Yloc; 151 CeedVector Xceed, Yceed; 152 CeedOperator op; 153 Ceed ceed; 154 PetscScalar loadIncrement; 155 }; 156 157 // Data for Jacobian setup routine 158 typedef struct FormJacobCtx_private *FormJacobCtx; 159 struct FormJacobCtx_private { 160 UserMult *jacobCtx; 161 PetscInt numLevels; 162 SNES snesCoarse; 163 Mat *jacobMat, jacobMatCoarse; 164 Vec Ucoarse; 165 }; 166 167 // Data for PETSc Prolongation/Restriction Matshell 168 typedef struct UserMultProlongRestr_private *UserMultProlongRestr; 169 struct UserMultProlongRestr_private { 170 MPI_Comm comm; 171 DM dmC, dmF; 172 Vec locVecC, locVecF, multVec; 173 CeedVector ceedVecC, ceedVecF; 174 CeedOperator opProlong, opRestrict; 175 Ceed ceed; 176 }; 177 178 // libCEED data struct for level 179 typedef struct CeedData_private *CeedData; 180 struct CeedData_private { 181 Ceed ceed; 182 CeedBasis basisx, basisu, basisCtoF, basisEnergy; 183 CeedElemRestriction Erestrictx, Erestrictu, Erestrictqdi, 184 ErestrictGradui, ErestrictEnergy; 185 CeedQFunction qfApply, qfJacob, qfEnergy; 186 CeedOperator opApply, opJacob, opRestrict, opProlong, opEnergy; 187 CeedVector qdata, gradu, xceed, yceed, truesoln, energy; 188 }; 189 190 // ----------------------------------------------------------------------------- 191 // Process command line options 192 // ----------------------------------------------------------------------------- 193 // Process general command line options 194 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx appCtx); 195 196 // Process physics options 197 PetscErrorCode ProcessPhysics(MPI_Comm comm, Physics phys, Units units); 198 199 // ----------------------------------------------------------------------------- 200 // Setup DM 201 // ----------------------------------------------------------------------------- 202 PetscErrorCode CreateBCLabel(DM dm, const char name[]); 203 204 // Create FE by degree 205 PetscErrorCode PetscFECreateByDegree(DM dm, PetscInt dim, PetscInt Nc, 206 PetscBool isSimplex, const char prefix[], 207 PetscInt order, PetscFE *fem); 208 209 // Read mesh and distribute DM in parallel 210 PetscErrorCode CreateDistributedDM(MPI_Comm comm, AppCtx appCtx, DM *dm); 211 212 // Setup DM with FE space of appropriate degree 213 PetscErrorCode SetupDMByDegree(DM dm, AppCtx appCtx, PetscInt order, 214 PetscInt ncompu); 215 216 // ----------------------------------------------------------------------------- 217 // libCEED Functions 218 // ----------------------------------------------------------------------------- 219 // Destroy libCEED objects 220 PetscErrorCode CeedDataDestroy(CeedInt level, CeedData data); 221 222 // Get libCEED restriction data from DMPlex 223 PetscErrorCode CreateRestrictionPlex(Ceed ceed, CeedInterlaceMode imode, 224 CeedInt P, CeedInt ncomp, 225 CeedElemRestriction *Erestrict, DM dm); 226 227 // Set up libCEED for a given degree 228 PetscErrorCode SetupLibceedFineLevel(DM dm, Ceed ceed, AppCtx appCtx, 229 Physics phys, CeedData *data, 230 PetscInt fineLevel, PetscInt ncompu, 231 PetscInt Ugsz, PetscInt Ulocsz, 232 CeedVector forceCeed, 233 CeedQFunction qfRestrict, 234 CeedQFunction qfProlong); 235 236 // Set up libCEED for a given degree 237 PetscErrorCode SetupLibceedLevel(DM dm, Ceed ceed, AppCtx appCtx, Physics phys, 238 CeedData *data, PetscInt level, 239 PetscInt ncompu, PetscInt Ugsz, 240 PetscInt Ulocsz, CeedVector forceCeed, 241 CeedQFunction qfRestrict, 242 CeedQFunction qfProlong); 243 244 // Setup context data for Jacobian evaluation 245 PetscErrorCode SetupJacobianCtx(MPI_Comm comm, AppCtx appCtx, DM dm, Vec V, 246 Vec Vloc, CeedData ceedData, Ceed ceed, 247 UserMult jacobianCtx); 248 249 // Setup context data for prolongation and restriction operators 250 PetscErrorCode SetupProlongRestrictCtx(MPI_Comm comm, DM dmC, DM dmF, Vec VF, 251 Vec VlocC, Vec VlocF, CeedData ceedDataC, 252 CeedData ceedDataF, Ceed ceed, 253 UserMultProlongRestr prolongRestrCtx); 254 255 // ----------------------------------------------------------------------------- 256 // Jacobian setup 257 // ----------------------------------------------------------------------------- 258 PetscErrorCode FormJacobian(SNES snes, Vec U, Mat J, Mat Jpre, void *ctx); 259 260 // ----------------------------------------------------------------------------- 261 // SNES Monitor 262 // ----------------------------------------------------------------------------- 263 PetscErrorCode ViewSolution(MPI_Comm comm, Vec U, PetscInt increment, 264 PetscScalar loadIncrement); 265 266 // ----------------------------------------------------------------------------- 267 // libCEED Operators for MatShell 268 // ----------------------------------------------------------------------------- 269 // This function uses libCEED to compute the local action of an operator 270 PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user); 271 272 // This function uses libCEED to compute the non-linear residual 273 PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 274 275 // This function uses libCEED to apply the Jacobian for assembly via a SNES 276 PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx); 277 278 // This function uses libCEED to compute the action of the Jacobian 279 PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y); 280 281 // This function uses libCEED to compute the action of the prolongation operator 282 PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y); 283 284 // This function uses libCEED to compute the action of the restriction operator 285 PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y); 286 287 // This function returns the computed diagonal of the operator 288 PetscErrorCode GetDiag_Ceed(Mat A, Vec D); 289 290 // This function calculates the strain energy in the final solution 291 PetscErrorCode ComputeStrainEnergy(UserMult user, CeedOperator opEnergy, Vec X, 292 CeedVector energyLoc, PetscReal *energy); 293 294 // ----------------------------------------------------------------------------- 295 // Boundary Functions 296 // ----------------------------------------------------------------------------- 297 // Note: If additional boundary conditions are added, an update is needed in 298 // elasticity.h for the boundaryOptions variable. 299 300 // BCMMS - boundary function 301 // Values on all points of the mesh is set based on given solution below 302 // for u[0], u[1], u[2] 303 PetscErrorCode BCMMS(PetscInt dim, PetscReal loadIncrement, 304 const PetscReal coords[], PetscInt ncompu, 305 PetscScalar *u, void *ctx); 306 307 // BCZero - fix boundary values at zero 308 PetscErrorCode BCZero(PetscInt dim, PetscReal loadIncrement, 309 const PetscReal coords[], PetscInt ncompu, 310 PetscScalar *u, void *ctx); 311 312 // BCClampTranslate - translate boundary values at fraction of load increment 313 PetscErrorCode BCClampTranslate(PetscInt dim, PetscReal loadIncrement, 314 const PetscReal coords[], PetscInt ncompu, 315 PetscScalar *u, void *ctx); 316 317 // BCClampRotate - rotate boundary values at fraction of load increment 318 PetscErrorCode BCClampRotate(PetscInt dim, PetscReal loadIncrement, 319 const PetscReal coords[], PetscInt ncompu, 320 PetscScalar *u, void *ctx); 321 322 #endif //setup_h 323