// #include "grid/multigrid3D.h"
#include "grid/octree.h"
#include "navier-stokes/centered.h"
#include "vof.h"
#include "tension.h"
scalar f[], * interfaces = {f};
#define rho1 900.
#define rho2 1000.
#define mu1 0.1
#define mu2 0.001
#define rho(f) (clamp(f,0,1)*(rho1 - rho2) + rho2)
#define mu(f) (1./(clamp(f,0,1)*(1./(mu1) - 1./(mu2)) + 1./(mu2)))
#define H 0.1
#define G 9.81
#define Ug ((rho2 - rho1)/rho1*sqrt(H*G/2.))
#define tc (H/(2.*Ug))
face vector alphav[], muv[], av[];
scalar rhov[];
int maxlevel = 6;
#if 0
uf.n[left] = 0.;
uf.n[right] = 0.;
p[left] = neumann(0);
p[right] = neumann(0);
f[top] = 0.;
f[right] = 0.;
f[left] = 0.;
f[bottom] = 0.;
#if dimension == 3
f[front] = 0.;
f[back] = 0.;
#endif
#endif
timer tt;
int main (int argc, char * argv[]) {
maxlevel = argc > 1 ? atoi(argv[1]) : 7;
size (H);
origin (-H/2., -H/2., -H/2.);
#if !TREE
N = 1 << maxlevel;
#endif
a = av;
mu = muv;
alpha = alphav;
rho = rhov;
f.sigma = 0.045;
DT = 2e-2;
tt = timer_start();
run();
}
event init (i = 0) {
#if TREE
scalar f1[];
foreach()
f1[] = (x <= 0 && y <= 0 && z <= 0);
astats s;
do {
s = adapt_wavelet ({f1}, (double[]){0.0}, maxlevel, list = NULL);
foreach()
f1[] = (x <= 0 && y <= 0 && z <= 0);
} while (s.nf);
foreach()
f[] = (x <= 0 && y <= 0 && z <= 0);
#else
foreach()
f[] = (x <= 0 && y <= 0 && z <= 0);
#endif
}
event acceleration (i++) {
foreach_face(y)
av.y[] -= G;
}
event properties (i++) {
#if TREE
f.prolongation = refine_bilinear;
f.dirty = true;
#endif
foreach_face() {
double ff = (f[] + f[-1])/2.;
alphav.x[] = fm.x[]/rho(ff);
muv.x[] = fm.x[]*mu(ff);
}
foreach()
rhov[] = cm[]*rho(f[]);
#if TREE
f.prolongation = fraction_refine;
f.dirty = true;
#endif
}
event logfile (i++; t <= 20) {
double ke1 = 0., ke2 = 0., vd = 0., vol1 = 0.;
double ep1 = 0., ep2 = 0.;
double er1 = 0., er2 = 0.;
double area = 0.;
int nc = 0;
static long tnc = 0;
foreach(reduction(+:ke1) reduction(+:ke2) reduction(+:vd)
reduction(+:vol1) reduction(+:ep1) reduction(+:ep2)
reduction(+:er1) reduction(+:er2) reduction(+:area)
reduction(+:nc)) {
if (y > H/2. - H/8.)
vol1 += f[]*dv();
ep1 += rho1*f[]*G*(y + H/2.)*dv();
ep2 += rho2*(1. - f[])*G*(y + H/2.)*dv();
// interfacial area
if (f[] > 1e-4 && f[] < 1. - 1e-4) {
coord m = mycs (point, f);
double alpha = plane_alpha (f[], m);
coord p;
area += sq(Delta)*plane_area_center (m, alpha, &p);
}
double w2 = 0.;
foreach_dimension() {
// kinetic energy
ke1 += dv()*f[]*rho1*sq(u.x[]);
ke2 += dv()*(1. - f[])*rho2*sq(u.x[]);
// viscous dissipation
vd += dv()*(sq(u.x[1] - u.x[-1]) +
sq(u.x[0,1] - u.x[0,-1]) +
sq(u.x[0,0,1] - u.x[0,0,-1]))/sq(2.*Delta);
// enstrophy
w2 += sq(u.x[0,1] - u.x[0,-1] - u.y[1,0] + u.y[-1,0]);
}
w2 /= sq(2.*Delta);
er1 += dv()*f[]*w2;
er2 += dv()*(1. - f[])*w2;
nc++;
}
ke1 /= 2.;
ke2 /= 2.;
er1 /= 2.;
er2 /= 2.;
// vd *= MU/vol;
if (i == 0)
fprintf (stderr,
"t ke1 ke2 ep1 ep2 er1 er2 R2 area mgp.i mgu.i nc time speed\n");
double elapsed = timer_elapsed (tt);
tnc += nc;
fprintf (stderr, "%g %g %g %g %g %g %g %g %g %d %d %d %g %g\n",
t/tc,
ke1/(1./16.*rho1*sq(Ug)*cube(H)),
ke2/(1./16.*rho2*sq(Ug)*cube(H)),
ep1/(rho1*G*15.*sq(H)*sq(H)/128.),
ep2/(rho2*G*49.*sq(H)*sq(H)/128.),
er1/0.0733,
er2/1.3759,
8.*vol1/cube(H),
area,
mgp.i, mgu.i, nc, elapsed, tnc/elapsed);
#if 0
nc = 0;
foreach()
nc++;
fprintf (stderr, "nc: %d\n", nc);
fflush (stderr);
#endif
}
#if 0
event movies (t += 0.1*tc) {
char name[80];
sprintf (name, "f-%d.ppm", maxlevel);
static FILE * fp = fopen (name, "w");
output_ppm (f, fp, min = 0, max = 1, n = 256);
}
#endif
#if !_MPI
event gfsview (i += 10) {
scalar pid[];
foreach()
pid[] = tid();
#if dimension == 3
static FILE * fp = popen ("gfsview3D inversion.gfv", "w");
#else
static FILE * fp = popen ("gfsview2D inversion2D.gfv", "w");
#endif
output_gfs (fp);
}
#endif
#if 0 //TREE && !_MPI
event gfsview (t += 0.1*tc) {
@if _MPI
char name[80];
sprintf (name, "output-%g.gfs", t);
FILE * fp = fopen (name, "w");
@else
static FILE * fp =
popen ("gfsview3D ../inversion.gfv", "w");
@endif
output_gfs (fp, translate = true);
@if _MPI
fclose (fp);
@endif
// fprintf (fp, "Save stdout { format = PPM width = 512 height = 512 }\n");
}
#endif
#if 0
event snapshot (i = 100; i += 100) {
dump (file = "snapshot", t = t);
char name[80];
sprintf (name, "snapshot-%d.gfs", i);
scalar pid[];
foreach()
pid[] = tid();
output_gfs (file = name);
}
#endif
#if TREE
event adapt (i++) {
adapt_wavelet ({f,u}, (double[]){0.005,0.005,0.005,0.005}, maxlevel);
}
#endif
