OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

werkin

Browse Source
This commit is contained in:
James McClure
2021-01-05 00:15:36 -05:00
parent 7e51ba962b
commit 77c2949b07
5 changed files with 81 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
models/GreyscaleColorModel.cpp
View File

@@ -273,7 +273,7 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx] = sumReduce( Dm->Comm, label_count[idx]);
label_count_global[idx] = Dm->Comm.sumReduce( label_count[idx]);
if (rank==0){
printf("Number of component labels: %lu \n",NLABELS);
@@ -286,8 +286,7 @@ void ScaLBL_GreyscaleColorModel::AssignComponentLabels()
}
ScaLBL_CopyToDevice(Phi, phase, N*sizeof(double));
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
ScaLBL_Comm->Barrier();
delete [] phase;
}
@@ -447,7 +446,7 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//Model-4
ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] SolidPotential_host;
delete [] GreySolidGrad_host;
delete [] Dst;
@@ -545,7 +544,7 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
// Set Dm to match Mask
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
//Initialize a weighted porosity after considering grey voxels
GreyPorosity=0.0;
@@ -571,7 +570,7 @@ void ScaLBL_GreyscaleColorModel::AssignGreyPoroPermLabels()
ScaLBL_CopyToDevice(Porosity_dvc, Porosity, Np*sizeof(double));
ScaLBL_CopyToDevice(Permeability_dvc, Permeability, Np*sizeof(double));
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] Porosity;
delete [] Permeability;
}
@@ -658,7 +657,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] TmpMap;
// copy the neighbor list
@@ -739,7 +738,7 @@ void ScaLBL_GreyscaleColorModel::Initialize(){
ScaLBL_CopyToDevice(Den,cDen,2*Np*sizeof(double));
ScaLBL_CopyToDevice(fq,cDist,19*Np*sizeof(double));
ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double));
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
MPI_Barrier(comm);
@@ -912,7 +911,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
MPI_Barrier(comm);
starttime = MPI_Wtime();
//.........................................
@@ -934,7 +933,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
@@ -955,7 +954,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
// Set BCs
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
@@ -978,16 +977,15 @@ void ScaLBL_GreyscaleColorModel::Run(){
//ScaLBL_D3Q19_AAodd_GreyscaleColor(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
ScaLBL_Comm->Barrier();
// *************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->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np);
// Perform the collision operation
@@ -1008,7 +1006,7 @@ void ScaLBL_GreyscaleColorModel::Run(){
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm_Regular->RecvHalo(Phi);
ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
// Set boundary conditions
if (BoundaryCondition == 3){
ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep);
@@ -1031,9 +1029,8 @@ void ScaLBL_GreyscaleColorModel::Run(){
//ScaLBL_D3Q19_AAeven_GreyscaleColor(dvcMap, fq, Aq, Bq, Den, Phi,GreySolidPhi,Porosity_dvc,Permeability_dvc,Velocity,
// rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff,
// alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
//************************************************************************
ScaLBL_Comm->Barrier();
//************************************************************************
PROFILE_STOP("Update");
//TODO For temporary use - writing Restart and Vis files should be included in the analysis framework in the future
@@ -1315,15 +1312,14 @@ void ScaLBL_GreyscaleColorModel::Run(){
}
morph_timesteps += analysis_interval;
}
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
ScaLBL_Comm->Barrier();
}
//analysis.finish();
PROFILE_STOP("Loop");
PROFILE_SAVE("lbpm_color_simulator",1);
//************************************************************************
ScaLBL_DeviceBarrier();
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
stoptime = MPI_Wtime();
ScaLBL_Comm->Barrier();
stoptime = MPI_Wtime();
if (rank==0) printf("-------------------------------------------------------------------\n");
// Compute the walltime per timestep
cputime = (stoptime - starttime)/timestep;
@@ -1369,14 +1365,14 @@ void ScaLBL_GreyscaleColorModel::ImageInit(std::string Filename){
// }
// }
// }
// Count=sumReduce( Dm->Comm, Count);
// PoreCount=sumReduce( Dm->Comm, PoreCount);
// Count=Dm->Comm.sumReduce( Count);
// PoreCount=Dm->Comm.sumReduce( PoreCount);
// if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount);
ScaLBL_D3Q19_Init(fq, Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
MPI_Barrier(ScaLBL_Comm->MPI_COMM_SCALBL);
ScaLBL_Comm->Barrier();
//ScaLBL_CopyToHost(Averages->Phi.data(),Phi,Nx*Ny*Nz*sizeof(double));
@@ -1447,8 +1443,8 @@ double ScaLBL_GreyscaleColorModel::SeedPhaseField(const double seed_water_in_oil
mass_loss += random_value*seed_water_in_oil;
}
count= sumReduce( Dm->Comm, count);
mass_loss= sumReduce( Dm->Comm, mass_loss);
count= Dm->Comm.sumReduce( count);
mass_loss= Dm->Comm.sumReduce( mass_loss);
if (rank == 0) printf("Remove mass %.5g from %.5g voxels \n",mass_loss,count);
// Need to initialize Aq, Bq, Den, Phi directly
@@ -1828,7 +1824,7 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
//
//
// ScaLBL_CopyToDevice(GreySolidGrad, GreySolidGrad_host, 3*Np*sizeof(double));
// ScaLBL_DeviceBarrier();
// ScaLBL_Comm->Barrier();
// delete [] SolidPotential_host;
// delete [] GreySolidGrad_host;
// delete [] Dst;
@@ -1976,7 +1972,7 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_DeviceBarrier();
// ScaLBL_Comm->Barrier();
//
// //debug
// //FILE *OUTFILE;
@@ -2047,7 +2043,7 @@ void ScaLBL_GreyscaleColorModel::WriteDebug(){
// }
//
// ScaLBL_CopyToDevice(GreySolidPhi, GreySolidPhi_host, Nx*Ny*Nz*sizeof(double));
// ScaLBL_DeviceBarrier();
// ScaLBL_Comm->Barrier();
//
// //debug
// FILE *OUTFILE;

12
models/GreyscaleModel.cpp
View File

@@ -267,7 +267,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity, double *Perm
// Set Dm to match Mask
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = Mask->id[i];
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
//Initialize a weighted porosity after considering grey voxels
GreyPorosity=0.0;
for (unsigned int idx=0; idx<NLABELS; idx++){
@@ -331,7 +331,7 @@ void ScaLBL_GreyscaleModel::AssignComponentLabels(double *Porosity,double *Perme
}
}
}
GreyPorosity = sumReduce( Dm->Comm, GreyPorosity_loc);
GreyPorosity = Dm->Comm.sumReduce( GreyPorosity_loc);
GreyPorosity = GreyPorosity/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (rank==0){
@@ -683,10 +683,10 @@ void ScaLBL_GreyscaleModel::Run(){
double As = Morphology.A();
double Hs = Morphology.H();
double Xs = Morphology.X();
Vs = sumReduce( Dm->Comm, Vs);
As = sumReduce( Dm->Comm, As);
Hs = sumReduce( Dm->Comm, Hs);
Xs = sumReduce( Dm->Comm, Xs);
Vs = Dm->Comm.sumReduce( Vs);
As = Dm->Comm.sumReduce( As);
Hs = Dm->Comm.sumReduce( Hs);
Xs = Dm->Comm.sumReduce( Xs);
double h = Dm->voxel_length;
//double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;

24
models/IonModel.cpp
View File

@@ -588,7 +588,7 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
}
for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
if (rank==0){
printf("LB Ion Solver: number of ion solid labels: %lu \n",NLABELS);
@@ -684,7 +684,7 @@ void ScaLBL_IonModel::Create(){
IonSolid_host = new double[Nx*Ny*Nz];
AssignSolidBoundary(IonSolid_host);
ScaLBL_CopyToDevice(IonSolid, IonSolid_host, Nx*Ny*Nz*sizeof(double));
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] IonSolid_host;
}
}
@@ -782,7 +782,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//.......create and start timer............
//double starttime,stoptime,cputime;
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//starttime = MPI_Wtime();
for (int ic=0; ic<number_ion_species; ic++){
@@ -795,7 +795,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_IonConcentration(NeighborList, &fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
//--------------------------------------- Set boundary conditions -------------------------------------//
if (BoundaryConditionInlet[ic]>0){
switch (BoundaryConditionInlet[ic]){
@@ -830,7 +830,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
if (BoundaryConditionSolid==1){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
}
// *************EVEN TIMESTEP*************//
@@ -839,7 +839,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
ScaLBL_Comm->SendD3Q7AA(fq, ic); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_IonConcentration(&fq[ic*Np*7],&Ci[ic*Np],ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, ic); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
//--------------------------------------- Set boundary conditions -------------------------------------//
if (BoundaryConditionInlet[ic]>0){
switch (BoundaryConditionInlet[ic]){
@@ -874,7 +874,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
if (BoundaryConditionSolid==1){
//TODO IonSolid may also be species-dependent
ScaLBL_Comm->SolidDirichletD3Q7(&fq[ic*Np*7], IonSolid);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
}
}
}
@@ -904,7 +904,7 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const i
//This function wirte out the data in a normal layout (by aggregating all decomposed domains)
ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
IonConcentration_LB_to_Phys(IonConcentration);
}
@@ -914,7 +914,7 @@ void ScaLBL_IonModel::getIonConcentration_debug(int timestep){
DoubleArray PhaseField(Nx,Ny,Nz);
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
IonConcentration_LB_to_Phys(PhaseField);
FILE *OUTFILE;
@@ -955,7 +955,7 @@ void ScaLBL_IonModel::DummyFluidVelocity(){
}
ScaLBL_AllocateDeviceMemory((void **) &FluidVelocityDummy, sizeof(double)*3*Np);
ScaLBL_CopyToDevice(FluidVelocityDummy, FluidVelocity_host, sizeof(double)*3*Np);
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] FluidVelocity_host;
}
@@ -976,7 +976,7 @@ void ScaLBL_IonModel::DummyElectricField(){
}
ScaLBL_AllocateDeviceMemory((void **) &ElectricFieldDummy, sizeof(double)*3*Np);
ScaLBL_CopyToDevice(ElectricFieldDummy, ElectricField_host, sizeof(double)*3*Np);
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] ElectricField_host;
}
@@ -1028,7 +1028,7 @@ double ScaLBL_IonModel::CalIonDenConvergence(vector<double> &ci_avg_previous){
// DoubleArray PhaseField(Nx,Ny,Nz);
// for (int ic=0; ic<number_ion_species; ic++){
// ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],PhaseField);
// ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"Ion%02i.%05i.raw",ic+1,rank);

38
models/PoissonSolver.cpp
View File

@@ -231,7 +231,7 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid)
}
for (size_t idx=0; idx<NLABELS; idx++)
label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
if (rank==0){
printf("LB-Poisson Solver: number of Poisson solid labels: %lu \n",NLABELS);
@@ -324,16 +324,16 @@ void ScaLBL_Poisson::Create(){
}
}
ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np);
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] TmpMap;
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
MPI_Barrier(comm);
delete [] neighborList;
// copy node ID
//ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz);
//ScaLBL_DeviceBarrier();
//ScaLBL_Comm->Barrier();
//Initialize solid boundary for electric potential
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np);
@@ -388,7 +388,7 @@ void ScaLBL_Poisson::Initialize(){
AssignSolidBoundary(psi_host);//step1
Potential_Init(psi_host);//step2
ScaLBL_CopyToDevice(Psi, psi_host, Nx*Ny*Nz*sizeof(double));
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np);
delete [] psi_host;
@@ -408,7 +408,7 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
//.......create and start timer............
//double starttime,stoptime,cputime;
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//starttime = MPI_Wtime();
timestep=0;
@@ -421,13 +421,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
SolveElectricPotentialAAodd();//update electric potential
SolvePoissonAAodd(ChargeDensity);//perform collision
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
// *************EVEN TIMESTEP*************//
timestep++;
SolveElectricPotentialAAeven();//update electric potential
SolvePoissonAAeven(ChargeDensity);//perform collision
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//************************************************************************/
// Check convergence of steady-state solution
@@ -509,7 +509,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAodd(){
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FROM NORMAL
ScaLBL_D3Q7_AAodd_Poisson_ElectricPotential(NeighborList, dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
// Set boundary conditions
if (BoundaryCondition == 1){
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
@@ -523,7 +523,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(){
ScaLBL_Comm->SendD3Q7AA(fq, 0); //READ FORM NORMAL
ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
ScaLBL_Comm->RecvD3Q7AA(fq, 0); //WRITE INTO OPPOSITE
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
// Set boundary conditions
if (BoundaryCondition == 1){
ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep);
@@ -570,7 +570,7 @@ void ScaLBL_Poisson::DummyChargeDensity(){
}
ScaLBL_AllocateDeviceMemory((void **) &ChargeDensityDummy, sizeof(double)*Np);
ScaLBL_CopyToDevice(ChargeDensityDummy, ChargeDensity_host, sizeof(double)*Np);
ScaLBL_DeviceBarrier();
ScaLBL_Comm->Barrier();
delete [] ChargeDensity_host;
}
@@ -579,7 +579,7 @@ void ScaLBL_Poisson::getElectricPotential_debug(int timestep){
DoubleArray PhaseField(Nx,Ny,Nz);
//ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
ScaLBL_CopyToHost(PhaseField.data(),Psi,sizeof(double)*Nx*Ny*Nz);
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
FILE *OUTFILE;
sprintf(LocalRankFilename,"Electric_Potential_Time_%i.%05i.raw",timestep,rank);
OUTFILE = fopen(LocalRankFilename,"wb");
@@ -597,22 +597,22 @@ void ScaLBL_Poisson::getElectricField(DoubleArray &Values_x, DoubleArray &Values
ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],Values_x);
ElectricField_LB_to_Phys(Values_x);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],Values_y);
ElectricField_LB_to_Phys(Values_y);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],Values_z);
ElectricField_LB_to_Phys(Values_z);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
}
void ScaLBL_Poisson::getElectricField_debug(int timestep){
//ScaLBL_D3Q7_Poisson_getElectricField(fq,ElectricField,tau,Np);
//ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],PhaseField);
@@ -658,7 +658,7 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
// ScaLBL_D3Q7_Poisson_ElectricField(NeighborList, dvcMap, dvcID, Psi, ElectricField, BoundaryConditionSolid,
// Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
// ScaLBL_Comm_Regular->RecvHalo(Psi);
// ScaLBL_DeviceBarrier();
// ScaLBL_Comm->Barrier();
// if (BoundaryCondition == 1){
// ScaLBL_Comm->Poisson_D3Q7_BC_z(dvcMap,Psi,Vin);
// ScaLBL_Comm->Poisson_D3Q7_BC_Z(dvcMap,Psi,Vout);
@@ -671,7 +671,7 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
//
// DoubleArray PhaseField(Nx,Ny,Nz);
// ScaLBL_Comm->RegularLayout(Map,Psi,PhaseField);
// //ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
// //ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"Electric_Potential.%05i.raw",rank);
// OUTFILE = fopen(LocalRankFilename,"wb");
@@ -728,7 +728,7 @@ void ScaLBL_Poisson::ElectricField_LB_to_Phys(DoubleArray &Efield_reg){
// }
//
// for (size_t idx=0; idx<NLABELS; idx++)
// label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
// label_count_global[idx]=Dm->Comm.sumReduce( label_count[idx]);
//
// if (rank==0){
// printf("LB-Poisson Solver: number of Poisson solid labels: %lu \n",NLABELS);

34
models/StokesModel.cpp
View File

@@ -343,7 +343,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
}
ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
@@ -365,7 +365,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
}
ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx, Fy, Fz,rho0,den_scale,h,time_conv,
0, ScaLBL_Comm->LastExterior(), Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//************************************************************************/
}
}
@@ -373,25 +373,25 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity, double *ElectricField){
void ScaLBL_StokesModel::getVelocity(DoubleArray &Vel_x, DoubleArray &Vel_y, DoubleArray &Vel_z){
//get velocity in physical unit [m/sec]
ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Vel_x);
Velocity_LB_to_Phys(Vel_x);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Vel_y);
Velocity_LB_to_Phys(Vel_y);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Vel_z);
Velocity_LB_to_Phys(Vel_z);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
}
void ScaLBL_StokesModel::getVelocity_debug(int timestep){
//get velocity in physical unit [m/sec]
ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
@@ -492,7 +492,7 @@ double ScaLBL_StokesModel::CalVelocityConvergence(double& flow_rate_previous,dou
//-----------------------------------------------------
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
@@ -574,7 +574,7 @@ void ScaLBL_StokesModel::Run(){
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
starttime = MPI_Wtime();
if (rank==0) printf("****************************************************************\n");
if (rank==0) printf("LB Single-Fluid Navier-Stokes Solver: timestepMax = %i\n", timestepMax);
@@ -602,7 +602,7 @@ void ScaLBL_StokesModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAodd_MRT(NeighborList, fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
timestep++;
ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
ScaLBL_D3Q19_AAeven_MRT(fq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
@@ -621,12 +621,12 @@ void ScaLBL_StokesModel::Run(){
ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq);
}
ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//************************************************************************/
if (timestep%1000==0){
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Velocity_x);
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Velocity_y);
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Velocity_z);
@@ -681,10 +681,10 @@ void ScaLBL_StokesModel::Run(){
double As = Morphology.A();
double Hs = Morphology.H();
double Xs = Morphology.X();
Vs=sumReduce( Dm->Comm, Vs);
As=sumReduce( Dm->Comm, As);
Hs=sumReduce( Dm->Comm, Hs);
Xs=sumReduce( Dm->Comm, Xs);
Vs=Dm->Comm.sumReduce( Vs);
As=Dm->Comm.sumReduce( As);
Hs=Dm->Comm.sumReduce( Hs);
Xs=Dm->Comm.sumReduce( Xs);
double h = Dm->voxel_length;
double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag;
if (rank==0) {
@@ -718,7 +718,7 @@ void ScaLBL_StokesModel::VelocityField(){
/* Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));