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