#include "poisson.h"
struct Viscosity {
vector u;
face vector mu;
scalar rho;
double dt;
int nrelax;
#if COMPRESSIBLE
face vector lambdav;
#endif
scalar * res;
};
#if AXI
# define lambda ((coord){1., 1. + dt/rho[]*(mu.x[] + mu.x[1] + \
mu.y[] + mu.y[0,1])/2./sq(y)})
#if COMPRESSIBLE
# define gamma ((coord) {1., 0})
# define gamma1 ((coord) {0., \
- (lambdav.x[] + lambdav.x[1] + \
lambdav.y[] + lambdav.y[0,1])/4./sq(y)})
#endif
#else // not AXI
# if dimension == 1
# define lambda ((coord){1.})
# elif dimension == 2
# define lambda ((coord){1.,1.})
# elif dimension == 3
# define lambda ((coord){1.,1.,1.})
#endif
#endif
static void relax_viscosity (scalar * a, scalar * b, int l, void * data)
{
struct Viscosity * p = (struct Viscosity *) data;
(const) face vector mu = p->mu;
(const) scalar rho = p->rho;
double dt = p->dt;
vector u = vector(a[0]), r = vector(b[0]);
#if COMPRESSIBLE
(const) face vector lambdav = p->lambdav;
#if AXI
scalar ur = u.y;
#endif
#endif
#if JACOBI
vector w[];
#else
vector w = u;
#endif
foreach_level_or_leaf (l) {
foreach_dimension() {
w.x[] = (dt/rho[]*(2.*mu.x[1]*u.x[1] + 2.*mu.x[]*u.x[-1]
#if COMPRESSIBLE
+ (lambdav.x[1]*u.x[1] + lambdav.x[]*u.x[-1]
#if dimension > 1
+ lambdav.x[1,0]*((u.y[1,1] + u.y[0,1])/4 -
(u.y[1,-1] + u.y[0,-1])/4.)
- lambdav.x[]*((u.y[0,1] + u.y[-1,1])/4 -
(u.y[0,-1] + u.y[-1,-1])/4.)
#endif
#if dimension > 2
+ lambdav.x[1,0,0]*((u.z[1,0,1] + u.z[0,0,1])/4 -
(u.z[1,0,-1] + u.z[0,0,-1])/4.)
- lambdav.x[]*((u.z[0,0,1] + u.z[-1,0,1])/4 -
(u.z[0,0,-1] + u.z[-1,0,-1])/4.)
#endif
#if AXI
+ (lambdav.x[1]*(ur[1] + gamma.x*ur[]) -
lambdav.x[]*(ur[-1] + gamma.x*ur[]))/2.*Delta/y)*y
#else
)
#endif
#endif
#if dimension > 1
+ mu.y[0,1]*(u.x[0,1] +
(u.y[1,0] + u.y[1,1])/4. -
(u.y[-1,0] + u.y[-1,1])/4.)
- mu.y[]*(- u.x[0,-1] +
(u.y[1,-1] + u.y[1,0])/4. -
(u.y[-1,-1] + u.y[-1,0])/4.)
#endif
#if dimension > 2
+ mu.z[0,0,1]*(u.x[0,0,1] +
(u.z[1,0,0] + u.z[1,0,1])/4. -
(u.z[-1,0,0] + u.z[-1,0,1])/4.)
- mu.z[]*(- u.x[0,0,-1] +
(u.z[1,0,-1] + u.z[1,0,0])/4. -
(u.z[-1,0,-1] + u.z[-1,0,0])/4.)
#endif
) + r.x[]*sq(Delta))/
(sq(Delta)*lambda.x + dt/rho[]*(2.*mu.x[1] + 2.*mu.x[]
#if COMPRESSIBLE
+ (lambdav.x[1] + lambdav.x[]
#if AXI
- (lambdav.x[1] - lambdav.x[])*Delta/2/y*gamma.y
- gamma1.x*sq(Delta))*y
#else
)
#endif
#endif
#if dimension > 1
+ mu.y[0,1] + mu.y[]
#endif
#if dimension > 2
+ mu.z[0,0,1] + mu.z[]
#endif
));
}
}
#if JACOBI
foreach_level_or_leaf (l)
foreach_dimension()
u.x[] = (u.x[] + 2.*w.x[])/3.;
#endif
#if TRASH
vector u1[];
foreach_level_or_leaf (l)
foreach_dimension()
u1.x[] = u.x[];
trash ({u});
foreach_level_or_leaf (l)
foreach_dimension()
u.x[] = u1.x[];
#endif
}
static double residual_viscosity (scalar * a, scalar * b, scalar * resl,
void * data)
{
struct Viscosity * p = (struct Viscosity *) data;
(const) face vector mu = p->mu;
(const) scalar rho = p->rho;
double dt = p->dt;
vector u = vector(a[0]), r = vector(b[0]), res = vector(resl[0]);
double maxres = 0.;
#if COMPRESSIBLE
(const) face vector lambdav = p->lambdav;
#if AXI
scalar ur = u.y;
#endif
#endif
#if TREE
/* conservative coarse/fine discretisation (2nd order) */
foreach_dimension() {
face vector taux[];
#if COMPRESSIBLE
face vector tauc[];
#if AXI
face vector axic[];
#endif
#endif
foreach_face(x) {
#if COMPRESSIBLE
#if AXI
axic.x[] = lambdav.x[]*(ur[]+ur[-1])/2.;
#endif
tauc.x[] = lambdav.x[]*(u.x[] - u.x[-1]
#if dimension > 1
+ (u.y[0,1] + u.y[-1,1])/4
- (u.y[0,-1] + u.y[-1,-1])/4.
#endif
#if dimension > 2
+ (u.z[0,0,1] + u.z[-1,0,1])/4
- (u.z[0,0,-1] + u.z[-1,0,-1])/4.
#endif
)/Delta;
#endif
taux.x[] = 2.*mu.x[]*(u.x[] - u.x[-1])/Delta;
}
#if dimension > 1
foreach_face(y)
taux.y[] = mu.y[]*(u.x[] - u.x[0,-1] +
(u.y[1,-1] + u.y[1,0])/4. -
(u.y[-1,-1] + u.y[-1,0])/4.)/Delta;
#endif
#if dimension > 2
foreach_face(z)
taux.z[] = mu.z[]*(u.x[] - u.x[0,0,-1] +
(u.z[1,0,-1] + u.z[1,0,0])/4. -
(u.z[-1,0,-1] + u.z[-1,0,0])/4.)/Delta;
#endif
boundary_flux ({taux});
foreach (reduction(max:maxres)) {
double d = 0.;
foreach_dimension()
d += taux.x[1] - taux.x[];
res.x[] = r.x[] - lambda.x*u.x[] + dt/rho[]*d/Delta;
#if COMPRESSIBLE
#if AXI
res.x[] += dt/rho[]*((tauc.x[1] - tauc.x[] +
(axic.x[1] - axic.x[])/y)/Delta
+ gamma1.x*ur[])*y;
#else
res.x[] += dt/rho[]*(tauc.x[1]-tauc.x[])/Delta;
#endif
#endif
if (fabs (res.x[]) > maxres)
maxres = fabs (res.x[]);
}
}
boundary (resl);
#else
/* "naive" discretisation (only 1st order on trees) */
foreach (reduction(max:maxres))
foreach_dimension() {
res.x[] = r.x[] - lambda.x*u.x[] +
dt/rho[]*(2.*mu.x[1,0]*(u.x[1] - u.x[])
- 2.*mu.x[]*(u.x[] - u.x[-1])
#if COMPRESSIBLE
+ (lambdav.x[1]*(u.x[1] - u.x[])
- lambdav.x[]*(u.x[] - u.x[-1])
#if dimension > 1
+ lambdav.x[1]*((u.y[1,1] + u.y[0,1])/4 -
(u.y[1,-1] + u.y[0,-1])/4.)
- lambdav.x[]*((u.y[0,1] + u.y[-1,1])/4 -
(u.y[0,-1] + u.y[-1,-1])/4.)
#endif
#if dimension > 2
+ lambdav.x[1]*((u.z[1,0,1] + u.z[0,0,1])/4 -
(u.z[1,0,-1] + u.z[0,0,-1])/4.)
- lambdav.x[]*((u.z[0,0,1] + u.z[-1,0,1])/4 -
(u.z[0,0,-1] + u.z[-1,0,-1])/4.)
#endif
#if AXI
+ (lambdav.x[1]*(ur[1] + ur[])/2.
- lambdav.x[]*(ur[-1] + ur[])/2.)*Delta/y
+ gamma1.x*ur[]*sq(Delta))*y
#else
)
#endif
#endif
#if dimension > 1
+ mu.y[0,1]*(u.x[0,1] - u.x[] +
(u.y[1,0] + u.y[1,1])/4. -
(u.y[-1,0] + u.y[-1,1])/4.)
- mu.y[]*(u.x[] - u.x[0,-1] +
(u.y[1,-1] + u.y[1,0])/4. -
(u.y[-1,-1] + u.y[-1,0])/4.)
#endif
#if dimension > 2
+ mu.z[0,0,1]*(u.x[0,0,1] - u.x[] +
(u.z[1,0,0] + u.z[1,0,1])/4. -
(u.z[-1,0,0] + u.z[-1,0,1])/4.)
- mu.z[]*(u.x[] - u.x[0,0,-1] +
(u.z[1,0,-1] + u.z[1,0,0])/4. -
(u.z[-1,0,-1] + u.z[-1,0,0])/4.)
#endif
)/sq(Delta);
if (fabs (res.x[]) > maxres)
maxres = fabs (res.x[]);
}
#endif
return maxres;
}
#if COMPRESSIBLE & AXI
#undef gamma
#undef gamma1
#endif
#undef lambda
mgstats viscosity (struct Viscosity p)
{
vector u = p.u, r[];
foreach()
foreach_dimension()
r.x[] = u.x[];
face vector mu = p.mu;
scalar rho = p.rho;
#if COMPRESSIBLE
face vector lambdav = p.lambdav;
restriction ({mu, lambdav, rho});
#else
restriction ({mu,rho});
#endif
return mg_solve ((scalar *){u}, (scalar *){r},
residual_viscosity, relax_viscosity, &p, p.nrelax, p.res);
}
mgstats viscosity_explicit (struct Viscosity p)
{
vector u = p.u, r[];
mgstats mg = {0};
mg.resb = residual_viscosity ((scalar *){u}, (scalar *){u}, (scalar *){r}, &p);
foreach()
foreach_dimension()
u.x[] += r.x[];
boundary ((scalar *){u});
return mg;
}
