simulationCases/asyBubbleinSheet_03.c

// Asymmetric Bubble inside a draining sheet. || Uses CLSVOF // Id 1 is liquid pool, and Id 2 is Newtonian gas.

#include “axi.h” #include “navier-stokes/centered.h” #define FILTERED // Smear density and viscosity jumps

#include “two-phase-clsvof.h” #include “integral.h” scalar sigmaf[]; #include “distance.h”

#define tsnap (1e-1) // 0.001 only for some cases.

// Error tolerancs #define fErr (1e-3) // error tolerance in f1 VOF #define KErr (1e-6) // error tolerance in VoF curvature calculated using height function method (see adapt event) #define VelErr (1e-3) // error tolerances in velocity – Use 1e-2 for low Oh and 1e-3 to 5e-3 for high Oh/moderate to high J

#define Ldomain 4

int MAXlevel; // Oh -> Solvent Ohnesorge number // Oha -> air Ohnesorge number // Bo -> Bond number = rhoR^3w2/gamma double Oh, Oha, Bo, offset, tmax; char nameOut[80], dumpFile[80];

// Boundary conditions u.n[left] = dirichlet(-2 * pow(Bo,0.5) * x); u.t[left] = dirichlet(pow(Bo,0.5) * y); p[left] = dirichlet(0.); pf[left] = dirichlet(0.);

u.n[right] = dirichlet(-2 * pow(Bo,0.5) * x); u.t[right] = dirichlet(pow(Bo,0.5) * y); p[right] = dirichlet(0.); pf[right] = dirichlet(0.);

u.n[top] = dirichlet(pow(Bo,0.5) * y); u.t[top] = dirichlet(-2 * pow(Bo,0.5) * x); p[top] = dirichlet(0.); pf[top] = dirichlet(0.);

int main(int argc, char const argv[]) { dtmax = 1e-5; // BEWARE of this for stability issues. L0 = Ldomain; X0 = -0.5 Ldomain;

// Values taken from the terminal MAXlevel = atoi(argv[1]); Oh = atof(argv[2]); Bo = atof(argv[3]); offset = atof(argv[4]); tmax = atof(argv[5]);

// Ensure that all the variables were transferred properly from the terminal or job script. if (argc < 6) { fprintf(ferr, “Lack of command line arguments. Check! Need %d more arguments”, 6 - argc); return 1; } fprintf(ferr, “Level %d, Oh %2.1e, Bo %2.1e, offset %4.3f”, MAXlevel, Oh, Bo, offset); init_grid(1 << 8);

// Create intermediate for all snapshots. char comm[80]; sprintf(comm, “mkdir -p intermediate”); system(comm);

// Name of the restart file, used in writingFiles event. sprintf(dumpFile, “dump”);

rho1 = 1., rho2 = 1e-3; Oha = 2e-5; mu1 = Oh, mu2 = Oha;

// f.sigma = 1.0; d.sigmaf = sigmaf;

run(); }

event init(t = 0) { if (!restore(file = dumpFile)) { // fraction(f, difference(2.25 - x * x, 1 - (x - offset) * (x - offset) - y * y)); foreach () { d[] = -(1 - (x - offset) * (x - offset) - y * y); d[] = fabs(x) > 1.2 ? -d[] : d[]; sigmaf[] = 1.; u.x[] = -2 * pow(Bo,0.5) * x; u.y[] = pow(Bo,0.5) * y; } boundary({f, u}); } // return 1; }

event adapt(i++) { // scalar KAPPA[]; // curvature(f, KAPPA); adapt_wavelet((scalar *){f, u.x, u.y}, (double[]){fErr, VelErr, VelErr}, MAXlevel, MAXlevel - 3); }

// Dumping snapshots event writingFiles(t = 0; t += tsnap; t <= tmax) { dump(file = dumpFile); sprintf(nameOut, “intermediate/snapshot-%5.4f”, t); dump(file = nameOut); }

// Ending Simulation event end(t = end) { }

// Log writing event logWriting(i++) { double ke = 0.; foreach (reduction(+ : ke)) { ke += (2 * pi * y) * (0.5 * rho(f[]) * (sq(u.x[]) + sq(u.y[]))) * sq(Delta); } static FILE *fp; if (pid() == 0) { if (i == 0) { fprintf(ferr, “i dt t ke”); fp = fopen(“log”, “w”); fprintf(fp, “Level %d, Oh %2.1e, Bo %2.1e, offset %4.3f”, MAXlevel, Oh, Bo, offset); fprintf(fp, “i dt t ke”); } else { fp = fopen(“log”, “a”); fprintf(fp, “%d %g %g %g”, i, dt, t, ke); } fprintf(fp, “%d %g %g %g”, i, dt, t, ke); fclose(fp); fprintf(ferr, “%d %g %g %g”, i, dt, t, ke); }

assert(ke > -1e-10); assert(ke < 1e4);

if ((ke > 1e4 || ke < 1e-6) && i > 1e1 && pid() == 0) { const char *message = ke > 1e4 ? “The kinetic energy blew up. Stopping simulation” : “kinetic energy too small now! Stopping!”; fprintf(ferr, “%s”, message); fp = fopen(“log”, “a”); fprintf(fp, “%s”, message); fclose(fp); dump(file = dumpFile); return 1; } }