14 double f, d2, d2p, d2m;
16 double scrh1,scrh2, fmin, fmax;
17 double fAV, fMD, fLC, fUL, fMP;
18 static double alpha = 4.0, epsm = 1.e-12;
20 f = 2.0*F[j-2] - 13.0*F[j-1] + 47.0*F[
j] + 27.0*F[j+1] - 3.0*F[j+2];
23 fMP = F[
j] +
MINMOD(F[j+1]-F[j],alpha*(F[j]-F[j-1]));
25 if ((f - F[j])*(f - fMP) <= epsm)
return f;
27 d2m = F[j-2] + F[
j ] - 2.0*F[j-1];
28 d2 = F[j-1] + F[j+1] - 2.0*F[
j];
29 d2p = F[
j ] + F[j+2] - 2.0*F[j+1];
31 scrh1 =
MINMOD(4.0*d2 - d2p, 4.0*d2p - d2);
33 dMMp =
MINMOD(scrh1,scrh2);
35 scrh1 =
MINMOD(4.0*d2m - d2, 4.0*d2 - d2m);
37 dMMm =
MINMOD(scrh1,scrh2);
39 fUL = F[
j] + alpha*(F[
j] - F[j-1]);
40 fAV = 0.5*(F[
j] + F[j+1]);
42 fLC = 0.5*(3.0*F[
j] - F[j-1]) + 4.0/3.0*dMMm;
44 scrh1 =
MIN(F[j], F[j+1]); scrh1 =
MIN(scrh1, fMD);
45 scrh2 =
MIN(F[j], fUL); scrh2 =
MIN(scrh2, fLC);
46 fmin =
MAX(scrh1, scrh2);
48 scrh1 =
MAX(F[j], F[j+1]); scrh1 =
MAX(scrh1, fMD);
49 scrh2 =
MAX(F[j], fUL); scrh2 =
MAX(scrh2, fLC);
50 fmax =
MIN(scrh1, scrh2);
64 double dvm2, dvm1, dvp1, dvp2;
65 double dvc, Sm1, Sp1, Sp2, SM;
66 double ap, am, vp, vm;
68 dvm2 = v[i-1] - v[i-2];
71 dvp2 = v[i+2] - v[i+1];
73 dvc = 0.5*(dvm2 + dvm1);
74 SM = 2.0*
MINMOD(dvm2, dvm1);
77 dvc = 0.5*(dvm1 + dvp1);
78 SM = 2.0*
MINMOD(dvm1, dvp1);
81 dvc = 0.5*(dvp2 + dvp1);
82 SM = 2.0*
MINMOD(dvp2, dvp1);
85 vp = 0.5*(v[
i] + v[i+1]) - (Sp2 - Sp1)/6.0;
86 vm = 0.5*(v[
i] + v[i-1]) - (Sp1 - Sm1)/6.0;
94 if (fabs(ap) >= 2.0*fabs(am)) ap = -2.0*am;
95 if (fabs(am) >= 2.0*fabs(ap)) am = -2.0*ap;
110 double dvp, dvm, r = 1.;
111 double a,b,
c,
q, th, lim;
112 double eta, psi,
eps = 1.e-12;
117 th = dvm/(dvp + 1.e-16);
127 eta = (dvm*dvm + dvp*dvp)/(eta*eta);
128 if ( eta <= 1.0 - eps) {
130 }
else if (eta >= 1.0 + eps){
133 psi = (1.0 - (eta - 1.0)/
eps)*q
134 + (1.0 + (eta - 1.0)/
eps)*psi;
137 return (v[i] + 0.5*dvp*lim);
152 double IS0, IS1, IS2, IS3, IS4, IS5;
157 static double thirteen_12 = 13.0/12.0;
158 double b0, b1, b2,t5;
160 a0 = F[j-2] - 2.0*F[j-1] + F[
j];
161 a1 = F[j-2] - 4.0*F[j-1] + 3.0*F[
j];
162 b0 = thirteen_12*a0*a0 + 0.25*a1*a1;
164 a0 = F[j-1] - 2.0*F[
j] + F[j+1];
165 a1 = F[j-1] - F[j+1];
166 b1 = thirteen_12*a0*a0 + 0.25*a1*a1;
168 a0 = F[
j] - 2.0*F[j+1] + F[j+2];
169 a1 = 3.0*F[
j] - 4.0*F[j+1] + F[j+2];
170 b2 = thirteen_12*a0*a0 + 0.25*a1*a1;
174 a0 = 1.0*(1.0 + t5/(b0 + 1.e-40));
175 a1 = 6.0*(1.0 + t5/(b1 + 1.e-40));
176 a2 = 3.0*(1.0 + t5/(b2 + 1.e-40));
178 sum_a = 1.0/(a0 + a1 + a2);
183 f0 = (2.0*F[j-2] - 7.0*F[j-1] + 11.0*F[
j]);
184 f1 = ( -F[j-1] + 5.0*F[
j] + 2.0*F[j+1]);
185 f2 = (2.0*F[
j] + 5.0*F[j+1] - F[j+2]);
187 f = (w0*f0 + w1*f1 + w2*f2)/6.0;
203 double a0, b0,
w0, f0, a1, b1, w1,
f1;
225 tau = (F[j+1] - 2.0*F[
j] + F[j-1]);
228 a0 = 2.0*(1.0 + tau/(dx2 + b0));
229 a1 = 1.0*(1.0 + tau/(dx2 + b1));
234 f0 = 0.5*( F[j+1] + F[
j]);
235 f1 = 0.5*(-F[j-1] + 3.0*F[
j]);
236 return (w0*f0 + w1*f1);
259 return (F[j] + 0.5*dF);
271 return (a +
MINMOD(b-a,c-a));
double WENO3_Reconstruct(double *F, double dx, int j)
double WENOZ_Reconstruct(double *F, double dx, int j)
double MP5_Reconstruct(double *F, double dx, int j)
double Median(double a, double b, double c)
double PPM_Reconstruct(double *v, double dx, int i)
double LIMO3_Reconstruct(double *v, double dx, int i)
double LIN_Reconstruct(double *F, double dx, int j)