Refactoring color simulator
This commit is contained in:
James E McClure
parent
bcd54408d2
commit
c5477b9b39
@ -16,83 +16,14 @@
|
||||
|
||||
/*
|
||||
* Simulator for two-phase flow in porous media
|
||||
* James E. McClure 2013-2016
|
||||
* James E. McClure 2013-2018
|
||||
*/
|
||||
|
||||
using namespace std;
|
||||
|
||||
//*************************************************************************
|
||||
// Implementation of Two-Phase Immiscible LBM using CUDA
|
||||
// Implementation of Two-Phase Immiscible LBM
|
||||
//*************************************************************************
|
||||
inline void PackID(int *list, int count, char *sendbuf, char *ID){
|
||||
// Fill in the phase ID values from neighboring processors
|
||||
// This packs up the values that need to be sent from one processor to another
|
||||
int idx,n;
|
||||
|
||||
for (idx=0; idx<count; idx++){
|
||||
n = list[idx];
|
||||
sendbuf[idx] = ID[n];
|
||||
}
|
||||
}
|
||||
//***************************************************************************************
|
||||
|
||||
inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
|
||||
// Fill in the phase ID values from neighboring processors
|
||||
// This unpacks the values once they have been recieved from neighbors
|
||||
int idx,n;
|
||||
|
||||
for (idx=0; idx<count; idx++){
|
||||
n = list[idx];
|
||||
ID[n] = recvbuf[idx];
|
||||
}
|
||||
}
|
||||
|
||||
//***************************************************************************************
|
||||
|
||||
inline void ZeroHalo(double *Data, int Nx, int Ny, int Nz)
|
||||
{
|
||||
int i,j,k,n;
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
i=0;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
i=Nx-1;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
for (k=0;k<Nz;k++){
|
||||
for (i=0;i<Nx;i++){
|
||||
j=0;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
j=Ny-1;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
k=0;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
k=Nz-1;
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
Data[2*n] = 0.0;
|
||||
Data[2*n+1] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
//***************************************************************************************
|
||||
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
@ -129,9 +60,7 @@ int main(int argc, char **argv)
|
||||
|
||||
// Variables that specify the computational domain
|
||||
string FILENAME;
|
||||
unsigned int nBlocks, nthreads;
|
||||
int Nx,Ny,Nz; // local sub-domain size
|
||||
int nspheres; // number of spheres in the packing
|
||||
double Lx,Ly,Lz; // Domain length
|
||||
double D = 1.0; // reference length for non-dimensionalization
|
||||
// Color Model parameters
|
||||
@ -139,18 +68,19 @@ int main(int argc, char **argv)
|
||||
double tau1, tau2, rho1,rho2;
|
||||
double Fx,Fy,Fz,tol,err;
|
||||
double alpha, beta;
|
||||
double das, dbs, phi_s;
|
||||
double din,dout;
|
||||
double wp_saturation;
|
||||
int BoundaryCondition;
|
||||
int InitialCondition;
|
||||
// bool pBC,Restart;
|
||||
// bool pBC,Restart;
|
||||
int i,j,k;
|
||||
|
||||
// pmmc threshold values
|
||||
//double fluid_isovalue,solid_isovalue;
|
||||
//fluid_isovalue = 0.0;
|
||||
//solid_isovalue = 0.0;
|
||||
double din, dout, flux;
|
||||
double inletA,inletB,outletA,outletB;
|
||||
inletA=1.f;
|
||||
inletB=0.f;
|
||||
outletA=0.f;
|
||||
outletB=1.f;
|
||||
typeBC=4;
|
||||
flux = 10.f;
|
||||
dout=1.f;
|
||||
|
||||
int RESTART_INTERVAL=20000;
|
||||
//int ANALYSIS_)INTERVAL=1000;
|
||||
@ -170,10 +100,7 @@ int main(int argc, char **argv)
|
||||
input >> rho2; // density wetting
|
||||
input >> alpha; // Surface Tension parameter
|
||||
input >> beta; // Width of the interface
|
||||
input >> phi_s; // value of phi at the solid surface
|
||||
// Line 2: wetting phase saturation to initialize
|
||||
input >> wp_saturation;
|
||||
// Line 3: External force components (Fx,Fy, Fz)
|
||||
// Line 2: External force components (Fx,Fy, Fz)
|
||||
input >> Fx;
|
||||
input >> Fy;
|
||||
input >> Fz;
|
||||
@ -210,7 +137,6 @@ int main(int argc, char **argv)
|
||||
//.......................................................................
|
||||
ifstream domain("Domain.in");
|
||||
if (input.is_open()){
|
||||
|
||||
domain >> nprocx;
|
||||
domain >> nprocy;
|
||||
domain >> nprocz;
|
||||
@ -243,10 +169,6 @@ int main(int argc, char **argv)
|
||||
MPI_Bcast(&rho2,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&alpha,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&beta,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&das,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&dbs,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&phi_s,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&wp_saturation,1,MPI_DOUBLE,0,comm);
|
||||
MPI_Bcast(&BoundaryCondition,1,MPI_INT,0,comm);
|
||||
MPI_Bcast(&InitialCondition,1,MPI_INT,0,comm);
|
||||
MPI_Bcast(&din,1,MPI_DOUBLE,0,comm);
|
||||
@ -278,15 +200,6 @@ int main(int argc, char **argv)
|
||||
|
||||
MPI_Barrier(comm);
|
||||
|
||||
// **************************************************************
|
||||
// **************************************************************
|
||||
double Ps = -(das-dbs)/(das+dbs);
|
||||
//double xIntPos = log((1.0+phi_s)/(1.0-phi_s))/(2.0*beta);
|
||||
|
||||
// Set the density values inside the solid based on the input value phi_s
|
||||
das = (phi_s+1.0)*0.5;
|
||||
dbs = 1.0 - das;
|
||||
|
||||
if (nprocs != nprocx*nprocy*nprocz){
|
||||
printf("nprocx = %i \n",nprocx);
|
||||
printf("nprocy = %i \n",nprocy);
|
||||
@ -302,8 +215,6 @@ int main(int argc, char **argv)
|
||||
printf("density (wetting) = %f \n", rho2);
|
||||
printf("alpha = %f \n", alpha);
|
||||
printf("beta = %f \n", beta);
|
||||
printf("das = %f \n", das);
|
||||
printf("dbs = %f \n", dbs);
|
||||
printf("gamma_{wn} = %f \n", 5.796*alpha);
|
||||
printf("Force(x) = %f \n", Fx);
|
||||
printf("Force(y) = %f \n", Fy);
|
||||
@ -342,12 +253,9 @@ int main(int argc, char **argv)
|
||||
|
||||
// Mask that excludes the solid phase
|
||||
Domain Mask(Nx,Ny,Nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
|
||||
|
||||
MPI_Barrier(comm);
|
||||
|
||||
Nx+=2; Ny+=2; Nz += 2;
|
||||
//Nx = Ny = Nz; // Cubic domain
|
||||
|
||||
int N = Nx*Ny*Nz;
|
||||
int dist_mem_size = N*sizeof(double);
|
||||
|
||||
@ -365,9 +273,9 @@ int main(int argc, char **argv)
|
||||
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
|
||||
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
|
||||
|
||||
// printf("Local File Name = %s \n",LocalRankFilename);
|
||||
// printf("Local File Name = %s \n",LocalRankFilename);
|
||||
// .......... READ THE INPUT FILE .......................................
|
||||
// char value;
|
||||
// char value;
|
||||
char *id;
|
||||
id = new char[N];
|
||||
double sum_local;
|
||||
@ -378,11 +286,9 @@ int main(int argc, char **argv)
|
||||
if (rank == 0) cout << "Reading in domain from signed distance function..." << endl;
|
||||
|
||||
//.......................................................................
|
||||
// Read the signed distance
|
||||
sprintf(LocalRankString,"%05d",rank);
|
||||
// sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
|
||||
// WriteLocalSolidID(LocalRankFilename, id, N);
|
||||
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
|
||||
// ReadBinaryFile(LocalRankFilename, Averages.SDs.data(), N);
|
||||
ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
|
||||
MPI_Barrier(comm);
|
||||
if (rank == 0) cout << "Domain set." << endl;
|
||||
@ -425,140 +331,55 @@ int main(int argc, char **argv)
|
||||
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank);
|
||||
fclose(IDFILE);
|
||||
|
||||
// Calculate the actual inlet cross-sectional area
|
||||
double inlet_area_local=0.f;
|
||||
double InletArea;
|
||||
if (Dm.kproc==0)
|
||||
{
|
||||
for ( k=1;k<2;k++){
|
||||
for ( j=1;j<Ny-1;j++){
|
||||
for ( i=1;i<Nx-1;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
//.......................................................................
|
||||
// Compute the media porosity, assign phase labels and solid composition
|
||||
//.......................................................................
|
||||
double sum;
|
||||
double sum_local=0.0, porosity;
|
||||
int Np=0; // number of local pore nodes
|
||||
double *PhaseLabel;
|
||||
PhaseLabel = new double[N];
|
||||
Dm.AssignComponentLabels(PhaseLabel);
|
||||
//.......................................................................
|
||||
for (k=1;k<Nz-1;k++){
|
||||
for (j=1;j<Ny-1;j++){
|
||||
for (i=1;i<Nx-1;i++){
|
||||
n = k*Nx*Ny+j*Nx+i;
|
||||
if (id[n] > 0){
|
||||
inlet_area_local+=1.f;
|
||||
sum_local+=1.0;
|
||||
Np++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
MPI_Allreduce(&inlet_area_local,&InletArea,1,MPI_DOUBLE,MPI_SUM,comm);
|
||||
|
||||
// Set up kstart, kfinish so that the reservoirs are excluded from averaging
|
||||
int kstart,kfinish;
|
||||
kstart = 1;
|
||||
kfinish = Nz-1;
|
||||
if (BoundaryCondition > 0 && Dm.kproc==0) kstart = 4;
|
||||
if (BoundaryCondition > 0 && Dm.kproc==nprocz-1) kfinish = Nz-4;
|
||||
|
||||
// Compute the pore volume
|
||||
sum_local = 0.0;
|
||||
for ( k=kstart;k<kfinish;k++){
|
||||
for ( j=1;j<Ny-1;j++){
|
||||
for ( i=1;i<Nx-1;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
if (id[n] > 0){
|
||||
sum_local += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
MPI_Allreduce(&sum_local,&pore_vol,1,MPI_DOUBLE,MPI_SUM,comm);
|
||||
porosity = pore_vol*iVol_global;
|
||||
MPI_Allreduce(&sum_local,&sum,1,MPI_DOUBLE,MPI_SUM,comm);
|
||||
porosity = sum*iVol_global;
|
||||
if (rank==0) printf("Media porosity = %f \n",porosity);
|
||||
//.........................................................
|
||||
/* // If external boundary conditions are applied remove solid
|
||||
if (BoundaryCondition > 0 && Dm.kproc == 0){
|
||||
// zero out the first layer so that no comms are performed on that boundary
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
int n = j*Nx+i;
|
||||
id[n] = 0;
|
||||
}
|
||||
}
|
||||
for (k=0; k<3; k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
//id[n] = 1;
|
||||
//Averages.SDs(n) = max(Averages.SDs(n),1.0*(2.5-k));
|
||||
Averages->SDs(n) = max(Averages->SDs(n),1.0*(2.5-k));
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if (BoundaryCondition > 0 && Dm.kproc == nprocz-1){
|
||||
// zero out the last layer so that no comms are performed on that boundary
|
||||
for (j=0; j<Ny; j++){
|
||||
for (i=0; i<Nx; i++){
|
||||
int n = (Nz-1)*Nx*Ny+j*Nx+i;
|
||||
id[n] = 0;
|
||||
}
|
||||
}
|
||||
for (k=Nz-3; k<Nz; k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
//id[n] = 2;
|
||||
//Averages.SDs(n) = max(Averages.SDs(n),1.0*(k-Nz+2.5));
|
||||
Averages->SDs(n) = max(Averages->SDs(n),1.0*(k-Nz+2.5));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
*/
|
||||
//.........................................................
|
||||
// don't perform computations at the eight corners
|
||||
id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0;
|
||||
id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0;
|
||||
//.........................................................
|
||||
|
||||
// set reservoirs
|
||||
if (BoundaryCondition > 0){
|
||||
for ( k=0;k<Nz;k++){
|
||||
for ( j=0;j<Ny;j++){
|
||||
for ( i=0;i<Nx;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
if (id[n] > 0){
|
||||
if (Dm.kproc==0 && k==0) id[n]=1;
|
||||
if (Dm.kproc==0 && k==1) id[n]=1;
|
||||
if (Dm.kproc==nprocz-1 && k==Nz-2) id[n]=2;
|
||||
if (Dm.kproc==nprocz-1 && k==Nz-1) id[n]=2;
|
||||
}
|
||||
Mask.id[n] = id[n];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize communication structures in averaging domain
|
||||
for (i=0; i<Mask.Nx*Mask.Ny*Mask.Nz; i++) Mask.id[i] = id[i];
|
||||
Mask.CommInit(comm);
|
||||
|
||||
//...........................................................................
|
||||
if (rank==0) printf ("Create ScaLBL_Communicator \n");
|
||||
// Create a communicator for the device
|
||||
// Create a communicator for the device (will use optimized layout)
|
||||
ScaLBL_Communicator ScaLBL_Comm(Mask);
|
||||
//Create a second communicator based on the regular data layout
|
||||
ScaLBL_Communicator ScaLBL_Comm_Regular(Mask);
|
||||
|
||||
//...........device phase ID.................................................
|
||||
if (rank==0) printf ("Copy phase ID to device \n");
|
||||
char *ID;
|
||||
ScaLBL_AllocateDeviceMemory((void **) &ID, N); // Allocate device memory
|
||||
// Don't compute in the halo
|
||||
for (k=0;k<Nz;k++){
|
||||
for (j=0;j<Ny;j++){
|
||||
for (i=0;i<Nx;i++){
|
||||
int n = k*Nx*Ny+j*Nx+i;
|
||||
if (i==0 || i==Nx-1 || j==0 || j==Ny-1 || k==0 || k==Nz-1) id[n] = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Copy to the device
|
||||
ScaLBL_CopyToDevice(ID, id, N);
|
||||
ScaLBL_DeviceBarrier();
|
||||
//...........................................................................
|
||||
if (rank==0) printf ("Set up memory efficient layout \n");
|
||||
int neighborSize=18*Np*sizeof(int);
|
||||
int *neighborList;
|
||||
IntArray Map(Nx,Ny,Nz);
|
||||
neighborList= new int[18*Np];
|
||||
ScaLBL_Comm.MemoryOptimizedLayoutAA(Map,neighborList,Dm.id,Np);
|
||||
MPI_Barrier(comm);
|
||||
|
||||
//...........................................................................
|
||||
// MAIN VARIABLES ALLOCATED HERE
|
||||
@ -566,50 +387,75 @@ int main(int argc, char **argv)
|
||||
// LBM variables
|
||||
if (rank==0) printf ("Allocating distributions \n");
|
||||
//......................device distributions.................................
|
||||
double *f_even,*f_odd;
|
||||
double *A_even,*A_odd,*B_even,*B_odd;
|
||||
int dist_mem_size = Np*sizeof(double);
|
||||
if (rank==0) printf ("Allocating distributions \n");
|
||||
|
||||
int *NeighborList;
|
||||
int *dvcMap;
|
||||
// double *f_even,*f_odd;
|
||||
double *fq, *Aq, *Bq;
|
||||
double *Den, *Phi;
|
||||
double *ColorGrad;
|
||||
double *Vel;
|
||||
double *Pressure;
|
||||
|
||||
//...........................................................................
|
||||
ScaLBL_AllocateDeviceMemory((void **) &f_even, 10*dist_mem_size); // Allocate device memory
|
||||
ScaLBL_AllocateDeviceMemory((void **) &f_odd, 9*dist_mem_size); // Allocate device memory
|
||||
ScaLBL_AllocateDeviceMemory((void **) &A_even, 4*dist_mem_size); // Allocate device memory
|
||||
ScaLBL_AllocateDeviceMemory((void **) &A_odd, 3*dist_mem_size); // Allocate device memory
|
||||
ScaLBL_AllocateDeviceMemory((void **) &B_even, 4*dist_mem_size); // Allocate device memory
|
||||
ScaLBL_AllocateDeviceMemory((void **) &B_odd, 3*dist_mem_size); // Allocate device memory
|
||||
//...........................................................................
|
||||
double *Phi,*Den;
|
||||
double *ColorGrad, *Velocity, *Pressure, *dvcSignDist;
|
||||
//...........................................................................
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Phi, dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Pressure, dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &dvcSignDist, dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np);
|
||||
|
||||
ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Aq, 7*dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*dist_mem_size);
|
||||
//...........................................................................
|
||||
|
||||
// Copy signed distance for device initialization
|
||||
//ScaLBL_CopyToDevice(dvcSignDist, Averages.SDs.data(), dist_mem_size);
|
||||
ScaLBL_CopyToDevice(dvcSignDist, Averages->SDs.data(), dist_mem_size);
|
||||
//...........................................................................
|
||||
|
||||
int logcount = 0; // number of surface write-outs
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Vel, 3*sizeof(double)*Np);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np);
|
||||
|
||||
//...........................................................................
|
||||
// MAIN VARIABLES INITIALIZED HERE
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
if (rank==0) printf("Setting the distributions, size = %i\n", N);
|
||||
//...........................................................................
|
||||
// Update GPU data structures
|
||||
if (rank==0) printf ("Setting up device map and neighbor list \n");
|
||||
int *TmpMap;
|
||||
TmpMap=new int[Np];
|
||||
for (k=1; k<Nz-1; k++){
|
||||
for (j=1; j<Ny-1; j++){
|
||||
for (i=1; i<Nx-1; i++){
|
||||
int idx=Map(i,j,k);
|
||||
if (!(idx < 0))
|
||||
TmpMap[idx] = k*Nx*Ny+j*Nx+i;
|
||||
}
|
||||
}
|
||||
}
|
||||
//for (int idx=0; idx<Np; idx++) printf("Map=%i\n",TmpMap[idx]);
|
||||
|
||||
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_D3Q19_Init(ID, f_even, f_odd, Nx, Ny, Nz);
|
||||
ScaLBL_Color_Init(ID, Den, Phi, das, dbs, Nx, Ny, Nz);
|
||||
ScaLBL_DeviceBarrier();
|
||||
//......................................................................
|
||||
delete [] TmpMap;
|
||||
|
||||
// copy the neighbor list
|
||||
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
|
||||
// initialize phi based on PhaseLabel (include solid component labels)
|
||||
ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double));
|
||||
//...........................................................................
|
||||
|
||||
if (rank==0) printf ("Initializing distributions \n");
|
||||
// Initialize the phase field and variables
|
||||
ScaLBL_D3Q19_Init(fq, Np);
|
||||
if (rank==0) printf ("Initializing phase field \n");
|
||||
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, Np);
|
||||
if (Dm.kproc==0){
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
|
||||
}
|
||||
if (Dm.kproc == nprocz-1){
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
|
||||
}
|
||||
|
||||
|
||||
if (Restart == true){
|
||||
|
||||
|
||||
if (rank==0){
|
||||
printf("Reading restart file! \n");
|
||||
ifstream restart("Restart.txt");
|
||||
@ -640,11 +486,6 @@ int main(int argc, char **argv)
|
||||
MPI_Barrier(comm);
|
||||
}
|
||||
|
||||
//......................................................................
|
||||
ScaLBL_D3Q7_Init(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
|
||||
ScaLBL_D3Q7_Init(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
//.......................................................................
|
||||
// Once phase has been initialized, map solid to account for 'smeared' interface
|
||||
//for (i=0; i<N; i++) Averages.SDs(i) -= (1.0);
|
||||
@ -652,147 +493,30 @@ int main(int argc, char **argv)
|
||||
for (i=0; i<N; i++) Dm.id[i] = Mask.id[i];
|
||||
//.......................................................................
|
||||
// Finalize setup for averaging domain
|
||||
//Averages.SetupCubes(Mask);
|
||||
// Averages.UpdateSolid();
|
||||
Averages->UpdateSolid();
|
||||
|
||||
//.......................................................................
|
||||
|
||||
//*************************************************************************
|
||||
// Compute the phase indicator field and reset Copy, Den
|
||||
//*************************************************************************
|
||||
ScaLBL_ComputePhaseField(ID, Phi, Den, N);
|
||||
//*************************************************************************
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_Comm.SendHalo(Phi);
|
||||
ScaLBL_Comm.RecvHalo(Phi);
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
//*************************************************************************
|
||||
|
||||
if (rank==0 && BoundaryCondition==1){
|
||||
printf("Setting inlet pressure = %f \n", din);
|
||||
printf("Setting outlet pressure = %f \n", dout);
|
||||
}
|
||||
if (BoundaryCondition==1 && Mask.kproc == 0) {
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
}
|
||||
|
||||
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
}
|
||||
|
||||
if (rank==0 && BoundaryCondition==2){
|
||||
printf("Setting inlet velocity = %f \n", din);
|
||||
printf("Setting outlet velocity = %f \n", dout);
|
||||
}
|
||||
if (BoundaryCondition==2 && Mask.kproc == 0) {
|
||||
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
}
|
||||
|
||||
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
}
|
||||
|
||||
// Set dynamic pressure boundary conditions
|
||||
double dp, slope;
|
||||
dp = slope = 0.0;
|
||||
if (BoundaryCondition==3){
|
||||
slope = abs(dout-din)/timestepMax;
|
||||
dp = abs(timestep)*slope;
|
||||
if (rank==0) printf("Change in pressure / time =%.3e \n",slope);
|
||||
// set the initial value
|
||||
din = 1.0+0.5*dp;
|
||||
dout = 1.0-0.5*dp;
|
||||
// set the initial boundary conditions
|
||||
if (Mask.kproc == 0) {
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
|
||||
}
|
||||
if (Mask.kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
}
|
||||
}
|
||||
|
||||
// Set flux boundary condition
|
||||
double SumLocal,SumGlobal;
|
||||
if (BoundaryCondition==4){
|
||||
double SumLocal=0.f;
|
||||
if (rank==0){
|
||||
printf("Using flux boundary condition \n");
|
||||
printf(" Volumetric flux: Q = %f \n",flux);
|
||||
printf(" Inlet area = %f \n",InletArea);
|
||||
printf(" outlet pressure: %f \n",dout);
|
||||
|
||||
}
|
||||
SumLocal=0.0;
|
||||
if (pBC && Dm.kproc == 0){
|
||||
SumLocal = ScaLBL_D3Q19_Flux_BC_z(ID,f_even,f_odd,flux,InletArea,Nx,Ny,Nz);
|
||||
}
|
||||
MPI_Allreduce(&SumLocal,&SumGlobal,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
|
||||
din = flux / InletArea + SumGlobal;
|
||||
|
||||
if (pBC && Dm.kproc == 0){
|
||||
if (rank==0) printf("Flux = %.3e, Computed inlet pressure: %f \n",flux,din);
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
|
||||
}
|
||||
|
||||
if (pBC && Dm.kproc == nprocz-1){
|
||||
// if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
|
||||
}
|
||||
}
|
||||
|
||||
ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
|
||||
ScaLBL_D3Q19_Pressure(fq,Pressure,Np);
|
||||
ScaLBL_D3Q19_Momentum(fq,Vel,Np);
|
||||
//...........................................................................
|
||||
// Copy the phase indicator field for the earlier timestep
|
||||
ScaLBL_DeviceBarrier();
|
||||
//ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Phase_tplus.data(),Phi,N*sizeof(double));
|
||||
|
||||
//...........................................................................
|
||||
//...........................................................................
|
||||
// Copy the data for for the analysis timestep
|
||||
//...........................................................................
|
||||
// Copy the phase from the GPU -> CPU
|
||||
//...........................................................................
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_D3Q19_Pressure(ID,f_even,f_odd,Pressure,Nx,Ny,Nz);
|
||||
ScaLBL_CopyToHost(Averages->Phase.data(),Phi,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Press.data(),Pressure,N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_x.data(),&Velocity[0],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_y.data(),&Velocity[N],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],N*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Press.data(),Pressure,Np*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_x.data(),&Velocity[0],Np*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_y.data(),&Velocity[N],Np*sizeof(double));
|
||||
ScaLBL_CopyToHost(Averages->Vel_z.data(),&Velocity[2*N],Np*sizeof(double));
|
||||
|
||||
double TmpDat;
|
||||
TmpDat = new double [Np];
|
||||
ScaLBL_CopyToHost(&TmpDat[0],&Pressure[0],SIZE);
|
||||
ScaLBL_Comm.RegularLayout(Map,TmpDat,Averages->Press.data());
|
||||
//...........................................................................
|
||||
|
||||
if (rank==0) printf("********************************************************\n");
|
||||
@ -861,164 +585,84 @@ int main(int argc, char **argv)
|
||||
while (timestep < timestepMax && err > tol ) {
|
||||
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
|
||||
PROFILE_START("Update");
|
||||
// *************ODD TIMESTEP*************
|
||||
timestep++;
|
||||
// Compute the Phase indicator field
|
||||
// Read for Aq, Bq happens in this routine (requires communication)
|
||||
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
|
||||
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm.next, Np, Np);
|
||||
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
|
||||
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
|
||||
|
||||
//*************************************************************************
|
||||
// Fused Color Gradient and Collision
|
||||
//*************************************************************************
|
||||
ScaLBL_D3Q19_ColorCollide_gen( ID,f_even,f_odd,Phi,ColorGrad,
|
||||
Velocity,Nx,Ny,Nz,tau1,tau2,rho1,rho2,alpha,beta,Fx,Fy,Fz);
|
||||
//*************************************************************************
|
||||
// Halo exchange for phase field
|
||||
ScaLBL_Comm_Regular.SendHalo(Phi);
|
||||
ScaLBL_Comm_Regular.RecvHalo(Phi);
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
//*************************************************************************
|
||||
// Pack and send the D3Q19 distributions
|
||||
ScaLBL_Comm.SendD3Q19(f_even, f_odd);
|
||||
//*************************************************************************
|
||||
|
||||
//*************************************************************************
|
||||
// Carry out the density streaming step for mass transport
|
||||
//*************************************************************************
|
||||
ScaLBL_D3Q7_ColorCollideMass(ID, A_even, A_odd, B_even, B_odd, Den, Phi,
|
||||
ColorGrad, Velocity, beta, N, pBC);
|
||||
//*************************************************************************
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
//*************************************************************************
|
||||
// Swap the distributions for momentum transport
|
||||
//*************************************************************************
|
||||
ScaLBL_D3Q19_Swap(ID, f_even, f_odd, Nx, Ny, Nz);
|
||||
//*************************************************************************
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
//*************************************************************************
|
||||
// Wait for communications to complete and unpack the distributions
|
||||
ScaLBL_Comm.RecvD3Q19(f_even, f_odd);
|
||||
//*************************************************************************
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
//*************************************************************************
|
||||
// Pack and send the D3Q7 distributions
|
||||
ScaLBL_Comm.BiSendD3Q7(A_even, A_odd, B_even, B_odd);
|
||||
//*************************************************************************
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_D3Q7_Swap(ID, A_even, A_odd, Nx, Ny, Nz);
|
||||
ScaLBL_D3Q7_Swap(ID, B_even, B_odd, Nx, Ny, Nz);
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
|
||||
//*************************************************************************
|
||||
// Wait for communication and unpack the D3Q7 distributions
|
||||
ScaLBL_Comm.BiRecvD3Q7(A_even, A_odd, B_even, B_odd);
|
||||
//*************************************************************************
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
//..................................................................................
|
||||
ScaLBL_D3Q7_Density(ID, A_even, A_odd, &Den[0], Nx, Ny, Nz);
|
||||
ScaLBL_D3Q7_Density(ID, B_even, B_odd, &Den[N], Nx, Ny, Nz);
|
||||
//*************************************************************************
|
||||
// Compute the phase indicator field
|
||||
//*************************************************************************
|
||||
ScaLBL_DeviceBarrier();
|
||||
MPI_Barrier(comm);
|
||||
|
||||
ScaLBL_ComputePhaseField(ID, Phi, Den, N);
|
||||
//*************************************************************************
|
||||
ScaLBL_Comm.SendHalo(Phi);
|
||||
ScaLBL_DeviceBarrier();
|
||||
ScaLBL_Comm.RecvHalo(Phi);
|
||||
//*************************************************************************
|
||||
|
||||
ScaLBL_DeviceBarrier();
|
||||
|
||||
// Pressure boundary conditions
|
||||
if (BoundaryCondition==1 && Mask.kproc == 0) {
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
// Perform the collision operation
|
||||
ScaLBL_Comm.SendD3Q19AA(fq); //READ FROM NORMAL
|
||||
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
|
||||
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm.next, Np, Np);
|
||||
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
|
||||
// Set BCs
|
||||
if (typeBC > 0){
|
||||
ScaLBL_Comm.Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
|
||||
ScaLBL_Comm.Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
|
||||
}
|
||||
if (BoundaryCondition==1 && Mask.kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
if (typeBC == 3){
|
||||
ScaLBL_Comm.D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
|
||||
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
|
||||
}
|
||||
if (typeBC == 4){
|
||||
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
|
||||
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
|
||||
}
|
||||
ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
|
||||
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm.next, Np);
|
||||
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
|
||||
|
||||
// Velocity boundary conditions
|
||||
if (BoundaryCondition==2 && Mask.kproc == 0) {
|
||||
ScaLBL_D3Q19_Velocity_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
}
|
||||
if (BoundaryCondition==2 && Mask.kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Velocity_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
}
|
||||
// *************EVEN TIMESTEP*************
|
||||
timestep++;
|
||||
// Compute the Phase indicator field
|
||||
ScaLBL_Comm.BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
|
||||
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm.next, Np, Np);
|
||||
ScaLBL_Comm.BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
|
||||
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm.next, Np);
|
||||
|
||||
if (BoundaryCondition==3){
|
||||
// Increase the pressure difference
|
||||
dp += slope;
|
||||
din = 1.0+0.5*dp;
|
||||
dout = 1.0-0.5*dp;
|
||||
// set the initial boundary conditions
|
||||
if (Mask.kproc == 0) {
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
}
|
||||
if (Mask.kproc == nprocz-1){
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
}
|
||||
}
|
||||
// Halo exchange for phase field
|
||||
ScaLBL_Comm_Regular.SendHalo(Phi);
|
||||
ScaLBL_Comm_Regular.RecvHalo(Phi);
|
||||
|
||||
// Set flux boundary condition
|
||||
if (BoundaryCondition==4){
|
||||
SumLocal=0.0;
|
||||
if (pBC && Dm.kproc == 0){
|
||||
SumLocal = ScaLBL_D3Q19_Flux_BC_z(ID,f_even,f_odd,flux,InletArea,Nx,Ny,Nz);
|
||||
// Perform the collision operation
|
||||
ScaLBL_Comm.SendD3Q19AA(fq); //READ FORM NORMAL
|
||||
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
|
||||
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm.next, Np, Np);
|
||||
ScaLBL_Comm.RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
|
||||
// Set boundary conditions
|
||||
if (typeBC > 0){
|
||||
ScaLBL_Comm.Color_BC_z(dvcMap, Phi, Den, inletA, inletB);
|
||||
ScaLBL_Comm.Color_BC_Z(dvcMap, Phi, Den, outletA, outletB);
|
||||
}
|
||||
MPI_Allreduce(&SumLocal,&SumGlobal,1,MPI_DOUBLE,MPI_SUM,Dm.Comm);
|
||||
din = flux / InletArea + SumGlobal;
|
||||
|
||||
if (pBC && Dm.kproc == 0){
|
||||
ScaLBL_D3Q19_Pressure_BC_z(f_even,f_odd,din,Nx,Ny,Nz);
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
|
||||
}
|
||||
|
||||
if (pBC && Dm.kproc == nprocz-1){
|
||||
// if (rank==nprocx*nprocy*nprocz-1) printf("Flux = %.3e, Computed outlet pressure: %f \n",flux,dout);
|
||||
ScaLBL_D3Q19_Pressure_BC_Z(f_even,f_odd,dout,Nx,Ny,Nz,Nx*Ny*(Nz-2));
|
||||
ScaLBL_D3Q19_Velocity(ID,f_even,f_odd,Velocity,Nx,Ny,Nz);
|
||||
ScaLBL_Color_BC_Z(Phi,Den,Velocity,A_even,A_odd,B_even,B_odd,Nx,Ny,Nz);
|
||||
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
|
||||
if (typeBC == 3){
|
||||
ScaLBL_Comm.D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
|
||||
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
|
||||
}
|
||||
else if (typeBC == 4){
|
||||
din = ScaLBL_Comm.D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep);
|
||||
ScaLBL_Comm.D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep);
|
||||
}
|
||||
ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Vel, rhoA, rhoB, tauA, tauB,
|
||||
alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm.next, Np);
|
||||
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
|
||||
//************************************************************************
|
||||
|
||||
MPI_Barrier(comm);
|
||||
PROFILE_STOP("Update");
|
||||
|
||||
// Timestep completed!
|
||||
timestep++;
|
||||
|
||||
// Run the analysis
|
||||
run_analysis(timestep,RESTART_INTERVAL,rank_info,*Averages,last_ids,last_index,last_id_map,
|
||||
Nx,Ny,Nz,pBC,beta,err,Phi,Pressure,Velocity,ID,f_even,f_odd,Den,
|
||||
LocalRestartFile,meshData,fillData,tpool,work_ids);
|
||||
|
||||
|
||||
// Save the timers
|
||||
if ( timestep%50==0 )
|
||||
PROFILE_SAVE("lbpm_color_simulator",1);
|
||||
|
||||
Loading…
Reference in New Issue
Block a user