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