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

Merge branch 'morphLBM'

Browse Source
This commit is contained in:
JamesEMcclure 2019-10-10 09:12:38 -04:00
commit 6cc4b1deb2
24 changed files with 357 additions and 222 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
analysis/SubPhase.cpp
View File

@ -688,31 +688,12 @@ void SubPhase::AggregateLabels(char *FILENAME){
int nx = Dm->Nx;
int ny = Dm->Ny;
int nz = Dm->Nz;
int npx = Dm->nprocx();
int npy = Dm->nprocy();
int npz = Dm->nprocz();
int ipx = Dm->iproc();
int ipy = Dm->jproc();
int ipz = Dm->kproc();
int nprocs = Dm->nprocx()*Dm->nprocy()*Dm->nprocz();
int full_nx = npx*(nx-2);
int full_ny = npy*(ny-2);
int full_nz = npz*(nz-2);
int local_size = (nx-2)*(ny-2)*(nz-2);
long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
signed char *LocalID;
LocalID = new signed char [local_size];
//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
// assign the ID for the local sub-region
for (int k=1; k<nz-1; k++){
for (int j=1; j<ny-1; j++){
for (int i=1; i<nx-1; i++){
// assign the ID from the phase indicator field
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;
signed char local_id_val = Dm->id[n];
if (local_id_val > 0){
@ -720,63 +701,13 @@ void SubPhase::AggregateLabels(char *FILENAME){
if (value > 0.0) local_id_val = 1;
else local_id_val = 2;
}
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
Dm->id[n] = local_id_val;
}
}
}
MPI_Barrier(Dm->Comm);
// populate the FullID
if (Dm->rank() == 0){
signed char *FullID;
FullID = new signed char [full_size];
// first handle local ID for rank 0
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 x = i-1;
int y = j-1;
int z = k-1;
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
// next get the local ID from the other ranks
for (int rnk = 1; rnk<nprocs; rnk++){
ipz = rnk / (npx*npy);
ipy = (rnk - ipz*npx*npy) / npx;
ipx = (rnk - ipz*npx*npy - ipy*npx);
//printf("ipx=%i ipy=%i ipz=%i\n", ipx, ipy, ipz);
int tag = 15+rnk;
MPI_Recv(LocalID,local_size,MPI_CHAR,rnk,tag,Dm->Comm,MPI_STATUS_IGNORE);
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 x = i-1 + ipx*(nx-2);
int y = j-1 + ipy*(ny-2);
int z = k-1 + ipz*(nz-2);
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
}
// write the output
FILE *OUTFILE;
OUTFILE = fopen(FILENAME,"wb");
fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE);
}
else{
// send LocalID to rank=0
int tag = 15+ Dm->rank();
int dstrank = 0;
MPI_Send(LocalID,local_size,MPI_CHAR,dstrank,tag,Dm->Comm);
}
MPI_Barrier(Dm->Comm);
Dm->AggregateLabels(FILENAME);
}

27
analysis/runAnalysis.cpp
View File

@ -223,16 +223,16 @@ private:
class IOWorkItem: public ThreadPool::WorkItemRet<void>
{
public:
IOWorkItem( std::shared_ptr<Database> input_db_, std::vector<IO::MeshDataStruct>& visData_,
IOWorkItem(int timestep_, std::shared_ptr<Database> input_db_, std::vector<IO::MeshDataStruct>& visData_,
SubPhase& Averages_, fillHalo<double>& fillData_, runAnalysis::commWrapper&& comm_ ):
input_db(input_db_), visData(visData_), Averages(Averages_), fillData(fillData_), comm(std::move(comm_))
timestep(timestep_), input_db(input_db_), visData(visData_), Averages(Averages_), fillData(fillData_), comm(std::move(comm_))
{
}
~IOWorkItem() { }
virtual void run() {
auto color_db = input_db->getDatabase( "Color" );
auto vis_db = input_db->getDatabase( "Visualization" );
int timestep = color_db->getWithDefault<int>( "timestep", 0 );
// int timestep = color_db->getWithDefault<int>( "timestep", 0 );
PROFILE_START("Save Vis",1);
@ -573,7 +573,7 @@ runAnalysis::runAnalysis(std::shared_ptr<Database> input_db, const RankInfoStruc
IO::initialize("","silo","false");
// Create the MeshDataStruct
d_meshData.resize(1);
d_meshData[0].meshName = "domain";
d_meshData[0].mesh = std::make_shared<IO::DomainMesh>( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz );
auto PhaseVar = std::make_shared<IO::Variable>();
@ -778,13 +778,13 @@ AnalysisType runAnalysis::computeAnalysisType( int timestep )
/******************************************************************
* Run the analysis *
******************************************************************/
void runAnalysis::run( std::shared_ptr<Database> input_db, TwoPhase& Averages, const double *Phi,
void runAnalysis::run(int timestep, std::shared_ptr<Database> input_db, TwoPhase& Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den)
{
int N = d_N[0]*d_N[1]*d_N[2];
auto db = input_db->getDatabase( "Analysis" );
int timestep = db->getWithDefault<int>( "timestep", 0 );
//int timestep = db->getWithDefault<int>( "timestep", 0 );
// Check which analysis steps we need to perform
auto type = computeAnalysisType( timestep );
@ -950,15 +950,15 @@ void runAnalysis::run( std::shared_ptr<Database> input_db, TwoPhase& Averages, c
/******************************************************************
* Run the analysis *
******************************************************************/
void runAnalysis::basic( std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
void runAnalysis::basic(int timestep, std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
{
int N = d_N[0]*d_N[1]*d_N[2];
// Check which analysis steps we need to perform
auto color_db = input_db->getDatabase( "Color" );
//auto vis_db = input_db->getDatabase( "Visualization" );
auto vis_db = input_db->getDatabase( "Visualization" );
int timestep = color_db->getWithDefault<int>( "timestep", 0 );
//int timestep = color_db->getWithDefault<int>( "timestep", 0 );
auto type = computeAnalysisType( timestep );
if ( type == AnalysisType::AnalyzeNone )
return;
@ -1027,6 +1027,7 @@ void runAnalysis::basic( std::shared_ptr<Database> input_db, SubPhase &Averages,
ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double));
if (d_rank==0) {
color_db->putScalar<int>("timestep",timestep);
color_db->putScalar<bool>( "Restart", true );
input_db->putDatabase("Color", color_db);
std::ofstream OutStream("Restart.db");
@ -1043,7 +1044,7 @@ void runAnalysis::basic( std::shared_ptr<Database> input_db, SubPhase &Averages,
if (timestep%d_visualization_interval==0){
// Write the vis files
auto work = new IOWorkItem( input_db, d_meshData, Averages, d_fillData, getComm() );
auto work = new IOWorkItem( timestep, input_db, d_meshData, Averages, d_fillData, getComm() );
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_subphase);
work->add_dependency(d_wait_vis);
@ -1053,12 +1054,12 @@ void runAnalysis::basic( std::shared_ptr<Database> input_db, SubPhase &Averages,
PROFILE_STOP("run");
}
void runAnalysis::WriteVisData( std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
void runAnalysis::WriteVisData(int timestep, std::shared_ptr<Database> input_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den)
{
int N = d_N[0]*d_N[1]*d_N[2];
auto color_db = input_db->getDatabase( "Color" );
auto vis_db = input_db->getDatabase( "Visualization" );
int timestep = color_db->getWithDefault<int>( "timestep", 0 );
//int timestep = color_db->getWithDefault<int>( "timestep", 0 );
// Check which analysis steps we need to perform
auto type = computeAnalysisType( timestep );
@ -1077,7 +1078,7 @@ void runAnalysis::WriteVisData( std::shared_ptr<Database> input_db, SubPhase &Av
PROFILE_START("write vis",1);
// if (Averages.WriteVis == true){
auto work2 = new IOWorkItem( input_db, d_meshData, Averages, d_fillData, getComm() );
auto work2 = new IOWorkItem(timestep, input_db, d_meshData, Averages, d_fillData, getComm() );
work2->add_dependency(d_wait_vis);
d_wait_vis = d_tpool.add_work(work2);

8
analysis/runAnalysis.h
View File

@ -39,18 +39,18 @@ class runAnalysis
public:
//! Constructor
runAnalysis( std::shared_ptr<Database> db, const RankInfoStruct& rank_info,
runAnalysis(std::shared_ptr<Database> db, const RankInfoStruct& rank_info,
std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm, std::shared_ptr <Domain> dm, int Np, bool Regular, IntArray Map );
//! Destructor
~runAnalysis();
//! Run the next analysis
void run( std::shared_ptr<Database> db, TwoPhase &Averages, const double *Phi,
void run(int timestep, std::shared_ptr<Database> db, TwoPhase &Averages, const double *Phi,
double *Pressure, double *Velocity, double *fq, double *Den );
void basic( std::shared_ptr<Database> db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den );
void WriteVisData( std::shared_ptr<Database> vis_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den);
void basic( int timestep, std::shared_ptr<Database> db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den );
void WriteVisData(int timestep, std::shared_ptr<Database> vis_db, SubPhase &Averages, const double *Phi, double *Pressure, double *Velocity, double *fq, double *Den);
//! Finish all active analysis
void finish();

122
common/Domain.cpp
View File

@ -46,6 +46,7 @@ Domain::Domain( int nx, int ny, int nz, int rnk, int npx, int npy, int npz,
Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0), voxel_length(1),
Comm(MPI_COMM_WORLD),
inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
inlet_layers_phase(1),outlet_layers_phase(2),
sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0),
sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0),
sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0),
@ -93,6 +94,7 @@ Domain::Domain( std::shared_ptr<Database> db, MPI_Comm Communicator):
Comm(MPI_COMM_NULL),
inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
outlet_layers_x(0), outlet_layers_y(0), outlet_layers_z(0),
inlet_layers_phase(1),outlet_layers_phase(2),
sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0),
sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0),
sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0),
@ -199,6 +201,12 @@ void Domain::initialize( std::shared_ptr<Database> db )
outlet_layers_y = OutletCount[1];
outlet_layers_z = OutletCount[2];
}
if (d_db->keyExists( "InletLayersPhase" )){
inlet_layers_phase = d_db->getScalar<int>( "InletLayersPhase" );
}
if (d_db->keyExists( "OutletLayersPhase" )){
outlet_layers_phase = d_db->getScalar<int>( "OutletLayersPhase" );
}
voxel_length = 1.0;
if (d_db->keyExists( "voxel_length" )){
voxel_length = d_db->getScalar<double>("voxel_length");
@ -263,6 +271,8 @@ void Domain::Decomp(std::string Filename)
outlet_layers_x = 0;
outlet_layers_y = 0;
outlet_layers_z = 0;
inlet_layers_phase=1;
outlet_layers_phase=2;
checkerSize = 32;
// Read domain parameters
@ -295,6 +305,12 @@ void Domain::Decomp(std::string Filename)
else {
checkerSize = SIZE[0];
}
if (database->keyExists( "InletLayersPhase" )){
inlet_layers_phase = database->getScalar<int>( "InletLayersPhase" );
}
if (database->keyExists( "OutletLayersPhase" )){
outlet_layers_phase = database->getScalar<int>( "OutletLayersPhase" );
}
auto ReadValues = database->getVector<int>( "ReadValues" );
auto WriteValues = database->getVector<int>( "WriteValues" );
auto ReadType = database->getScalar<std::string>( "ReadType" );
@ -361,10 +377,10 @@ void Domain::Decomp(std::string Filename)
// relabel the data
std::vector<long int> LabelCount(ReadValues.size(),0);
for (int k = 0; k<Nz; k++){
for (int j = 0; j<Ny; j++){
for (int i = 0; i<Nx; i++){
n = k*Nx*Ny+j*Nx+i;
for (int k = 0; k<global_Nz; k++){
for (int j = 0; j<global_Ny; j++){
for (int i = 0; i<global_Nx; i++){
n = k*global_Nx*global_Ny+j*global_Nx+i;
//char locval = loc_id[n];
char locval = SegData[n];
for (int idx=0; idx<ReadValues.size(); idx++){
@ -433,7 +449,8 @@ void Domain::Decomp(std::string Filename)
for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize+j/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
//SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = inlet_layers_phase;
}
else{
// solid checkers
@ -490,7 +507,8 @@ void Domain::Decomp(std::string Filename)
for (int i = 0; i<global_Nx; i++){
if ( (i/checkerSize+j/checkerSize)%2 == 0){
// void checkers
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
//SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 2;
SegData[k*global_Nx*global_Ny+j*global_Nx+i] = outlet_layers_phase;
}
else{
// solid checkers
@ -577,6 +595,98 @@ void Domain::Decomp(std::string Filename)
MPI_Barrier(Comm);
}
void Domain::AggregateLabels(char *FILENAME){
int nx = Nx;
int ny = Ny;
int nz = Nz;
int npx = nprocx();
int npy = nprocy();
int npz = nprocz();
int ipx = iproc();
int ipy = jproc();
int ipz = kproc();
int nprocs = nprocx()*nprocy()*nprocz();
int full_nx = npx*(nx-2);
int full_ny = npy*(ny-2);
int full_nz = npz*(nz-2);
int local_size = (nx-2)*(ny-2)*(nz-2);
long int full_size = long(full_nx)*long(full_ny)*long(full_nz);
signed char *LocalID;
LocalID = new signed char [local_size];
//printf("aggregate labels: local size=%i, global size = %i",local_size, full_size);
// assign the ID for the local sub-region
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 n = k*nx*ny+j*nx+i;
signed char local_id_val = id[n];
LocalID[(k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1] = local_id_val;
}
}
}
MPI_Barrier(Comm);
// populate the FullID
if (rank() == 0){
signed char *FullID;
FullID = new signed char [full_size];
// first handle local ID for rank 0
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 x = i-1;
int y = j-1;
int z = k-1;
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
// next get the local ID from the other ranks
for (int rnk = 1; rnk<nprocs; rnk++){
ipz = rnk / (npx*npy);
ipy = (rnk - ipz*npx*npy) / npx;
ipx = (rnk - ipz*npx*npy - ipy*npx);
//printf("ipx=%i ipy=%i ipz=%i\n", ipx, ipy, ipz);
int tag = 15+rnk;
MPI_Recv(LocalID,local_size,MPI_CHAR,rnk,tag,Comm,MPI_STATUS_IGNORE);
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 x = i-1 + ipx*(nx-2);
int y = j-1 + ipy*(ny-2);
int z = k-1 + ipz*(nz-2);
int n_local = (k-1)*(nx-2)*(ny-2) + (j-1)*(nx-2) + i-1;
int n_full = z*full_nx*full_ny + y*full_nx + x;
FullID[n_full] = LocalID[n_local];
}
}
}
}
// write the output
FILE *OUTFILE;
OUTFILE = fopen(FILENAME,"wb");
fwrite(FullID,1,full_size,OUTFILE);
fclose(OUTFILE);
}
else{
// send LocalID to rank=0
int tag = 15+ rank();
int dstrank = 0;
MPI_Send(LocalID,local_size,MPI_CHAR,dstrank,tag,Comm);
}
MPI_Barrier(Comm);
}
/********************************************************
* Initialize communication *
********************************************************/

3
common/Domain.h
View File

@ -127,6 +127,8 @@ public: // Public variables (need to create accessors instead)
int Nx,Ny,Nz,N;
int inlet_layers_x, inlet_layers_y, inlet_layers_z;
int outlet_layers_x, outlet_layers_y, outlet_layers_z;
int inlet_layers_phase; //as usual: 1->n, 2->w
int outlet_layers_phase;
double porosity;
RankInfoStruct rank_info;
@ -195,6 +197,7 @@ public: // Public variables (need to create accessors instead)
void CommunicateMeshHalo(DoubleArray &Mesh);
void CommInit();
int PoreCount();
void AggregateLabels(char *FILENAME);
private:

3
example/Bubble/input.db
View File

@ -30,4 +30,5 @@ Analysis {
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
}

4
example/CMakeLists.txt
View File

@ -1,5 +1,5 @@
# Copy the examples to the install folder
INSTALL_EXAMPLE (NonNewtonianChannelFlow )
#INSTALL_EXAMPLE (NonNewtonianChannelFlow )
INSTALL_EXAMPLE( Bubble )
INSTALL_EXAMPLE( ConstrainedBubble )
INSTALL_EXAMPLE( Piston )
@ -10,7 +10,7 @@ INSTALL_EXAMPLE( Sph1896 )
INSTALL_EXAMPLE( drainage )
INSTALL_EXAMPLE( imbibition )
INSTALL_EXAMPLE( relperm )
INSTALL_EXAMPLE( Poiseuille )
#INSTALL_EXAMPLE( Poiseuille )
INSTALL_EXAMPLE( Juanes )
INSTALL_EXAMPLE( MicroModel )
INSTALL_EXAMPLE( CornerFlow )

2
example/CornerFlow/input.db
View File

@ -47,3 +47,5 @@ Analysis {
}
Visualization {
}

6
example/InkBottle/Color.in
View File

@ -1,6 +0,0 @@
1.0
1.0e-2 0.95 -1.0
0.7
0.0 0.0 0.0
0 1 1.01 0.99
50000 20000 1e-5

3
example/InkBottle/Domain.in
View File

@ -1,3 +0,0 @@
1 1 1
80 80 80
1.0 1.0 1.0

2
example/Piston/poiseuille.db
View File

@ -18,3 +18,5 @@ Domain {
BC = 0 // Boundary condition type
}
Visualization {
}

11
example/Poiseuille/README.md
View File

@ -1,11 +0,0 @@
# 1. Edit the Domain.in file based on the desired system size
# 2. Run the pre-processor
export TUBEWIDTH=20
export LAYERWIDTH=0
mpirun -np 1 ../../tests/lbpm_plates_pp $TUBEWIDTH $LAYERWIDTH
# 3. Run single-phase simulation within the domain
mpirun -np 1 ../../tests/lbpm_permeability_simulator

2
example/Sph125/input.db
View File

@ -37,3 +37,5 @@ Analysis {
}
Visualization {
}

3
example/Sph1896/input.db
View File

@ -51,3 +51,6 @@ Analysis {
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}
Visualization {
}

2
example/SquareTube/input.db
View File

@ -43,3 +43,5 @@ Analysis {
}
Visualization {
}

53
example/systems/summit/input.db Normal file
View File

@ -0,0 +1,53 @@
B1;95;0cMRT {
tau = 1.0 // relaxation time
F = 0, 0, 1e-4 // external body force applied to system
timestepMax = 1000 // max number of timesteps
din = 1.0
dout = 1.0
Restart = false
flux = 0.0
}
Color {
tauA = 0.7; // relaxation time for fluid A
tauB = 0.7; // relaxation time for fluid B
rhoA = 1.0; // mass density for fluid A
rhoB = 1.0; // mass density for fluid B
alpha = 1e-3; // controls interfacial tension between fluids
beta = 0.95; // controls interface width
F = 0, 0, 1.0e-5 // external body force applied to the system
Restart = false // restart from checkpoint file?
din = 1.0 // density at inlet (if external BC is applied)
dout = 1.0 // density at outlet (if external BC is applied )
timestepMax = 3000 // maximum number of timesteps to simulate
flux = 0.0 // volumetric flux in voxels per timestep (if flux BC is applied)
ComponentLabels = 0 // comma separated list of solid mineral labels
ComponentAffinity = -1.0 // comma separated list of phase indicato field value to assign for each mineral label
}
Domain {
nproc = 1, 1, 1 // Number of processors (Npx,Npy,Npz)
n = 400, 400, 400 // Size of local domain (Nx,Ny,Nz)
nspheres = 1896 // Number of spheres (only needed if using a sphere packing)
L = 1, 1, 1 // Length of domain (x,y,z)
BC = 0 // Boundary condition type
// BC = 0 for periodic BC
// BC = 1 for pressure BC (applied in z direction)
// BC = 4 for flux BC (applied in z direction
Filename = "mineralmodel.raw" // name of raw binary file to read digital image from
ReadType = "8bit"
ReadValues = 0, 1, 2, 3, 4, 5, 6, 7, 8 // list of labels within the binary file (read)
WriteValues = 1, -1, -2, -2, -3, -3, -4, -4, -4 // list of labels within the output files (write)
}
Analysis {
blobid_interval = 1000 // Frequency to perform blob identification
analysis_interval = 1000 // Frequency to perform analysis
restart_interval = 2000 // Frequency to write restart data
visualization_interval = 2000 // Frequency to write visualization data
restart_file = "Restart" // Filename to use for restart file (will append rank)
N_threads = 4 // Number of threads to use
load_balance = "independent" // Load balance method to use: "none", "default", "independent"
}

8
example/systems/summit/summit-spheres-1g.lsf
View File

@ -1,5 +1,5 @@
#!/bin/bash
#BSUB -P CSC275MCCLURE
#BSUB -P CSC380
#BSUB -J spheres
#BSUB -o spheres.o%J
#BSUB -W 10
@ -10,10 +10,9 @@ date
module load gcc cuda
export SCALBL_DIR=$HOME/summit/build/LBPM-WIA/tests
jsrun -n1 -r1 -g1 -c1 -brs $SCALBL_DIR/GenerateSphereTest 1896
export LBPM_DIR=/ccs/proj/csc380/mcclurej/install/LBPM/tests
jsrun -n1 -r1 -g1 -c1 -brs $LBPM_DIR/GenerateSphereTest input.db
# Create the 1200 GPU case
@ -45,4 +44,3 @@ for i in `seq -w 0 $NRANKS`; do distfile="$BASEDIST$i"; echo $distfile; cp SignD
exit;

20
models/ColorModel.cpp
View File

@ -622,7 +622,7 @@ void ScaLBL_ColorModel::Run(){
printf(" morph_delta = %f \n",morph_delta);
}
else if (protocol == "shell aggregation"){
printf(" using protocol = shell aggregation \n");
printf(" using protocol = shell aggregation \n");
printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS);
printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS);
printf(" tolerance = %f \n",tolerance);
@ -724,17 +724,17 @@ void ScaLBL_ColorModel::Run(){
MPI_Barrier(comm);
PROFILE_STOP("Update");
// Run the analysis
//analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
color_db->putScalar<int>("timestep",timestep);
current_db->putDatabase("Color", color_db);
analysis.basic( current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
printf("....inlet pressure=%f \n",din);
printf("%i %f \n",timestep,din);
}
// Run the analysis
//analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
//color_db->putScalar<int>("timestep",timestep);
//current_db->putDatabase("Color", color_db);
analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
// allow initial ramp-up to get closer to steady state
if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
@ -784,7 +784,7 @@ void ScaLBL_ColorModel::Run(){
delta_volume_target = Dm->Volume*volA *morph_delta; // set target volume change
Averages->Full();
Averages->Write(timestep);
analysis.WriteVisData( current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.finish();
if (rank==0){

32
models/DFHModel.cpp
View File

@ -60,20 +60,22 @@ void ScaLBL_DFHModel::ReadParams(string filename){
analysis_db = db->getDatabase( "Analysis" );
// Color Model parameters
timestepMax = color_db->getScalar<int>( "timestepMax" );
tauA = color_db->getScalar<double>( "tauA" );
tauB = color_db->getScalar<double>( "tauB" );
rhoA = color_db->getScalar<double>( "rhoA" );
rhoB = color_db->getScalar<double>( "rhoB" );
Fx = color_db->getVector<double>( "F" )[0];
Fy = color_db->getVector<double>( "F" )[1];
Fz = color_db->getVector<double>( "F" )[2];
alpha = color_db->getScalar<double>( "alpha" );
beta = color_db->getScalar<double>( "beta" );
Restart = color_db->getScalar<bool>( "Restart" );
din = color_db->getScalar<double>( "din" );
dout = color_db->getScalar<double>( "dout" );
flux = color_db->getScalar<double>( "flux" );
timestepMax = color_db->getWithDefault<int>( "timestepMax", 100 );
tauA = color_db->getWithDefault<double>( "tauA", 1.0 );
tauB = color_db->getWithDefault<double>( "tauB", 1.0 );
rhoA = color_db->getWithDefault<double>( "rhoA", 1.0 );
rhoB = color_db->getWithDefault<double>( "rhoB", 1.0 );
alpha = color_db->getWithDefault<double>( "alpha", 0.001 );
beta = color_db->getWithDefault<double>( "beta", 0.95 );
Restart = color_db->getWithDefault<bool>( "Restart", true );
din = color_db->getWithDefault<double>( "din", 1.0 );
dout = color_db->getWithDefault<double>( "dout", 1.0 );
flux = color_db->getWithDefault<double>( "flux", 0.0 );
if (color_db->keyExists( "F" )){
Fx = color_db->getVector<double>( "F" )[0];
Fy = color_db->getVector<double>( "F" )[1];
Fz = color_db->getVector<double>( "F" )[2];
}
inletA=1.f;
inletB=0.f;
outletA=0.f;
@ -573,7 +575,7 @@ void ScaLBL_DFHModel::Run(){
PROFILE_STOP("Update");
// Run the analysis
analysis.run( analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.run(timestep, analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den );
}
analysis.finish();
PROFILE_STOP("Loop");

11
sample_scripts/configure_summit
View File

@ -5,6 +5,9 @@
module load cmake gcc
module load cuda
export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
# configure
rm -rf CMake*
cmake \
@ -20,11 +23,11 @@ cmake \
-D MPIEXEC=mpirun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_HDF5=1 \
-D HDF5_DIRECTORY="/ccs/proj/bip176/TPLs/hdf5/" \
-D HDF5_LIB="/ccs/proj/bip176/TPLs/hdf5/lib/libhdf5.a" \
-D HDF5_DIRECTORY="$HDF5_DIR" \
-D HDF5_LIB="$HDF5_DIR/lib/libhdf5.a" \
-D USE_SILO=1 \
-D SILO_LIB="/ccs/proj/bip176/TPLs/silo/lib/libsiloh5.a" \
-D SILO_DIRECTORY="/ccs/proj/bip176/TPLs/silo/" \
-D SILO_LIB="$SILO_DIR/lib/libsiloh5.a" \
-D SILO_DIRECTORY="$SILO_DIR" \
-D USE_DOXYGEN:BOOL=false \
-D USE_TIMER=0 \
~/LBPM-WIA

24
tests/TestBubbleDFH.cpp
View File

@ -86,18 +86,18 @@ int main(int argc, char **argv)
// Color Model parameters
int timestepMax = color_db->getScalar<int>( "timestepMax" );
double tauA = color_db->getScalar<double>( "tauA" );
double tauB = color_db->getScalar<double>( "tauB" );
double rhoA = color_db->getScalar<double>( "rhoA" );
double rhoB = color_db->getScalar<double>( "rhoB" );
double tauA = color_db->getWithDefault<double>( "tauA", 1.0 );
double tauB = color_db->getWithDefault<double>( "tauB", 1.0 );
double rhoA = color_db->getWithDefault<double>( "rhoA", 1.0 );
double rhoB = color_db->getWithDefault<double>( "rhoB", 1.0 );
double Fx = color_db->getVector<double>( "F" )[0];
double Fy = color_db->getVector<double>( "F" )[1];
double Fz = color_db->getVector<double>( "F" )[2];
double alpha = color_db->getScalar<double>( "alpha" );
double beta = color_db->getScalar<double>( "beta" );
bool Restart = color_db->getScalar<bool>( "Restart" );
double din = color_db->getScalar<double>( "din" );
double dout = color_db->getScalar<double>( "dout" );;
double alpha = color_db->getWithDefault<double>( "alpha", 0.001 );
double beta = color_db->getWithDefault<double>( "beta", 0.95 );
bool Restart = color_db->getWithDefault<bool>( "Restart", false );
double din = color_db->getWithDefault<double>( "din", 1.0 );
double dout = color_db->getWithDefault<double>( "dout", 1.0 );;
double inletA=1.f;
double inletB=0.f;
double outletA=0.f;
@ -108,7 +108,7 @@ int main(int argc, char **argv)
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
int BoundaryCondition = domain_db->getScalar<int>( "BC" );
int BoundaryCondition = domain_db->getWithDefault<int>( "BC", 0 );
int Nx = size[0];
int Ny = size[1];
int Nz = size[2];
@ -405,7 +405,7 @@ int main(int argc, char **argv)
//************ MAIN ITERATION LOOP ***************************************/
PROFILE_START("Loop");
//std::shared_ptr<Database> analysis_db;
runAnalysis analysis( analysis_db, rank_info, ScaLBL_Comm, Dm, Np, pBC, Map );
runAnalysis analysis(db, rank_info, ScaLBL_Comm, Dm, Np, pBC, Map );
//analysis.createThreads( analysis_method, 4 );
while (timestep < timestepMax && err > tol ) {
//if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); }
@ -487,7 +487,7 @@ int main(int argc, char **argv)
PROFILE_STOP("Update");
// Run the analysis
analysis.run( analysis_db, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.run( timestep, db, *Averages, Phi, Pressure, Velocity, fq, Den );
}
analysis.finish();

56
tests/lbpm_morph_pp.cpp
View File

@ -32,22 +32,15 @@ int main(int argc, char **argv)
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int i,j,k,n;
int BC=0;
// char fluidValue,solidValue;
int MAXTIME=1000;
int READ_FROM_BLOCK=0;
char LocalRankString[8];
int n, nx, ny, nz;
char LocalRankFilename[40];
char FILENAME[128];
string filename;
double Rcrit_new, SW;
double SW;
if (argc > 1){
filename=argv[1];
Rcrit_new=0.f;
//Rcrit_new=0.f;
//SW=strtod(argv[2],NULL);
}
else ERROR("No input database provided\n");
@ -61,6 +54,7 @@ int main(int argc, char **argv)
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
SW = domain_db->getScalar<double>("Sw");
auto READFILE = domain_db->getScalar<std::string>( "Filename" );
// Generate the NWP configuration
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
@ -70,28 +64,21 @@ int main(int argc, char **argv)
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
int N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit();
signed char *id;
signed char *id_connected;
id = new signed char [N];
signed char *id_connected;
id_connected = new signed char [N];
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_morph_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
Mask->Decomp(READFILE);
Mask->CommInit();
nx+=2; ny+=2; nz+=2;
// Generate the signed distance map
@ -105,6 +92,7 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n=k*nx*ny+j*nx+i;
id[n] = Mask->id[n];
if (id[n] == 1){
phase(i,j,k) = 1.0;
}
@ -120,8 +108,11 @@ int main(int argc, char **argv)
for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
// Initialize the solid phase
if (id[n] > 0) id_solid(i,j,k) = 1;
else id_solid(i,j,k) = 0;
if (Mask->id[n] > 0){
id_solid(i,j,k) = 1;
}
else
id_solid(i,j,k) = 0;
}
}
}
@ -129,7 +120,6 @@ int main(int argc, char **argv)
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;
// Initialize distance to +/- 1
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
@ -237,6 +227,22 @@ int main(int argc, char **argv)
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);
fclose(ID);
// write the geometry to a single file
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;
Mask->id[n] = id[n];
}
}
}
MPI_Barrier(comm);
sprintf(FILENAME,READFILE.c_str());
sprintf(FILENAME+strlen(FILENAME),".morph.raw");
if (rank==0) printf("Writing file to: %s \n", FILENAME);
Mask->AggregateLabels(FILENAME);
}
MPI_Barrier(comm);

44
tests/lbpm_morphdrain_pp.cpp
View File

@ -32,19 +32,13 @@ int main(int argc, char **argv)
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int i,j,k,n;
int BC=0;
// char fluidValue,solidValue;
int MAXTIME=1000;
int READ_FROM_BLOCK=0;
int n, nprocx, nprocy, nprocz, nx, ny, nz;
char LocalRankString[8];
char LocalRankFilename[40];
char FILENAME[128];
string filename;
double Rcrit_new, SW;
double SW,Rcrit_new;
if (argc > 1){
filename=argv[1];
Rcrit_new=0.f;
@ -61,6 +55,7 @@ int main(int argc, char **argv)
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
SW = domain_db->getScalar<double>("Sw");
auto READFILE = domain_db->getScalar<std::string>( "Filename" );
// Generate the NWP configuration
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
@ -77,19 +72,24 @@ int main(int argc, char **argv)
int N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit();
signed char *id;
id = new signed char [N];
sprintf(LocalRankFilename,"ID.%05i",rank);
Mask->Decomp(READFILE);
Mask->CommInit();
/*sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_morphopen_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
*/
nx+=2; ny+=2; nz+=2;
// Generate the signed distance map
@ -102,9 +102,13 @@ int main(int argc, char **argv)
for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
id[n] = Mask->id[n];
// Initialize the solid phase
if (id[n] > 0) id_solid(i,j,k) = 1;
else id_solid(i,j,k) = 0;
if (Mask->id[n] > 0){
id_solid(i,j,k) = 1;
}
else
id_solid(i,j,k) = 0;
}
}
}
@ -188,6 +192,22 @@ int main(int argc, char **argv)
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);
fclose(ID);
// write the geometry to a single file
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;
Mask->id[n] = id[n];
}
}
}
MPI_Barrier(comm);
sprintf(FILENAME,READFILE.c_str());
sprintf(FILENAME+strlen(FILENAME),".morphdrain.raw");
if (rank==0) printf("Writing file to: %s \n", FILENAME);
Mask->AggregateLabels(FILENAME);
}
MPI_Barrier(comm);

52
tests/lbpm_morphopen_pp.cpp
View File

@ -32,19 +32,13 @@ int main(int argc, char **argv)
//.......................................................................
// Reading the domain information file
//.......................................................................
int nprocx, nprocy, nprocz, nx, ny, nz, nspheres;
double Lx, Ly, Lz;
int i,j,k,n;
int BC=0;
// char fluidValue,solidValue;
int MAXTIME=1000;
int READ_FROM_BLOCK=0;
int n, nprocx, nprocy, nprocz, nx, ny, nz;
char LocalRankString[8];
char LocalRankFilename[40];
char FILENAME[128];
string filename;
double Rcrit_new, SW;
double SW,Rcrit_new;
if (argc > 1){
filename=argv[1];
Rcrit_new=0.f;
@ -56,6 +50,7 @@ int main(int argc, char **argv)
auto domain_db = db->getDatabase( "Domain" );
// Read domain parameters
auto READFILE = domain_db->getScalar<std::string>( "Filename" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
@ -84,19 +79,20 @@ int main(int argc, char **argv)
int N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
std::shared_ptr<Domain> Mask (new Domain(domain_db,comm));
// std::shared_ptr<Domain> Dm (new Domain(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BC));
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit();
signed char *id;
id = new signed char [N];
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("lbpm_morphopen_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
Mask->Decomp(READFILE);
Mask->CommInit();
// Generate the NWP configuration
//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);
// GenerateResidual(id,nx,ny,nz,Saturation);
nx+=2; ny+=2; nz+=2;
// Generate the signed distance map
@ -109,9 +105,13 @@ int main(int argc, char **argv)
for (int j=0;j<ny;j++){
for (int i=0;i<nx;i++){
int n = k*nx*ny+j*nx+i;
id[n] = Mask->id[n];
// Initialize the solid phase
if (id[n] > 0) id_solid(i,j,k) = 1;
else id_solid(i,j,k) = 0;
if (Mask->id[n] > 0){
id_solid(i,j,k) = 1;
}
else
id_solid(i,j,k) = 0;
}
}
}
@ -195,6 +195,22 @@ int main(int argc, char **argv)
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);
fclose(ID);
// write the geometry to a single file
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;
Mask->id[n] = id[n];
}
}
}
MPI_Barrier(comm);
sprintf(FILENAME,READFILE.c_str());
sprintf(FILENAME+strlen(FILENAME),".morphopen.raw");
if (rank==0) printf("Writing file to: %s \n", FILENAME);
Mask->AggregateLabels(FILENAME);
}
MPI_Barrier(comm);