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 {FVBC_PERIODIC, FVBC_OUTFLOW} FVBCType; 34 extern const char *FVBCTypes[]; 35 /* 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 typedef PetscErrorCode (*RiemannFunction)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*,PetscReal,PetscReal,PetscReal); 37 typedef PetscErrorCode (*ReconstructFunction)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*,PetscReal); 38 39 PetscErrorCode RiemannListAdd(PetscFunctionList*,const char*,RiemannFunction); 40 PetscErrorCode RiemannListFind(PetscFunctionList,const char*,RiemannFunction*); 41 PetscErrorCode ReconstructListAdd(PetscFunctionList*,const char*,ReconstructFunction); 42 PetscErrorCode ReconstructListFind(PetscFunctionList,const char*,ReconstructFunction*); 43 PetscErrorCode PhysicsDestroy_SimpleFree(void*); 44 45 typedef struct { 46 PetscErrorCode (*sample)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*); 47 RiemannFunction riemann; 48 ReconstructFunction characteristic; 49 PetscErrorCode (*destroy)(void*); 50 void *user; 51 PetscInt dof; 52 char *fieldname[16]; 53 } PhysicsCtx; 54 55 void Limit2_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 56 void Limit2_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 57 void Limit2_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 58 void Limit2_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 59 void Limit2_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 60 void Limit2_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 61 void Limit2_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 62 void Limit2_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt len_slow,const PetscInt len_fast,PetscInt n,PetscScalar *lmt); 63 64 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 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 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 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 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 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 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 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 73 typedef PetscErrorCode (*RiemannFunction_2WaySplit)(void*,PetscInt,const PetscScalar*,const PetscScalar*,PetscScalar*,PetscReal*); 74 typedef PetscErrorCode (*ReconstructFunction_2WaySplit)(void*,PetscInt,const PetscScalar*,PetscScalar*,PetscScalar*,PetscReal*); 75 76 PetscErrorCode RiemannListAdd_2WaySplit(PetscFunctionList*,const char*,RiemannFunction_2WaySplit); 77 PetscErrorCode RiemannListFind_2WaySplit(PetscFunctionList,const char*,RiemannFunction_2WaySplit*); 78 PetscErrorCode ReconstructListAdd_2WaySplit(PetscFunctionList*,const char*,ReconstructFunction_2WaySplit); 79 PetscErrorCode ReconstructListFind_2WaySplit(PetscFunctionList,const char*,ReconstructFunction_2WaySplit*); 80 81 typedef struct { 82 PetscErrorCode (*sample2)(void*,PetscInt,FVBCType,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal*); 83 RiemannFunction_2WaySplit riemann2; 84 ReconstructFunction_2WaySplit characteristic2; 85 PetscErrorCode (*destroy)(void*); 86 void *user; 87 PetscInt dof; 88 char *fieldname[16]; 89 } PhysicsCtx2; 90 91 typedef struct { 92 void (*limit)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscScalar*); 93 PhysicsCtx physics; 94 MPI_Comm comm; 95 char prefix[256]; 96 97 /* Local work arrays */ 98 PetscScalar *R,*Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */ 99 PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */ 100 PetscScalar *cslope; /* Limited slope, written in characteristic basis */ 101 PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */ 102 PetscScalar *flux; /* Flux across interface */ 103 PetscReal *speeds; /* Speeds of each wave */ 104 105 PetscReal cfl_idt; /* Max allowable value of 1/Delta t */ 106 PetscReal cfl; 107 PetscReal xmin,xmax; 108 PetscInt initial; 109 PetscBool simulation; 110 FVBCType bctype; 111 PetscBool exact; 112 113 /* context for partitioned system */ 114 void (*limit3)(LimitInfo,const PetscScalar*,const PetscScalar*,const PetscInt,const PetscInt,const PetscInt,const PetscInt,PetscInt,PetscScalar*); 115 void (*limit2)(LimitInfo,const PetscScalar*,const PetscScalar*,PetscInt,PetscInt,PetscInt,PetscScalar*); 116 PhysicsCtx2 physics2; 117 PetscInt hratio; /* hratio = hslow/hfast */ 118 IS isf,iss,isf2,iss2,ism,issb,ismb; 119 PetscBool recursive; 120 PetscInt sm,mf,fm,ms; /* positions (array index) for slow-medium, medium-fast, fast-medium, medium-slow interfaces */ 121 PetscInt sf,fs; /* slow-fast and fast-slow interfaces */ 122 PetscInt lsbwidth,rsbwidth; /* left slow buffer width and right slow buffer width */ 123 PetscInt lmbwidth,rmbwidth; /* left medium buffer width and right medium buffer width */ 124 } FVCtx; 125 126 /* --------------------------------- Finite Volume Solver ----------------------------------- */ 127 PetscErrorCode FVRHSFunction(TS,PetscReal,Vec,Vec,void*); 128 PetscErrorCode FVSample(FVCtx*,DM,PetscReal,Vec); 129 PetscErrorCode SolutionStatsView(DM,Vec,PetscViewer); 130