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

merged to master -- fixed runAnalysis

Browse Source
This commit is contained in:
James E McClure
2018-09-07 14:10:57 -04:00
parent b9b07bc4db 40cc07d508
commit a990bd2f75
25 changed files with 865 additions and 374 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
analysis/TwoPhase.cpp
View File

@@ -394,7 +394,7 @@ void TwoPhase::UpdateMeshValues()
for (j=0; j<Ny; j++){
for (i=0; i<Nx; i++){
n = k*Nx*Ny+j*Nx+i;
if (Dm->id[n] == 0){
if (!(Dm->id[n] > 0)){
// Solid phase
PhaseID(i,j,k) = 0;
}
@@ -431,7 +431,7 @@ void TwoPhase::ComputeLocal()
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm->id[n] != 0 ){
if ( Dm->id[n] > 0 ){
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution
@@ -624,7 +624,7 @@ void TwoPhase::ComponentAverages()
// Compute volume averages
for (int p=0;p<8;p++){
n = i+cube[p][0] + (j+cube[p][1])*Nx + (k+cube[p][2])*Nx*Ny;
if ( Dm->id[n] != 0 ){
if ( Dm->id[n] > 0 ){
// 1-D index for this cube corner
// compute the norm of the gradient of the phase indicator field
// Compute the non-wetting phase volume contribution

2
analysis/distance.cpp
View File

@@ -153,7 +153,7 @@ void CalcVecDist( Array<Vec> &d, const Array<int> &ID0, const Domain &Dm,
for (int i=0; i<N[0]; i++) {
for (int j=0; j<N[1]; j++) {
ID(i,j,0) = ID(i,j,1);
ID(i,j,N[1]-1) = ID(i,j,N[2]-2);
ID(i,j,N[2]-1) = ID(i,j,N[2]-2);
}
}
// Communicate ghosts

305
analysis/runAnalysis.cpp
View File

@@ -13,7 +13,6 @@
#include "ProfilerApp.h"
AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)
{
lhs = static_cast<AnalysisType>(
@@ -40,16 +39,21 @@ void DeleteArray( const TYPE *p )
class WriteRestartWorkItem: public ThreadPool::WorkItemRet<void>
{
public:
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cphi_, std::shared_ptr<double> cfq_, int N_ ):
filename(filename_), cphi(cphi_), cfq(cfq_), N(N_) {}
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_, std::shared_ptr<double> cfq_, int N_ ):
filename(filename_), cDen(cDen_), cfq(cfq_), N(N_) {}
virtual void run() {
PROFILE_START("Save Checkpoint",1);
double value;
ofstream File(filename,ios::binary);
for (int n=0; n<N; n++){
// Write the two density values
value = cphi.get()[n];
value = cDen.get()[n];
File.write((char*) &value, sizeof(value));
value = cDen.get()[N+n];
File.write((char*) &value, sizeof(value));
}
for (int n=0; n<N; n++){
// Write the distributions
for (int q=0; q<19; q++){
value = cfq.get()[q*N+n];
@@ -60,12 +64,12 @@ public:
PROFILE_STOP("Save Checkpoint",1);
};
private:
WriteRestartWorkItem();
const char* filename;
std::shared_ptr<double> cfq,cphi;
// const DoubleArray& phase;
//const DoubleArray& dist;
const int N;
WriteRestartWorkItem();
const char* filename;
std::shared_ptr<double> cfq,cDen;
// const DoubleArray& phase;
//const DoubleArray& dist;
const int N;
};
@@ -153,29 +157,38 @@ private:
class WriteVisWorkItem: public ThreadPool::WorkItemRet<void>
{
public:
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
TwoPhase& Avgerages_, fillHalo<double>& fillData_, runAnalysis::commWrapper&& comm_ ):
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_))
{
}
~WriteVisWorkItem() { }
virtual void run() {
PROFILE_START("Save Vis",1);
ASSERT(visData[0].vars[0]->name=="phase");
ASSERT(visData[0].vars[1]->name=="Pressure");
ASSERT(visData[0].vars[2]->name=="SignDist");
ASSERT(visData[0].vars[3]->name=="BlobID");
Array<double>& PhaseData = visData[0].vars[0]->data;
Array<double>& PressData = visData[0].vars[1]->data;
Array<double>& SignData = visData[0].vars[2]->data;
Array<double>& BlobData = visData[0].vars[3]->data;
fillData.copy(Averages.SDn,PhaseData);
fillData.copy(Averages.Press,PressData);
fillData.copy(Averages.SDs,SignData);
fillData.copy(Averages.Label_NWP,BlobData);
IO::writeData( timestep, visData, comm.comm );
PROFILE_STOP("Save Vis",1);
};
WriteVisWorkItem( int timestep_, std::vector<IO::MeshDataStruct>& visData_,
TwoPhase& Avgerages_, fillHalo<double>& fillData_, runAnalysis::commWrapper&& comm_ ):
timestep(timestep_), visData(visData_), Averages(Avgerages_), fillData(fillData_), comm(std::move(comm_))
{
}
~WriteVisWorkItem() { }
virtual void run() {
PROFILE_START("Save Vis",1);
ASSERT(visData[0].vars[0]->name=="phase");
ASSERT(visData[0].vars[1]->name=="Pressure");
ASSERT(visData[0].vars[2]->name=="Velocity_x");
ASSERT(visData[0].vars[3]->name=="Velocity_y");
ASSERT(visData[0].vars[4]->name=="Velocity_z");
ASSERT(visData[0].vars[5]->name=="SignDist");
ASSERT(visData[0].vars[6]->name=="BlobID");
Array<double>& PhaseData = visData[0].vars[0]->data;
Array<double>& PressData = visData[0].vars[1]->data;
Array<double>& VelxData = visData[0].vars[2]->data;
Array<double>& VelyData = visData[0].vars[3]->data;
Array<double>& VelzData = visData[0].vars[4]->data;
Array<double>& SignData = visData[0].vars[5]->data;
Array<double>& BlobData = visData[0].vars[6]->data;
fillData.copy(Averages.SDn,PhaseData);
fillData.copy(Averages.Press,PressData);
fillData.copy(Averages.SDs,SignData);
fillData.copy(Averages.Vel_x,VelxData);
fillData.copy(Averages.Vel_y,VelyData);
fillData.copy(Averages.Vel_z,VelzData);
fillData.copy(Averages.Label_NWP,BlobData);
IO::writeData( timestep, visData, comm.comm );
PROFILE_STOP("Save Vis",1);
};
private:
WriteVisWorkItem();
int timestep;
@@ -315,8 +328,12 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
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>();
auto PressVar = std::make_shared<IO::Variable>();
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
auto BlobIDVar = std::make_shared<IO::Variable>();
PhaseVar->name = "phase";
PhaseVar->type = IO::VariableType::VolumeVariable;
PhaseVar->dim = 1;
@@ -327,6 +344,23 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
PressVar->dim = 1;
PressVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(PressVar);
VxVar->name = "Velocity_x";
VxVar->type = IO::VariableType::VolumeVariable;
VxVar->dim = 1;
VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(VxVar);
VyVar->name = "Velocity_y";
VyVar->type = IO::VariableType::VolumeVariable;
VyVar->dim = 1;
VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(VyVar);
VzVar->name = "Velocity_z";
VzVar->type = IO::VariableType::VolumeVariable;
VzVar->dim = 1;
VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
d_meshData[0].vars.push_back(VzVar);
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VariableType::VolumeVariable;
SignDistVar->dim = 1;
@@ -531,113 +565,116 @@ void runAnalysis::run( int timestep, TwoPhase& Averages, const double *Phi,
}
delete [] TmpDat;
}
*/
//if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus);
else
ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
//memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
}
if ( timestep%d_analysis_interval == 0 ) {
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus);
else
ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double));
//memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
}
//if ( matches(type,AnalysisType::CopySimState) ) {
if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
// Copy the members of Averages to the cpu (phase was copied above)
PROFILE_START("Copy-Pressure",1);
ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
ScaLBL_DeviceBarrier();
PROFILE_STOP("Copy-Pressure",1);
PROFILE_START("Copy-Wait",1);
PROFILE_STOP("Copy-Wait",1);
PROFILE_START("Copy-State",1);
//memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
else
ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
// copy other variables
d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
PROFILE_STOP("Copy-State",1);
}
std::shared_ptr<double> cfq,cPhi;
//if ( matches(type,AnalysisType::CreateRestart) ) {
if (timestep%d_restart_interval==0){
// Copy restart data to the CPU
cPhi = std::shared_ptr<double>(new double[d_Np],DeleteArray<double>);
cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
ScaLBL_CopyToHost(cPhi.get(),Phi,d_Np*sizeof(double));
}
PROFILE_STOP("Copy data to host",1);
*/
//if ( matches(type,AnalysisType::CopyPhaseIndicator) ) {
if ( timestep%d_analysis_interval + 8 == d_analysis_interval ) {
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tplus);
else
ScaLBL_CopyToHost(Averages.Phase_tplus.data(),Phi,N*sizeof(double));
//memcpy(Averages.Phase_tplus.data(),phase->data(),N*sizeof(double));
}
if ( timestep%d_analysis_interval == 0 ) {
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase_tminus);
else
ScaLBL_CopyToHost(Averages.Phase_tminus.data(),Phi,N*sizeof(double));
//memcpy(Averages.Phase_tminus.data(),phase->data(),N*sizeof(double));
}
//if ( matches(type,AnalysisType::CopySimState) ) {
if ( timestep%d_analysis_interval + 4 == d_analysis_interval ) {
// Copy the members of Averages to the cpu (phase was copied above)
PROFILE_START("Copy-Pressure",1);
ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
//ScaLBL_D3Q19_Momentum(fq,Velocity,d_Np);
ScaLBL_DeviceBarrier();
PROFILE_STOP("Copy-Pressure",1);
PROFILE_START("Copy-Wait",1);
PROFILE_STOP("Copy-Wait",1);
PROFILE_START("Copy-State",1);
//memcpy(Averages.Phase.data(),phase->data(),N*sizeof(double));
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,Averages.Phase);
else
ScaLBL_CopyToHost(Averages.Phase.data(),Phi,N*sizeof(double));
// copy other variables
d_ScaLBL_Comm->RegularLayout(d_Map,Pressure,Averages.Press);
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[0],Averages.Vel_x);
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[d_Np],Averages.Vel_y);
d_ScaLBL_Comm->RegularLayout(d_Map,&Velocity[2*d_Np],Averages.Vel_z);
PROFILE_STOP("Copy-State",1);
}
std::shared_ptr<double> cfq,cDen;
//if ( matches(type,AnalysisType::CreateRestart) ) {
if (timestep%d_restart_interval==0){
// Copy restart data to the CPU
cDen = std::shared_ptr<double>(new double[2*d_Np],DeleteArray<double>);
cfq = std::shared_ptr<double>(new double[19*d_Np],DeleteArray<double>);
ScaLBL_CopyToHost(cfq.get(),fq,19*d_Np*sizeof(double));
ScaLBL_CopyToHost(cDen.get(),Den,2*d_Np*sizeof(double));
}
PROFILE_STOP("Copy data to host",1);
// Spawn threads to do blob identification work
if ( matches(type,AnalysisType::IdentifyBlobs) ) {
phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase);
// Spawn threads to do blob identification work
if ( matches(type,AnalysisType::IdentifyBlobs) ) {
phase = std::shared_ptr<DoubleArray>(new DoubleArray(d_N[0],d_N[1],d_N[2]));
if (d_regular)
d_ScaLBL_Comm->RegularLayout(d_Map,Phi,*phase);
else
ScaLBL_CopyToHost(phase->data(),Phi,N*sizeof(double));
BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
BlobIDList new_list(new std::vector<BlobIDType>());
auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
work1->add_dependency(d_wait_blobID);
work2->add_dependency(d_tpool.add_work(work1));
d_wait_blobID = d_tpool.add_work(work2);
d_last_index = new_index;
d_last_ids = new_ids;
d_last_id_map = new_list;
}
BlobIDstruct new_index(new std::pair<int,IntArray>(0,IntArray()));
BlobIDstruct new_ids(new std::pair<int,IntArray>(0,IntArray()));
BlobIDList new_list(new std::vector<BlobIDType>());
auto work1 = new BlobIdentificationWorkItem1(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
auto work2 = new BlobIdentificationWorkItem2(timestep,d_N[0],d_N[1],d_N[2],d_rank_info,
phase,Averages.SDs,d_last_ids,new_index,new_ids,new_list,getComm());
work1->add_dependency(d_wait_blobID);
work2->add_dependency(d_tpool.add_work(work1));
d_wait_blobID = d_tpool.add_work(work2);
d_last_index = new_index;
d_last_ids = new_ids;
d_last_id_map = new_list;
}
// Spawn threads to do the analysis work
//if (timestep%d_restart_interval==0){
// if ( matches(type,AnalysisType::ComputeAverages) ) {
if ( timestep%d_analysis_interval == 0 ) {
auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta);
work->add_dependency(d_wait_blobID);
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying
d_wait_analysis = d_tpool.add_work(work);
}
// Spawn threads to do the analysis work
//if (timestep%d_restart_interval==0){
// if ( matches(type,AnalysisType::ComputeAverages) ) {
if ( timestep%d_analysis_interval == 0 ) {
auto work = new AnalysisWorkItem(type,timestep,Averages,d_last_index,d_last_id_map,d_beta);
work->add_dependency(d_wait_blobID);
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_vis); // Make sure we are done using analysis before modifying
d_wait_analysis = d_tpool.add_work(work);
}
// Spawn a thread to write the restart file
// if ( matches(type,AnalysisType::CreateRestart) ) {
if (timestep%d_restart_interval==0){
// Spawn a thread to write the restart file
// if ( matches(type,AnalysisType::CreateRestart) ) {
if (timestep%d_restart_interval==0){
if (d_rank==0) {
FILE *Rst = fopen("Restart.txt","w");
fprintf(Rst,"%i\n",timestep+4);
fclose(Rst);
}
// Write the restart file (using a seperate thread)
auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cPhi,cfq,d_Np);
work->add_dependency(d_wait_restart);
d_wait_restart = d_tpool.add_work(work);
}
if (d_rank==0) {
FILE *Rst = fopen("Restart.txt","w");
fprintf(Rst,"%i\n",timestep+4);
fclose(Rst);
}
// Write the restart file (using a seperate thread)
auto work = new WriteRestartWorkItem(d_restartFile.c_str(),cDen,cfq,d_Np);
work->add_dependency(d_wait_restart);
d_wait_restart = d_tpool.add_work(work);
}
// Save the results for visualization
// if ( matches(type,AnalysisType::CreateRestart) ) {
if (timestep%d_restart_interval==0){
// Write the vis files
auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
work->add_dependency(d_wait_blobID);
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_vis);
d_wait_vis = d_tpool.add_work(work);
}
PROFILE_STOP("run");
// Save the results for visualization
// if ( matches(type,AnalysisType::CreateRestart) ) {
if (timestep%d_restart_interval==0){
// Write the vis files
auto work = new WriteVisWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
work->add_dependency(d_wait_blobID);
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_vis);
d_wait_vis = d_tpool.add_work(work);
}
PROFILE_STOP("run");
}

2
common/ScaLBL.cpp
View File

@@ -1529,7 +1529,7 @@ void ScaLBL_Communicator::RegularLayout(IntArray map, const double *data, Double
int Nz = map.size(2);
// initialize the array
regdata.fill(-1.f);
regdata.fill(0.f);
double *TmpDat;
double value;

24
cpu/Color.cpp
View File

@@ -1252,7 +1252,7 @@ extern "C" void ScaLBL_SetSlice_z(double *Phi, double value, int Nx, int Ny, int
// double *ColorGrad, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
// double Fx, double Fy, double Fz, int start, int finish, int Np){
extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, double *Bq, double *Den, double *Phi,
double *Velocity, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
double *Vel, double rhoA, double rhoB, double tauA, double tauB, double alpha, double beta,
double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){
int ijk,nn,n;
@@ -1770,9 +1770,9 @@ extern "C" void ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *Aq, do
ux = jx / rho0;
uy = jy / rho0;
uz = jz / rho0;
//Velocity[n] = ux;
//Velocity[Np+n] = uy;
//Velocity[2*Np+n] = uz;
Vel[n] = ux;
Vel[Np+n] = uy;
Vel[2*Np+n] = uz;
// Instantiate mass transport distributions
// Stationary value - distribution 0
@@ -2418,9 +2418,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double *di
ux = jx / rho0;
uy = jy / rho0;
uz = jz / rho0;
//Velocity[n] = ux;
//Velocity[Np+n] = uy;
//Velocity[2*Np+n] = uz;
Vel[n] = ux;
Vel[Np+n] = uy;
Vel[2*Np+n] = uz;
// Instantiate mass transport distributions
// Stationary value - distribution 0
@@ -2777,11 +2777,15 @@ extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, doubl
n = Map[idx];
phi = Phi[n];
if (phi > 0.f){
nA = 1.0; nB = 0.f;
if (phi > 1.f){
nA = 1.0; nB = 0.f;
}
else if (phi < -1.f){
nB = 1.0; nA = 0.f;
}
else{
nB = 1.0; nA = 0.f;
nA=0.5*(phi+1.f);
nB=0.5*(1.f-phi);
}
Den[idx] = nA;
Den[Np+idx] = nB;

1
example/Plates/MixedWet.sh
View File

@@ -4,5 +4,4 @@ LBPM_DIR=../../tests
python ParallelPlates.py
mpirun -np 1 $LBPM_DIR/lbpm_serial_decomp input.db
mpirun -np 4 $LBPM_DIR/lbpm_segmented_pp input.db
mpirun -np 4 $LBPM_DIR/lbpm_color_simulator input.db

6
example/Sph1896/Simulation.sh Normal file
View File

@@ -0,0 +1,6 @@
#!/bin/bash
LBPM_INSTALL_DIR=../../bin
mpirun -np 8 $LBPM_INSTALL_DIR/GenerateSphereTest input.db
mpirun -np 8 $LBPM_INSTALL_DIR/lbpm_color_simulator input.db

56
example/Sph1896/input.db
View File

@@ -1,7 +1,7 @@
MRT {
tau = 1.0
F = 0, 0, 1e-4
timestepMax = 1000
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
@@ -9,39 +9,45 @@ MRT {
}
Color {
tau = 1.0;
alpha = 1e-2;
beta = 0.95;
phi_s = 0.8;
wp_saturation = 0.7
F = 0, 0, 0
Restart = false
pBC = 0
din = 1.0
dout = 1.0
timestepMax = 200
interval = 1000
tol = 1e-5;
das = 0.1
dbs = 0.9
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)
nproc = 2, 2, 2 // Number of processors (Npx,Npy,Npz)
n = 120, 120, 120 // Size of local domain (Nx,Ny,Nz)
nspheres = 1896 // Number of spheres
L = 1, 1, 1 // Length of domain (x,y,z)
BC = 0 // Boundary condition type
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
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 = 20000 // Frequency to write restart data
vis_interval = 20000 // Frequency to write visualization data
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"
}

32
gpu/Color.cu
View File

@@ -1762,9 +1762,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_Color(int *Map, double *dist, double *A
ux = jx / rho0;
uy = jy / rho0;
uz = jz / rho0;
//Velocity[n] = ux;
//Velocity[Np+n] = uy;
//Velocity[2*Np+n] = uz;
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
// Instantiate mass transport distributions
// Stationary value - distribution 0
@@ -2411,9 +2411,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_Color(int *neighborList, int *Map, double
ux = jx / rho0;
uy = jy / rho0;
uz = jz / rho0;
//Velocity[n] = ux;
//Velocity[Np+n] = uy;
//Velocity[2*Np+n] = uz;
Velocity[n] = ux;
Velocity[Np+n] = uy;
Velocity[2*Np+n] = uz;
// Instantiate mass transport distributions
// Stationary value - distribution 0
@@ -3897,17 +3897,21 @@ __global__ void dvc_ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, d
int S = Np/NBLOCKS/NTHREADS + 1;
for (int s=0; s<S; s++){
//........Get 1-D index for this thread....................
idx = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x;
idx = S*blockIdx.x*blockDim.x + s*blockDim.x + threadIdx.x + start;
if (idx<finish) {
n = Map[idx];
phi = Phi[n];
if (phi > 0.f){
nA = 1.0; nB = 0.f;
}
else{
nB = 1.0; nA = 0.f;
}
if (phi > 1.f){
nA = 1.0; nB = 0.f;
}
else if (phi < -1.f){
nB = 1.0; nA = 0.f;
}
else{
nA=0.5*(phi+1.f);
nB=0.5*(1.f-phi);
}
Den[idx] = nA;
Den[Np+idx] = nB;
@@ -4057,7 +4061,7 @@ extern "C" void ScaLBL_PhaseField_Init(int *Map, double *Phi, double *Den, doubl
dvc_ScaLBL_PhaseField_Init<<<NBLOCKS,NTHREADS >>>(Map, Phi, Den, Aq, Bq, start, finish, Np);
cudaError_t err = cudaGetLastError();
if (cudaSuccess != err){
printf("CUDA error in ScaLBL_D3Q19_ColorGrad: %s \n",cudaGetErrorString(err));
printf("CUDA error in ScaLBL_PhaseField_Init: %s \n",cudaGetErrorString(err));
}
}

14
gpu/D3Q19.cu
View File

@@ -30,6 +30,20 @@ __constant__ __device__ double mrt_V12=0.04166666666666666;
//__shared__ double Transform[722]=
// {};
#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600
#else
__device__ double atomicAdd(double* address, double val) {
unsigned long long int* address_as_ull = (unsigned long long int*)address;
unsigned long long int old = *address_as_ull, assumed;
do {
assumed = old;
old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val+__longlong_as_double(assumed)));
} while (assumed != old);
return __longlong_as_double(old);
}
#endif
using namespace cooperative_groups;
__device__ double reduce_sum(thread_group g, double *temp, double val)
{

14
gpu/dfh.cu
View File

@@ -5,6 +5,20 @@
#define NBLOCKS 1024
#define NTHREADS 256
#if !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600
#else
__device__ double atomicAdd(double* address, double val) {
unsigned long long int* address_as_ull = (unsigned long long int*)address;
unsigned long long int old = *address_as_ull, assumed;
do {
assumed = old;
old = atomicCAS(address_as_ull, assumed, __double_as_longlong(val+__longlong_as_double(assumed)));
} while (assumed != old);
return __longlong_as_double(old);
}
#endif
__global__ void dvc_ScaLBL_Gradient_Unpack(double weight, double Cqx, double Cqy, double Cqz,
int *list, int start, int count, double *recvbuf, double *phi, double *grad, int N){
//....................................................................................

79
models/ColorModel.cpp
View File

@@ -2,6 +2,7 @@
color lattice boltzmann model
*/
#include "models/ColorModel.h"
#include "analysis/distance.h"
ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, MPI_Comm COMM):
rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rhoA(0),rhoB(0),alpha(0),beta(0),
@@ -113,24 +114,42 @@ void ScaLBL_ColorModel::SetDomain(){
void ScaLBL_ColorModel::ReadInput(){
size_t readID;
//.......................................................................
if (rank == 0) printf("Read input media... \n");
//.......................................................................
Mask->ReadIDs();
for (int i=0; i<Nx*Ny*Nz; i++) id[i] = Mask->id[i]; // save what was read
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString);
// .......... READ THE INPUT FILE .......................................
//...........................................................................
if (rank == 0) cout << "Reading in signed distance function..." << endl;
//.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
ReadBinaryFile(LocalRankFilename, Averages->SDs.data(), N);
MPI_Barrier(comm);
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(Nx,Ny,Nz);
int count = 0;
// Solve for the position of the solid phase
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 the solid phase
if (Mask->id[n] > 0) id_solid(i,j,k) = 1;
else id_solid(i,j,k) = 0;
}
}
}
// Initialize the signed distance function
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
Averages->SDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
}
}
// MeanFilter(Averages->SDs);
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(Averages->SDs,id_solid,*Mask);
if (rank == 0) cout << "Domain set." << endl;
}
@@ -272,6 +291,9 @@ void ScaLBL_ColorModel::Create(){
********************************************************/
void ScaLBL_ColorModel::Initialize(){
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
/*
* This function initializes model
*/
@@ -290,13 +312,29 @@ void ScaLBL_ColorModel::Initialize(){
}
MPI_Bcast(&timestep,1,MPI_INT,0,comm);
// Read in the restart file to CPU buffers
double *cPhi = new double[Np];
double *cDist = new double[19*Np];
int *TmpMap;
TmpMap = new int[Np];
double *cPhi, *cDist;
cPhi = new double[N];
cDist = new double[19*Np];
ifstream File(LocalRestartFile,ios::binary);
double value;
int idx;
double value,va,vb;
ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int));
ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double));
for (int n=0; n<Np; n++){
File.read((char*) &va, sizeof(va));
File.read((char*) &vb, sizeof(vb));
value = (va-vb)/(va+vb);
idx = TmpMap[n];
if (!(idx < 0) && idx<N)
cPhi[idx] = value;
}
for (int n=0; n<Np; n++){
File.read((char*) &value, sizeof(value));
cPhi[n] = value;
// Read the distributions
for (int q=0; q<19; q++){
File.read((char*) &value, sizeof(value));
@@ -306,15 +344,12 @@ void ScaLBL_ColorModel::Initialize(){
File.close();
// Copy the restart data to the GPU
ScaLBL_CopyToDevice(fq,cDist,19*Np*sizeof(double));
ScaLBL_CopyToDevice(Phi,cPhi,Np*sizeof(double));
ScaLBL_CopyToDevice(Phi,cPhi,N*sizeof(double));
ScaLBL_DeviceBarrier();
delete [] cPhi;
delete [] cDist;
MPI_Barrier(comm);
}
if (rank==0) printf ("Initializing distributions \n");
ScaLBL_D3Q19_Init(fq, Np);
if (rank==0) printf ("Initializing phase field \n");
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);

112
models/MRTModel.cpp
View File

@@ -53,6 +53,10 @@ void ScaLBL_MRTModel::SetDomain(){
Nx+=2; Ny+=2; Nz += 2;
N = Nx*Ny*Nz;
Distance.resize(Nx,Ny,Nz);
Velocity_x.resize(Nx,Ny,Nz);
Velocity_y.resize(Nx,Ny,Nz);
Velocity_z.resize(Nx,Ny,Nz);
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);
@@ -138,6 +142,13 @@ void ScaLBL_MRTModel::Initialize(){
void ScaLBL_MRTModel::Run(){
double rlx_setA=1.0/tau;
double rlx_setB = 8.f*(2.f-rlx_setA)/(8.f-rlx_setA);
Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
memcpy(Morphology.SDn.data(), Distance.data(), Nx*Ny*Nz*sizeof(double));
if (rank==0) printf("time Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
@@ -160,6 +171,49 @@ void ScaLBL_MRTModel::Run(){
ScaLBL_D3Q19_AAeven_MRT(fq, 0, ScaLBL_Comm->LastExterior(), Np, rlx_setA, rlx_setB, Fx, Fy, Fz);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
//************************************************************************/
if (timestep%1000==0){
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); 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);
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 k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
if (Distance(i,j,k) > 0){
vax_loc += Velocity_x(i,j,k);
vay_loc += Velocity_y(i,j,k);
vaz_loc += Velocity_z(i,j,k);
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("%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu,
Morphology.V(),Morphology.A(),Morphology.J(),Morphology.X(),vax,vay,vaz);
}
}
//************************************************************************/
stoptime = MPI_Wtime();
@@ -178,9 +232,9 @@ void ScaLBL_MRTModel::Run(){
}
void ScaLBL_MRTModel::VelocityField(double *VELOCITY){
void ScaLBL_MRTModel::VelocityField(){
Minkowski Morphology(Mask);
/* Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
ScaLBL_D3Q19_Momentum(fq,Velocity, Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
@@ -223,6 +277,58 @@ void ScaLBL_MRTModel::VelocityField(double *VELOCITY){
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);
auto VxVar = std::make_shared<IO::Variable>();
auto VyVar = std::make_shared<IO::Variable>();
auto VzVar = std::make_shared<IO::Variable>();
auto SignDistVar = std::make_shared<IO::Variable>();
IO::initialize("","silo","false");
// Create the MeshDataStruct
visData.resize(1);
visData[0].meshName = "domain";
visData[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 );
SignDistVar->name = "SignDist";
SignDistVar->type = IO::VariableType::VolumeVariable;
SignDistVar->dim = 1;
SignDistVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(SignDistVar);
VxVar->name = "Velocity_x";
VxVar->type = IO::VariableType::VolumeVariable;
VxVar->dim = 1;
VxVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VxVar);
VyVar->name = "Velocity_y";
VyVar->type = IO::VariableType::VolumeVariable;
VyVar->dim = 1;
VyVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VyVar);
VzVar->name = "Velocity_z";
VzVar->type = IO::VariableType::VolumeVariable;
VzVar->dim = 1;
VzVar->data.resize(Dm->Nx-2,Dm->Ny-2,Dm->Nz-2);
visData[0].vars.push_back(VzVar);
Array<double>& SignData = visData[0].vars[0]->data;
Array<double>& VelxData = visData[0].vars[1]->data;
Array<double>& VelyData = visData[0].vars[2]->data;
Array<double>& VelzData = visData[0].vars[3]->data;
ASSERT(visData[0].vars[0]->name=="SignDist");
ASSERT(visData[0].vars[1]->name=="Velocity_x");
ASSERT(visData[0].vars[2]->name=="Velocity_y");
ASSERT(visData[0].vars[3]->name=="Velocity_z");
fillData.copy(Distance,SignData);
fillData.copy(Velocity_x,VelxData);
fillData.copy(Velocity_y,VelyData);
fillData.copy(Velocity_z,VelzData);
IO::writeData( timestep, visData, Dm->Comm );
}

7
models/MRTModel.h
View File

@@ -28,7 +28,7 @@ public:
void Create();
void Initialize();
void Run();
void VelocityField(double *Vz);
void VelocityField();
bool Restart,pBC;
int timestep,timestepMax;
@@ -58,9 +58,12 @@ public:
//Minkowski Morphology;
DoubleArray Velocity_x;
DoubleArray Velocity_y;
DoubleArray Velocity_z;
private:
MPI_Comm comm;
// filenames
char LocalRankString[8];
char LocalRankFilename[40];

60
sample_scripts/config_build_rhea Executable file
View File

@@ -0,0 +1,60 @@
# Clear all modules (except modules)
#for var in ${LOADEDMODULES//:/ }; do if [ " ${var///*/}" != " modules" ]; then module unload " ${var///*/}" > /dev/null 2>&1; fi; done
source $MODULESHOME/init/bash
module unload PE-intel
module load PE-gnu/6.2.0-2.0.1
#module load gcc/6.2.0
module load cmake3/3.5.2
module load openmpi netcdf hdf5
export MPICH_RDMA_ENABLED_CUDA=0
echo $GNU_VERSION
module list
# Remove CMake files from previous configures
rm -rf CMake*
# Configure
cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CMAKE_C_COMPILER:PATH=mpicc \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_CXX_COMPILER:PATH=mpicxx \
-D CMAKE_CXX_STD=11 \
-D USE_TIMER=false \
-D TIMER_DIRECTORY=${HOME}/timerutility/build/opt \
-D MPI_COMPILER:BOOL=TRUE \
-D MPIEXEC=aprun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_CUDA=0 \
-D CUDA_FLAGS="-arch sm_35" \
-D USE_HDF5=1 \
-D USE_SILO=1 \
-D HDF5_DIRECTORY=/sw/rhea/hdf5/1.8.11/rhel6.6_gnu4.8.2/ \
-D HDF5_LIB=/sw/rhea/hdf5/1.8.11/rhel6.6_gnu4.8.2/lib/libhdf5.so \
-D SILO_DIRECTORY=/lustre/atlas1/geo106/proj-shared/rhea/silo\
~/LBPM-WIA
# -D PREFIX=$MEMBERWORK/geo106/eos-LBPM-WIA \
#-D CUDA_HOST_COMPILER="/usr/bin/gcc" \
# Build the code
make install -j 8
# Fix permissions
#chmod -R g+w $PROJWORK/geo106/eos-LBPM-WIA
# Run the fast tests
# ctest -E WEEKLY
# Run the slow tests
# ctest -R WEEKLY -VV

18
sample_scripts/config_build_titan
View File

@@ -1,11 +1,13 @@
# Set the modules and enviornmental variables
source $MODULESHOME/init/bash
module load cudatoolkit
module load cmake
module load cmake3/3.9.0
export MPICH_RDMA_ENABLED_CUDA=1
#module swap cray-mpich2 cray-mpich2/5.6.3
module swap PrgEnv-pgi PrgEnv-gnu
module load cray-hdf5 silo
#module load PrgEnv-pgi
#module swap gcc gcc/4.8.2
module load cudatoolkit
#module load cray-hdf5 silo
# Remove CMake files from previous configures
rm -rf CMake*
@@ -17,20 +19,24 @@ cmake \
-D CMAKE_C_COMPILER:PATH=cc \
-D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_FLAGS="-fPIC" \
-D CXX_STD=11 \
-D USE_CUDA=1 \
-D CMAKE_CUDA_FLAGS="-arch sm_35" \
-D CMAKE_CUDA_HOST_COMPILER="/usr/bin/gcc" \
-D CMAKE_CUDA_HOST_COMPILER="/opt/gcc/6.3.0/bin/gcc" \
-D USE_MPI=1 \
-D MPI_COMPILER:BOOL=TRUE \
-D MPIEXEC=aprun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_SILO=1 \
-D SILO_DIRECTORY=/sw/xk6/silo/4.8/sles11.1_gnu4.5.3 \
-D HDF5_DIRECTORY=/opt/cray/hdf5/1.8.16/GNU/4.9 \
-D SILO_DIRECTORY=/ccs/proj/geo106/titan/TPLS/silo \
-D HDF5_DIRECTORY=/ccs/proj/geo106/titan/TPLS/hdf5 \
-D USE_TIMER=0 \
${HOME}/LBPM-WIA
# -D SILO_DIRECTORY=/sw/xk6/silo/4.8/sles11.1_gnu4.5.3 \
# -D HDF5_DIRECTORY=/opt/cray/hdf5/1.8.16/GNU/4.9 \
# Build the code

57
sample_scripts/config_build_titan_silo
View File

@@ -1,42 +1,33 @@
# Set the modules and enviornmental variables
source $MODULESHOME/init/bash
module load cudatoolkit
module load cmake
module load cmake3/3.9.0
export MPICH_RDMA_ENABLED_CUDA=1
#module swap cray-mpich2 cray-mpich2/5.6.3
module swap PrgEnv-pgi PrgEnv-gnu
#module load PrgEnv-pgi
#module swap gcc gcc/4.8.2
module load cudatoolkit
module load cray-hdf5 silo
# Remove CMake files from previous configures
rm -rf CMake*
# Configure
cmake \
-D CMAKE_C_COMPILER:PATH=cc \
-D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_FLAGS="" \
-D MPI_COMPILER:BOOL=TRUE \
-D MPIEXEC=aprun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CXX_STD=11 \
-D CUDA_FLAGS="-arch sm_35" \
-D USE_SILO=1 \
-D SILO_DIRECTORY=/sw/xk6/silo/4.8/sles11.1_gnu4.5.3 \
-D HDF5_DIRECTORY=/opt/cray/hdf5/1.8.16/GNU/4.9 \
-D CUDA_HOST_COMPILER="/usr/bin/gcc" \
-D USE_CUDA=1 \
-D USE_TIMER=0 \
${HOME}/LBPM-WIA
# Build the code
make install -j 8
# Run the fast tests
#ctest -E WEEKLY
# Run the slow tests
# ctest -R WEEKLY -VV
cmake \
-D CMAKE_BUILD_TYPE:STRING=Release \
-D CMAKE_C_COMPILER:PATH=cc \
-D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_STANDARD=11 \
-D USE_CUDA=1 \
-D CMAKE_CUDA_FLAGS="-arch sm_35" \
-D CMAKE_CUDA_HOST_COMPILER="/opt/gcc/6.3.0/bin/gcc" \
-D USE_MPI=1 \
-D MPI_COMPILER:BOOL=TRUE \
-D MPIEXEC=aprun \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_NETCDF=0 \
-D USE_SILO=1 \
-D SILO_DIRECTORY=/sw/xk6/silo/4.8/sles11.1_gnu4.5.3 \
-D HDF5_DIRECTORY=/sw/xk6/hdf5/1.8.7/cle4.1_gnu4.7.2 \
-D USE_DOXYGEN=0 \
-D EXTERNAL_LIBS="/sw/xk6/szip/2.1/sles11.1_gnu4.7.2/lib/libsz.a" \
-D USE_TIMER=0 \
${HOME}/master_LBPM-WIA/LBPM-WIA/

1
tests/CMakeLists.txt
View File

@@ -33,6 +33,7 @@ ADD_LBPM_EXECUTABLE( BlobIdentifyParallel )
ADD_LBPM_EXECUTABLE( convertIO )
ADD_LBPM_EXECUTABLE( DataAggregator )
#ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )
ADD_LBPM_EXECUTABLE( TestMinkowski )
CONFIGURE_FILE( ${CMAKE_CURRENT_SOURCE_DIR}/cylindertest ${CMAKE_CURRENT_BINARY_DIR}/cylindertest COPYONLY )

18
tests/GenerateSphereTest.cpp
View File

@@ -360,9 +360,9 @@ int main(int argc, char **argv)
nprocx = Dm->nprocx();
nprocy = Dm->nprocy();
nprocz = Dm->nprocz();
nspheres = domain_db->getScalar<int>( "nspheres");
nspheres = domain_db->getScalar<int>( "nspheres");
printf("Set domain \n");
//printf("Set domain \n");
int BoundaryCondition=1;
//Nz += 2;
//Nx = Ny = Nz; // Cubic domain
@@ -379,14 +379,6 @@ int main(int argc, char **argv)
}
Dm->CommInit();
if (rank==0) printf("Number of nodes per side = %i \n", Nx);
if (rank==0) printf("Total Number of nodes = %i \n", N);
if (rank==0) printf("********************************************************\n");
//.......................................................................
if (rank == 0) printf("Read input media... \n");
//.......................................................................
//.......................................................................
// Filenames used
char LocalRankString[8];
@@ -410,7 +402,7 @@ int main(int argc, char **argv)
//.......................................................................
// Read in sphere pack
if (rank==1) printf("nspheres =%i \n",nspheres);
//if (rank==1) printf("nspheres =%i \n",nspheres);
//.......................................................................
double *cx,*cy,*cz,*rad;
cx = new double[nspheres];
@@ -489,7 +481,7 @@ int main(int argc, char **argv)
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;
//.........................................................
//.......................................................................
/* //.......................................................................
sprintf(LocalRankString,"%05d",rank);
sprintf(LocalRankFilename,"%s%s","SignDist.",LocalRankString);
FILE *DIST = fopen(LocalRankFilename,"wb");
@@ -497,7 +489,7 @@ int main(int argc, char **argv)
fwrite(SignDist.data(),1,N*sizeof(double),DIST);
fclose(DIST);
//......................................................................
*/
//.......................................................................
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
FILE *IDFILE = fopen(LocalRankFilename,"wb");

222
tests/TestMinkowski.cpp Normal file
View File

@@ -0,0 +1,222 @@
// Sequential blob analysis
// Reads parallel simulation data and performs connectivity analysis
// and averaging on a blob-by-blob basis
// James E. McClure 2014
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <functional>
#include "common/Array.h"
#include "common/Domain.h"
#include "common/Communication.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
#include "IO/Mesh.h"
#include "IO/Writer.h"
#include "IO/netcdf.h"
#include "analysis/analysis.h"
#include "analysis/filters.h"
#include "analysis/distance.h"
#include "analysis/Minkowski.h"
#include "ProfilerApp.h"
int main(int argc, char **argv)
{
// Initialize MPI
int rank, nprocs;
MPI_Init(&argc,&argv);
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
{
Utilities::setErrorHandlers();
PROFILE_START("Main");
//std::vector<std::string> filenames;
if ( argc<2 ) {
if ( rank == 0 ){
printf("At least one filename must be specified\n");
}
return 1;
}
std::string filename = std::string(argv[1]);
if ( rank == 0 ){
printf("Input data file: %s\n",filename.c_str());
}
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
// Read domain parameters
auto Filename = domain_db->getScalar<std::string>( "Filename" );
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto SIZE = domain_db->getVector<int>( "N" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
auto nx = size[0];
auto ny = size[1];
auto nz = size[2];
auto nprocx = nproc[0];
auto nprocy = nproc[1];
auto nprocz = nproc[2];
auto Nx = SIZE[0];
auto Ny = SIZE[1];
auto Nz = SIZE[2];
int i,j,k,n;
char *SegData = NULL;
// Rank=0 reads the entire segmented data and distributes to worker processes
if (rank==0){
printf("Dimensions of segmented image: %i x %i x %i \n",Nx,Ny,Nz);
SegData = new char[Nx*Ny*Nz];
FILE *SEGDAT = fopen(Filename.c_str(),"rb");
if (SEGDAT==NULL) ERROR("Error reading segmented data");
size_t ReadSeg;
ReadSeg=fread(SegData,1,Nx*Ny*Nz,SEGDAT);
if (ReadSeg != size_t(Nx*Ny*Nz)) printf("lbpm_segmented_decomp: Error reading segmented data (rank=%i)\n",rank);
fclose(SEGDAT);
printf("Read segmented data from %s \n",Filename.c_str());
}
MPI_Barrier(comm);
// Get the rank info
int N = (nx+2)*(ny+2)*(nz+2);
std::shared_ptr<Domain> Dm (new Domain(domain_db,comm));
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
n = k*(nx+2)*(ny+2)+j*(nx+2)+i;
Dm->id[n] = 1;
}
}
}
Dm->CommInit();
int z_transition_size = 0;
int xStart = 0;
int yStart = 0;
int zStart = 0;
// Set up the sub-domains
if (rank==0){
printf("Distributing subdomain across %i processors \n",nprocs);
printf("Process grid: %i x %i x %i \n",Dm->nprocx(),Dm->nprocy(),Dm->nprocz());
printf("Subdomain size: %i \n",N);
//printf("Size of transition region: %i \n", z_transition_size);
char *tmp;
tmp = new char[N];
for (int kp=0; kp<nprocz; kp++){
for (int jp=0; jp<nprocy; jp++){
for (int ip=0; ip<nprocx; ip++){
// rank of the process that gets this subdomain
int rnk = kp*Dm->nprocx()*Dm->nprocy() + jp*Dm->nprocx() + ip;
// Pack and send the subdomain for rnk
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int x = xStart + ip*nx + i-1;
int y = yStart + jp*ny + j-1;
// int z = zStart + kp*nz + k-1;
int z = zStart + kp*nz + k-1 - z_transition_size;
if (x<xStart) x=xStart;
if (!(x<Nx)) x=Nx-1;
if (y<yStart) y=yStart;
if (!(y<Ny)) y=Ny-1;
if (z<zStart) z=zStart;
if (!(z<Nz)) z=Nz-1;
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
int nglobal = z*Nx*Ny+y*Nx+x;
tmp[nlocal] = SegData[nglobal];
}
}
}
if (rnk==0){
for (k=0;k<nz+2;k++){
for (j=0;j<ny+2;j++){
for (i=0;i<nx+2;i++){
int nlocal = k*(nx+2)*(ny+2) + j*(nx+2) + i;
Dm->id[nlocal] = tmp[nlocal];
}
}
}
}
else{
printf("Sending data to process %i \n", rnk);
MPI_Send(tmp,N,MPI_CHAR,rnk,15,comm);
}
}
}
}
}
else{
// Recieve the subdomain from rank = 0
printf("Ready to recieve data %i at process %i \n", N,rank);
MPI_Recv(Dm->id,N,MPI_CHAR,0,15,comm,MPI_STATUS_IGNORE);
}
MPI_Barrier(comm);
// Compute the Minkowski functionals
MPI_Barrier(comm);
std::shared_ptr<Minkowski> Averages(new Minkowski(Dm));
// Calculate the distance
// Initialize the domain and communication
nx+=2; ny+=2; nz+=2;
Array<char> id(nx,ny,nz);
//if (rank==0){
//printf("ID: %i, %i, %i \n",Dm->Nx, Dm->Ny, Dm->Nz);
// printf("ID: %i, %i, %i \n",id.size(0),id.size(1),id.size(2));
// printf("SDn: %i, %i, %i \n",Averages->SDn.size(0),Averages->SDn.size(1),Averages->SDn.size(2));
//}
// Solve for the position of the solid phase
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n = k*nx*ny+j*nx+i;
// Initialize the object
if (Dm->id[n] == ReadValues[0]) id(i,j,k) = 1;
else id(i,j,k) = 0;
}
}
}
for (k=0;k<nz;k++){
for (j=0;j<ny;j++){
for (i=0;i<nx;i++){
n=k*nx*ny+j*nx+i;
// Initialize distance to +/- 1
Averages->SDn(i,j,k) = 2.0*double(id(i,j,k))-1.0;
}
}
}
//MeanFilter(Averages->SDn);
//std::array<bool> bc(3)={1,1,1};
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(Averages->SDn,id,*Dm);
if (rank==0) printf("Computing Minkowski functionals \n");
Averages->Initialize();
Averages->UpdateMeshValues();
Averages->ComputeLocal();
Averages->Reduce();
Averages->PrintAll();
}
PROFILE_STOP("Main");
PROFILE_SAVE("Minkowski",true);
MPI_Barrier(comm);
MPI_Finalize();
return 0;
}

6
tests/TestPoiseuille.cpp
View File

@@ -58,8 +58,12 @@ int main(int argc, char **argv)
MRT.Initialize(); // initializing the model will set initial conditions for variables
MRT.Run();
double *Vz; Vz= new double [3*MRT.Np];
MRT.VelocityField(Vz);
int SIZE=MRT.Np*sizeof(double);
ScaLBL_D3Q19_Momentum(MRT.fq,MRT.Velocity, MRT.Np);
ScaLBL_DeviceBarrier(); MPI_Barrier(comm);
ScaLBL_CopyToHost(&Vz[0],&MRT.Velocity[0],3*SIZE);
if (rank == 0) printf("Force: %f,%f,%f \n",MRT.Fx,MRT.Fy,MRT.Fz);
double mu = MRT.mu;
int Nx = MRT.Nx;

85
tests/lbpm_morphdrain_pp.cpp
View File

@@ -53,14 +53,6 @@ int main(int argc, char **argv)
MPI_Comm_rank(comm,&rank);
MPI_Comm_size(comm,&nprocs);
//double Rcrit;
double SW=strtod(argv[1],NULL);
if (rank==0){
//printf("Critical radius %f (voxels)\n",Rcrit);
printf("Target saturation %f \n",SW);
}
// }
//.......................................................................
// Reading the domain information file
//.......................................................................
@@ -69,34 +61,27 @@ int main(int argc, char **argv)
int i,j,k,n;
int BC=0;
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
string filename;
if (argc > 1){
filename=argv[1];
}
MPI_Barrier(comm);
// Computational domain
MPI_Bcast(&nx,1,MPI_INT,0,comm);
MPI_Bcast(&ny,1,MPI_INT,0,comm);
MPI_Bcast(&nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//.................................................
MPI_Barrier(comm);
else ERROR("No input database provided\n");
// read the input database
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
// Read domain parameters
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
double SW = domain_db->getScalar<double>( "Sw" );
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
@@ -128,28 +113,8 @@ int main(int argc, char **argv)
}
}
}
/*// Exclude the maximum / minimum z rows
if (rank/(nprocx*nprocy)==0){
int k=0;
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n = k*nx*ny+j*nx+i;
Dm.id[n] = 0;
}
}
}
if (rank/(nprocx*nprocy)==nprocz-1){
int k=nz-1;
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
n = k*nx*ny+j*nx+i;
Dm.id[n] = 0;
}
}
}
*/
Dm.CommInit();
Dm.CommInit();
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
@@ -176,7 +141,7 @@ int main(int argc, char **argv)
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
Dm.CommInit();
int iproc = Dm.iproc();
int jproc = Dm.jproc();
@@ -479,12 +444,6 @@ int main(int argc, char **argv)
}
}
// if (rank==0){
// printf("Final saturation=%f\n",sw);
// printf("Final critical radius=%f\n",Rcrit);
//
// }
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);

3
tests/lbpm_morphopen_pp.cpp
View File

@@ -74,7 +74,6 @@ int main(int argc, char **argv)
filename=argv[1];
Rcrit_new=0.f;
//SW=strtod(argv[2],NULL);
if (rank==0) printf("Target saturation %f \n",SW);
}
else ERROR("No input database provided\n");
// read the input database
@@ -89,6 +88,8 @@ int main(int argc, char **argv)
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
SW = domain_db->getScalar<double>("Sw");
if (rank==0) printf("Target saturation %f \n",SW);
nx = size[0];
ny = size[1];
nz = size[2];

3
tests/lbpm_permeability_simulator.cpp
View File

@@ -51,8 +51,7 @@ int main(int argc, char **argv)
MRT.Create(); // creating the model will create data structure to match the pore structure and allocate variables
MRT.Initialize(); // initializing the model will set initial conditions for variables
MRT.Run();
double *Velocity; Velocity= new double [3*MRT.Np];
MRT.VelocityField(Velocity);
MRT.VelocityField();
}
// ****************************************************
MPI_Barrier(comm);

106
tests/lbpm_refine_pp.cpp
View File

@@ -31,37 +31,28 @@ int main(int argc, char **argv)
int i,j,k,n;
int BC=0;
if ( rank==0 ) {
printf("lbpm_refine_pp: Refining signed distance mesh (x2) \n");
string filename;
if (argc > 1){
filename=argv[1];
}
if (rank==0){
ifstream domain("Domain.in");
domain >> nprocx;
domain >> nprocy;
domain >> nprocz;
domain >> nx;
domain >> ny;
domain >> nz;
domain >> nspheres;
domain >> Lx;
domain >> Ly;
domain >> Lz;
else ERROR("No input database provided\n");
// read the input database
auto db = std::make_shared<Database>( filename );
auto domain_db = db->getDatabase( "Domain" );
}
MPI_Barrier(comm);
// Computational domain
MPI_Bcast(&nx,1,MPI_INT,0,comm);
MPI_Bcast(&ny,1,MPI_INT,0,comm);
MPI_Bcast(&nz,1,MPI_INT,0,comm);
MPI_Bcast(&nprocx,1,MPI_INT,0,comm);
MPI_Bcast(&nprocy,1,MPI_INT,0,comm);
MPI_Bcast(&nprocz,1,MPI_INT,0,comm);
MPI_Bcast(&nspheres,1,MPI_INT,0,comm);
MPI_Bcast(&Lx,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Ly,1,MPI_DOUBLE,0,comm);
MPI_Bcast(&Lz,1,MPI_DOUBLE,0,comm);
//.................................................
MPI_Barrier(comm);
// Read domain parameters
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
nx = size[0];
ny = size[1];
nz = size[2];
nprocx = nproc[0];
nprocy = nproc[1];
nprocz = nproc[2];
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
@@ -112,7 +103,15 @@ int main(int argc, char **argv)
ReadSignDist=fread(SignDist.data(),8,N,DIST);
if (ReadSignDist != size_t(N)) printf("lbpm_refine_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);
char *Labels;
Labels = new char[N];
sprintf(LocalRankFilename,"ID.%05i",rank);
FILE *LABELS = fopen(LocalRankFilename,"rb");
size_t ReadLabels;
ReadLabels=fread(Labels,1,N,LABELS);
if (ReadLabels != size_t(N)) printf("lbpm_refine_pp: Error reading ID (rank=%i)\n",rank);
fclose(LABELS);
if ( rank==0 ) printf("Set up Domain, read input distance \n");
@@ -130,6 +129,9 @@ int main(int argc, char **argv)
}
int ri,rj,rk,rn; //refined mesh indices
//char *RefineLabel;
//RefineLabel = new char [rnx*rny*rnz];
Array <char> RefineLabel(rnx,rny,rnz);
for (rk=1; rk<rnz-1; rk++){
for (rj=1; rj<rny-1; rj++){
for (ri=1; ri<rnx-1; ri++){
@@ -146,6 +148,7 @@ int main(int argc, char **argv)
pt.y=0.5*(rj-1)+1.f;
pt.z=0.5*(rk-1)+1.f;
RefinedSignDist(ri,rj,rk) = LocalApprox.eval(pt);
RefineLabel(ri,rj,rk) = Labels[k*nx*ny+j*nx+i];
}
}
}
@@ -175,6 +178,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i,j,k);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i,j,k);
}
}
}
@@ -183,16 +188,17 @@ int main(int argc, char **argv)
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
for (int k=0; k<nz; k++){
/* for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
if (BlockDist(i,j,k) > 0.f)
id[k*nx*ny + j*nx + i]=2;
else
id[k*nx*ny + j*nx + i]=0;
id[k*nx*ny + j*nx + i]= 0;
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -203,6 +209,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j,k);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j,k);
}
}
}
@@ -211,7 +219,7 @@ int main(int argc, char **argv)
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
for (int k=0; k<nz; k++){
/* for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
if (BlockDist(i,j,k) > 0.f)
@@ -221,6 +229,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -232,6 +241,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j+ny-2,k);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j+ny-2,k);
}
}
}
@@ -240,7 +251,7 @@ int main(int argc, char **argv)
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
for (int k=0; k<nz; k++){
/* for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
if (BlockDist(i,j,k) > 0.f)
@@ -250,6 +261,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -260,6 +272,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i,j+ny-2,k);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i,j+ny-2,k);
}
}
}
@@ -267,7 +281,7 @@ int main(int argc, char **argv)
REFINEDIST = fopen(LocalRankFilename,"wb");
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
/*
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@@ -278,6 +292,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -288,6 +303,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i,j,k+nz-2);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i,j,k+nz-2);
}
}
}
@@ -295,7 +312,7 @@ int main(int argc, char **argv)
REFINEDIST = fopen(LocalRankFilename,"wb");
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
/*
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@@ -306,6 +323,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -316,6 +334,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j,k+nz-2);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j,k+nz-2);
}
}
}
@@ -324,7 +344,7 @@ int main(int argc, char **argv)
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
for (int k=0; k<nz; k++){
/* for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
if (BlockDist(i,j,k) > 0.f)
@@ -334,6 +354,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -344,6 +365,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i,j+ny-2,k+nz-2);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i,j+ny-2,k+nz-2);
}
}
}
@@ -351,7 +374,7 @@ int main(int argc, char **argv)
REFINEDIST = fopen(LocalRankFilename,"wb");
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
/*
for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
@@ -362,6 +385,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);
@@ -372,6 +396,8 @@ int main(int argc, char **argv)
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
BlockDist(i,j,k) = RefinedSignDist(i+nx-2,j+ny-2,k+nz-2);
if (BlockDist(i,j,k) > 0) id[k*nx*ny + j*nx + i]=2;
else id[k*nx*ny + j*nx + i] = RefineLabel(i+nx-2,j+ny-2,k+nz-2);
}
}
}
@@ -380,7 +406,8 @@ int main(int argc, char **argv)
REFINEDIST = fopen(LocalRankFilename,"wb");
fwrite(BlockDist.data(),8,nx*ny*nz,REFINEDIST);
fclose(REFINEDIST);
for (int k=0; k<nz; k++){
/* for (int k=0; k<nz; k++){
for (int j=0; j<ny; j++){
for (int i=0; i<nx; i++){
if (BlockDist(i,j,k) > 0.f)
@@ -390,6 +417,7 @@ int main(int argc, char **argv)
}
}
}
*/
sprintf(LocalRankFilename,"RefineID.%05i",writerank);
WRITEID = fopen(LocalRankFilename,"wb");
fwrite(id,1,nx*ny*nz,WRITEID);