Merge branch 'VariableDensity' of https://github.com/JamesEMcClure/LBPM-WIA into VariableDensity

This commit is contained in:
James E McClure
2017-11-27 07:35:25 -05:00
4 changed files with 152 additions and 35 deletions

View File

@@ -103,10 +103,12 @@ extern "C" void ScaLBL_Color_InitDistance(char *ID, double *Den, double *Phi, do
//*************************************************************************
extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, double *A_odd,
extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz)
{
int i,j,k,n,N;
double na,nb,ux,uy,uz;
N = Nx*Ny*Nz;
// Fill the inlet with component a
for (k=0; k<1; k++){
@@ -126,6 +128,7 @@ extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, doub
Den[n] = 1.0;
Den[N+n] = 0.0;
/*
A_even[n] = 0.3333333333333333;
A_odd[n] = 0.1111111111111111;
A_even[N+n] = 0.1111111111111111;
@@ -135,21 +138,80 @@ extern "C" void ScaLBL_Color_BC_z(double *Phi, double *Den, double *A_even, doub
A_even[3*N+n] = 0.1111111111111111;
B_even[n] = 0.0;
B_odd[n] = 0.0;
B_odd[n] = 0.0;
B_even[N+n] = 0.0;
B_odd[N+n] = 0.0;
B_odd[N+n] = 0.0;
B_even[2*N+n] = 0.0;
B_odd[2*N+n] = 0.0;
B_even[3*N+n] = 0.0;
*/
na = 1.0;
nb = 0.0;
//....Load the flow velocity...........
ux = Velocity[n];
uy = Velocity[N+n];
uz = Velocity[2*N+n];
A_even[n] = 0.3333333333333333*na;
B_even[n] = 0.3333333333333333*nb;
// Non-Stationary equilibrium distributions
//feq[0] = 0.1111111111111111*(1+4.5*ux);
//feq[1] = 0.1111111111111111*(1-4.5*ux);
//feq[2] = 0.1111111111111111*(1+4.5*uy);
//feq[3] = 0.1111111111111111*(1-4.5*uy);
//feq[4] = 0.1111111111111111*(1+4.5*uz);
//feq[5] = 0.1111111111111111*(1-4.5*uz);
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*ux));
b1 = nb*(0.1111111111111111*(1+4.5*ux));
a2 = na*(0.1111111111111111*(1-4.5*ux));
b2 = nb*(0.1111111111111111*(1-4.5*ux));
A_odd[n] = a1;
A_even[N+n] = a2;
B_odd[n] = b1;
B_even[N+n] = b2;
//...............................................
// q = 2
// Cq = {0,1,0}
a1 = na*(0.1111111111111111*(1+4.5*uy));
b1 = nb*(0.1111111111111111*(1+4.5*uy));
a2 = na*(0.1111111111111111*(1-4.5*uy));
b2 = nb*(0.1111111111111111*(1-4.5*uy));
A_odd[N+n] = a1;
A_even[2*N+n] = a2;
B_odd[N+n] = b1;
B_even[2*N+n] = b2;
//...............................................
// q = 4
// Cq = {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*uz));
b1 = nb*(0.1111111111111111*(1+4.5*uz));
a2 = na*(0.1111111111111111*(1-4.5*uz));
b2 = nb*(0.1111111111111111*(1-4.5*uz));
A_odd[2*N+n] = a1;
A_even[3*N+n] = a2;
B_odd[2*N+n] = b1;
B_even[3*N+n] = b2;
}
}
}
}
//*************************************************************************
extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, double *A_odd,
extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *Velocity, double *A_even, double *A_odd,
double *B_even, double *B_odd, int Nx, int Ny, int Nz)
{
int i,j,k,n,N;
double na,nb,ux,uy,uz;
N = Nx*Ny*Nz;
// Fill the outlet with component b
for (k=Nz-3; k<Nz-1; k++){
@@ -161,7 +223,7 @@ extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, doub
Den[n] = 0.0;
Den[N+n] = 1.0;
A_even[n] = 0.0;
/* A_even[n] = 0.0;
A_odd[n] = 0.0;
A_even[N+n] = 0.0;
A_odd[N+n] = 0.0;
@@ -176,6 +238,61 @@ extern "C" void ScaLBL_Color_BC_Z(double *Phi, double *Den, double *A_even, doub
B_even[2*N+n] = 0.1111111111111111;
B_odd[2*N+n] = 0.1111111111111111;
B_even[3*N+n] = 0.1111111111111111;
*/
na = 0.0;
nb = 1.0;
//....Load the flow velocity...........
ux = Velocity[n];
uy = Velocity[N+n];
uz = Velocity[2*N+n];
A_even[n] = 0.3333333333333333*na;
B_even[n] = 0.3333333333333333*nb;
// Non-Stationary equilibrium distributions
//feq[0] = 0.1111111111111111*(1+4.5*ux);
//feq[1] = 0.1111111111111111*(1-4.5*ux);
//feq[2] = 0.1111111111111111*(1+4.5*uy);
//feq[3] = 0.1111111111111111*(1-4.5*uy);
//feq[4] = 0.1111111111111111*(1+4.5*uz);
//feq[5] = 0.1111111111111111*(1-4.5*uz);
//...............................................
// q = 0,2,4
// Cq = {1,0,0}, {0,1,0}, {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*ux));
b1 = nb*(0.1111111111111111*(1+4.5*ux));
a2 = na*(0.1111111111111111*(1-4.5*ux));
b2 = nb*(0.1111111111111111*(1-4.5*ux));
A_odd[n] = a1;
A_even[N+n] = a2;
B_odd[n] = b1;
B_even[N+n] = b2;
//...............................................
// q = 2
// Cq = {0,1,0}
a1 = na*(0.1111111111111111*(1+4.5*uy));
b1 = nb*(0.1111111111111111*(1+4.5*uy));
a2 = na*(0.1111111111111111*(1-4.5*uy));
b2 = nb*(0.1111111111111111*(1-4.5*uy));
A_odd[N+n] = a1;
A_even[2*N+n] = a2;
B_odd[N+n] = b1;
B_even[2*N+n] = b2;
//...............................................
// q = 4
// Cq = {0,0,1}
a1 = na*(0.1111111111111111*(1+4.5*uz));
b1 = nb*(0.1111111111111111*(1+4.5*uz));
a2 = na*(0.1111111111111111*(1-4.5*uz));
b2 = nb*(0.1111111111111111*(1-4.5*uz));
A_odd[2*N+n] = a1;
A_even[3*N+n] = a2;
B_odd[2*N+n] = b1;
B_even[3*N+n] = b2;
}
}