1c4762a1bSJed Brown #include <petscts.h> 2c4762a1bSJed Brown #include <petscdm.h> 3c4762a1bSJed Brown 4c4762a1bSJed Brown typedef struct _LimitInfo { 5c4762a1bSJed Brown PetscReal hx; 6c4762a1bSJed Brown PetscInt m; 7c4762a1bSJed Brown 8c4762a1bSJed Brown /* context for partitioned system */ 9c4762a1bSJed Brown PetscReal hxs,hxm,hxf; 10c4762a1bSJed Brown } *LimitInfo; 11c4762a1bSJed Brown void Limit_Upwind(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 12c4762a1bSJed Brown void Limit_LaxWendroff(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 13c4762a1bSJed Brown void Limit_BeamWarming(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 14c4762a1bSJed Brown void Limit_Fromm(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 15c4762a1bSJed Brown void Limit_Minmod(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 16c4762a1bSJed Brown void Limit_Superbee(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 17c4762a1bSJed Brown void Limit_MC(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 18c4762a1bSJed Brown void Limit_VanLeer(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 19c4762a1bSJed Brown void Limit_VanAlbada(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); /* differentiable */ 20c4762a1bSJed Brown void Limit_VanAlbadaTVD(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 21c4762a1bSJed Brown void Limit_Koren(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); /* differentiable */ 22c4762a1bSJed Brown void Limit_KorenSym(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); /* differentiable */ 23c4762a1bSJed Brown void Limit_Koren3(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 24c4762a1bSJed Brown void Limit_CadaTorrilhon2(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 25c4762a1bSJed Brown void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 26c4762a1bSJed Brown void Limit_CadaTorrilhon3R0p1(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 27c4762a1bSJed Brown void Limit_CadaTorrilhon3R1(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 28c4762a1bSJed Brown void Limit_CadaTorrilhon3R10(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 29c4762a1bSJed Brown void Limit_CadaTorrilhon3R100(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 30c4762a1bSJed Brown 31c4762a1bSJed Brown /* --------------------------------- Finite Volume data structures ----------------------------------- */ 32c4762a1bSJed Brown 33*84ff8c8bSHong Zhang typedef enum {FVBC_PERIODIC, FVBC_OUTFLOW, FVBC_INFLOW} FVBCType; 34c4762a1bSJed Brown extern const char *FVBCTypes[]; 35*84ff8c8bSHong Zhang /* we add three new variables at the end of input parameters of function to be position of cell center, left bounday of domain, right boundary of domain */ 36c4762a1bSJed Brown typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*,PetscReal,PetscReal,PetscReal); 37c4762a1bSJed Brown typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*,PetscReal); 38c4762a1bSJed Brown 39c4762a1bSJed Brown PetscErrorCode RiemannListAdd(PetscFunctionList*,const char*,RiemannFunction); 40c4762a1bSJed Brown PetscErrorCode RiemannListFind(PetscFunctionList,const char*,RiemannFunction*); 41c4762a1bSJed Brown PetscErrorCode ReconstructListAdd(PetscFunctionList*,const char*,ReconstructFunction); 42c4762a1bSJed Brown PetscErrorCode ReconstructListFind(PetscFunctionList,const char*,ReconstructFunction*); 43c4762a1bSJed Brown PetscErrorCode PhysicsDestroy_SimpleFree(void*); 44c4762a1bSJed Brown 45c4762a1bSJed Brown typedef struct { 46c4762a1bSJed Brown PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*); 47c4762a1bSJed Brown RiemannFunction riemann; 48c4762a1bSJed Brown ReconstructFunction characteristic; 49c4762a1bSJed Brown PetscErrorCode (*destroy)(void*); 50c4762a1bSJed Brown void *user; 51c4762a1bSJed Brown PetscInt dof; 52c4762a1bSJed Brown char *fieldname[16]; 53c4762a1bSJed Brown } PhysicsCtx; 54c4762a1bSJed Brown 55c4762a1bSJed Brown void Limit2_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 56c4762a1bSJed Brown void Limit2_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 57c4762a1bSJed Brown void Limit2_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 58c4762a1bSJed Brown void Limit2_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 59c4762a1bSJed Brown void Limit2_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 60c4762a1bSJed Brown void Limit2_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 61c4762a1bSJed Brown void Limit2_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 62c4762a1bSJed Brown void Limit2_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 63c4762a1bSJed Brown 64c4762a1bSJed Brown void Limit3_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 65c4762a1bSJed Brown void Limit3_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 66c4762a1bSJed Brown void Limit3_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 67c4762a1bSJed Brown void Limit3_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 68c4762a1bSJed Brown void Limit3_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 69c4762a1bSJed Brown void Limit3_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 70c4762a1bSJed Brown void Limit3_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 71c4762a1bSJed Brown void Limit3_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sm,const PetscInt mf,const PetscInt fm,const PetscInt ms,PetscInt n,PetscScalar *lmt); 72c4762a1bSJed Brown 73c4762a1bSJed Brown typedef PetscErrorCode (*RiemannFunction_2WaySplit)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*); 74c4762a1bSJed Brown typedef PetscErrorCode (*ReconstructFunction_2WaySplit)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*); 75c4762a1bSJed Brown 76c4762a1bSJed Brown PetscErrorCode RiemannListAdd_2WaySplit(PetscFunctionList*,const char*,RiemannFunction_2WaySplit); 77c4762a1bSJed Brown PetscErrorCode RiemannListFind_2WaySplit(PetscFunctionList,const char*,RiemannFunction_2WaySplit*); 78c4762a1bSJed Brown PetscErrorCode ReconstructListAdd_2WaySplit(PetscFunctionList*,const char*,ReconstructFunction_2WaySplit); 79c4762a1bSJed Brown PetscErrorCode ReconstructListFind_2WaySplit(PetscFunctionList,const char*,ReconstructFunction_2WaySplit*); 80c4762a1bSJed Brown 81c4762a1bSJed Brown typedef struct { 82c4762a1bSJed Brown PetscErrorCode (*sample2)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*); 83*84ff8c8bSHong Zhang PetscErrorCode (*inflow)(void*,PetscReal,PetscReal,PetscReal*); 84c4762a1bSJed Brown RiemannFunction_2WaySplit riemann2; 85c4762a1bSJed Brown ReconstructFunction_2WaySplit characteristic2; 86c4762a1bSJed Brown PetscErrorCode (*destroy)(void*); 87c4762a1bSJed Brown void *user; 88c4762a1bSJed Brown PetscInt dof; 89c4762a1bSJed Brown char *fieldname[16]; 90*84ff8c8bSHong Zhang PetscBool *bcinflowindex; /* Boolean array where bcinflowindex[dof*i+j] = TRUE indicates that the jth component of the solution 91*84ff8c8bSHong Zhang is an inflow boundary condition and i = 0 is left bc, i = 1 is right bc. FALSE implies outflow 92*84ff8c8bSHong Zhang outflow boundary condition. */ 93c4762a1bSJed Brown } PhysicsCtx2; 94c4762a1bSJed Brown 95c4762a1bSJed Brown typedef struct { 96c4762a1bSJed Brown void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 97c4762a1bSJed Brown PhysicsCtx physics; 98c4762a1bSJed Brown MPI_Comm comm; 99c4762a1bSJed Brown char prefix[256]; 100c4762a1bSJed Brown 101c4762a1bSJed Brown /* Local work arrays */ 102c4762a1bSJed Brown PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */ 103c4762a1bSJed Brown PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */ 104c4762a1bSJed Brown PetscScalar *cslope; /* Limited slope, written in characteristic basis */ 105c4762a1bSJed Brown PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */ 106c4762a1bSJed Brown PetscScalar *flux; /* Flux across interface */ 107c4762a1bSJed Brown PetscReal *speeds; /* Speeds of each wave */ 108*84ff8c8bSHong Zhang PetscReal *ub; /* Boundary data for inflow boundary conditions */ 109c4762a1bSJed Brown 110c4762a1bSJed Brown PetscReal cfl_idt; /* Max allowable value of 1/Delta t */ 111c4762a1bSJed Brown PetscReal cfl; 112c4762a1bSJed Brown PetscReal xmin,xmax; 113c4762a1bSJed Brown PetscInt initial; 114c4762a1bSJed Brown PetscBool simulation; 115c4762a1bSJed Brown FVBCType bctype; 116c4762a1bSJed Brown PetscBool exact; 117c4762a1bSJed Brown 118c4762a1bSJed Brown /* context for partitioned system */ 119c4762a1bSJed Brown void (*limit3)(LimitInfo,const PetscScalar*,const PetscScalar*,const PetscInt,const PetscInt,const PetscInt,const PetscInt,PetscInt,PetscScalar*); 120c4762a1bSJed Brown void (*limit2)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscInt,PetscInt,PetscInt,PetscScalar*); 121c4762a1bSJed Brown PhysicsCtx2 physics2; 122c4762a1bSJed Brown PetscInt hratio; /* hratio = hslow/hfast */ 123c4762a1bSJed Brown IS isf,iss,isf2,iss2,ism,issb,ismb; 124c4762a1bSJed Brown PetscBool recursive; 125c4762a1bSJed Brown PetscInt sm,mf,fm,ms; /* positions (array index) for slow-medium, medium-fast, fast-medium, medium-slow interfaces */ 126c4762a1bSJed Brown PetscInt sf,fs; /* slow-fast and fast-slow interfaces */ 127c4762a1bSJed Brown PetscInt lsbwidth,rsbwidth; /* left slow buffer width and right slow buffer width */ 128c4762a1bSJed Brown PetscInt lmbwidth,rmbwidth; /* left medium buffer width and right medium buffer width */ 129c4762a1bSJed Brown } FVCtx; 130c4762a1bSJed Brown 131c4762a1bSJed Brown /* --------------------------------- Finite Volume Solver ----------------------------------- */ 132c4762a1bSJed Brown PetscErrorCode FVRHSFunction(TS,PetscReal,Vec,Vec,void*); 133c4762a1bSJed Brown PetscErrorCode FVSample(FVCtx*,DM,PetscReal,Vec); 134c4762a1bSJed Brown PetscErrorCode SolutionStatsView(DM,Vec,PetscViewer); 135