/** Here we try to obtain a flow driven by a pressure drop on a tri-periodic cubic domain.
Solution: since the pressure drop cannot be imposed in a periodic domain, the solution is to replace it with an equivalent body force i.e. we have in the r.h.s. of the momentum equation
$$
-\partial_x p + \rho a_x
$$
so that adding a constant pressure gradient is equivalent to adding an acceleration
$$
a_x = - 1/\rho\partial_x p
$$
*/
#define LEVEL 6
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#define rhoval 1000.
#define muval 60.
#define mydt 0.00375
#define maxiter (20)
double deltau;
scalar un[];
int main() {
L0 = 1;
DT = mydt;
init_grid(1 << (LEVEL));
/* Tri-periodic flow driven by body force */
foreach_dimension()
periodic(top);
/*
p[left] = dirichlet(0);
p[right] = dirichlet(-10); */
const face vector dp[] = {10./L0/rhoval, 0, 0};
a = dp;
run();
}
event init (i = 0) {
origin (0, 0, 0);
/* set dynamic viscosity */
const face vector muc[] = {muval, muval, muval};
mu = muc;
/* set density of the flow */
const scalar rhoc[] = rhoval;
rho = rhoc;
/* The flow is at rest initially. */
foreach(){
u.x[] = 0.;
un[] = u.x[];
}
}
/* Logging the variable "deltau" that measures the difference between the current velocity to the initial one. */
event logfile (i++;i