Issue of convergence error for u.x in the case of rho_r = 0.001 and mu_r = 1
23 views
Skip to first unread message
Anirban Chakraborty
Aug 29, 2025, 4:04:06 AM (9 days ago) Aug 29
to basilisk-fr
/**
# Taylor-Culick Retraction of a Liquid drop */
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
//#include "output_mpi.h"
//#include "adapt.h"
#define LEVEL 10
#define INLEVEL 8
#define MINLEVEL 5
#define MAXLEVEL 10
/**
We are simulating only one half of the liquid sheet.*/
#define L 4
#define L0 (1.5*(L+1))
/**Parameters */
#define R 1.0
#define Oh 1
#define muRatio 1
#define rhoRatio 0.001
//double tC = R*R*Oh;
double tC = 1.0;
double tEnd = 5000;
double utc = 1./(1.*R*Oh);
/**
The boundary conditions are slip walls except right*/
u.n[right] = neumann(0.);
p[right] = dirichlet(0.);
pf[right] = dirichlet(0.);
int main() {
/**
The domain */
size (L0);
init_grid (1 << INLEVEL);
//mass conservation tolerance : very important parameter
TOLERANCE = 1e-5;
NITERMAX = 200;
/**
We define the parameters of the fluids. The internal fluid (the water) is
the fluid 1. The external fluid (the air) is the fluid 2.*/
double rhoInt = 1;
double muInt = 1.;
rho1 = rhoInt, mu1 = muInt*Oh;
rho2 = rhoInt*rhoRatio, mu2 = muInt*muRatio*Oh;
f.sigma = 1.0;
run();
}
/**
We define the geometry with boolean operations. We define a
rectangular area (the intersection of Line and Line_vert). We also
define a half circle. Then, our geometry will be the union of the
rectangle with the half circle. */
double geometry (double x, double y) {
double Line = y-R;
double Line_vert = x-L;
double Circle = sq(x-L)+sq(y)-sq(R);
double right_part = max(-Circle,-Line_vert);
return min(-Line, right_part);
}
scalar fdrop[];
event init (t = 0) {
if (!restore (file = "dump/dump-1118.000000")){
/**
We initialise the geometry of the interface.*/
refine (x < (L + 5*R) && y < (5*R) && level < LEVEL);
refine (x < (L + 6*R) && y > (5*R) && y < (6*R) && level < LEVEL-1.);
refine (x > (L + 5*R) && x < (L + 6*R) && y > 0. && y < (6*R) && level < LEVEL-1.);
refine (x < (L + 7*R) && y > (6*R) && y < (7*R) && level < LEVEL-2.);
refine (x > (L + 6*R) && x < (L + 7*R) && y > 0. && y < (7*R) && level < LEVEL-2.);
refine (x < (L + 8*R) && y > (7*R) && y < (8*R) && level < LEVEL-3.);
refine (x > (L + 7*R) && x < (L + 8*R) && y > 0. && y < (8*R) && level < LEVEL-3.);
refine (x < (L + 9*R) && y > (8*R) && y < (9*R) && level < LEVEL-4.);
refine (x > (L + 8*R) && x < (L + 9*R) && y > 0. && y < (9*R) && level < LEVEL-4.);
refine (x < (L + 10*R) && y > (9*R) && y < (10*R) && level < LEVEL-5.);
refine (x > (L + 9*R) && x < (L + 10*R) && y > 0. && y < (10*R) && level < LEVEL-5.);
refine (x < (L + 11*R) && y > (10*R) && y < (11*R) && level < LEVEL-6.);
refine (x > (L + 10*R) && x < (L + 11*R) && y > 0. && y < (11*R) && level < LEVEL-6.);
// refine (x < (L + 12*R) && y > (11*R) && y < (12*R) && level < LEVEL-7.);
// refine (x > (L + 11*R) && x < (L + 12*R) && y > 0. && y < (12*R) && level < LEVEL-7.);
fraction(f, geometry(x,y));
fraction(fdrop, geometry(x,y));
}
scalar le[];
foreach(){
le[] = level;
}
/* char name_grid[100];
sprintf (name_grid, "grid.png");
FILE* fgrid = fopen (name_grid, "w");
output_ppm (le, file = name_grid,min=MINLEVEL,max=LEVEL,n = 1000,box = {{0.,0.},{L0,L0}});
fclose (fgrid);
*/
}
event adapt (i++) {
double uemax = 1e-3;
double femax = 1e-3;
adapt_wavelet ({f,u.x, u.y, fdrop}, (double[]){femax,uemax,uemax,femax}, maxlevel=MAXLEVEL);
// adapt_wavelet ({f,u.x, u.y}, (double[]){femax,uemax,uemax}, maxlevel=MAXLEVEL);
// adapt_wavelet_leave_interface ({f,u.x, u.y},{f}, (double[]){femax,uemax, uemax}, maxlevel=MAXLEVEL, padding = 0);
}
/**
We output : the interface position and scalar / vector fields . */
event saveDatas (t += tC; t <= tEnd) {
/**
We only output the interface in this function. */
char name[100];
sprintf (name, "interface/interface_%f_%d.dat", t,pid());
FILE * fp = fopen (name, "w");
scalar pid[], ff[];
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1.0e-6 ? 0 : f[] > 1. - 1.0e-6 ? 1. : f[];
}
boundary ({pid,ff});
output_facets (ff, fp);
fclose (fp);
/*
char name[80];
sprintf (name, "interface-%f.txt", t);
FILE* fp = fopen (name, "w");
output_facets (f, fp);
fclose (fp);
*/
/* char frac[80];
sprintf (frac, "velo-%f.dat", t);
FILE * fpu = fopen (frac, "w");
foreach(){
fprintf(fpu, "%f %f %f %f \n",x,y,u.x[],u.y[]);
}
//output_field ({p}, fpf, linear = false);
fclose (fpu);
*/
/* char s[80];
sprintf (s, "field_%g_%d.dat", t,pid());
FILE * fpmu = fopen (s, "w");
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1.0e-6 ? 0 : f[] > 1. - 1.0e-6 ? 1. : f[];
}
boundary ({pid,ff});
foreach(){
// fprintf (fpmu," %f %f %g %g \n",x,y,u.x[],u.y[]);
// output_field ({f, muv}, fpf, linear = true);
}
fclose (fpmu); */
char name_u[80];
sprintf (name_u, "field/field_%f.dat",t);
FILE * fp_u = fopen (name_u, "w");
output_field ({f,p, u.x, u.y}, fp_u, n = 1000,linear = true, box = {{0,0},{1.2*(L+R),4*R}});
// output_matrix_part_mpi(u.x,fp_u, n = 1000,linear =true,box = {{0,0},{1.2*L,4*R}});
fclose (fp_u);
/* char frac2[80];
sprintf (frac2, "vof-%f.dat", t);
FILE * fpf = fopen (frac2, "w");
foreach(){
fprintf(fpf, "%f %f %g\n",x,y,f[]);
}
// output_field ({p}, fpf, linear = true);
fclose (fpf);
*/
/* char frac2[80];
sprintf (frac2, "vof-%f_%d.dat", t, pid());
FILE * fpf = fopen (frac2, "w");
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1.0e-6 ? 0 : f[] > 1. - 1.0e-6 ? 1. : f[];
}
boundary ({pid,ff});
foreach(){
fprintf(fpf, "%f %f %g\n",x,y,f[]);
}
fclose (fpf);*/
scalar vort[];
vorticity (u, vort);
char namev[80];
// sprintf (namev, "field/vort-%f.txt", t);
sprintf (namev, "field/vort-%f", t);
FILE * fpv = fopen (namev, "w");
output_field ({vort}, fpv, n = 1000,linear = true, box = {{0,0},{1.2*(L+R),3*R}});
// output_matrix_part_mpi(vort,fpv, n = 1000,linear =true,box = {{0.,0.},{1.2*L,3*R}});
fclose (fpv);
char name_vort[100];
char name_vof[100];
sprintf (name_vort, "image/vort-%f.png", t);
sprintf (name_vof, "image/vof-%f.png", t);
FILE* fvort = fopen (name_vort, "w");
FILE* fvof = fopen (name_vof, "w");
double om = 0.15 ;//(mu1 /Oh)*1./rho1/R/R;
//output_ppm (vort, file = name_vort, mask =m,min=-1.,max=1., n = 1000,box = {{0.,0.},{10.,2.}});
output_ppm (vort, fvort,min=-om,max=om, map=cool_warm, n = 1000,box = {{0,0},{1.2*(L+R),5*R}}, linear = false);
output_ppm (f, fvof,min=0,max=1, map=jet, n = 1000,box = {{0,0},{1.2*(L+R),5*R}}, linear = false);
fclose (fvort);
fclose (fvof);
char named[80];
sprintf (named, "dump/dump-%f", t);
dump (file = named);
// char d4[80];
// sprintf (d4, "vor-%g", t);
// FILE * g4 = fopen (d4, "w");
// output_matrix_mpi (vort, g4);
// fclose (g4);
}
/**
We output the maximum x position of the liquid finger. Indeed this
position should linearly evolve in time, with a few variations coming
from the capillary wave. */
double xPrev = -1, tPrev = -1, veloPrev = 1;
event extractPosition (i ++) {
/**
We define a new vector field, h. */
vector h[];
/**
We reconstruct the height function field and take the corresponding
maximum along the x axis. */
heights (f, h);
double xMax = -HUGE;;
foreach(reduction(max:xMax))
if (h.x[] != nodata) {
double xi = x + height(h.x[])*Delta;
if (xi > xMax)
xMax = xi;
}
/**
We also output the velocity of the end of the bulge.*/
//if (fabs(xPrev - xMax) >= 0.01){
double veloTip = tPrev >= 0 ? (xPrev - xMax)/(t - tPrev) : 0.;
char name[80];
sprintf (name, "velocity_pos.dat");
static FILE * fp = fopen (name, "w");
fprintf (fp," %g %g %g\n", t, xMax, veloTip);
fflush (fp);
//fprintf(stderr, "%g %g %g\n", t, xMax, veloTip);
//fflush (stderr);
tPrev = t, xPrev = xMax, veloPrev=abs(veloTip);
//}
/**
We output, in the standard output file, the step with the corresponding
time. */
}
event logfile(i+=10) {
printf ("i = %d t = %g\n", i,t);
//printf ("utc = %g", utc);
fflush(stdout);
}
/* Movies post processing */
event movies (t += 0.05*tC; t <= tEnd){
scalar vort[];
vorticity (u, vort);
scalar m[];
foreach() {
if (f[]<0.95 && f[] >0.05)
m[] = -10;
else
m[] = 0.;
}
boundary ({m});
// static FILE * fo = popen ("ppm2gif > vort.gif", "w");
// output_ppm (vort, fo, mask =m,min=-1.,max=1., n = 1000,box = {{0.,0.},{10.,2.}});
// char name_vort[100];
// sprintf (name_vort, "image/vort-%f.png", t);
// FILE* fvort = fopen (name_vort, "w");
double om = 0.15 ;//(mu1 /Oh)*1./rho1/R/R;
//output_ppm (vort, file = name_vort, mask =m,min=-1.,max=1., n = 1000,box = {{0.,0.},{10.,2.}});
// output_ppm (vort, fvort,min=-om,max=om, map=cool_warm, n = 1000,box = {{0,0},{1.2*L,5*R}}, linear = false);
output_ppm (f, file = "f.mp4", box = {{0,0},{(L+R),5*R}},
map = cool_warm, n=1000, linear = false, min = 0, max = 1);
output_ppm (vort, file = "vort.mp4", box = {{0,0},{1.1*(L+R),5*R}},
map = cool_warm, n=1000, linear = false, min = -om, max = om);
// output_ppm (f, file = "vort.mp4", box = {{0,0},{1.2*L,5*R}}, linear = false, min=-om, max=om,);
// output_ppm (vort, file = "vort.mp4",n = 1000, box = {{0,0.},{1.2*L,5*R}}, min = -om, max = om, linear = true, mask = m);
// fclose (fvort);
}
event end (t = tEnd) {}
# Taylor-Culick Retraction of a Liquid drop */
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
//#include "output_mpi.h"
//#include "adapt.h"
#define LEVEL 10
#define INLEVEL 8
#define MINLEVEL 5
#define MAXLEVEL 10
/**
We are simulating only one half of the liquid sheet.*/
#define L 4
#define L0 (1.5*(L+1))
/**Parameters */
#define R 1.0
#define Oh 1
#define muRatio 1
#define rhoRatio 0.001
//double tC = R*R*Oh;
double tC = 1.0;
double tEnd = 5000;
double utc = 1./(1.*R*Oh);
/**
The boundary conditions are slip walls except right*/
u.n[right] = neumann(0.);
p[right] = dirichlet(0.);
pf[right] = dirichlet(0.);
int main() {
/**
The domain */
size (L0);
init_grid (1 << INLEVEL);
//mass conservation tolerance : very important parameter
TOLERANCE = 1e-5;
NITERMAX = 200;
/**
We define the parameters of the fluids. The internal fluid (the water) is
the fluid 1. The external fluid (the air) is the fluid 2.*/
double rhoInt = 1;
double muInt = 1.;
rho1 = rhoInt, mu1 = muInt*Oh;
rho2 = rhoInt*rhoRatio, mu2 = muInt*muRatio*Oh;
f.sigma = 1.0;
run();
}
/**
We define the geometry with boolean operations. We define a
rectangular area (the intersection of Line and Line_vert). We also
define a half circle. Then, our geometry will be the union of the
rectangle with the half circle. */
double geometry (double x, double y) {
double Line = y-R;
double Line_vert = x-L;
double Circle = sq(x-L)+sq(y)-sq(R);
double right_part = max(-Circle,-Line_vert);
return min(-Line, right_part);
}
scalar fdrop[];
event init (t = 0) {
if (!restore (file = "dump/dump-1118.000000")){
/**
We initialise the geometry of the interface.*/
refine (x < (L + 5*R) && y < (5*R) && level < LEVEL);
refine (x < (L + 6*R) && y > (5*R) && y < (6*R) && level < LEVEL-1.);
refine (x > (L + 5*R) && x < (L + 6*R) && y > 0. && y < (6*R) && level < LEVEL-1.);
refine (x < (L + 7*R) && y > (6*R) && y < (7*R) && level < LEVEL-2.);
refine (x > (L + 6*R) && x < (L + 7*R) && y > 0. && y < (7*R) && level < LEVEL-2.);
refine (x < (L + 8*R) && y > (7*R) && y < (8*R) && level < LEVEL-3.);
refine (x > (L + 7*R) && x < (L + 8*R) && y > 0. && y < (8*R) && level < LEVEL-3.);
refine (x < (L + 9*R) && y > (8*R) && y < (9*R) && level < LEVEL-4.);
refine (x > (L + 8*R) && x < (L + 9*R) && y > 0. && y < (9*R) && level < LEVEL-4.);
refine (x < (L + 10*R) && y > (9*R) && y < (10*R) && level < LEVEL-5.);
refine (x > (L + 9*R) && x < (L + 10*R) && y > 0. && y < (10*R) && level < LEVEL-5.);
refine (x < (L + 11*R) && y > (10*R) && y < (11*R) && level < LEVEL-6.);
refine (x > (L + 10*R) && x < (L + 11*R) && y > 0. && y < (11*R) && level < LEVEL-6.);
// refine (x < (L + 12*R) && y > (11*R) && y < (12*R) && level < LEVEL-7.);
// refine (x > (L + 11*R) && x < (L + 12*R) && y > 0. && y < (12*R) && level < LEVEL-7.);
fraction(f, geometry(x,y));
fraction(fdrop, geometry(x,y));
}
scalar le[];
foreach(){
le[] = level;
}
/* char name_grid[100];
sprintf (name_grid, "grid.png");
FILE* fgrid = fopen (name_grid, "w");
output_ppm (le, file = name_grid,min=MINLEVEL,max=LEVEL,n = 1000,box = {{0.,0.},{L0,L0}});
fclose (fgrid);
*/
}
event adapt (i++) {
double uemax = 1e-3;
double femax = 1e-3;
adapt_wavelet ({f,u.x, u.y, fdrop}, (double[]){femax,uemax,uemax,femax}, maxlevel=MAXLEVEL);
// adapt_wavelet ({f,u.x, u.y}, (double[]){femax,uemax,uemax}, maxlevel=MAXLEVEL);
// adapt_wavelet_leave_interface ({f,u.x, u.y},{f}, (double[]){femax,uemax, uemax}, maxlevel=MAXLEVEL, padding = 0);
}
/**
We output : the interface position and scalar / vector fields . */
event saveDatas (t += tC; t <= tEnd) {
/**
We only output the interface in this function. */
char name[100];
sprintf (name, "interface/interface_%f_%d.dat", t,pid());
FILE * fp = fopen (name, "w");
scalar pid[], ff[];
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1.0e-6 ? 0 : f[] > 1. - 1.0e-6 ? 1. : f[];
}
boundary ({pid,ff});
output_facets (ff, fp);
fclose (fp);
/*
char name[80];
sprintf (name, "interface-%f.txt", t);
FILE* fp = fopen (name, "w");
output_facets (f, fp);
fclose (fp);
*/
/* char frac[80];
sprintf (frac, "velo-%f.dat", t);
FILE * fpu = fopen (frac, "w");
foreach(){
fprintf(fpu, "%f %f %f %f \n",x,y,u.x[],u.y[]);
}
//output_field ({p}, fpf, linear = false);
fclose (fpu);
*/
/* char s[80];
sprintf (s, "field_%g_%d.dat", t,pid());
FILE * fpmu = fopen (s, "w");
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1.0e-6 ? 0 : f[] > 1. - 1.0e-6 ? 1. : f[];
}
boundary ({pid,ff});
foreach(){
// fprintf (fpmu," %f %f %g %g \n",x,y,u.x[],u.y[]);
// output_field ({f, muv}, fpf, linear = true);
}
fclose (fpmu); */
char name_u[80];
sprintf (name_u, "field/field_%f.dat",t);
FILE * fp_u = fopen (name_u, "w");
output_field ({f,p, u.x, u.y}, fp_u, n = 1000,linear = true, box = {{0,0},{1.2*(L+R),4*R}});
// output_matrix_part_mpi(u.x,fp_u, n = 1000,linear =true,box = {{0,0},{1.2*L,4*R}});
fclose (fp_u);
/* char frac2[80];
sprintf (frac2, "vof-%f.dat", t);
FILE * fpf = fopen (frac2, "w");
foreach(){
fprintf(fpf, "%f %f %g\n",x,y,f[]);
}
// output_field ({p}, fpf, linear = true);
fclose (fpf);
*/
/* char frac2[80];
sprintf (frac2, "vof-%f_%d.dat", t, pid());
FILE * fpf = fopen (frac2, "w");
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1.0e-6 ? 0 : f[] > 1. - 1.0e-6 ? 1. : f[];
}
boundary ({pid,ff});
foreach(){
fprintf(fpf, "%f %f %g\n",x,y,f[]);
}
fclose (fpf);*/
scalar vort[];
vorticity (u, vort);
char namev[80];
// sprintf (namev, "field/vort-%f.txt", t);
sprintf (namev, "field/vort-%f", t);
FILE * fpv = fopen (namev, "w");
output_field ({vort}, fpv, n = 1000,linear = true, box = {{0,0},{1.2*(L+R),3*R}});
// output_matrix_part_mpi(vort,fpv, n = 1000,linear =true,box = {{0.,0.},{1.2*L,3*R}});
fclose (fpv);
char name_vort[100];
char name_vof[100];
sprintf (name_vort, "image/vort-%f.png", t);
sprintf (name_vof, "image/vof-%f.png", t);
FILE* fvort = fopen (name_vort, "w");
FILE* fvof = fopen (name_vof, "w");
double om = 0.15 ;//(mu1 /Oh)*1./rho1/R/R;
//output_ppm (vort, file = name_vort, mask =m,min=-1.,max=1., n = 1000,box = {{0.,0.},{10.,2.}});
output_ppm (vort, fvort,min=-om,max=om, map=cool_warm, n = 1000,box = {{0,0},{1.2*(L+R),5*R}}, linear = false);
output_ppm (f, fvof,min=0,max=1, map=jet, n = 1000,box = {{0,0},{1.2*(L+R),5*R}}, linear = false);
fclose (fvort);
fclose (fvof);
char named[80];
sprintf (named, "dump/dump-%f", t);
dump (file = named);
// char d4[80];
// sprintf (d4, "vor-%g", t);
// FILE * g4 = fopen (d4, "w");
// output_matrix_mpi (vort, g4);
// fclose (g4);
}
/**
We output the maximum x position of the liquid finger. Indeed this
position should linearly evolve in time, with a few variations coming
from the capillary wave. */
double xPrev = -1, tPrev = -1, veloPrev = 1;
event extractPosition (i ++) {
/**
We define a new vector field, h. */
vector h[];
/**
We reconstruct the height function field and take the corresponding
maximum along the x axis. */
heights (f, h);
double xMax = -HUGE;;
foreach(reduction(max:xMax))
if (h.x[] != nodata) {
double xi = x + height(h.x[])*Delta;
if (xi > xMax)
xMax = xi;
}
/**
We also output the velocity of the end of the bulge.*/
//if (fabs(xPrev - xMax) >= 0.01){
double veloTip = tPrev >= 0 ? (xPrev - xMax)/(t - tPrev) : 0.;
char name[80];
sprintf (name, "velocity_pos.dat");
static FILE * fp = fopen (name, "w");
fprintf (fp," %g %g %g\n", t, xMax, veloTip);
fflush (fp);
//fprintf(stderr, "%g %g %g\n", t, xMax, veloTip);
//fflush (stderr);
tPrev = t, xPrev = xMax, veloPrev=abs(veloTip);
//}
/**
We output, in the standard output file, the step with the corresponding
time. */
}
event logfile(i+=10) {
printf ("i = %d t = %g\n", i,t);
//printf ("utc = %g", utc);
fflush(stdout);
}
/* Movies post processing */
event movies (t += 0.05*tC; t <= tEnd){
scalar vort[];
vorticity (u, vort);
scalar m[];
foreach() {
if (f[]<0.95 && f[] >0.05)
m[] = -10;
else
m[] = 0.;
}
boundary ({m});
// static FILE * fo = popen ("ppm2gif > vort.gif", "w");
// output_ppm (vort, fo, mask =m,min=-1.,max=1., n = 1000,box = {{0.,0.},{10.,2.}});
// char name_vort[100];
// sprintf (name_vort, "image/vort-%f.png", t);
// FILE* fvort = fopen (name_vort, "w");
double om = 0.15 ;//(mu1 /Oh)*1./rho1/R/R;
//output_ppm (vort, file = name_vort, mask =m,min=-1.,max=1., n = 1000,box = {{0.,0.},{10.,2.}});
// output_ppm (vort, fvort,min=-om,max=om, map=cool_warm, n = 1000,box = {{0,0},{1.2*L,5*R}}, linear = false);
output_ppm (f, file = "f.mp4", box = {{0,0},{(L+R),5*R}},
map = cool_warm, n=1000, linear = false, min = 0, max = 1);
output_ppm (vort, file = "vort.mp4", box = {{0,0},{1.1*(L+R),5*R}},
map = cool_warm, n=1000, linear = false, min = -om, max = om);
// output_ppm (f, file = "vort.mp4", box = {{0,0},{1.2*L,5*R}}, linear = false, min=-om, max=om,);
// output_ppm (vort, file = "vort.mp4",n = 1000, box = {{0,0.},{1.2*L,5*R}}, min = -om, max = om, linear = true, mask = m);
// fclose (fvort);
}
event end (t = tEnd) {}
kindly help
Reply all
Reply to author
Forward
0 new messages