Finite difference solver with 5th-order WENO method for two-phase Rayleigh-Taylor fluid instability
Date:
The nonlinear convex splitting method for Cahn-Hilliard equation, projection method for Navier-Stokes equations, and finite difference method for spatial discretization can refer to: H.G. Lee, K. Kim, J. Kim, On the long time simulation of the Rayleigh-Taylor instability, International Journal for Numerical Methods in Engineering, Vol 85, pp. 1633-1647, 2011.
The 5th-order WENO finite difference scheme for treating velocity advection can refer to: G.-S. Jiang, C.-C. Wu, A High-Order WENO Finite Difference Scheme for the Equations of Ideal Magnetohydrodynamics, Journal of Computational Physics, Vol. 150, pp. 561-594, 1999.
Note: The periodic boundary condition is used for all variables along x-direction. On the top and bottom boundaries, the velocities are zero, the pressure and phase-fielsd variable satisfy the homogeneous-Neumann boundary condition.
/** The following C code is used to calculate phase-field Cahn-Hilliard equation and incompressible Navier-Stokes equations. ***/
#include
int nx, ny, n_level, p_relax, c_relax, count; double pi, xleft, xright, yleft, yright, h, dt, gam, Cahn, rho1, rho2, vis1, vis2, Pe, Re, We, Fr, **adv_u, **adv_v, **adv_c, **fx, **fy, **tu, **tv, **worku, **workv, **workp, **rho, **vis, **ct, **sc, **smu, **mu; char bufferu[100], bufferv[100],bufferp[100],bufferc[100]; int main() { extern int nx, ny, n_level, p_relax, c_relax, count; extern double pi, xleft, xright, yleft, yright, h, dt, gam, Cahn, rho1, rho2, vis1, vis2, Pe, Re, We, Fr, **adv_u, **adv_v, **adv_c, **fx, **fy, **tu, **tv, **worku, **workv, **workp, **rho, **vis, **ct, **sc, **smu, **mu; extern char bufferu[100], bufferv[100],bufferp[100],bufferc[100]; int it, max_it, ns; double **u, **nu, **v, **nv, **oc, **nc, **p;
FILE *fu, *fv, *fc, *fp;
nx = gnx;
ny = gny;
n_level = (int)(log(nx)/log(2)-0.9);
p_relax = 7;
c_relax = 5;
pi = 4.0*atan(1.0);
xleft = 0.0, xright = 2.0;
yleft = 0.0, yright = xright*(double)ny/(double)nx;
count = 0;
h = xright/(double)nx;
dt = 0.16*h;
gam = 4.0*h/(2.0*sqrt(2.0)*atanh(0.9));
Cahn = pow(gam,2);
max_it = 8000;
ns = max_it/40;
rho1 = 3.0; //oil
rho2 = 1.0; // water
vis1 = 1.0; //oil
vis2 = 1.0; // water
Pe = 10.0/gam;
Re = 3000.0;
We = 50.0;
Fr = 1.0;
printf("nx = %d, ny = %d\n", nx, ny);
printf("n_level = %d\n", n_level);
printf("dt = %f\n", dt);
printf("max_it = %d\n", max_it);
printf("ns = %d\n", ns);
printf("Reynolds number = %f\n", Re);
printf("Weber number = %f\n", We);
u = dmatrix(-1, nx+1, 0, ny+1);
nu = dmatrix(-1, nx+1, 0, ny+1);
v = dmatrix(0, nx+1, -1, ny+1);
nv = dmatrix(0, nx+1, -1, ny+1);
oc = dmatrix(0, nx+1, 0, ny+1);
nc = dmatrix(0, nx+1, 0, ny+1);
p = dmatrix(1, nx, 1, ny);
adv_u = dmatrix(0, nx, 1, ny);
adv_v = dmatrix(1, nx, 0, ny);
adv_c = dmatrix(1, nx, 1, ny);
fx = dmatrix(0, nx, 1, ny);
fy = dmatrix(1, nx, 0, ny);
tu = dmatrix(0, nx, 1, ny);
tv = dmatrix(1, nx, 0, ny);
worku = dmatrix(0, nx, 1, ny);
workv = dmatrix(1, nx, 0, ny);
workp = dmatrix(1, nx, 1, ny);
rho = dmatrix(0, nx+1, 0, ny+1);
vis = dmatrix(0, nx+1, 0, ny+1);
ct = dmatrix(1, nx, 1, ny);
sc = dmatrix(1, nx, 1, ny);
smu = dmatrix(1, nx, 1, ny);
mu = dmatrix(1, nx, 1, ny);
zero_matrix(mu, 1, nx, 1, ny);
sprintf(bufferu,"dataur.m");
sprintf(bufferv,"datavr.m");
sprintf(bufferp,"datapr.m");
sprintf(bufferc,"datacr.m");
fu = fopen(bufferu,"w");
fv = fopen(bufferv,"w");
fp = fopen(bufferp,"w");
fc = fopen(bufferc,"w");
fclose(fu);
fclose(fv);
fclose(fp);
fclose(fc);
initialization(u, v, oc, p);
mat_copy(nu, u, 0, nx, 1, ny);
mat_copy(nv, v, 1, nx, 0, ny);
mat_copy(nc, oc, 1, nx, 1, ny);
print_data(nu, nv, nc, p);
for (it=1; it<=max_it; it++) {
cal_den(oc, rho);
augmenc(rho, nx, ny);
cal_vis(oc, vis);
augmenc(vis, nx, ny);
full_step(u, v, oc, nu, nv, nc, p);
mat_copy(u, nu, 0, nx, 1, ny);
mat_copy(v, nv, 1, nx, 0, ny);
mat_copy(oc, nc, 1, nx, 1, ny);
printf("it = %d\n", it);
if (it % ns == 0) {
count++;
print_data(nu, nv, nc, p);
printf("print out counts %d\n", count);
}
}
return 0; }
void initialization(double **u, double **v, double **c, double **p) { extern double xright, yright, h, gam;
int i, j;
double x, y;
ijloop {
x = ((double)i-0.5)*h;
y = ((double)j-0.5)*h;
c[i][j] = tanh((y-0.5+0.01*(1.0-2.0*rand()/RAND_MAX))/(sqrt(2.0)*gam));
p[i][j] = 0.0;
}
zero_matrix(u, 0, nx, 1, ny);
zero_matrix(v, 1, nx, 0, ny);
}
void cal_den(double **c, double **rho) { extern double rho1, rho2;
int i, j;
ijloop {
rho[i][j] = 0.5*fabs(rho1*(1.0+c[i][j])+rho2*(1.0-c[i][j]))/rho2;
}
}
void cal_vis(double **c, double **vis) { extern double vis1, vis2;
int i, j;
ijloop {
vis[i][j] = 0.5*fabs(vis1*(1.0+c[i][j])+vis2*(1.0-c[i][j]))/vis2;
}
}
void full_step(double **u, double **v, double **c, double **nu, double **nv, double **nc, double **p) { extern int nx, ny; extern double dt, **adv_u, **adv_v, **adv_c, **fx, **fy, **tu, **tv, **worku, **workv, **rho;
int i, j;
advection_step(u, v, c, adv_u, adv_v, adv_c);
i0jloop {
fx[i][j] = 0.0;
}
ij0loop {
fy[i][j] = -1.0/Fr;
}
temp_uv(tu, tv, u, v, adv_u, adv_v, fx, fy);
Poisson(tu, tv, p);
grad_p(p, worku, workv, nx, ny);
i0jloop {
nu[i][j] = tu[i][j] - dt*worku[i][j]/(0.5*(rho[i+1][j]+rho[i][j]));
}
ij0loop {
nv[i][j] = tv[i][j] - dt*workv[i][j]/(0.5*(rho[i][j+1]+rho[i][j]));
}
iloop {
nv[i][0] = nv[i][ny] = 0.0;
}
cahn(c, adv_c, nc);
}
void advection_step(double **u, double **v, double **c, double **adv_u, double **adv_v, double **adv_c) { extern int nx, ny; extern double h;
int i, j;
augmenuv(u, v);
augmenc(c, nx, ny);
advection_uv(u,v,adv_u,adv_v);
ijloop { // conservative
adv_c[i][j] = (u[i][j]*(c[i+1][j]+c[i][j])-u[i-1][j]*(c[i][j]+c[i-1][j]))/(2.0*h)
+ (v[i][j]*(c[i][j+1]+c[i][j])-v[i][j-1]*(c[i][j]+c[i][j-1]))/(2.0*h);
}
}
void advection_uv(double **u, double **v, double **adv_u, double **adv_v){ extern int nx, ny; extern double h;
int i,j,k; double *temp_x, *temp_y, **fuc, **fplus, **fminus, **iv, **iu, max_vel, *temp_x2, *temp_y2, **fuc2, **fplus2, **fminus2, max_vel2;
augmenuv(u,v);
temp_x = dvector(0,nx+1); temp_y = dvector(0,ny);
temp_x2 = dvector(0,nx); temp_y2 = dvector(0,ny+1);
iv = dmatrix(0,nx,1,ny); iu = dmatrix(1,nx,0,ny);
fuc = dmatrix(-3,nx+3,-2,ny+3); fplus = dmatrix(-3,nx+3,-2,ny+3); fminus = dmatrix(-3,nx+3,-2,ny+3);
fuc2 = dmatrix(-2,nx+3,-3,ny+3); fplus2 = dmatrix(-2,nx+3,-3,ny+3); fminus2 = dmatrix(-2,nx+3,-3,ny+3);
//start computation of adv_u
// x - flux of adv_u
max_vel = mat_max(u, 0, nx, 1, ny);
matmult(fuc, u, u, 0, nx, 1, ny); mat_comb(fplus, 0.5, fuc, max_vel, u, 0, nx, 1, ny); mat_comb(fminus, 0.5, fuc, -max_vel, u, 0, nx, 1, ny);
for (j = 1; j<= ny; j++){
for (k = 1; k <= 3; k++){ fuc[-k][j] = fuc[nx-1][j]; fuc[nx+k][j] = fuc[1][j]; fplus[-k][j] = fplus[nx-1][j]; fplus[nx+k][j] = fplus[1][j]; fminus[-k][j] = fminus[nx-1][j]; fminus[nx+k][j] = fminus[1][j]; }
for (i = 0; i <= nx+1; i++){ temp_x[i] = ( -fuc[i-2][j] + 7.0fuc[i-1][j] + 7.0fuc[i][j] - fuc[i+1][j] )/12.0
- vpn( fplus[i-2][j] - fplus[i-3][j], fplus[i-1][j]-fplus[i-2][j], fplus[i][j]-fplus[i-1][j], fplus[i+1][j]-fplus[i][j] )
- vpn( fminus[i+2][j]-fminus[i+1][j], fminus[i+1][j]-fminus[i][j], fminus[i][j]-fminus[i-1][j], fminus[i-1][j]-fminus[i-2][j] ); }
for (i = 0; i<= nx; i++){ adv_u[i][j] = (temp_x[i+1] - temp_x[i])/h; }
}
// y - flux of adv_u
i0jloop{ iv[i][j] = 0.25*( v[i][j-1]+v[i+1][j-1]+v[i][j]+v[i+1][j] ); }
max_vel = mat_max(iv, 0, nx, 1, ny);
matmult(fuc, iv, u, 0, nx, 1, ny); mat_comb(fplus, 0.5, fuc, max_vel, u, 0, nx, 1, ny); mat_comb(fminus, 0.5, fuc, -max_vel, u, 0, nx, 1, ny);
for (i = 0; i <= nx; i++){
for(k = 1; k<=3; k++){
fuc[i][1-k] = -fuc[i][1];
fuc[i][ny+k] = -fuc[i][ny];
fplus[i][1-k] = -fplus[i][1];
fplus[i][ny+k] = -fplus[i][ny];
fminus[i][1-k] = -fminus[i][1];
fminus[i][ny+k] = -fminus[i][ny]; }
for (j = 0; j<= ny; j++){
temp_y[j] = (-fuc[i][j-1] + 7.0*fuc[i][j] + 7.0*fuc[i][j+1] - fuc[i][j+2] )/12.0
- vpn(fplus[i][j-1]-fplus[i][j-2], fplus[i][j]-fplus[i][j-1], fplus[i][j+1]-fplus[i][j], fplus[i][j+2]-fplus[i][j+1])
+ vpn(fminus[i][j+3]-fminus[i][j+2], fminus[i][j+2]-fminus[i][j+1], fminus[i][j+1]-fminus[i][j], fminus[i][j]-fminus[i][j-1]);
}
for (j = 1; j <= ny; j++){
adv_u[i][j] += (temp_y[j] - temp_y[j-1])/h; }
}
//end computation of adv_u
//start computation of adv_v
// x - flux of adv_v
ij0loop{ iu[i][j] = 0.25*( u[i-1][j]+u[i][j]+u[i-1][j+1]+u[i][j+1] ); }
max_vel2 = mat_max(iu, 1, nx, 0, ny);
matmult(fuc2, iu, v, 1, nx, 0, ny); mat_comb(fplus2, 0.5, fuc2, max_vel2, v, 1, nx, 0, ny); mat_comb(fminus2, 0.5, fuc2, -max_vel2, v, 1, nx, 0, ny);
for (j = 0; j<= ny; j++){
for (k = 1; k<=3; k++){
fuc2[1-k][j] = fuc2[nx][j];
fuc2[nx+k][j] = fuc2[1][j];
fplus2[1-k][j] = fplus2[nx][j];
fplus2[nx+k][j] = fplus2[1][j];
fminus2[1-k][j] = fminus2[nx][j];
fminus2[nx+k][j] = fminus2[1][j]; }
for( i = 0; i<= nx; i++){
temp_x2[i] = (-fuc2[i-1][j]+7.0*fuc2[i][j]+7.0*fuc2[i+1][j]-fuc2[i+2][j] )/12.0
- vpn(fplus2[i-1][j]-fplus2[i-2][j], fplus2[i][j]-fplus2[i-1][j], fplus2[i+1][j]-fplus2[i][j], fplus2[i+2][j]-fplus2[i+1][j]) + vpn(fminus2[i+3][j]-fminus2[i+2][j], fminus2[i+2][j]-fminus2[i+1][j], fminus2[i+1][j]-fminus2[i][j], fminus2[i][j]-fminus2[i-1][j]);
}
for (i = 1; i <= nx; i++){
adv_v[i][j] = (temp_x2[i] - temp_x2[i-1])/h; }
}
// v - flux of adv_v
max_vel2 = mat_max(v, 1, nx, 0, ny);
matmult(fuc2, v, v, 1, nx, 0, ny); mat_comb(fplus2, 0.5, fuc2, max_vel2, v, 1, nx, 0, ny); mat_comb(fminus2, 0.5, fuc2, -max_vel2, v, 1, nx, 0, ny);
for (i = 1; i<= nx; i++){
for(k = 1; k<=3; k++){
fuc2[i][-k] = -fuc2[i][1];
fuc2[i][ny+k] = -fuc2[i][ny-1];
fplus2[i][-k] = -fplus2[i][1];
fplus2[i][ny+k] = -fplus2[i][ny-1];
fminus2[i][-k] = -fminus2[i][1];
fminus2[i][ny+k] = -fminus2[i][ny-1]; }
for (j = 0; j <= ny+1; j++){
temp_y2[j] = ( -fuc2[i][j-2] + 7.0*fuc2[i][j-1] + 7.0*fuc2[i][j] - fuc2[i][j+1] )/12.0 - vpn( fplus2[i][j-2] - fplus2[i][j-3], fplus2[i][j-1]-fplus2[i][j-2], fplus2[i][j]-fplus2[i][j-1], fplus2[i][j+1]-fplus2[i][j] ) + vpn( fminus2[i][j+2]-fminus2[i][j+1], fminus2[i][j+1]-fminus2[i][j], fminus2[i][j]-fminus2[i][j-1], fminus2[i][j-1]-fminus2[i][j-2] );
}
for (j = 0; j<= ny; j++){ adv_v[i][j] += (temp_y2[j+1] - temp_y2[j])/h; }
}
free_dvector(temp_x, 0, nx+1);
free_dvector(temp_y, 0, ny);
free_dvector(temp_x2, 0, nx);
free_dvector(temp_y2, 0, ny+1);
free_dmatrix(fuc, -3, nx+3, -2, ny+3);
free_dmatrix(fplus, -3, nx+3, -2, ny+3);
free_dmatrix(fminus, -3, nx+3, -2, ny+3);
free_dmatrix(fuc2, -2, nx+3, -3, ny+3);
free_dmatrix(fplus2, -2, nx+3, -3, ny+3);
free_dmatrix(fminus2, -2, nx+3, -3, ny+3);
free_dmatrix(iv, 0, nx, 1, ny);
free_dmatrix(iu, 1, nx, 0, ny);
}
void augmenuv(double **u, double **v) { extern int nx, ny;
int i, j;
double aa;
for (j=1; j<=ny; j++) {
u[-1][j] = u[nx-1][j];
u[nx+1][j] = u[1][j];
}
for (i=-1; i<=nx+1; i++) {
u[i][0] = -u[i][1];
u[i][ny+1] = -u[i][ny];
}
for (j=0; j<=ny; j++) {
v[0][j] = v[nx][j];
v[nx+1][j] = v[1][j];
}
for (i=0; i<=nx+1; i++) {
v[i][0] = v[i][ny] = 0.0;
v[i][-1] = -v[i][1];
v[i][ny+1] = -v[i][ny-1];
}
}
void augmenc(double **c, int nxt, int nyt) { int i, j;
for (j=1; j<=nyt; j++) {
c[0][j] = c[nxt][j];
c[nxt+1][j] = c[1][j];
}
for (i=0; i<=nxt+1; i++) {
c[i][0] = c[i][1];
c[i][nyt+1] = c[i][nyt];
}
}
void temp_uv(double **tu, double **tv, double **u, double **v, double **adv_u, double **adv_v, double **fx, double **fy) { extern double h, Re, **rho, **vis;
int i, j;
double vis_u, vis_d, vis_r, vis_l;
i0jloop {
vis_u = 0.25*(vis[i][j]+vis[i+1][j]+vis[i][j+1]+vis[i+1][j+1]);
vis_d = 0.25*(vis[i][j-1]+vis[i+1][j-1]+vis[i][j]+vis[i+1][j]);
tu[i][j] = u[i][j] + dt*( ( 2.0*(vis[i+1][j]*(u[i+1][j]-u[i][j])-vis[i][j]*(u[i][j]-u[i-1][j]))
+ vis_u*(u[i][j+1]-u[i][j]) - vis_d*(u[i][j]-u[i][j-1])
+ vis_u*(v[i+1][j]-v[i][j]) - vis_d*(v[i+1][j-1]-v[i][j-1]) )/(h*h*(0.5*(rho[i+1][j]+rho[i][j]))*Re)
- adv_u[i][j] + fx[i][j] );
}
ij0loop {
vis_r = 0.25*(vis[i][j]+vis[i+1][j]+vis[i][j+1]+vis[i+1][j+1]);
vis_l = 0.25*(vis[i-1][j]+vis[i][j]+vis[i-1][j+1]+vis[i][j+1]);
tv[i][j] = v[i][j] + dt*( ( 2.0*vis[i][j+1]*(v[i][j+1]-v[i][j])-2.0*vis[i][j]*(v[i][j]-v[i][j-1])
+ vis_r*(v[i+1][j]-v[i][j]) - vis_l*(v[i][j]-v[i-1][j])
+ vis_r*(u[i][j+1]-u[i][j]) - vis_l*(u[i-1][j+1]-u[i-1][j]) )/(h*h*(0.5*(rho[i][j+1]+rho[i][j]))*Re)
- adv_v[i][j] + fy[i][j]);
}
}
void Poisson(double **tu, double **tv, double **p) { extern int nx, ny; extern double **workp;
source_uv(tu, tv, workp, nx, ny);
MG_Poisson(p, workp);
}
void source_uv(double **tu, double **tv, double **divuv, int nxt, int nyt) { extern double dt;
int i, j;
div_uv(tu, tv, divuv, nxt, nyt);
ijloopt {
divuv[i][j] /= dt;
}
}
void div_uv(double **tu, double **tv, double **divuv, int nxt, int nyt) { extern double xright;
int i, j;
double ht;
ht = xright/(double)nxt;
ijloopt {
divuv[i][j] = (tu[i][j] - tu[i-1][j] + tv[i][j] - tv[i][j-1])/ht;
}
}
void MG_Poisson(double **p, double **f) { extern int nx, ny; extern double **rho;
int it_mg = 1, max_it = 100;
double resid = 1.0, resid2 = 10.0, tol = 1.0e-5, **sor;
sor = dmatrix(1, nx, 1, ny);
mat_copy(sor, p, 1, nx, 1, ny);
while (it_mg <= max_it && resid >= tol) {
vcycle_uv(p, f, rho, nx, ny, 1);
pressure_update(p);
mat_sub(sor, sor, p, 1, nx, 1, ny);
resid = mat_max(sor, 1, nx, 1, ny);
mat_copy(sor, p, 1, nx, 1, ny);
if (resid > resid2)
it_mg = max_it;
else
resid2 = resid;
it_mg++;
}
printf("Mac pressure iteration = %d residual = %16.14f\n", it_mg-1, resid);
free_dmatrix(sor, 1, nx, 1, ny);
return; }
void vcycle_uv(double **uf, double **ff, double **wf, int nxf, int nyf, int ilevel) { extern int n_level;
relax_p(uf, ff, wf, ilevel, nxf, nyf);
if (ilevel < n_level) {
int nxc, nyc;
double **rf, **fc, **uc, **wc;
nxc = nxf/2, nyc = nyf/2;
rf = dmatrix(1, nxf, 1, nyf);
fc = dmatrix(1, nxc, 1, nyc);
uc = dmatrix(1, nxc, 1, nyc);
wc = dmatrix(0, nxc+1, 0, nyc+1);
residual_den(rf, uf, ff, wf, nxf, nyf);
restrict1(rf, fc, nxc, nyc);
restrict1(wf, wc, nxc, nyc);
augmenc(wc, nxc, nyc);
zero_matrix(uc, 1, nxc, 1, nyc);
vcycle_uv(uc, fc, wc, nxc, nyc, ilevel+1);
prolong(uc, rf, nxc, nyc);
mat_add(uf, uf, rf, 1, nxf, 1, nyf);
relax_p(uf, ff, wf, ilevel, nxf, nyf);
free_dmatrix(rf, 1, nxf, 1, nyf);
free_dmatrix(fc, 1, nxc, 1, nyc);
free_dmatrix(uc, 1, nxc, 1, nyc);
free_dmatrix(wc, 0, nxc+1, 0, nyc+1);
}
}
void relax_p(double **p, double **f, double **w, int ilevel, int nxt, int nyt) { extern int ny, p_relax; extern double xright, Fr;
int i, j, iter;
double ht, ht2, a[4], sorc, coef;
ht = xright/(double)nxt;
ht2 = pow(ht,2);
for (iter=1; iter<=p_relax; iter++) {
ijloopt {
a[0] = 2.0/(w[i+1][j]+w[i][j]);
a[1] = 2.0/(w[i][j]+w[i-1][j]);
a[2] = 2.0/(w[i][j+1]+w[i][j]);
a[3] = 2.0/(w[i][j]+w[i][j-1]);
sorc = f[i][j];
coef = -(a[0] + a[1])/ht2;
if (i==1)
sorc -= (a[0]*p[i+1][j] + a[1]*p[nxt][j])/ht2;
else if (i==nxt)
sorc -= (a[0]*p[1][j] + a[1]*p[i-1][j])/ht2;
else
sorc -= (a[0]*p[i+1][j] + a[1]*p[i-1][j])/ht2;
if (j==1) {
if (nyt==ny)
sorc -= a[2]*p[i][j+1]/ht2 + 1.0/(ht*Fr);
else
sorc -= a[2]*p[i][j+1]/ht2;
coef -= a[2]/ht2;
}
else if (j==nyt) {
if (nyt==ny)
sorc -= -1.0/(ht*Fr) + a[3]*p[i][j-1]/ht2;
else
sorc -= a[3]*p[i][j-1]/ht2;
coef -= a[3]/ht2;
}
else {
sorc -= (a[2]*p[i][j+1] + a[3]*p[i][j-1])/ht2;
coef -= (a[2] + a[3])/ht2;
}
p[i][j] = sorc/coef;
}
}
}
void residual_den(double **r, double **u, double **f, double **den, int nxt, int nyt) { int i, j; double **dpdx, **dpdy;
dpdx = dmatrix(0, nxt, 1, nyt);
dpdy = dmatrix(1, nxt, 0, nyt);
grad_p(u, dpdx, dpdy, nxt, nyt);
i0jloopt {
dpdx[i][j] = dpdx[i][j]/(0.5*(den[i+1][j]+den[i][j]));
}
ij0loopt {
dpdy[i][j] = dpdy[i][j]/(0.5*(den[i][j+1]+den[i][j]));
}
div_uv(dpdx, dpdy, r, nxt, nyt);
mat_sub(r, f, r, 1, nxt, 1, nyt);
free_dmatrix(dpdx, 0, nxt, 1, nyt);
free_dmatrix(dpdy, 1, nxt, 0, nyt); }
void grad_p(double **p, double **dpdx, double **dpdy, int nxt, int nyt) { extern double xright, Fr, **rho;
int i, j;
double ht;
ht = xright/(double)nxt;
i0jloopt {
if (i==0)
dpdx[0][j] = (p[1][j] - p[nxt][j])/ht;
else if (i==nxt)
dpdx[nxt][j] = (p[1][j] - p[nxt][j])/ht;
else
dpdx[i][j] = (p[i+1][j] - p[i][j])/ht;
}
ij0loopt {
if (j==0) {
if (nyt==ny)
dpdy[i][0] = -0.5*(rho[i][1]+rho[i][0])/Fr;
else
dpdy[i][0] = 0.0;
}
else if (j==nyt) {
if (nyt==ny)
dpdy[i][nyt] = -0.5*(rho[i][nyt+1]+rho[i][nyt])/Fr;
else
dpdy[i][nyt] = 0.0;
}
else
dpdy[i][j] = (p[i][j+1] - p[i][j])/ht;
}
}
void restrict1(double **u_fine, double **u_coarse, int nxt, int nyt) { int i, j;
ijloopt {
u_coarse[i][j] = 0.25*(u_fine[2*i-1][2*j-1] + u_fine[2*i-1][2*j]
+ u_fine[2*i][2*j-1] + u_fine[2*i][2*j]);
}
}
void prolong(double **u_coarse, double **u_fine, int nxt, int nyt) { int i, j;
ijloopt {
u_fine[2*i-1][2*j-1] = u_fine[2*i-1][2*j] =
u_fine[2*i][2*j-1] = u_fine[2*i][2*j] = u_coarse[i][j];
}
}
void cahn(double **c_old, double **adv_c, double **c_new) { extern int nx, ny; extern double **ct, **sc, **smu, **mu;
int it_mg = 1, max_it_CH = 50;
double resid = 1.0, tol = 1.0e-8;
mat_copy(ct, c_old, 1, nx, 1, ny);
source(c_old, adv_c, sc, smu);
while (it_mg <= max_it_CH && resid > tol) {
vcycle(c_new, mu, sc, smu, nx ,ny, 1);
resid = error(ct, c_new, nx, ny);
mat_copy(ct, c_new, 1, nx, 1, ny);
it_mg++;
}
printf("cahn %16.14f %d\n", resid, it_mg-1); }
void source(double **c_old, double **adv_c, double **src_c, double **src_mu) { extern int nx, ny; extern double dt, Pe;
int i, j;
double **lap_c;
lap_c = dmatrix(1, nx, 1, ny);
laplace_ch(c_old, lap_c, nx, ny);
ijloop {
src_c[i][j] = c_old[i][j]/dt - lap_c[i][j]/Pe - adv_c[i][j];
src_mu[i][j] = 0.0;
}
free_dmatrix(lap_c, 1, nx, 1, ny); }
void laplace_ch(double **a, double **lap_a, int nxt, int nyt) { extern double xright;
int i, j;
double ht2, dadx_L, dadx_R, dady_B, dady_T;
ht2 = pow(xright/(double)nxt,2);
ijloopt {
if (i > 1)
dadx_L = a[i][j] - a[i-1][j];
else
dadx_L = a[i][j] - a[nxt][j];
if (i < nxt)
dadx_R = a[i+1][j] - a[i][j];
else
dadx_R = a[1][j] - a[i][j];
if (j > 1)
dady_B = a[i][j] - a[i][j-1];
else
dady_B = 0.0;
if (j < nyt)
dady_T = a[i][j+1] - a[i][j];
else
dady_T = 0.0;
lap_a[i][j] = (dadx_R - dadx_L + dady_T - dady_B)/ht2;
}
}
void vcycle(double **uf_new, double **wf_new, double **su, double **sw, int nxf, int nyf, int ilevel) { extern int n_level;
relax(uf_new, wf_new, su, sw, ilevel, nxf, nyf);
if (ilevel < n_level) {
int nxc, nyc;
double **uc_new, **wc_new, **duc, **dwc,
**uc_def, **wc_def, **uf_def, **wf_def;
nxc = nxf/2, nyc = nyf/2;
uc_new = dmatrix(1, nxc, 1, nyc);
wc_new = dmatrix(1, nxc, 1, nyc);
duc = dmatrix(1, nxc, 1, nyc);
dwc = dmatrix(1, nxc, 1, nyc);
uc_def = dmatrix(1, nxc, 1, nyc);
wc_def = dmatrix(1, nxc, 1, nyc);
uf_def = dmatrix(1, nxf, 1, nyf);
wf_def = dmatrix(1, nxf, 1, nyf);
restrict2(uf_new, uc_new, wf_new, wc_new, nxc, nyc);
defect(duc, dwc, uf_new, wf_new, su, sw, nxf, nyf, uc_new, wc_new, nxc, nyc);
mat_copy2(uc_def, uc_new, wc_def, wc_new, 1, nxc, 1, nyc);
vcycle(uc_def, wc_def, duc, dwc, nxc, nyc, ilevel+1);
mat_sub2(uc_def, uc_def, uc_new, wc_def, wc_def, wc_new, 1, nxc, 1, nyc);
prolong_ch(uc_def, uf_def, wc_def, wf_def, nxc, nyc);
mat_add2(uf_new, uf_new, uf_def, wf_new, wf_new, wf_def, 1, nxf, 1, nyf);
relax(uf_new, wf_new, su, sw, ilevel, nxf, nyf);
free_dmatrix(uc_new, 1, nxc, 1, nyc);
free_dmatrix(wc_new, 1, nxc, 1, nyc);
free_dmatrix(duc, 1, nxc, 1, nyc);
free_dmatrix(dwc, 1, nxc, 1, nyc);
free_dmatrix(uc_def, 1, nxc, 1, nyc);
free_dmatrix(wc_def, 1, nxc, 1, nyc);
free_dmatrix(uf_def, 1, nxf, 1, nyf);
free_dmatrix(wf_def, 1, nxf, 1, nyf);
}
}
void relax(double **c_new, double **mu_new, double **su, double **sw, int ilevel, int nxt, int nyt) { extern int c_relax; extern double xright, dt, Cahn, Pe;
int i, j, iter;
double ht2, a[4], f[2], det;
ht2 = pow(xright/(double)nxt,2);
for (iter=1; iter<=c_relax; iter++) {
ijloopt {
a[0] = 1.0/dt;
a[1] = 2.0/(ht2*Pe);
if (j > 1) a[1] += 1.0/(ht2*Pe);
if (j < nyt) a[1] += 1.0/(ht2*Pe);
a[2] = -d2f(c_new[i][j]) - 2.0*Cahn/ht2;
if (j > 1) a[2] += -Cahn/ht2;
if (j < nyt) a[2] += -Cahn/ht2;
a[3] = 1.0;
f[0] = su[i][j];
if (i > 1) f[0] += mu_new[i-1][j]/(ht2*Pe);
else f[0] += mu_new[nxt][j]/(ht2*Pe);
if (i < nxt) f[0] += mu_new[i+1][j]/(ht2*Pe);
else f[0] += mu_new[1][j]/(ht2*Pe);
if (j > 1) f[0] += mu_new[i][j-1]/(ht2*Pe);
if (j < nyt) f[0] += mu_new[i][j+1]/(ht2*Pe);
f[1] = sw[i][j] + df(c_new[i][j]) - d2f(c_new[i][j])*c_new[i][j];
if (i > 1) f[1] -= Cahn*c_new[i-1][j]/ht2;
else f[1] -= Cahn*c_new[nxt][j]/ht2;
if (i < nxt) f[1] -= Cahn*c_new[i+1][j]/ht2;
else f[1] -= Cahn*c_new[1][j]/ht2;
if (j > 1) f[1] -= Cahn*c_new[i][j-1]/ht2;
if (j < nyt) f[1] -= Cahn*c_new[i][j+1]/ht2;
det = a[0]*a[3] - a[1]*a[2];
c_new[i][j] = (a[3]*f[0] - a[1]*f[1])/det;
mu_new[i][j] = (-a[2]*f[0] + a[0]*f[1])/det;
}
}
}
void defect(double **duc, double **dwc, double **uf_new, double **wf_new, double **suf, double **swf, int nxf, int nyf, double **uc_new, double **wc_new, int nxc, int nyc) { double **ruf, **rwf, **ruc, **rwc, **rruf, **rrwf;
ruf = dmatrix(1, nxf, 1, nyf);
rwf = dmatrix(1, nxf, 1, nyf);
ruc = dmatrix(1, nxc, 1, nyc);
rwc = dmatrix(1, nxc, 1, nyc);
rruf = dmatrix(1, nxc, 1, nyc);
rrwf = dmatrix(1, nxc, 1, nyc);
nonL(ruf, rwf, uf_new, wf_new, nxf, nyf);
nonL(ruc, rwc, uc_new, wc_new, nxc, nyc);
mat_sub2(ruf, suf, ruf, rwf, swf, rwf, 1, nxf, 1, nyf);
restrict2(ruf, rruf, rwf, rrwf, nxc, nyc);
mat_add2(duc, ruc, rruf, dwc, rwc, rrwf, 1, nxc, 1, nyc);
free_dmatrix(ruf, 1, nxf, 1, nyf);
free_dmatrix(rwf, 1, nxf, 1, nyf);
free_dmatrix(ruc, 1, nxc, 1, nyc);
free_dmatrix(rwc, 1, nxc, 1, nyc);
free_dmatrix(rruf, 1, nxc, 1, nyc);
free_dmatrix(rrwf, 1, nxc, 1, nyc); }
void nonL(double **ru, double **rw, double **c_new, double **mu_new, int nxt, int nyt) { extern double dt, Cahn, Pe;
int i, j;
double **lap_c, **lap_mu;
lap_c = dmatrix(1, nxt, 1, nyt);
lap_mu = dmatrix(1, nxt, 1, nyt);
laplace_ch(c_new, lap_c, nxt, nyt);
laplace_ch(mu_new, lap_mu, nxt, nyt);
ijloopt {
ru[i][j] = c_new[i][j]/dt - lap_mu[i][j]/Pe;
rw[i][j] = mu_new[i][j] - df(c_new[i][j]) + Cahn*lap_c[i][j];
}
free_dmatrix(lap_c, 1, nxt, 1, nyt);
free_dmatrix(lap_mu, 1, nxt, 1, nyt); }
void restrict2(double **uf, double **uc, double **vf, double **vc, int nxt, int nyt) { int i, j;
ijloopt {
uc[i][j] = 0.25*(uf[2*i-1][2*j-1] + uf[2*i-1][2*j] + uf[2*i][2*j-1] + uf[2*i][2*j]);
vc[i][j] = 0.25*(vf[2*i-1][2*j-1] + vf[2*i-1][2*j] + vf[2*i][2*j-1] + vf[2*i][2*j]);
}
}
void prolong_ch(double **uc, double **uf, double **vc, double **vf, int nxt, int nyt) { int i, j;
ijloopt {
uf[2*i-1][2*j-1] = uf[2*i-1][2*j] = uf[2*i][2*j-1] = uf[2*i][2*j] = uc[i][j];
vf[2*i-1][2*j-1] = vf[2*i-1][2*j] = vf[2*i][2*j-1] = vf[2*i][2*j] = vc[i][j];
}
}
double error(double **c_old, double **c_new, int nxt, int nyt) { double **r, res;
r = dmatrix(1, nxt, 1, nyt);
mat_sub(r, c_new, c_old, 1, nxt, 1, nyt);
res = mat_max(r, 1, nxt, 1, nyt);
free_dmatrix(r, 1, nxt, 1, nyt);
return res; }
double df(double c) { double value;
value = pow(c,3);
return value; }
double d2f(double c) { double value;
value = 3.0*c*c;
return value; }
/***** util ******/ double * dvector (long nl,long nh) { double * v; v=(double * ) malloc((nh-nl+2) * sizeof(double)); return v-nl+1; }
double **dmatrix(long nrl, long nrh, long ncl, long nch) { double **m; long i, nrow=nrh-nrl+1+NR_END, ncol=nch-ncl+1+NR_END;
m=(double **) malloc((nrow)*sizeof(double*));
m+=NR_END;
m-=nrl;
m[nrl]=(double *) malloc((nrow*ncol)*sizeof(double));
m[nrl]+=NR_END;
m[nrl]-=ncl;
for (i=nrl+1; i<=nrh; i++) m[i]=m[i-1]+ncol;
return m; }
void free_dvector(double * v,long nl,long nh) { free (v+nl-1); }
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch) { free(m[nrl]+ncl-NR_END); free(m+nrl-NR_END);
return; }
double vpn(double a, double b, double c, double d) { double ep, is0, is1, is2, a0, a1, a2, w0, w2; ep = 1.0e-6;
is0 = 13.0pow(a-b,2) + 3.0pow(a-3.0b,2); is1 = 13.0pow(b-c,2) + 3.0pow(b+c,2); is2 = 13.0pow(c-d,2) + 3.0pow(3.0c-d,2);
a0 = 1.0/pow(ep+is0,2); a1 = 6.0/pow(ep+is1,2); a2 = 3.0/pow(ep+is2,2);
w0 = a0/(a0+a1+a2); w2 = a2/(a0+a1+a2);
return w0(a-2.0b+c)/3.0+(w2-0.5)(b-2.0c+d)/6.0;
}
void matmult(double **a, double **b, double **c, int nrl, int nrh, int ncl, int nch) { int i, j;
for (i=nrl; i<=nrh; i++) for (j=ncl; j<=nch; j++) a[i][j]=b[i][j] * c[i][j]; }
void mat_comb(double **a, double s1, double **b, double s2, double **c, int nrl, int nrh, int ncl, int nch) { int i,j;
for (i=nrl; i<=nrh; i++) for (j=ncl; j<=nch; j++)
a[i][j] = s1 * b[i][j] + s2 * c[i][j];
return; }
void zero_matrix(double **a, int xl, int xr, int yl, int yr) { int i, j;
for (i=xl; i<=xr; i++) {
for (j=yl; j<=yr; j++) {
a[i][j] = 0.0;
}
}
return; }
void mat_copy(double **a, double **b, int xl, int xr, int yl, int yr) { int i, j;
for (i=xl; i<=xr; i++) {
for (j=yl; j<=yr; j++) {
a[i][j] = b[i][j];
}
}
return; }
void mat_copy2(double **a, double **b, double **a2, double **b2, int xl, int xr, int yl, int yr) { int i, j;
for (i=xl; i<=xr; i++) {
for (j=yl; j<=yr; j++) {
a[i][j] = b[i][j];
a2[i][j] = b2[i][j];
}
}
return; }
void mat_add(double **a, double **b, double **c, int xl, int xr, int yl, int yr) { int i, j;
for (i=xl; i<=xr; i++) {
for (j=yl; j<=yr; j++) {
a[i][j] = b[i][j]+c[i][j];
}
}
return; }
void mat_add2(double **a, double **b, double **c, double **a2, double **b2, double **c2, int xl, int xr, int yl, int yr) { int i, j;
for (i=xl; i<=xr; i++) {
for (j=yl; j<=yr; j++) {
a[i][j] = b[i][j]+c[i][j];
a2[i][j] = b2[i][j]+c2[i][j];
}
}
return; }
void mat_sub(double **a, double **b, double **c, int nrl, int nrh, int ncl, int nch) { int i, j;
for (i=nrl; i<=nrh; i++) {
for (j=ncl; j<=nch; j++) {
a[i][j] = b[i][j]-c[i][j];
}
}
return; }
void mat_sub2(double **a, double **b, double **c, double **a2, double **b2, double **c2, int nrl, int nrh, int ncl, int nch) { int i, j;
for (i=nrl; i<=nrh; i++) {
for (j=ncl; j<=nch; j++) {
a[i][j] = b[i][j]-c[i][j];
a2[i][j] = b2[i][j]-c2[i][j];
}
}
return; }
double mat_max(double **a, int nrl, int nrh, int ncl, int nch) { int i, j; double x = 0.0;
for (i=nrl; i<=nrh; i++) {
for (j=ncl; j<=nch; j++) {
if (fabs(a[i][j]) > x)
x = fabs(a[i][j]);
}
}
return x; }
void print_mat(FILE *fptr, double **a, int nrl, int nrh, int ncl, int nch) { int i, j;
for (i=nrl; i<=nrh; i++) {
for (j=ncl; j<=nch; j++)
fprintf(fptr, " %16.14f", a[i][j]);
fprintf(fptr, "\n");
}
return; }
void print_data(double **u, double **v, double **c, double **p) { extern char bufferu[100], bufferv[100],bufferp[100],bufferc[100]; int i, j; FILE *fu, *fv, *fp, *fc; fu = fopen(bufferu,”a”); fv = fopen(bufferv,”a”); fp = fopen(bufferp,”a”); fc = fopen(bufferc,”a”);
iloop {
jloop {
fprintf(fu, " %16.14f", 0.5*(u[i][j]+u[i-1][j]));
fprintf(fv, " %16.14f", 0.5*(v[i][j]+v[i][j-1]));
fprintf(fp, " %16.14f", p[i][j]);
fprintf(fc, " %16.14f", c[i][j]);
}
fprintf(fv, "\n");
fprintf(fp, "\n");
fprintf(fc, "\n");
}
fclose(fu);
fclose(fv);
fclose(fp);
fclose(fc);
return; }
void pressure_update(double **a) { extern int nx, ny;
int i, j;
double ave = 0.0;
ijloop {
ave = ave + a[i][j];
}
ave /= (nx+0.0)*(ny+0.0);
ijloop {
a[i][j] -= ave;
}
return; }
/** The following util.h file is associated with C code. ***/
#define gnx 512 #define gny 256 #define iloop for(i=1;i<=gnx;i++) #define i0loop for(i=0;i<=gnx;i++) #define jloop for(j=1;j<=gny;j++) #define j0loop for(j=0;j<=gny;j++) #define ijloop iloop jloop #define i0jloop i0loop jloop #define ij0loop iloop j0loop #define iloopt for(i=1;i<=nxt;i++) #define i0loopt for(i=0;i<=nxt;i++) #define jloopt for(j=1;j<=nyt;j++) #define j0loopt for(j=0;j<=nyt;j++) #define ijloopt iloopt jloopt #define i0jloopt i0loopt jloopt #define ij0loopt iloopt j0loopt double * dvector (long nl,long nh); double **dmatrix(long nrl, long nrh, long ncl, long nch); void free_dvector(double * v,long nl,long nh); void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch); double vpn(double a, double b, double c, double d); void matmult(double **a, double **b, double **c, int nrl, int nrh, int ncl, int nch); void mat_comb(double **a, double s1, double **b, double s2, double **c, int nrl, int nrh, int ncl, int nch); void zero_matrix(double **a, int xl, int xr, int yl, int yr); void mat_copy(double **a, double **b, int xl, int xr, int yl, int yr); void mat_copy2(double **a, double **b, double **a2, double **b2, int xl, int xr, int yl, int yr); void mat_add(double **a, double **b, double **c, int xl, int xr, int yl, int yr); void mat_add2(double **a, double **b, double **c, double **a2, double **b2, double **c2, int xl, int xr, int yl, int yr); void mat_sub(double **a, double **b, double **c, int nrl, int nrh, int ncl, int nch); void mat_sub2(double **a, double **b, double **c, double **a2, double **b2, double **c2, int nrl, int nrh, int ncl, int nch); double mat_max(double **a, int nrl, int nrh, int ncl, int nch); void print_mat(FILE *fptr, double **a, int nrl, int nrh, int ncl, int nch); void print_data(double **u, double **v, double **c, double **p); void pressure_update(double **a);
/** The following bnsch.h file is associated with C code. **/
void initialization(double **u, double **v, double **c, double **p); void cal_den(double **c, double **rho); void cal_vis(double **c, double **vis); void full_step(double **u, double **v, double **c, double **nu, double **nv, double **nc, double **p); void advection_step(double **u, double **v, double **c, double **adv_u, double **adv_v, double **adv_c); void advection_uv(double **u, double **v, double **adv_u, double **adv_v); void augmenuv(double **u, double **v); void augmenc(double **c, int nxt, int nyt); void temp_uv(double **tu, double **tv, double **u, double **v, double **adv_u, double **adv_v, double **fx, double **fy); void Poisson(double **tu, double **tv, double **p); void source_uv(double **tu, double **tv, double **divuv, int nxt, int nyt); void div_uv(double **tu, double **tv, double **divuv, int nxt, int nyt); void MG_Poisson(double **p, double **f); void vcycle_uv(double **uf, double **ff, double **wf, int nxf, int nyf, int ilevel); void relax_p(double **p, double **f, double **w, int ilevel, int nxt, int nyt); void residual_den(double **r, double **u, double **f, double **den, int nxt, int nyt); void grad_p(double **p, double **dpdx, double **dpdy, int nxt, int nyt); void restrict1(double **u_fine, double **u_coarse, int nxt, int nyt); void prolong(double **u_coarse, double **u_fine, int nxt, int nyt); void cahn(double **c_old, double **adv_c, double **c_new); void source(double **c_old, double **adv_c, double **src_c, double **src_mu); void laplace_ch(double **a, double **lap_a, int nxt, int nyt); void vcycle(double **uf_new, double **wf_new, double **su, double **sw, int nxf, int nyf, int ilevel); void relax(double **c_new, double **mu_new, double **su, double **sw, int ilevel, int nxt, int nyt); void defect(double **duc, double **dwc, double **uf_new, double **wf_new, double **suf, double **swf, int nxf, int nyf, double **uc_new, double **wc_new, int nxc, int nyc); void nonL(double **ru, double **rw, double **c_new, double **mu_new, int nxt, int nyt); void restrict2(double **uf, double **uc, double **vf, double **vc, int nxt, int nyt); void prolong_ch(double **uc, double **uf, double **vc, double **vf, int nxt, int nyt); double error(double **c_old, double **c_new, int nxt, int nyt); double df(double c); double d2f(double c);
%The following MATLAB code is used to show the results.
clear all; xright=2.0; yright=1.0; xleft = 0; yleft = 0; nx=512; ny = 256;h= (xright-xleft)/nx; x=linspace(xleft+0.5h,xright-0.5h,nx); y=linspace(yleft+0.5h,yright-0.5h,ny); [xx,yy]=meshgrid(x,y); ss=sprintf(‘…/datacr.m’); phi = load(ss); for i = 1:41 clf; A = phi((i-1)nx+1:inx,:); surf(xx,yy,-A’); shading interp; colormap bone; view(0,90); axis image; axis([0 xright 0 yright]); axis off; pause(0.3); end