save the work;to be built and tested

This commit is contained in:
Rex Zhe Li 2021-01-31 19:06:07 -05:00
parent 4085deb5e3
commit e22de8ae7e
3 changed files with 1408 additions and 878 deletions

File diff suppressed because it is too large Load Diff

View File

@ -10,7 +10,7 @@ color lattice boltzmann model
#include <time.h>
ScaLBL_FreeLeeModel::ScaLBL_FreeLeeModel(int RANK, int NP, const Utilities::MPI& COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),W(0),gamma(0),
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),tauM(0),rhoA(0),rhoB(0),W(0),gamma(0),kappa(0),beta(0),
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
@ -30,10 +30,13 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
// set defaults
timestepMax = 100000;
tauA = tauB = 1.0;
tauM = 1.0;//relaxation time for phase field
rhoA = rhoB = 1.0;
Fx = Fy = Fz = 0.0;
gamma=1e-3;//surface tension
W=5.0;//interfacial thickness
beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
Restart=false;
din=dout=1.0;
flux=0.0;
@ -81,6 +84,9 @@ void ScaLBL_FreeLeeModel::ReadParams(string filename){
inletB=0.f;
outletA=0.f;
outletB=1.f;
//update secondary parameters
beta = 12.0*gamma/W;
kappa = 3.0*gamma*W/2.0;//beta and kappa are related to surface tension \gamma
//if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
BoundaryCondition = 0;
@ -258,15 +264,16 @@ void ScaLBL_FreeLeeModel::Create(){
TmpMap[idx] = Nxh*Nyh*Nzh-1;
}
}
// copy the device map
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
delete [] TmpMap;
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
comm.barrier();
delete [] TmpMap;
delete [] neighborList;
}
void ScaLBL_FreeLeeModel::AssignComponentLabels()
void ScaLBL_FreeLeeModel::AssignComponentLabels_ChemPotential_ColorGrad()
{
double *phase;
phase = new double[Nh];
@ -288,12 +295,26 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels()
// Assign the labels
for (size_t idx=0; idx<NLABELS; idx++) label_count[idx]=0;
for (int k=1;k<Nzh-1;k++){
for (int j=1;j<Nyh-1;j++){
for (int i=1;i<Nxh-1;i++){
int n = (k-1)*Nx*Ny+(j-1)*Nx+(i-1);
for (int k=0;k<Nzh;k++){
for (int j=0;j<Nyh;j++){
for (int i=0;i<Nxh;i++){
//idx for double-halo array 'phase'
int nh = k*Nxh*Nyh+j*Nxh+i;
//idx for single-halo array Mask->id[n]
int x=i-1;
int y=j-1;
int z=k-1;
if (x<0) x=0;
if (y<0) y=0;
if (z<0) z=0;
if (x>=Nx) x=Nx-1;
if (y>=Ny) y=Ny-1;
if (z>=Nz) z=Nz-1;
int n = z*Nx*Ny+y*Nx+x;
VALUE=id[n];
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
@ -307,7 +328,7 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels()
// fluid labels are reserved
if (VALUE == 1) AFFINITY=1.0;
else if (VALUE == 2) AFFINITY=-1.0;
phase[n] = AFFINITY;
phase[nh] = AFFINITY;
}
}
}
@ -329,56 +350,10 @@ void ScaLBL_FreeLeeModel::AssignComponentLabels()
}
//compute color gradient and laplacian of phase field
double *ColorGrad_host, mu_phi_host;
ColorGrad_host = new double[3*Np];
mu_phi_host = new double[Np];
//copy all data to device
ScaLBL_CopyToDevice(Phi, phase, N*sizeof(double));
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
delete [] phase;
}
void ScaLBL_FreeLeeModel::AssignChemPotential_ColorGrad()
{
double *SolidPotential_host = new double [Nx*Ny*Nz];
double *GreySolidGrad_host = new double [3*Np];
size_t NLABELS=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = greyscaleColor_db->getVector<int>( "GreySolidLabels" );
auto AffinityList = greyscaleColor_db->getVector<double>( "GreySolidAffinity" );
NLABELS=LabelList.size();
if (NLABELS != AffinityList.size()){
ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n");
}
for (int k=0;k<Nz;k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
VALUE=id[n];
AFFINITY=0.f;//all nodes except the specified grey nodes have grey-solid affinity = 0.0
// Assign the affinity from the paired list
for (unsigned int idx=0; idx < NLABELS; idx++){
//printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]);
if (VALUE == LabelList[idx]){
AFFINITY=AffinityList[idx];
idx = NLABELS;
//Mask->id[n] = 0; // set mask to zero since this is an immobile component
}
}
SolidPotential_host[n] = AFFINITY;
}
}
}
// Calculate grey-solid color-gradient
double *Dst;
Dst = new double [3*3*3];
for (int kk=0; kk<3; kk++){
@ -389,8 +364,8 @@ void ScaLBL_FreeLeeModel::AssignChemPotential_ColorGrad()
}
}
}
double w_face = 1.f;
double w_edge = 0.5;
double w_face = 1.0/18.0;
double w_edge = 1.0/36.0;
double w_corner = 0.f;
//local
Dst[13] = 0.f;
@ -424,14 +399,21 @@ void ScaLBL_FreeLeeModel::AssignChemPotential_ColorGrad()
Dst[23] = w_edge;
Dst[25] = w_edge;
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
int idx=Map(i,j,k);
double cs2_inv = 3.0;//inverse of c_s^2 for D3Q19 lattice
int width = 2;//For better readability: make halo width explicity wherever possible
for (int k=width; k<Nzh-width; k++){
for (int j=width; j<Nyh-width; j++){
for (int i=width; i<Nxh-width; i++){
//idx for double-halo array 'phase'
int nh = k*Nxh*Nyh+j*Nxh+i;
int idx=Map(i-width+1,j-width+1,k-width+1);
if (!(idx < 0)){
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
double phi_x = 0.f;
double phi_y = 0.f;
double phi_z = 0.f;
double phi_Lap = 0.f;//Laplacian of the phase field
for (int kk=0; kk<3; kk++){
for (int jj=0; jj<3; jj++){
for (int ii=0; ii<3; ii++){
@ -446,51 +428,43 @@ void ScaLBL_FreeLeeModel::AssignChemPotential_ColorGrad()
if (idi < 0) idi=0;
if (idj < 0) idj=0;
if (idk < 0) idk=0;
if (!(idi < Nx)) idi=Nx-1;
if (!(idj < Ny)) idj=Ny-1;
if (!(idk < Nz)) idk=Nz-1;
if (!(idi < Nxh)) idi=Nxh-1;
if (!(idj < Nyh)) idj=Nyh-1;
if (!(idk < Nzh)) idk=Nzh-1;
int nn = idk*Nx*Ny + idj*Nx + idi;
int nn = idk*Nxh*Nyh + idj*Nxh + idi;
double vec_x = double(ii-1);
double vec_y = double(jj-1);
double vec_z = double(kk-1);
double GWNS=SolidPotential_host[nn];
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
double GWNS=phase[nn];
double GWNS_local=phase[nh];
phi_x += GWNS*weight*vec_x;
phi_y += GWNS*weight*vec_y;
phi_z += GWNS*weight*vec_z;
phi_Lap += weight*(GWNS-GWNS_local);//Laplacian of the phase field
}
}
}
if (Averages->SDs(i,j,k)<2.0){
GreySolidGrad_host[idx+0*Np] = phi_x;
GreySolidGrad_host[idx+1*Np] = phi_y;
GreySolidGrad_host[idx+2*Np] = phi_z;
}
else{
GreySolidGrad_host[idx+0*Np] = 0.0;
GreySolidGrad_host[idx+1*Np] = 0.0;
GreySolidGrad_host[idx+2*Np] = 0.0;
}
//store color gradient
ColorGrad_host[idx+0*Np] = cs2_inv*phi_x;
ColorGrad_host[idx+1*Np] = cs2_inv*phi_y;
ColorGrad_host[idx+2*Np] = cs2_inv*phi_z;
//compute chemical potential
phi_Lap = 2.0*cs2_inv*phi_Lap;
mu_phi_host[idx] = 4.0*beta*phase[nh]*(phase[nh]+1.0)*(phase[nh]-1.0) - kappa*phi_Lap;
}
}
}
}
if (rank==0){
printf("Number of Grey-solid labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
printf(" grey-solid label=%d, grey-solid affinity=%f\n",VALUE,AFFINITY);
}
}
ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
ScaLBL_DeviceBarrier();
delete [] SolidPotential_host;
delete [] GreySolidGrad_host;
//copy all data to device
ScaLBL_CopyToDevice(Phi, phase, Nh*sizeof(double));
ScaLBL_CopyToDevice(ColorGrad, ColorGrad_host, 3*Np*sizeof(double));
ScaLBL_CopyToDevice(mu_phi, mu_phi_host, Np*sizeof(double));
comm.barrier();
delete [] phase;
delete [] ColorGrad_host;
delete [] mu_phi_host;
delete [] Dst;
}
@ -501,15 +475,12 @@ void ScaLBL_FreeLeeModel::Initialize(){
if (rank==0) printf ("Initializing phase field, chemical potential and color gradient\n");
AssignComponentLabels_ChemPotential_ColorGrad();//initialize phase field Phi
//if (rank==0) printf ("Initializing chemical potential and color gradient \n");
//AssignChemPotential_ColorGrad();
if (rank==0) printf ("Initializing distributions for momentum transport\n");
ScaLBL_D3Q19_FreeLeeModel_Init(gqbar, mu_phi, ColorGrad, Fx, Fy, Fz, Np);
if (rank==0) printf ("Initializing density field and distributions for phase-field transport\n");
ScaLBL_FreeLeeModel_PhaseField_Init(dvcMap, Phi, Den, hq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_FreeLeeModel_PhaseField_Init(dvcMap, Phi, Den, hq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_FreeLeeModel_PhaseField_Init(dvcMap, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_FreeLeeModel_PhaseField_Init(dvcMap, Phi, Den, hq, ColorGrad, rhoA, rhoB, tauM, W, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
if (Restart == true){
//TODO need to revise this function
@ -576,6 +547,7 @@ void ScaLBL_FreeLeeModel::Initialize(){
}
// establish reservoirs for external bC
// TODO to be revised
if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4 ){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
@ -617,13 +589,12 @@ void ScaLBL_FreeLeeModel::Run(){
timestep++;
//-------------------------------------------------------------------------------------------------------------------
// Compute the Phase indicator field
// Read for hq, Bq happens in this routine (requires communication)
//ScaLBL_Comm->SendD3Q7AA(hq); //READ FROM NORMAL
ScaLBL_Comm->SendD3Q7AA(hq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq); //WRITE INTO OPPOSITE
// Read for hq happens in this routine (requires communication)
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, hq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAodd_FreeLeeModel_PhaseField(NeighborList, dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FROM NORMAL
@ -635,8 +606,8 @@ void ScaLBL_FreeLeeModel::Run(){
// Halo exchange for phase field
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
@ -653,19 +624,19 @@ void ScaLBL_FreeLeeModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAodd_FreeLeeModel(NeighborList, dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
// *************EVEN TIMESTEP*************
timestep++;
// Compute the Phase indicator field
ScaLBL_Comm->SendD3Q7AA(hq); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, hq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq); //WRITE INTO OPPOSITE
ScaLBL_Comm->SendD3Q7AA(hq,0); //READ FROM NORMAL
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(hq,0); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, hq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q7_AAeven_FreeLeeModel_PhaseField(dvcMap, hq, Den, Phi, rhoA, rhoB, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
ScaLBL_Comm->SendD3Q19AA(gqbar); //READ FORM NORMAL
@ -675,8 +646,8 @@ void ScaLBL_FreeLeeModel::Run(){
ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
}
ScaLBL_Comm_WideHalo->Send(Phi);
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_WideHalo->Recv(Phi);
ScaLBL_Comm->RecvD3Q19AA(gqbar); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
@ -693,8 +664,8 @@ void ScaLBL_FreeLeeModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_z(fq);
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, hq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB,
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_D3Q19_AAeven_FreeLeeModel(dvcMap, gqbar, hq, Den, Phi, mu_phi, Velocity, Pressure, rhoA, rhoB, tauA, tauB, tauM,
kappa, beta, W, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_Comm->Barrier();
//************************************************************************
PROFILE_STOP("Update");
@ -722,30 +693,24 @@ void ScaLBL_FreeLeeModel::Run(){
void ScaLBL_FreeLeeModel::WriteDebug(){
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);
DoubleArray PhaseData(Nxh,Nyh,Nzh);
//ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField);
ScaLBL_CopyToHost(PhaseField.data(), Phi, sizeof(double)*N);
ScaLBL_CopyToHost(PhaseData.data(), Phi, sizeof(double)*Nh);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Phase.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,OUTFILE);
fwrite(PhaseData.data(),8,Nh,OUTFILE);
fclose(OUTFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField);
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,Den,PhaseField);
FILE *AFILE;
sprintf(LocalRankFilename,"A.%05i.raw",rank);
sprintf(LocalRankFilename,"Density.%05i.raw",rank);
AFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,AFILE);
fclose(AFILE);
ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField);
FILE *BFILE;
sprintf(LocalRankFilename,"B.%05i.raw",rank);
BFILE = fopen(LocalRankFilename,"wb");
fwrite(PhaseField.data(),8,N,BFILE);
fclose(BFILE);
ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField);
FILE *PFILE;
sprintf(LocalRankFilename,"Pressure.%05i.raw",rank);

View File

@ -35,7 +35,8 @@ public:
int timestep,timestepMax;
int BoundaryCondition;
double tauA,tauB,rhoA,rhoB;
double W,gamma;
double tauM;//relaxation time for phase field (or mass)
double W,gamma,kappa,beta;
double Fx,Fy,Fz,flux;
double din,dout,inletA,inletB,outletA,outletB;