1 #include "finitevolume1d.h"
2 #include <petscdmda.h>
3 #include <petscdraw.h>
4 #include <petsc/private/tsimpl.h>
5
6 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */
7 const char *FVBCTypes[] = {"PERIODIC","OUTFLOW","FVBCType","FVBC_",0};
8
Sgn(PetscReal a)9 PETSC_STATIC_INLINE PetscReal Sgn(PetscReal a) { return (a<0) ? -1 : 1; }
Abs(PetscReal a)10 PETSC_STATIC_INLINE PetscReal Abs(PetscReal a) { return (a<0) ? 0 : a; }
Sqr(PetscReal a)11 PETSC_STATIC_INLINE PetscReal Sqr(PetscReal a) { return a*a; }
12 //PETSC_STATIC_INLINE PetscReal MaxAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) > PetscAbs(b)) ? a : b; }
MinAbs(PetscReal a,PetscReal b)13 PETSC_UNUSED PETSC_STATIC_INLINE PetscReal MinAbs(PetscReal a,PetscReal b) { return (PetscAbs(a) < PetscAbs(b)) ? a : b; }
MinMod2(PetscReal a,PetscReal b)14 PETSC_STATIC_INLINE PetscReal MinMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscAbs(b)); }
MaxMod2(PetscReal a,PetscReal b)15 PETSC_STATIC_INLINE PetscReal MaxMod2(PetscReal a,PetscReal b) { return (a*b<0) ? 0 : Sgn(a)*PetscMax(PetscAbs(a),PetscAbs(b)); }
MinMod3(PetscReal a,PetscReal b,PetscReal c)16 PETSC_STATIC_INLINE PetscReal MinMod3(PetscReal a,PetscReal b,PetscReal c) {return (a*b<0 || a*c<0) ? 0 : Sgn(a)*PetscMin(PetscAbs(a),PetscMin(PetscAbs(b),PetscAbs(c))); }
17
18 //PETSC_STATIC_INLINE PetscReal RangeMod(PetscReal a,PetscReal xmin,PetscReal xmax) { PetscReal range = xmax-xmin; return xmin +PetscFmodReal(range+PetscFmodReal(a,range),range); }
19
20
21 /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
Limit_Upwind(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)22 void Limit_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
23 {
24 PetscInt i;
25 for (i=0; i<info->m; i++) lmt[i] = 0;
26 }
Limit_LaxWendroff(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)27 void Limit_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
28 {
29 PetscInt i;
30 for (i=0; i<info->m; i++) lmt[i] = jR[i];
31 }
Limit_BeamWarming(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)32 void Limit_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
33 {
34 PetscInt i;
35 for (i=0; i<info->m; i++) lmt[i] = jL[i];
36 }
Limit_Fromm(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)37 void Limit_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
38 {
39 PetscInt i;
40 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i]);
41 }
Limit_Minmod(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)42 void Limit_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
43 {
44 PetscInt i;
45 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i]);
46 }
Limit_Superbee(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)47 void Limit_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
48 {
49 PetscInt i;
50 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]));
51 }
Limit_MC(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)52 void Limit_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
53 {
54 PetscInt i;
55 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i]);
56 }
Limit_VanLeer(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)57 void Limit_VanLeer(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
58 { /* phi = (t + abs(t)) / (1 + abs(t)) */
59 PetscInt i;
60 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Abs(jR[i])+Abs(jL[i])*jR[i])/(Abs(jL[i])+Abs(jR[i])+1e-15);
61 }
Limit_VanAlbada(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)62 void Limit_VanAlbada(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
63 { /* phi = (t + t^2) / (1 + t^2) */
64 PetscInt i;
65 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*Sqr(jR[i])+Sqr(jL[i])*jR[i])/(Sqr(jL[i])+Sqr(jR[i])+1e-15);
66 }
Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)67 void Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
68 { /* phi = (t + t^2) / (1 + t^2) */
69 PetscInt i;
70 for (i=0; i<info->m; i++) lmt[i] = (jL[i]*jR[i]<0) ? 0 : (jL[i]*Sqr(jR[i])+Sqr(jL[i])*jR[i])/(Sqr(jL[i])+Sqr(jR[i])+1e-15);
71 }
Limit_Koren(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)72 void Limit_Koren(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
73 { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
74 PetscInt i;
75 for (i=0; i<info->m; i++) lmt[i] = ((jL[i]*Sqr(jR[i])+2*Sqr(jL[i])*jR[i])/(2*Sqr(jL[i])-jL[i]*jR[i]+2*Sqr(jR[i])+1e-15));
76 }
Limit_KorenSym(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)77 void Limit_KorenSym(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt) /* differentiable */
78 { /* Symmetric version of above */
79 PetscInt i;
80 for (i=0; i<info->m; i++) lmt[i] = (1.5*(jL[i]*Sqr(jR[i])+Sqr(jL[i])*jR[i])/(2*Sqr(jL[i])-jL[i]*jR[i]+2*Sqr(jR[i])+1e-15));
81 }
Limit_Koren3(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)82 void Limit_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
83 { /* Eq 11 of Cada-Torrilhon 2009 */
84 PetscInt i;
85 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i]);
86 }
CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)87 static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)
88 {
89 return PetscMax(0,PetscMin((L+2*R)/3,PetscMax(-0.5*L,PetscMin(2*L,PetscMin((L+2*R)/3,1.6*R)))));
90 }
Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)91 void Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
92 { /* Cada-Torrilhon 2009, Eq 13 */
93 PetscInt i;
94 for (i=0; i<info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]);
95 }
Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)96 void Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
97 { /* Cada-Torrilhon 2009, Eq 22 */
98 /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
99 const PetscReal eps = 1e-7,hx = info->hx;
100 PetscInt i;
101 for (i=0; i<info->m; i++) {
102 const PetscReal eta = (Sqr(jL[i])+Sqr(jR[i]))/Sqr(r*hx);
103 lmt[i] = ((eta < 1-eps) ? (jL[i]+2*jR[i])/3 : ((eta>1+eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]) : 0.5*((1-(eta-1)/eps)*(jL[i]+2*jR[i])/3+(1+(eta+1)/eps)*CadaTorrilhonPhiHatR_Eq13(jL[i],jR[i]))));
104 }
105 }
Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)106 void Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
107 {
108 Limit_CadaTorrilhon3R(0.1,info,jL,jR,lmt);
109 }
Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)110 void Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
111 {
112 Limit_CadaTorrilhon3R(1,info,jL,jR,lmt);
113 }
Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)114 void Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
115 {
116 Limit_CadaTorrilhon3R(10,info,jL,jR,lmt);
117 }
Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)118 void Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,PetscScalar *lmt)
119 {
120 Limit_CadaTorrilhon3R(100,info,jL,jR,lmt);
121 }
122
123 /* ----------------------- Limiters for split systems ------------------------- */
124
Limit2_Upwind(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)125 void Limit2_Upwind(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
126 {
127 PetscInt i;
128 for (i=0; i<info->m; i++) lmt[i] = 0;
129 }
Limit2_LaxWendroff(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)130 void Limit2_LaxWendroff(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
131 {
132 PetscInt i;
133 if (n < sf-1) { /* slow components */
134 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs;
135 } else if (n == sf-1) { /* slow component which is next to fast components */
136 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxs/2.0+info->hxf/2.0);
137 } else if (n == sf) { /* fast component which is next to slow components */
138 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxf;
139 } else if (n > sf && n < fs-1) { /* fast components */
140 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxf;
141 } else if (n == fs-1) { /* fast component next to slow components */
142 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxf/2.0+info->hxs/2.0);
143 } else if (n == fs) { /* slow component next to fast components */
144 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs;
145 } else { /* slow components */
146 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs;
147 }
148 }
Limit2_BeamWarming(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)149 void Limit2_BeamWarming(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
150 {
151 PetscInt i;
152 if (n < sf-1) {
153 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs;
154 } else if (n == sf-1) {
155 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs;
156 } else if (n == sf) {
157 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxs/2.0+info->hxf/2.0);
158 } else if (n > sf && n < fs-1) {
159 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxf;
160 } else if (n == fs-1) {
161 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxf;
162 } else if (n == fs) {
163 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxf/2.0+info->hxs/2.0);
164 } else {
165 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs;
166 }
167 }
Limit2_Fromm(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)168 void Limit2_Fromm(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
169 {
170 PetscInt i;
171 if (n < sf-1) {
172 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxs;
173 } else if (n == sf-1) {
174 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxs+jR[i]/(info->hxs/2.0+info->hxf/2.0));
175 } else if (n == sf) {
176 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxs/2.0+info->hxf/2.0)+jR[i]/info->hxf);
177 } else if (n > sf && n < fs-1) {
178 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxf;
179 } else if (n == fs-1) {
180 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxs/2.0));
181 } else if (n == fs) {
182 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxf/2.0+info->hxs/2.0)+jR[i]/info->hxs);
183 } else {
184 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxs;
185 }
186 }
Limit2_Minmod(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)187 void Limit2_Minmod(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
188 {
189 PetscInt i;
190 if (n < sf-1) {
191 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxs;
192 } else if (n == sf-1) {
193 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxf/2.0));
194 } else if (n == sf) {
195 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxs/2.0+info->hxf/2.0),jR[i]/info->hxf);
196 } else if (n > sf && n < fs-1) {
197 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxf;
198 } else if (n == fs-1) {
199 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxs/2.0));
200 } else if (n == fs) {
201 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxf/2.0+info->hxs/2.0),jR[i]/info->hxs);
202 } else {
203 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxs;
204 }
205 }
Limit2_Superbee(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)206 void Limit2_Superbee(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
207 {
208 PetscInt i;
209 if (n < sf-1) {
210 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxs;
211 } else if (n == sf-1) {
212 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/info->hxs,2*jR[i]/(info->hxs/2.0+info->hxf/2.0)),MinMod2(2*jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxf/2.0)));
213 } else if (n == sf) {
214 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/(info->hxs/2.0+info->hxf/2.0),2*jR[i]/info->hxf),MinMod2(2*jL[i]/(info->hxs/2.0+info->hxf/2.0),jR[i]/info->hxf));
215 } else if (n > sf && n < fs-1) {
216 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxf;
217 } else if (n == fs-1) {
218 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/info->hxf,2*jR[i]/(info->hxf/2.0+info->hxs/2.0)),MinMod2(2*jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxs/2.0)));
219 } else if (n == fs) {
220 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/(info->hxf/2.0+info->hxs/2.0),2*jR[i]/info->hxs),MinMod2(2*jL[i]/(info->hxf/2.0+info->hxs/2.0),jR[i]/info->hxs));
221 } else {
222 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxs;
223 }
224 }
Limit2_MC(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)225 void Limit2_MC(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
226 {
227 PetscInt i;
228 if (n < sf-1) {
229 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxs;
230 } else if (n == sf-1) {
231 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,0.5*(jL[i]/info->hxs+jR[i]/(info->hxf/2.0+info->hxs/2.0)),2*jR[i]/(info->hxf/2.0+info->hxs/2.0));
232 } else if (n == sf) {
233 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxf/2.0),0.5*(jL[i]/(info->hxs/2.0+info->hxf/2.0)+jR[i]/info->hxf),2*jR[i]/info->hxf);
234 } else if (n > sf && n < fs-1) {
235 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxf;
236 } else if (n == fs-1) {
237 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxs/2.0)),2*jR[i]/(info->hxf/2.0+info->hxs/2.0));
238 } else if (n == fs) {
239 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxs/2.0),0.5*(jL[i]/(info->hxf/2.0+info->hxs/2.0)+jR[i]/info->hxs),2*jR[i]/info->hxs);
240 } else {
241 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxs;
242 }
243 }
Limit2_Koren3(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)244 void Limit2_Koren3(LimitInfo info,const PetscScalar *jL,const PetscScalar *jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar *lmt)
245 { /* Eq 11 of Cada-Torrilhon 2009 */
246 PetscInt i;
247 if (n < sf-1) {
248 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs;
249 } else if (n == sf-1) {
250 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,(jL[i]/info->hxs+2*jR[i]/(info->hxf/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxf/2.0+info->hxs/2.0));
251 } else if (n == sf) {
252 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxf/2.0),(jL[i]/(info->hxs/2.0+info->hxf/2.0)+2*jR[i]/info->hxf)/3,2*jR[i]/info->hxf);
253 } else if (n > sf && n < fs-1) {
254 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxf;
255 } else if (n == fs-1) {
256 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,(jL[i]/info->hxf+2*jR[i]/(info->hxf/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxf/2.0+info->hxs/2.0));
257 } else if (n == fs) {
258 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxs/2.0),(jL[i]/(info->hxf/2.0+info->hxs/2.0)+2*jR[i]/info->hxs)/3,2*jR[i]/info->hxs);
259 } else {
260 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs;
261 }
262 }
263
264 /* ---- Three-way splitting ---- */
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)265 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)
266 {
267 PetscInt i;
268 for (i=0; i<info->m; i++) lmt[i] = 0;
269 }
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)270 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)
271 {
272 PetscInt i;
273 if (n < sm-1 || n > ms) { /* slow components */
274 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxs;
275 } else if (n == sm-1 || n == ms-1) { /* slow component which is next to medium components */
276 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxs/2.0+info->hxm/2.0);
277 } else if (n < mf-1 || n > fm) { /* medium components */
278 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxm;
279 } else if (n == mf-1 || n == fm) { /* medium component next to fast components */
280 for (i=0; i<info->m; i++) lmt[i] = jR[i]/(info->hxm/2.0+info->hxf/2.0);
281 } else { /* fast components */
282 for (i=0; i<info->m; i++) lmt[i] = jR[i]/info->hxf;
283 }
284 }
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)285 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)
286 {
287 PetscInt i;
288 if (n < sm || n > ms) {
289 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxs;
290 } else if (n == sm || n == ms) {
291 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxs/2.0+info->hxf/2.0);
292 }else if (n < mf || n > fm) {
293 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxm;
294 } else if (n == mf || n == fm) {
295 for (i=0; i<info->m; i++) lmt[i] = jL[i]/(info->hxm/2.0+info->hxf/2.0);
296 } else {
297 for (i=0; i<info->m; i++) lmt[i] = jL[i]/info->hxf;
298 }
299 }
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)300 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)
301 {
302 PetscInt i;
303 if (n < sm-1 || n > ms) {
304 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxs;
305 } else if (n == sm-1) {
306 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxs+jR[i]/(info->hxs/2.0+info->hxf/2.0));
307 } else if (n == sm) {
308 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxs/2.0+info->hxm/2.0)+jR[i]/info->hxm);
309 } else if (n == ms-1) {
310 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxm+jR[i]/(info->hxs/2.0+info->hxf/2.0));
311 } else if (n == ms) {
312 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxm/2.0+info->hxs/2.0)+jR[i]/info->hxs);
313 } else if (n < mf-1 || n > fm) {
314 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxm;
315 } else if (n == mf-1) {
316 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxm+jR[i]/(info->hxm/2.0+info->hxf/2.0));
317 } else if (n == mf) {
318 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxm/2.0+info->hxf/2.0)+jR[i]/info->hxf);
319 } else if (n == fm-1) {
320 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxm/2.0));
321 } else if (n == fm) {
322 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]/(info->hxf/2.0+info->hxm/2.0)+jR[i]/info->hxm);
323 } else {
324 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxf;
325 }
326 }
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)327 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)
328 {
329 PetscInt i;
330 if (n < sm-1 || n > ms) {
331 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxs;
332 } else if (n == sm-1) {
333 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxf/2.0));
334 } else if (n == sm) {
335 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxs/2.0+info->hxf/2.0),jR[i]/info->hxf);
336 } else if (n == ms-1) {
337 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxs/2.0));
338 } else if (n == ms) {
339 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxm/2.0+info->hxs/2.0),jR[i]/info->hxs);
340 } else if (n < mf-1 || n > fm) {
341 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i],jR[i])/info->hxm;
342 } else if (n == mf-1) {
343 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxf/2.0));
344 } else if (n == mf) {
345 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxm/2.0+info->hxf/2.0),jR[i]/info->hxf);
346 } else if (n == fm-1) {
347 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxm/2.0));
348 } else if (n == fm) {
349 for (i=0; i<info->m; i++) lmt[i] = MinMod2(jL[i]/(info->hxf/2.0+info->hxm/2.0),jR[i]/info->hxm);
350 } else {
351 for (i=0; i<info->m; i++) lmt[i] = 0.5*(jL[i]+jR[i])/info->hxf;
352 }
353 }
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)354 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)
355 {
356 PetscInt i;
357 if (n < sm-1 || n > ms) {
358 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxs;
359 } else if (n == sm-1) {
360 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/info->hxs,2*jR[i]/(info->hxs/2.0+info->hxm/2.0)),MinMod2(2*jL[i]/info->hxs,jR[i]/(info->hxs/2.0+info->hxm/2.0)));
361 } else if (n == sm) {
362 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/(info->hxs/2.0+info->hxm/2.0),2*jR[i]/info->hxm),MinMod2(2*jL[i]/(info->hxs/2.0+info->hxm/2.0),jR[i]/info->hxm));
363 } else if (n == ms-1) {
364 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/info->hxm,2*jR[i]/(info->hxm/2.0+info->hxs/2.0)),MinMod2(2*jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxs/2.0)));
365 } else if (n == ms) {
366 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/(info->hxm/2.0+info->hxs/2.0),2*jR[i]/info->hxs),MinMod2(2*jL[i]/(info->hxm/2.0+info->hxs/2.0),jR[i]/info->hxs));
367 } else if (n < mf-1 || n > fm) {
368 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxm;
369 } else if (n == mf-1) {
370 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/info->hxm,2*jR[i]/(info->hxm/2.0+info->hxf/2.0)),MinMod2(2*jL[i]/info->hxm,jR[i]/(info->hxm/2.0+info->hxf/2.0)));
371 } else if (n == mf) {
372 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i]/(info->hxm/2.0+info->hxf/2.0),2*jR[i]/info->hxf),MinMod2(2*jL[i]/(info->hxm/2.0+info->hxf/2.0),jR[i]/info->hxf));
373 } else if (n == fm-1) {
374 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/info->hxf,2*jR[i]/(info->hxf/2.0+info->hxm/2.0)),MinMod2(2*jL[i]/info->hxf,jR[i]/(info->hxf/2.0+info->hxm/2.0)));
375 } else if (n == fm) {
376 for (i=0; i<info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i]/(info->hxf/2.0+info->hxm/2.0),2*jR[i]/info->hxm),MinMod2(2*jL[i]/(info->hxf/2.0+info->hxm/2.0),jR[i]/info->hxm));
377 } else {
378 for (i=0; i<info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i],2*jR[i]),MinMod2(2*jL[i],jR[i]))/info->hxf;
379 }
380 }
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)381 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)
382 {
383 PetscInt i;
384 if (n < sm-1 || n > ms) {
385 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxs;
386 } else if (n == sm-1) {
387 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,0.5*(jL[i]/info->hxs+jR[i]/(info->hxs/2.0+info->hxm/2.0)),2*jR[i]/(info->hxs/2.0+info->hxm/2.0));
388 } else if (n == sm) {
389 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxm/2.0),0.5*(jL[i]/(info->hxs/2.0+info->hxm/2.0)+jR[i]/info->hxm),2*jR[i]/info->hxm);
390 } else if (n == ms-1) {
391 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,0.5*(jL[i]/info->hxm+jR[i]/(info->hxm/2.0+info->hxs/2.0)),2*jR[i]/(info->hxm/2.0+info->hxs/2.0));
392 } else if (n == ms) {
393 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxs/2.0),0.5*(jL[i]/(info->hxm/2.0+info->hxs/2.0)+jR[i]/info->hxs),2*jR[i]/info->hxs);
394 } else if (n < mf-1 || n > fm) {
395 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxm;
396 } else if (n == mf-1) {
397 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,0.5*(jL[i]/info->hxm+jR[i]/(info->hxm/2.0+info->hxf/2.0)),2*jR[i]/(info->hxm/2.0+info->hxf/2.0));
398 } else if (n == mf) {
399 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxf/2.0),0.5*(jL[i]/(info->hxm/2.0+info->hxf/2.0)+jR[i]/info->hxf),2*jR[i]/info->hxf);
400 } else if (n == fm-1) {
401 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,0.5*(jL[i]/info->hxf+jR[i]/(info->hxf/2.0+info->hxm/2.0)),2*jR[i]/(info->hxf/2.0+info->hxm/2.0));
402 } else if (n == fm) {
403 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxm/2.0),0.5*(jL[i]/(info->hxf/2.0+info->hxm/2.0)+jR[i]/info->hxm),2*jR[i]/info->hxm);
404 } else {
405 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],0.5*(jL[i]+jR[i]),2*jR[i])/info->hxf;
406 }
407 }
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)408 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)
409 { /* Eq 11 of Cada-Torrilhon 2009 */
410 PetscInt i;
411 if (n < sm-1 || n > ms) {
412 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs;
413 } else if (n == sm-1) {
414 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxs,(jL[i]/info->hxs+2*jR[i]/(info->hxm/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxm/2.0+info->hxs/2.0));
415 } else if (n == sm) {
416 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxs/2.0+info->hxm/2.0),(jL[i]/(info->hxs/2.0+info->hxm/2.0)+2*jR[i]/info->hxm)/3,2*jR[i]/info->hxm);
417 } else if (n == ms-1) {
418 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,(jL[i]/info->hxm+2*jR[i]/(info->hxm/2.0+info->hxs/2.0))/3,2*jR[i]/(info->hxm/2.0+info->hxs/2.0));
419 } else if (n == ms) {
420 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxs/2.0),(jL[i]/(info->hxm/2.0+info->hxs/2.0)+2*jR[i]/info->hxs)/3,2*jR[i]/info->hxs);
421 } else if (n < mf-1 || n > fm) {
422 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxm;
423 } else if (n == mf-1) {
424 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxm,(jL[i]/info->hxm+2*jR[i]/(info->hxm/2.0+info->hxf/2.0))/3,2*jR[i]/(info->hxm/2.0+info->hxf/2.0));
425 } else if (n == mf) {
426 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxm/2.0+info->hxf/2.0),(jL[i]/(info->hxm/2.0+info->hxf/2.0)+2*jR[i]/info->hxf)/3,2*jR[i]/info->hxf);
427 } else if (n == fm-1) {
428 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/info->hxf,(jL[i]/info->hxf+2*jR[i]/(info->hxf/2.0+info->hxm/2.0))/3,2*jR[i]/(info->hxf/2.0+info->hxm/2.0));
429 } else if (n == fm) {
430 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i]/(info->hxf/2.0+info->hxm/2.0),(jL[i]/(info->hxf/2.0+info->hxm/2.0)+2*jR[i]/info->hxm)/3,2*jR[i]/info->hxm);
431 } else {
432 for (i=0; i<info->m; i++) lmt[i] = MinMod3(2*jL[i],(jL[i]+2*jR[i])/3,2*jR[i])/info->hxs;
433 }
434 }
RiemannListAdd(PetscFunctionList * flist,const char * name,RiemannFunction rsolve)435 PetscErrorCode RiemannListAdd(PetscFunctionList *flist,const char *name,RiemannFunction rsolve)
436 {
437 PetscErrorCode ierr;
438
439 PetscFunctionBeginUser;
440 ierr = PetscFunctionListAdd(flist,name,rsolve);CHKERRQ(ierr);
441 PetscFunctionReturn(0);
442 }
443
RiemannListFind(PetscFunctionList flist,const char * name,RiemannFunction * rsolve)444 PetscErrorCode RiemannListFind(PetscFunctionList flist,const char *name,RiemannFunction *rsolve)
445 {
446 PetscErrorCode ierr;
447
448 PetscFunctionBeginUser;
449 ierr = PetscFunctionListFind(flist,name,rsolve);CHKERRQ(ierr);
450 if (!*rsolve) SETERRQ1(PETSC_COMM_SELF,1,"Riemann solver \"%s\" could not be found",name);
451 PetscFunctionReturn(0);
452 }
453
ReconstructListAdd(PetscFunctionList * flist,const char * name,ReconstructFunction r)454 PetscErrorCode ReconstructListAdd(PetscFunctionList *flist,const char *name,ReconstructFunction r)
455 {
456 PetscErrorCode ierr;
457
458 PetscFunctionBeginUser;
459 ierr = PetscFunctionListAdd(flist,name,r);CHKERRQ(ierr);
460 PetscFunctionReturn(0);
461 }
462
ReconstructListFind(PetscFunctionList flist,const char * name,ReconstructFunction * r)463 PetscErrorCode ReconstructListFind(PetscFunctionList flist,const char *name,ReconstructFunction *r)
464 {
465 PetscErrorCode ierr;
466
467 PetscFunctionBeginUser;
468 ierr = PetscFunctionListFind(flist,name,r);CHKERRQ(ierr);
469 if (!*r) SETERRQ1(PETSC_COMM_SELF,1,"Reconstruction \"%s\" could not be found",name);
470 PetscFunctionReturn(0);
471 }
472
473 /* --------------------------------- Physics ----------------------------------- */
474
PhysicsDestroy_SimpleFree(void * vctx)475 PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
476 {
477 PetscErrorCode ierr;
478
479 PetscFunctionBeginUser;
480 ierr = PetscFree(vctx);CHKERRQ(ierr);
481 PetscFunctionReturn(0);
482 }
483
484 /* --------------------------------- Finite Volume Solver ----------------------------------- */
485
FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void * vctx)486 PetscErrorCode FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void *vctx)
487 {
488 FVCtx *ctx = (FVCtx*)vctx;
489 PetscErrorCode ierr;
490 PetscInt i,j,k,Mx,dof,xs,xm;
491 PetscReal hx,cfl_idt = 0;
492 PetscScalar *x,*f,*slope;
493 Vec Xloc;
494 DM da;
495
496 PetscFunctionBeginUser;
497 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
498 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr); /* Xloc contains ghost points */
499 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr); /* Mx is the number of center points */
500 hx = (ctx->xmax-ctx->xmin)/Mx;
501 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr); /* X is solution vector which does not contain ghost points */
502 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
503 ierr = VecZeroEntries(F);CHKERRQ(ierr); /* F is the right hand side function corresponds to center points */
504 ierr = DMDAVecGetArray(da,Xloc,&x);CHKERRQ(ierr);
505 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
506 ierr = DMDAGetArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr); /* contains ghost points */
507 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
508
509 if (ctx->bctype == FVBC_OUTFLOW) {
510 for (i=xs-2; i<0; i++) {
511 for (j=0; j<dof; j++) x[i*dof+j] = x[j];
512 }
513 for (i=Mx; i<xs+xm+2; i++) {
514 for (j=0; j<dof; j++) x[i*dof+j] = x[(xs+xm-1)*dof+j];
515 }
516 }
517 for (i=xs-1; i<xs+xm+1; i++) {
518 struct _LimitInfo info;
519 PetscScalar *cjmpL,*cjmpR;
520 /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
521 ierr = (*ctx->physics.characteristic)(ctx->physics.user,dof,&x[i*dof],ctx->R,ctx->Rinv,ctx->speeds,ctx->xmin+hx*i);CHKERRQ(ierr);
522 /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
523 ierr = PetscArrayzero(ctx->cjmpLR,2*dof);CHKERRQ(ierr);
524 cjmpL = &ctx->cjmpLR[0];
525 cjmpR = &ctx->cjmpLR[dof];
526 for (j=0; j<dof; j++) {
527 PetscScalar jmpL,jmpR;
528 jmpL = x[(i+0)*dof+j]-x[(i-1)*dof+j];
529 jmpR = x[(i+1)*dof+j]-x[(i+0)*dof+j];
530 for (k=0; k<dof; k++) {
531 cjmpL[k] += ctx->Rinv[k+j*dof]*jmpL;
532 cjmpR[k] += ctx->Rinv[k+j*dof]*jmpR;
533 }
534 }
535 /* Apply limiter to the left and right characteristic jumps */
536 info.m = dof;
537 info.hx = hx;
538 (*ctx->limit)(&info,cjmpL,cjmpR,ctx->cslope);
539 for (j=0; j<dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
540 for (j=0; j<dof; j++) {
541 PetscScalar tmp = 0;
542 for (k=0; k<dof; k++) tmp += ctx->R[j+k*dof]*ctx->cslope[k];
543 slope[i*dof+j] = tmp;
544 }
545 }
546
547 for (i=xs; i<xs+xm+1; i++) {
548 PetscReal maxspeed;
549 PetscScalar *uL,*uR;
550 uL = &ctx->uLR[0];
551 uR = &ctx->uLR[dof];
552 for (j=0; j<dof; j++) {
553 uL[j] = x[(i-1)*dof+j]+slope[(i-1)*dof+j]*hx/2;
554 uR[j] = x[(i-0)*dof+j]-slope[(i-0)*dof+j]*hx/2;
555 }
556 ierr = (*ctx->physics.riemann)(ctx->physics.user,dof,uL,uR,ctx->flux,&maxspeed,ctx->xmin+hx*i,ctx->xmin,ctx->xmax);CHKERRQ(ierr);
557 cfl_idt = PetscMax(cfl_idt,PetscAbsScalar(maxspeed/hx)); /* Max allowable value of 1/Delta t */
558 if (i > xs) {
559 for (j=0; j<dof; j++) f[(i-1)*dof+j] -= ctx->flux[j]/hx;
560 }
561 if (i < xs+xm) {
562 for (j=0; j<dof; j++) f[i*dof+j] += ctx->flux[j]/hx;
563 }
564 }
565 ierr = DMDAVecRestoreArray(da,Xloc,&x);CHKERRQ(ierr);
566 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
567 ierr = DMDARestoreArray(da,PETSC_TRUE,&slope);CHKERRQ(ierr);
568 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
569 ierr = MPI_Allreduce(&cfl_idt,&ctx->cfl_idt,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);
570 if (0) {
571 /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
572 PetscReal dt,tnow;
573 ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
574 ierr = TSGetTime(ts,&tnow);CHKERRQ(ierr);
575 if (dt > 0.5/ctx->cfl_idt) {
576 if (1) {
577 ierr = PetscPrintf(ctx->comm,"Stability constraint exceeded at t=%g, dt %g > %g\n",(double)tnow,(double)dt,(double)(0.5/ctx->cfl_idt));CHKERRQ(ierr);
578 } else SETERRQ2(PETSC_COMM_SELF,1,"Stability constraint exceeded, %g > %g",(double)dt,(double)(ctx->cfl/ctx->cfl_idt));
579 }
580 }
581 PetscFunctionReturn(0);
582 }
583
FVSample(FVCtx * ctx,DM da,PetscReal time,Vec U)584 PetscErrorCode FVSample(FVCtx *ctx,DM da,PetscReal time,Vec U)
585 {
586 PetscErrorCode ierr;
587 PetscScalar *u,*uj;
588 PetscInt i,j,k,dof,xs,xm,Mx;
589
590 PetscFunctionBeginUser;
591 if (!ctx->physics.sample) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Physics has not provided a sampling function");
592 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
593 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
594 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
595 ierr = PetscMalloc1(dof,&uj);CHKERRQ(ierr);
596 for (i=xs; i<xs+xm; i++) {
597 const PetscReal h = (ctx->xmax-ctx->xmin)/Mx,xi = ctx->xmin+h/2+i*h;
598 const PetscInt N = 200;
599 /* Integrate over cell i using trapezoid rule with N points. */
600 for (k=0; k<dof; k++) u[i*dof+k] = 0;
601 for (j=0; j<N+1; j++) {
602 PetscScalar xj = xi+h*(j-N/2)/(PetscReal)N;
603 ierr = (*ctx->physics.sample)(ctx->physics.user,ctx->initial,ctx->bctype,ctx->xmin,ctx->xmax,time,xj,uj);CHKERRQ(ierr);
604 for (k=0; k<dof; k++) u[i*dof+k] += ((j==0 || j==N) ? 0.5 : 1.0)*uj[k]/N;
605 }
606 }
607 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
608 ierr = PetscFree(uj);CHKERRQ(ierr);
609 PetscFunctionReturn(0);
610 }
611
SolutionStatsView(DM da,Vec X,PetscViewer viewer)612 PetscErrorCode SolutionStatsView(DM da,Vec X,PetscViewer viewer)
613 {
614 PetscErrorCode ierr;
615 PetscReal xmin,xmax;
616 PetscScalar sum,tvsum,tvgsum;
617 const PetscScalar *x;
618 PetscInt imin,imax,Mx,i,j,xs,xm,dof;
619 Vec Xloc;
620 PetscBool iascii;
621
622 PetscFunctionBeginUser;
623 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
624 if (iascii) {
625 /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
626 ierr = DMGetLocalVector(da,&Xloc);CHKERRQ(ierr);
627 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
628 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);CHKERRQ(ierr);
629 ierr = DMDAVecGetArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr);
630 ierr = DMDAGetCorners(da,&xs,0,0,&xm,0,0);CHKERRQ(ierr);
631 ierr = DMDAGetInfo(da,0, &Mx,0,0, 0,0,0, &dof,0,0,0,0,0);CHKERRQ(ierr);
632 tvsum = 0;
633 for (i=xs; i<xs+xm; i++) {
634 for (j=0; j<dof; j++) tvsum += PetscAbsScalar(x[i*dof+j]-x[(i-1)*dof+j]);
635 }
636 ierr = MPI_Allreduce(&tvsum,&tvgsum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)da));CHKERRQ(ierr);
637 ierr = DMDAVecRestoreArrayRead(da,Xloc,(void*)&x);CHKERRQ(ierr);
638 ierr = DMRestoreLocalVector(da,&Xloc);CHKERRQ(ierr);
639 ierr = VecMin(X,&imin,&xmin);CHKERRQ(ierr);
640 ierr = VecMax(X,&imax,&xmax);CHKERRQ(ierr);
641 ierr = VecSum(X,&sum);CHKERRQ(ierr);
642 ierr = PetscViewerASCIIPrintf(viewer,"Solution range [%8.5f,%8.5f] with minimum at %D, mean %8.5f, ||x||_TV %8.5f\n",(double)xmin,(double)xmax,imin,(double)(sum/Mx),(double)(tvgsum/Mx));CHKERRQ(ierr);
643 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type not supported");
644 PetscFunctionReturn(0);
645 }
646