1 #include <petscts.h> 2 #include <petscdm.h> 3 4 typedef struct _LimitInfo { 5 PetscReal hx; 6 PetscInt m; 7 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 *); 30 31 /* --------------------------------- Finite Volume data structures ----------------------------------- */ 32 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 bounday 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); 42 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 *); 48 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; 58 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); 67 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); 76 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 *); 79 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 *); 84 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; 98 99 typedef struct { 100 void (*limit)(LimitInfo, const PetscScalar *, const PetscScalar *, PetscScalar *); 101 PhysicsCtx physics; 102 MPI_Comm comm; 103 char prefix[256]; 104 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 */ 113 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; 121 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; 134 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); 139