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

refactor new models to FOM

Browse Source
This commit is contained in:
James McClure 2021-01-05 13:51:32 -05:00
parent 7644a8288e
commit 5bcc10c5fb
10 changed files with 105 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
models/GreyscaleColorModel.cpp
View File

@ -167,7 +167,7 @@ void ScaLBL_GreyscaleColorModel::ReadInput(){
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
id_view.viewRaw( size1, Mask->id.data());
fill.copy( input_id, id_view );
fill.fill( id_view );
}
@ -600,7 +600,7 @@ void ScaLBL_GreyscaleColorModel::Create(){
if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N);
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np);
MPI_Barrier(comm);
//...........................................................................

10
models/GreyscaleModel.cpp
View File

@ -366,7 +366,7 @@ void ScaLBL_GreyscaleModel::Create(){
if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N);
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np);
MPI_Barrier(comm);
//...........................................................................
@ -649,10 +649,10 @@ void ScaLBL_GreyscaleModel::Run(){
}
}
}
vax = sumReduce( Mask->Comm, vax_loc);
vay = sumReduce( Mask->Comm, vay_loc);
vaz = sumReduce( Mask->Comm, vaz_loc);
count = sumReduce( Mask->Comm, count_loc);
vax = Dm->Comm.sumReduce( vax_loc);
vay = Dm->Comm.sumReduce( vay_loc);
vaz = Dm->Comm.sumReduce( vaz_loc);
count = Dm->Comm.sumReduce( count_loc);
vax /= count;
vay /= count;

38
models/IonModel.cpp
View File

@ -461,7 +461,7 @@ void ScaLBL_IonModel::SetDomain(){
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
MPI_Barrier(comm);
comm.barrier();
unsigned short int BC_inlet_min = *min_element(BoundaryConditionInlet.begin(),BoundaryConditionInlet.end());
unsigned short int BC_outlet_min = *min_element(BoundaryConditionOutlet.begin(),BoundaryConditionOutlet.end());
@ -478,7 +478,7 @@ void ScaLBL_IonModel::SetDomain(){
}
Dm->CommInit();
MPI_Barrier(comm);
comm.barrier();
rank = Dm->rank();
nprocx = Dm->nprocx();
@ -506,7 +506,7 @@ void ScaLBL_IonModel::ReadInput(){
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
id_view.viewRaw( size1, Mask->id.data() );
fill.copy( input_id, id_view );
fill.fill( id_view );
}
@ -650,8 +650,8 @@ void ScaLBL_IonModel::Create(){
if (rank==0) printf ("LB Ion Solver: Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
MPI_Barrier(comm);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
@ -670,15 +670,15 @@ void ScaLBL_IonModel::Create(){
if (rank==0) printf ("LB Ion Solver: Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
comm.barrier();
//Initialize solid boundary for electrical potential
//if ion concentration at solid surface is specified
if (BoundaryConditionSolid==1){
ScaLBL_AllocateDeviceMemory((void **) &IonSolid, sizeof(double)*Nx*Ny*Nz);
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np);
MPI_Barrier(comm);
ScaLBL_AllocateDeviceMemory((void **) &IonSolid, sizeof(double)*Nx*Ny*Nz);
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id.data(), Np);
comm.barrier();
double *IonSolid_host;
IonSolid_host = new double[Nx*Ny*Nz];
@ -704,7 +704,7 @@ void ScaLBL_IonModel::Initialize(){
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
}
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
MPI_Barrier(comm);
comm.barrier();
for (int ic=0; ic<number_ion_species; ic++){
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
}
@ -782,7 +782,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){
//.......create and start timer............
//double starttime,stoptime,cputime;
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); comm.barrier();
//starttime = MPI_Wtime();
for (int ic=0; ic<number_ion_species; 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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
}
// *************EVEN TIMESTEP*************//
@ -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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
}
}
}
@ -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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
IonConcentration_LB_to_Phys(PhaseField);
FILE *OUTFILE;
@ -1001,10 +1001,8 @@ double ScaLBL_IonModel::CalIonDenConvergence(vector<double> &ci_avg_previous){
ci_loc +=Ci_host[idx];
count_loc+=1.0;
}
MPI_Allreduce(&ci_loc,&ci_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
ci_avg = Mask->Comm.sumReduce( ci_loc);
count = Mask->Comm.sumReduce( count_loc);
ci_avg /= count;
double ci_avg_mag=ci_avg;
if (ci_avg==0.0) ci_avg_mag=1.0;
@ -1028,7 +1026,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_Comm->Barrier(); MPI_Barrier(comm);
// ScaLBL_Comm->Barrier(); comm.barrier();
//
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"Ion%02i.%05i.raw",ic+1,rank);

4
models/IonModel.h
View File

@ -16,7 +16,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
@ -85,7 +85,7 @@ public:
double *ElectricFieldDummy;
private:
MPI_Comm comm;
Utilities::MPI comm;
// filenames
char LocalRankString[8];

37
models/MRTModel.cpp
View File

@ -86,9 +86,9 @@ void ScaLBL_MRTModel::SetDomain(){
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
MPI_Barrier(comm);
comm.barrier();
Dm->CommInit();
MPI_Barrier(comm);
comm.barrier();
rank = Dm->rank();
nprocx = Dm->nprocx();
@ -116,7 +116,7 @@ void ScaLBL_MRTModel::ReadInput(){
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
id_view.viewRaw( size1, Mask->id.data() );
fill.copy( input_id, id_view );
fill.fill( id_view );
}
@ -193,7 +193,7 @@ void ScaLBL_MRTModel::Create(){
if (rank==0) printf ("Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
comm.barrier();
}
@ -228,7 +228,7 @@ void ScaLBL_MRTModel::Run(){
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_DeviceBarrier(); comm.barrier();
starttime = MPI_Wtime();
if (rank==0) printf("Beginning AA timesteps, timestepMax = %i \n", timestepMax);
if (rank==0) printf("********************************************************\n");
@ -255,7 +255,7 @@ void ScaLBL_MRTModel::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_DeviceBarrier(); comm.barrier();
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);
@ -274,12 +274,12 @@ void ScaLBL_MRTModel::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_DeviceBarrier(); comm.barrier();
//************************************************************************/
if (timestep%1000==0){
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_DeviceBarrier(); comm.barrier();
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);
@ -300,11 +300,11 @@ void ScaLBL_MRTModel::Run(){
}
}
}
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
}
vax=Dm->Comm.sumReduce( vax_loc);
vay=Dm->Comm.sumReduce( vay_loc);
vaz=Dm->Comm.sumReduce( vaz_loc);
count=Dm->Comm.sumReduce( count_loc);
vax /= count;
vay /= count;
@ -334,10 +334,11 @@ void ScaLBL_MRTModel::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) {
@ -371,7 +372,7 @@ void ScaLBL_MRTModel::VelocityField(){
/* Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_DeviceBarrier(); comm.barrier();
ScaLBL_CopyToHost(&VELOCITY[0],&Velocity[0],3*SIZE);
memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));

4
models/MRTModel.h
View File

@ -11,7 +11,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
@ -63,7 +63,7 @@ public:
DoubleArray Velocity_y;
DoubleArray Velocity_z;
private:
MPI_Comm comm;
Utilities::MPI comm;
// filenames
char LocalRankString[8];

40
models/PoissonSolver.cpp
View File

@ -116,11 +116,11 @@ void ScaLBL_Poisson::SetDomain(){
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
MPI_Barrier(comm);
comm.barrier();
Dm->BoundaryCondition = BoundaryCondition;
Mask->BoundaryCondition = BoundaryCondition;
Dm->CommInit();
MPI_Barrier(comm);
comm.barrier();
rank = Dm->rank();
nprocx = Dm->nprocx();
@ -148,7 +148,7 @@ void ScaLBL_Poisson::ReadInput(){
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
id_view.viewRaw( size1, Mask->id.data() );
fill.copy( input_id, id_view );
fill.fill( id_view );
}
@ -275,8 +275,8 @@ void ScaLBL_Poisson::Create(){
if (rank==0) printf ("LB-Poisson Solver: Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
MPI_Barrier(comm);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
@ -329,15 +329,15 @@ void ScaLBL_Poisson::Create(){
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
ScaLBL_Comm->Barrier();
MPI_Barrier(comm);
comm.barrier();
delete [] neighborList;
// copy node ID
//ScaLBL_CopyToDevice(dvcID, Mask->id, sizeof(signed char)*Nx*Ny*Nz);
//ScaLBL_Comm->Barrier();
//Initialize solid boundary for electric potential
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id, Np);
MPI_Barrier(comm);
ScaLBL_Comm->SetupBounceBackList(Map, Mask->id.data(), Np);
comm.barrier();
}
void ScaLBL_Poisson::Potential_Init(double *psi_init){
@ -408,7 +408,7 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
//.......create and start timer............
//double starttime,stoptime,cputime;
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); comm.barrier();
//starttime = MPI_Wtime();
timestep=0;
@ -421,13 +421,13 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
SolveElectricPotentialAAodd();//update electric potential
SolvePoissonAAodd(ChargeDensity);//perform collision
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
// *************EVEN TIMESTEP*************//
timestep++;
SolveElectricPotentialAAeven();//update electric potential
SolvePoissonAAeven(ChargeDensity);//perform collision
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
// Check convergence of steady-state solution
@ -450,9 +450,9 @@ void ScaLBL_Poisson::Run(double *ChargeDensity){
}
}
}
MPI_Allreduce(&psi_loc,&psi_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
psi_avg=Dm->Comm.sumReduce( psi_loc);
count=Dm->Comm.sumReduce( count_loc);
psi_avg /= count;
double psi_avg_mag=psi_avg;
if (psi_avg==0.0) psi_avg_mag=1.0;
@ -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_Comm->Barrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); comm.barrier();
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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&ElectricField[1*Np],Values_y);
ElectricField_LB_to_Phys(Values_y);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&ElectricField[2*Np],Values_z);
ElectricField_LB_to_Phys(Values_z);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_Poisson::getElectricField_debug(int timestep){
//ScaLBL_D3Q7_Poisson_getElectricField(fq,ElectricField,tau,Np);
//ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
//ScaLBL_Comm->Barrier(); comm.barrier();
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,&ElectricField[0*Np],PhaseField);
@ -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_Comm->Barrier(); MPI_Barrier(comm);
// //ScaLBL_Comm->Barrier(); comm.barrier();
// FILE *OUTFILE;
// sprintf(LocalRankFilename,"Electric_Potential.%05i.raw",rank);
// OUTFILE = fopen(LocalRankFilename,"wb");

4
models/PoissonSolver.h
View File

@ -12,7 +12,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
@ -77,7 +77,7 @@ public:
double *ChargeDensityDummy;// for debugging
private:
MPI_Comm comm;
Utilities::MPI comm;
// filenames
char LocalRankString[8];

109
models/StokesModel.cpp
View File

@ -189,11 +189,11 @@ void ScaLBL_StokesModel::SetDomain(){
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
MPI_Barrier(comm);
comm.barrier();
Dm->BoundaryCondition = BoundaryCondition;
Mask->BoundaryCondition = BoundaryCondition;
Dm->CommInit();
MPI_Barrier(comm);
comm.barrier();
rank = Dm->rank();
nprocx = Dm->nprocx();
@ -221,7 +221,7 @@ void ScaLBL_StokesModel::ReadInput(){
ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 );
fillHalo<signed char> fill( comm, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 );
Array<signed char> id_view;
id_view.viewRaw( size1, Mask->id );
id_view.viewRaw( size1, Mask->id.data() );
fill.copy( input_id, id_view );
fill.fill( id_view );
}
@ -278,8 +278,8 @@ void ScaLBL_StokesModel::Create(){
if (rank==0) printf ("LB Single-Fluid Solver: Set up memory efficient layout \n");
Map.resize(Nx,Ny,Nz); Map.fill(-2);
auto neighborList= new int[18*Npad];
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id,Np);
MPI_Barrier(comm);
Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np);
comm.barrier();
//...........................................................................
// MAIN VARIABLES ALLOCATED HERE
//...........................................................................
@ -298,7 +298,7 @@ void ScaLBL_StokesModel::Create(){
if (rank==0) printf ("LB Single-Fluid Solver: Setting up device map and neighbor list \n");
// copy the neighbor list
ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize);
MPI_Barrier(comm);
comm.barrier();
}
@ -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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
}
}
@ -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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],Vel_x);
Velocity_LB_to_Phys(Vel_x);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],Vel_y);
Velocity_LB_to_Phys(Vel_y);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],Vel_z);
Velocity_LB_to_Phys(Vel_z);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
}
void ScaLBL_StokesModel::getVelocity_debug(int timestep){
//get velocity in physical unit [m/sec]
ScaLBL_D3Q19_Momentum(fq, Velocity, Np);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
DoubleArray PhaseField(Nx,Ny,Nz);
ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField);
@ -469,10 +469,10 @@ vector<double> ScaLBL_StokesModel::computeElectricForceAvg(double *ChargeDensity
count_loc+=1.0;
}
MPI_Allreduce(&Fx_loc,&Fx_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&Fy_loc,&Fy_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&Fz_loc,&Fz_avg,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
Fx_avg=Dm->Comm.sumReduce( Fx_loc);
Fy_avg=Dm->Comm.sumReduce( Fy_loc);
Fz_avg=Dm->Comm.sumReduce( Fz_loc);
count=Dm->Comm.sumReduce( count_loc);
Fx_avg /= count;
Fy_avg /= count;
@ -492,7 +492,7 @@ double ScaLBL_StokesModel::CalVelocityConvergence(double& flow_rate_previous,dou
//-----------------------------------------------------
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
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);
@ -514,11 +514,11 @@ double ScaLBL_StokesModel::CalVelocityConvergence(double& flow_rate_previous,dou
}
}
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax=Dm->Comm.sumReduce( vax_loc);
vay=Dm->Comm.sumReduce( vay_loc);
vaz=Dm->Comm.sumReduce( vaz_loc);
count=Dm->Comm.sumReduce( count_loc);
vax /= count;
vay /= count;
vaz /= count;
@ -574,7 +574,7 @@ void ScaLBL_StokesModel::Run(){
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
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_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
//************************************************************************/
if (timestep%1000==0){
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_Comm->Barrier(); MPI_Barrier(comm);
ScaLBL_Comm->Barrier(); comm.barrier();
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);
@ -648,10 +648,12 @@ void ScaLBL_StokesModel::Run(){
}
}
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax=Dm->Comm.sumReduce( vax_loc);
vay=Dm->Comm.sumReduce( vay_loc);
vaz=Dm->Comm.sumReduce( vaz_loc);
count=Dm->Comm.sumReduce( count_loc);
vax /= count;
vay /= count;
@ -714,51 +716,6 @@ void ScaLBL_StokesModel::Run(){
}
void ScaLBL_StokesModel::VelocityField(){
/* Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
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));
Morphology.Initialize();
Morphology.UpdateMeshValues();
Morphology.ComputeLocal();
Morphology.Reduce();
double count_loc=0;
double count;
double vax,vay,vaz;
double vax_loc,vay_loc,vaz_loc;
vax_loc = vay_loc = vaz_loc = 0.f;
for (int n=0; n<ScaLBL_Comm->LastExterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
for (int n=ScaLBL_Comm->FirstInterior(); n<ScaLBL_Comm->LastInterior(); n++){
vax_loc += VELOCITY[n];
vay_loc += VELOCITY[Np+n];
vaz_loc += VELOCITY[2*Np+n];
count_loc+=1.0;
}
MPI_Allreduce(&vax_loc,&vax,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vay_loc,&vay,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&vaz_loc,&vaz,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
MPI_Allreduce(&count_loc,&count,1,MPI_DOUBLE,MPI_SUM,Mask->Comm);
vax /= count;
vay /= count;
vaz /= count;
double mu = (tau-0.5)/3.f;
if (rank==0) printf("Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
if (rank==0) printf("%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",Fx, Fy, Fz, mu,
Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
*/
std::vector<IO::MeshDataStruct> visData;
fillHalo<double> fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1);

4
models/StokesModel.h
View File

@ -14,7 +14,7 @@
#include "common/ScaLBL.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "common/MPI.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
@ -76,7 +76,7 @@ public:
DoubleArray Velocity_y;
DoubleArray Velocity_z;
private:
MPI_Comm comm;
Utilities::MPI comm;
// filenames
char LocalRankString[8];