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

Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA

Browse Source
This commit is contained in:
James E McClure
2018-11-25 20:55:58 -05:00
parent 8e636b808c 414abf0a55
commit c7566c3d4f
14 changed files with 1305 additions and 953 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
IO/Writer.h
View File

@@ -39,7 +39,7 @@ namespace IO {
* silo - Silo
* @param[in] append Append any existing data (default is false)
*/
void initialize( const std::string& path="", const std::string& format="new", bool append=false );
void initialize( const std::string& path="", const std::string& format="silo", bool append=false );
/*!

33
analysis/TwoPhase.cpp
View File

@@ -174,7 +174,7 @@ TwoPhase::TwoPhase(std::shared_ptr <Domain> dm):
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"time dEs "); // Timestep, Change in Surface Energy
fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn "); // Timestep, Change in Surface Energy
fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages
fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages
fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz ");
@@ -205,7 +205,7 @@ TwoPhase::TwoPhase(std::shared_ptr <Domain> dm):
sprintf(LocalRankFilename,"%s%s","timelog.tcat.",LocalRankString);
TIMELOG = fopen(LocalRankFilename,"a+");
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"time "); // Timestep, Change in Surface Energy
fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn ");; // Timestep,
fprintf(TIMELOG,"sw pw pn awn ans aws Jwn Kwn lwns cwns KNwns KGwns "); // Scalar averages
fprintf(TIMELOG,"vawx vawy vawz vanx vany vanz "); // Velocity averages
fprintf(TIMELOG,"vawnx vawny vawnz vawnsx vawnsy vawnsz ");
@@ -338,6 +338,19 @@ void TwoPhase::Initialize()
wwndnw = 0.0; wwnsdnwn = 0.0; Jwnwwndnw=0.0;
}
void TwoPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha)
{
Fx = force_x;
Fy = force_y;
Fz = force_z;
rho_n = rhoA;
rho_w = rhoB;
nu_n = (tauA-0.5)/3.f;
nu_w = (tauB-0.5)/3.f;
gamma_wn = 5.796*alpha;
}
/*
void TwoPhase::SetupCubes(Domain &Dm)
{
@@ -589,15 +602,15 @@ void TwoPhase::ComputeLocal()
n = k*Nx*Ny+j*Nx+i;
if (!(Dm->id[n] > 0)){
// Solid phase
phase_label(i,j,k) = 0;
phase_label(i,j,k) = 1;
}
else if (SDn(i,j,k) < 0.0){
// wetting phase
phase_label(i,j,k) = 1;
phase_label(i,j,k) = 0;
}
else {
// non-wetting phase
phase_label(i,j,k) = 0;
phase_label(i,j,k) = 1;
}
phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0;
}
@@ -613,15 +626,15 @@ void TwoPhase::ComputeLocal()
n = k*Nx*Ny+j*Nx+i;
if (!(Dm->id[n] > 0)){
// Solid phase
phase_label(i,j,k) = 0;
phase_label(i,j,k) = 1;
}
else if (SDn(i,j,k) < 0.0){
// wetting phase
phase_label(i,j,k) = 0;
phase_label(i,j,k) = 1;
}
else {
// non-wetting phase
phase_label(i,j,k) = 1;
phase_label(i,j,k) = 0;
}
phase_distance(i,j,k) =2.0*double(phase_label(i,j,k))-1.0;
}
@@ -1192,7 +1205,7 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT)
void TwoPhase::PrintAll(int timestep)
{
if (Dm->rank()==0){
fprintf(TIMELOG,"%i %.5g ",timestep,dEs); // change in surface energy
fprintf(TIMELOG,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw_global,pan_global); // saturation and pressure
fprintf(TIMELOG,"%.5g %.5g %.5g ",awn_global,ans_global,aws_global); // interfacial areas
fprintf(TIMELOG,"%.5g %.5g ",Jwn_global, Kwn_global); // curvature of wn interface
@@ -1219,7 +1232,7 @@ void TwoPhase::PrintAll(int timestep)
else{
sat_w = 1.0 - nwp_volume/(nwp_volume+wp_volume);
fprintf(TIMELOG,"%i ",timestep); // change in surface energy
fprintf(TIMELOG,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(TIMELOG,"%.5g %.5g %.5g ",sat_w,paw,pan); // saturation and pressure
fprintf(TIMELOG,"%.5g %.5g %.5g ",awn,ans,aws); // interfacial areas
fprintf(TIMELOG,"%.5g %.5g ",Jwn, Kwn); // curvature of wn interface

12
analysis/TwoPhase.h
View File

@@ -104,6 +104,11 @@ public:
double euler,Kn,Jn,An;
double euler_global,Kn_global,Jn_global,An_global;
double rho_n, rho_w;
double nu_n, nu_w;
double gamma_wn;
double Fx, Fy, Fz;
double Jwn,Jwn_global; // average mean curavture - wn interface
double Kwn,Kwn_global; // average Gaussian curavture - wn interface
double KNwns,KNwns_global; // wns common curve normal curavture
@@ -184,6 +189,13 @@ public:
int GetCubeLabel(int i, int j, int k, IntArray &BlobLabel);
void SortBlobs();
void PrintComponents(int timestep);
void SetParams(double rhoA, double rhoB, double tauA, double tauB, double force_x, double force_y, double force_z, double alpha);
double Volume_w(){
return wp_volume_global;
}
double Volume_n(){
return nwp_volume_global;
}
};

16
analysis/runAnalysis.cpp
View File

@@ -38,8 +38,8 @@ AnalysisType& operator |=(AnalysisType &lhs, AnalysisType rhs)
}
bool matches( AnalysisType x, AnalysisType y )
{
return static_cast<std::underlying_type<AnalysisType>::type>(x) &
static_cast<std::underlying_type<AnalysisType>::type>(y) != 0;
return ( static_cast<std::underlying_type<AnalysisType>::type>(x) &
static_cast<std::underlying_type<AnalysisType>::type>(y) ) != 0;
}
@@ -55,7 +55,7 @@ class WriteRestartWorkItem: public ThreadPool::WorkItemRet<void>
{
public:
WriteRestartWorkItem( const char* filename_, std::shared_ptr<double> cDen_, std::shared_ptr<double> cfq_, int N_ ):
filename(filename_), cDen(cDen_), cfq(cfq_), N(N_) {}
filename(filename_), cfq(cfq_), cDen(cDen_), N(N_) {}
virtual void run() {
PROFILE_START("Save Checkpoint",1);
double value;
@@ -318,10 +318,16 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
d_regular ( Regular),
d_rank_info( rank_info ),
d_Map( Map ),
d_ScaLBL_Comm( ScaLBL_Comm),
d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1)
d_fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1),
d_ScaLBL_Comm( ScaLBL_Comm)
{
// Ids of work items to use for dependencies
ThreadPool::thread_id_t d_wait_blobID;
ThreadPool::thread_id_t d_wait_analysis;
ThreadPool::thread_id_t d_wait_vis;
ThreadPool::thread_id_t d_wait_restart;
char rankString[20];
sprintf(rankString,"%05d",Dm->rank());
d_N[0] = Dm->Nx;

2
analysis/uCT.cpp
View File

@@ -163,7 +163,7 @@ void solve( const Array<float>& VOL, Array<float>& Mean, Array<char>& ID,
// Compute the median filter on the sparse array
Med3D( VOL, Mean );
fillFloat.fill( Mean );
segment( Mean, ID, 0.01 );
segment( Mean, ID, threshold );
// Compute the distance using the segmented volume
CalcDist( Dist, ID, Dm );
fillFloat.fill(Dist);

8
common/Array.hpp
View File

@@ -1032,8 +1032,8 @@ Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::coarsen(
throw std::invalid_argument( "Array must be multiple of filter size" );
}
Array<TYPE, FUN, Allocator> y( S2 );
if ( d_size.ndim() <= 3 )
throw std::logic_error( "Function programmed for more than 3 dimensions" );
if ( d_size.ndim() > 3 )
throw std::logic_error( "Function not programmed for more than 3 dimensions" );
const auto &Nh = filter.d_size;
for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) {
for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) {
@@ -1067,8 +1067,8 @@ Array<TYPE, FUN, Allocator> Array<TYPE, FUN, Allocator>::coarsen( const std::vec
}
Array<TYPE, FUN, Allocator> tmp( ratio );
Array<TYPE, FUN, Allocator> y( S2 );
if ( d_size.ndim() <= 3 )
throw std::logic_error( "Function programmed for more than 3 dimensions" );
if ( d_size.ndim() > 3 )
throw std::logic_error( "Function not programmed for more than 3 dimensions" );
for ( size_t k1 = 0; k1 < y.d_size[2]; k1++ ) {
for ( size_t j1 = 0; j1 < y.d_size[1]; j1++ ) {
for ( size_t i1 = 0; i1 < y.d_size[0]; i1++ ) {

292
models/ColorModel.cpp
View File

@@ -95,8 +95,6 @@ void ScaLBL_ColorModel::ReadParams(string filename){
outletA=0.f;
outletB=1.f;
if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
// Read domain parameters
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
@@ -112,6 +110,8 @@ void ScaLBL_ColorModel::ReadParams(string filename){
nprocy = nproc[1];
nprocz = nproc[2];
if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details)
}
void ScaLBL_ColorModel::SetDomain(){
Dm = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // full domain for analysis
@@ -166,6 +166,8 @@ void ScaLBL_ColorModel::ReadInput(){
CalcDist(Averages->SDs,id_solid,*Mask);
if (rank == 0) cout << "Domain set." << endl;
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha);
}
void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
@@ -410,12 +412,62 @@ void ScaLBL_ColorModel::Run(){
int nprocs=nprocx*nprocy*nprocz;
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
bool SET_CAPILLARY_NUMBER = false;
bool MORPH_ADAPT = false;
bool USE_MORPH = false;
int morph_interval;
double morph_delta;
int morph_timesteps = 0;
int ramp_timesteps = 50000;
double capillary_number;
double tolerance = 1.f;
double Ca_previous = 0.f;
int target_saturation_index=0;
std::vector<double> target_saturation;
double TARGET_SATURATION = 0.f;
if (color_db->keyExists( "target_saturation" )){
target_saturation = color_db->getVector<double>( "target_saturation" );
TARGET_SATURATION = target_saturation[0];
}
if (color_db->keyExists( "capillary_number" )){
capillary_number = color_db->getScalar<double>( "capillary_number" );
SET_CAPILLARY_NUMBER=true;
}
else{
capillary_number=0;
}
if (BoundaryCondition != 0 && SET_CAPILLARY_NUMBER==true){
if (rank == 0) printf("WARINING: capillary number target only supported for BC = 0 \n");
SET_CAPILLARY_NUMBER=false;
}
if (analysis_db->keyExists( "morph_delta" )){
morph_delta = analysis_db->getScalar<double>( "morph_delta" );
}
else{
morph_delta=0.5;
}
if (analysis_db->keyExists( "morph_interval" )){
morph_interval = analysis_db->getScalar<int>( "morph_interval" );
USE_MORPH = true;
}
else{
morph_interval=1000000;
USE_MORPH = false;
}
if (analysis_db->keyExists( "tolerance" )){
tolerance = analysis_db->getScalar<double>( "tolerance" );
}
else{
tolerance = 0.02;
}
int analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
if (rank==0){
printf("********************************************************\n");
printf("No. of timesteps: %i \n", timestepMax);
fflush(stdout);
}
//.......create and start timer............
double starttime,stoptime,cputime;
ScaLBL_DeviceBarrier();
@@ -511,6 +563,113 @@ void ScaLBL_ColorModel::Run(){
// Run the analysis
analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
// allow initial ramp-up to get closer to steady state
if (timestep > ramp_timesteps && timestep%analysis_interval == analysis_interval-20 && USE_MORPH){
if ( morph_timesteps > morph_interval ){
double volB = Averages->Volume_w();
double volA = Averages->Volume_n();
double vA_x = Averages->van_global(0);
double vA_y = Averages->van_global(1);
double vA_z = Averages->van_global(2);
double vB_x = Averages->vaw_global(0);
double vB_y = Averages->vaw_global(1);
double vB_z = Averages->vaw_global(2);
double muA = rhoA*(tauA-0.5)/3.f;
double muB = rhoB*(tauB-0.5)/3.f;
double flow_rate_A = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z);
double flow_rate_B = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z);
double current_saturation = volB/(volA+volB);
double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double(Nx*Ny*Nz*nprocs));
double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz);
//double krA = muA*volA*flow_rate_A/force_magnitude/double(Nx*Ny*Nz*nprocs);
//double krB = muB*volB*flow_rate_B/force_magnitude/double(Nx*Ny*Nz*nprocs);
if (fabs((Ca - Ca_previous)/Ca) < tolerance ){
MORPH_ADAPT = true;
if (rank==0){
printf("** WRITE STEADY POINT *** ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
volA /= double(Nx*Ny*Nz*nprocs);
volB /= double(Nx*Ny*Nz*nprocs);
FILE * kr_log_file = fopen("relperm.csv","a");
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g ",timestep-analysis_interval+20,muA,muB,5.796*alpha,Fx,Fy,Fz);
fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z);
fclose(kr_log_file);
printf(" Measured capillary number %f \n ",Ca);
}
if (SET_CAPILLARY_NUMBER ){
Fx *= capillary_number / Ca;
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
if (force_magnitude > 1e-3){
Fx *= 1e-3/force_magnitude; // impose ceiling for stability
Fy *= 1e-3/force_magnitude;
Fz *= 1e-3/force_magnitude;
}
if (force_magnitude < 1e-6){
Fx *= 1e-6/force_magnitude; // impose floor
Fy *= 1e-6/force_magnitude;
Fz *= 1e-6/force_magnitude;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha);
}
if (morph_delta > 0.f){
// wetting phase saturation will decrease
while (current_saturation < TARGET_SATURATION && target_saturation_index < target_saturation.size() ){
TARGET_SATURATION = target_saturation[target_saturation_index++];
}
}
else{
// wetting phase saturation will increase
while (current_saturation > TARGET_SATURATION && target_saturation_index < target_saturation.size() ){
TARGET_SATURATION = target_saturation[target_saturation_index++];
if (rank==0) printf(" Set target saturation as %f (currently %f)\n",TARGET_SATURATION,current_saturation);
}
}
}
else{
if (rank==0){
printf("** Continue to simulate steady *** \n ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
}
morph_timesteps=0;
}
Ca_previous = Ca;
}
if (MORPH_ADAPT ){
if (rank==0) printf("***Morphological step with target saturation %f ***\n",TARGET_SATURATION);
double volB = Averages->Volume_w();
double volA = Averages->Volume_n();
double delta_volume = MorphInit(beta,morph_delta);
double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A
// update the volume
volA += delta_volume;
volB -= delta_volume;
if ((delta_volume_target - delta_volume) / delta_volume > 0.f){
morph_delta *= 1.01*min((delta_volume_target - delta_volume) / delta_volume, 2.0);
if (morph_delta > 1.f) morph_delta = 1.f;
if (morph_delta < -1.f) morph_delta = -1.f;
if (fabs(morph_delta) < 0.05 ) morph_delta = 0.05*(morph_delta)/fabs(morph_delta); // set minimum
if (rank==0) printf(" Adjust morph delta: %f \n", morph_delta);
}
//MORPH_ADAPT = false;
if (volB/(volA + volB) > TARGET_SATURATION){
MORPH_ADAPT = false;
TARGET_SATURATION = target_saturation[target_saturation_index++];
}
MPI_Barrier(comm);
morph_timesteps = 0;
}
morph_timesteps += analysis_interval;
}
}
analysis.finish();
PROFILE_STOP("Loop");
@@ -535,6 +694,133 @@ void ScaLBL_ColorModel::Run(){
// ************************************************************************
}
double ScaLBL_ColorModel::MorphInit(const double beta, const double morph_delta){
const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz);
double vF = 0.f;
double vS = 0.f;
DoubleArray phase(Nx,Ny,Nz);
IntArray phase_label(Nx,Ny,Nz);;
DoubleArray phase_distance(Nx,Ny,Nz);
Array<char> phase_id(Nx,Ny,Nz);
// Basic algorithm to
// 1. Copy phase field to CPU
ScaLBL_CopyToHost(phase.data(), Phi, N*sizeof(double));
double count,count_global,volume_initial,volume_final;
count = 0.f;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm);
volume_initial = count_global;
// 2. Identify connected components of phase field -> phase_label
BlobIDstruct new_index;
ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,Averages->SDs,vF,vS,phase_label,comm);
MPI_Barrier(comm);
// only operate on component "0"
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
int label = phase_label(i,j,k);
if (label == 0 ) phase_id(i,j,k) = 0;
else phase_id(i,j,k) = 1;
}
}
}
// 3. Generate a distance map to the largest object -> phase_distance
CalcDist(phase_distance,phase_id,*Dm);
double temp,value;
double factor=0.5/beta;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase_distance(i,j,k) < 3.f ){
value = phase(i,j,k);
if (value > 1.f) value=1.f;
if (value < -1.f) value=-1.f;
// temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm.
temp = -factor*log((1.0+value)/(1.0-value));
/// use this approximation close to the object
if (fabs(value) < 0.8 && Averages->SDs(i,j,k) > 1.f ){
phase_distance(i,j,k) = temp;
}
}
}
}
}
// 4. Apply erosion / dilation operation to phase_distance
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
double walldist=Averages->SDs(i,j,k);
double wallweight = 1.f / (1+exp(-5.f*(walldist-1.f)));
phase_distance(i,j,k) -= wallweight*morph_delta;
}
}
}
// 5. Update phase indicator field based on new distnace
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;
double d = phase_distance(i,j,k);
if (Averages->SDs(i,j,k) > 0.f){
// only update phase field in immediate proximity of largest component
if (d < 3.f){
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
}
}
count = 0.f;
for (int k=0; k<Nz; k++){
for (int j=0; j<Ny; j++){
for (int i=0; i<Nx; i++){
if (phase(i,j,k) > 0.f && Averages->SDs(i,j,k) > 0.f) count+=1.f;
}
}
}
MPI_Allreduce(&count,&count_global,1,MPI_DOUBLE,MPI_SUM,comm);
volume_final=count_global;
double delta_volume = (volume_final-volume_initial);
if (rank == 0) printf("MorphInit: change fluid volume fraction by %f \n", delta_volume/volume_initial);
// 6. copy back to the device
//if (rank==0) printf("MorphInit: copy data back to device\n");
ScaLBL_CopyToDevice(Phi,phase.data(),N*sizeof(double));
// 7. Re-initialize phase field and density
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);
if (BoundaryCondition >0 ){
if (Dm->kproc()==0){
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1);
ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2);
}
if (Dm->kproc() == nprocz-1){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2);
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
return delta_volume;
}
void ScaLBL_ColorModel::WriteDebug(){
// Copy back final phase indicator field and convert to regular layout
DoubleArray PhaseField(Nx,Ny,Nz);

1
models/ColorModel.h
View File

@@ -92,6 +92,7 @@ private:
//int rank,nprocs;
void LoadParams(std::shared_ptr<Database> db0);
void AssignComponentLabels(double *phase);
double MorphInit(const double beta, const double morph_delta);
};

22
models/MRTModel.cpp
View File

@@ -161,7 +161,11 @@ void ScaLBL_MRTModel::Run(){
Minkowski Morphology(Mask);
int SIZE=Np*sizeof(double);
if (rank==0) printf("time Fx Fy Fz mu Vs As Js Xs vx vy vz\n");
if (rank==0){
FILE * log_file = fopen("Permeability.csv","a");
fprintf(log_file,"time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n");
fclose(log_file);
}
//.......create and start timer............
double starttime,stoptime,cputime;
@@ -219,13 +223,19 @@ void ScaLBL_MRTModel::Run(){
vay /= count;
vaz /= count;
if (rank==0) printf("Computing Minkowski functionals \n");
//if (rank==0) printf("Computing Minkowski functionals \n");
Morphology.ComputeScalar(Distance,0.f);
Morphology.PrintAll();
//Morphology.PrintAll();
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);
if (rank==0) {
double h = Lz/double(Nz);
double absperm = h*h*mu*sqrt(vax*vax+vay*vay+vaz*vaz)/sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
printf(" %f\n",absperm);
FILE * log_file = fopen("Permeability.csv","a");
fprintf(log_file,"%i %.8g %.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, absperm);
fclose(log_file);
}
}
}
//************************************************************************/

15
sample_scripts/config_build_eos
View File

@@ -7,11 +7,13 @@ module unload PrgEnv-intel
module load PrgEnv-gnu/6.0.4
module unload gcc cmake
module load gcc/6.3.0
module load cray-hdf5-parallel/1.10.0.1
module load cray-hdf5-parallel/1.10.2.0
module load cray-netcdf-hdf5parallel
module load mercurial git
module load cmake3/3.6.1
export LD_LIBRARY_PATH=/usr/lib64:/lustre/atlas/proj-shared/geo106/eos/netcdf/lib:/lustre/atlas/proj-shared/geo106/eos/zlib/lib:/lustre/atlas/proj-shared/geo106/eos/hdf5/lib:$LD_LIBRARY_PATH
export MPICH_RDMA_ENABLED_CUDA=0
echo $GNU_VERSION
@@ -27,7 +29,8 @@ cmake \
-D CMAKE_C_COMPILER:PATH=cc \
-D CMAKE_CXX_COMPILER:PATH=CC \
-D CMAKE_CXX_COMPILER:PATH=CC \
-D CXX_STD=11 \
-D CMAKE_CXX_FLAGS="-fPIC" \
-D CMAKE_CXX_STANDARD=14 \
-D USE_TIMER=false \
-D TIMER_DIRECTORY=${HOME}/timerutility/build/opt \
-D MPI_COMPILER:BOOL=TRUE \
@@ -35,10 +38,12 @@ cmake \
-D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
-D USE_CUDA=0 \
-D CUDA_FLAGS="-arch sm_35" \
-DBUILD_SHARED_LIBS=OFF \
-D USE_SILO=1 \
-D SILO_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/silo \
-D HDF5_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/hdf5 \
-D HDF5_LIB=/lustre/atlas/proj-shared/geo106/eos/hdf5/lib/libhdf5.a \
-D USE_NETCDF=1 \
-D NETCDF_DIRECTORY=$NETCDF_DIR \
-D HDF5_DIRECTORY=$HDF5_DIR \
-D NETCDF_DIRECTORY=/lustre/atlas/proj-shared/geo106/eos/netcdf \
-D CMAKE_SKIP_RPATH=true \
~/LBPM-WIA

2
tests/lbpm_color_simulator.cpp
View File

@@ -54,7 +54,7 @@ int main(int argc, char **argv)
ColorModel.Create(); // creating the model will create data structure to match the pore structure and allocate variables
ColorModel.Initialize(); // initializing the model will set initial conditions for variables
ColorModel.Run();
ColorModel.WriteDebug();
//ColorModel.WriteDebug();
PROFILE_STOP("Main");
PROFILE_SAVE("lbpm_color_simulator",1);

54
tests/lbpm_morphdrain_pp.cpp
View File

@@ -11,7 +11,7 @@
#include <sstream>
#include "common/Array.h"
#include "common/Domain.h"
#include "analysis/distance.h"
//*************************************************************************
// Morpohological drainage pre-processor
@@ -96,7 +96,7 @@ int main(int argc, char **argv)
char LocalRankFilename[40];
int BoundaryCondition=1;
Domain Dm(nx,ny,nz,rank,nprocx,nprocy,nprocz,Lx,Ly,Lz,BoundaryCondition);
Domain Dm(domain_db,comm);
nx+=2; ny+=2; nz+=2;
int N = nx*ny*nz;
@@ -116,32 +116,60 @@ int main(int argc, char **argv)
Dm.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_morphdrain_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
int xdim,ydim,zdim;
xdim=Dm.Nx-2;
ydim=Dm.Ny-2;
zdim=Dm.Nz-2;
fillHalo<double> fillData(Dm.Comm, Dm.rank_info,{xdim,ydim,zdim},{1,1,1},0,1);
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(nx,ny,nz);
DoubleArray SignDist(nx,ny,nz);
// Read the signed distance from file
// 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 (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
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
}
}
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(SignDist,id_solid,Dm);
/* // Read the signed distance from file
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"rb");
size_t ReadSignDist;
ReadSignDist=fread(SignDist.data(),8,N,DIST);
if (ReadSignDist != size_t(N)) printf("lbpm_morphdrain_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);
*/
fillData.fill(SignDist);
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_segmented_pp: Error reading ID (rank=%i) \n",rank);
fclose(IDFILE);
Dm.CommInit();
int iproc = Dm.iproc();
int jproc = Dm.jproc();

65
tests/lbpm_morphopen_pp.cpp
View File

@@ -11,6 +11,7 @@
#include <sstream>
#include "common/Array.h"
#include "common/Domain.h"
#include "analysis/distance.h"
//*************************************************************************
@@ -106,16 +107,45 @@ int main(int argc, char **argv)
char *id;
id = new 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);
nx+=2; ny+=2; nz+=2;
// Generate the signed distance map
// Initialize the domain and communication
Array<char> id_solid(nx,ny,nz);
DoubleArray SignDist(nx,ny,nz);
// Read the signed distance from file
sprintf(LocalRankFilename,"SignDist.%05i",rank);
FILE *DIST = fopen(LocalRankFilename,"rb");
size_t ReadSignDist;
ReadSignDist=fread(SignDist.data(),8,N,DIST);
if (ReadSignDist != size_t(N)) printf("lbpm_morphopen_pp: Error reading signed distance function (rank=%i)\n",rank);
fclose(DIST);
// 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 (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
SignDist(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0;
}
}
}
if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
CalcDist(SignDist,id_solid,*Dm);
MPI_Barrier(comm);
double count,countGlobal,totalGlobal;
@@ -135,7 +165,9 @@ 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;
if (SignDist(i,j,k) < 0.0) id[n] = 0;
if (SignDist(i,j,k) < 0.f){
// don't do anything
}
else{
// initially saturated with wetting phase
id[n] = 2;
@@ -392,21 +424,8 @@ int main(int argc, char **argv)
if (rank==0) printf("Writing ID file \n");
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *ID = fopen(LocalRankFilename,"rb");
readID=fread(Dm->id,1,N,ID);
if (readID != size_t(N)) printf("lbpm_segmented_pp: Error reading ID \n");
fclose(ID);
// Preserve mineral labels
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;
if (SignDist(i,j,k) < 0.0) id[n] = Dm->id[n];
}
}
}
ID = fopen(LocalRankFilename,"wb");
FILE *ID = fopen(LocalRankFilename,"wb");
fwrite(id,1,N,ID);
fclose(ID);
}

226
tests/lbpm_uCT_pp.cpp
View File

@@ -73,29 +73,28 @@ int main(int argc, char **argv)
int nprocz = nproc[2];
auto InputFile = uct_db->getScalar<std::string>( "InputFile" );
auto target=uct_db->getScalar<int>("target");
auto background=uct_db->getScalar<int>("background");
auto target = uct_db->getScalar<float>("target");
auto background = uct_db->getScalar<float>("background");
auto rough_cutoff = uct_db->getScalar<float>( "rough_cutoff" );
auto lamda = uct_db->getScalar<float>( "lamda" );
auto nlm_sigsq = uct_db->getScalar<float>( "nlm_sigsq" );
auto nlm_depth = uct_db->getScalar<int>( "nlm_depth" );
auto cx=uct_db->getScalar<int>( "center_x" );
auto cy=uct_db->getScalar<int>( "center_y" );
auto cz=uct_db->getScalar<int>( "center_z" );
auto center = uct_db->getVector<int>( "center" );
auto CylRad = uct_db->getScalar<float>( "cylinder_radius" );
//.......................................................................
// Reading the domain information file
//.......................................................................
// std::shared_ptr<Domain> Dm ();
//for (int i=0; i<Dm->Nx*Dm->Ny*Dm->Nz; i++) Dm->id[i] = 1;
//Dm->CommInit();
auto maxLevels = uct_db->getScalar<int>( "max_levels" );
std::vector<int> offset( 3, 0 );
if ( uct_db->keyExists( "offset" ) )
offset = uct_db->getVector<int>( "offset" );
// Check that the number of processors >= the number of ranks
if ( rank==0 ) {
printf("Number of MPI ranks required: %i \n", nprocx*nprocy*nprocz);
printf("Number of MPI ranks used: %i \n", nprocs);
printf("Full domain size: %i x %i x %i \n",nx*nprocx,ny*nprocy,nz*nprocz);
printf("target value = %f \n",target);
printf("background value = %f \n",background);
printf("cylinder center = %i, %i, %i \n",center[0],center[1],center[2]);
printf("cylinder radius = %f \n",CylRad);
}
if ( nprocs < nprocx*nprocy*nprocz ){
ERROR("Insufficient number of processors");
@@ -105,18 +104,20 @@ int main(int argc, char **argv)
int ratio[3] = {2,2,2};
//std::vector<size_t> ratio = {4,4,4};
// need to set up databases for each level of the mesh
std:vector<Database> multidomain_db;
std::vector<std::shared_ptr<Database>> multidomain_db(1,domain_db);
std::vector<int> Nx(1,nx), Ny(1,ny), Nz(1,nz);
while ( Nx.back()%ratio[0]==0 && Nx.back()>8 &&
Ny.back()%ratio[1]==0 && Ny.back()>8 &&
Nz.back()%ratio[2]==0 && Nz.back()>8 )
Nz.back()%ratio[2]==0 && Nz.back()>8 &&
(int) Nx.size() < maxLevels )
{
Nx.push_back( Nx.back()/ratio[0] );
Ny.push_back( Ny.back()/ratio[1] );
Nz.push_back( Nz.back()/ratio[2] );
// clone the domain and create coarse version based on Nx,Ny,Nz
//multidomain_db.push_back();
auto db2 = domain_db->cloneDatabase();
db2->putVector<int>( "n", { Nx.back(), Ny.back(), Nz.back() } );
multidomain_db.push_back(db2);
}
int N_levels = Nx.size();
@@ -125,7 +126,7 @@ int main(int argc, char **argv)
for (int i=0; i<N_levels; i++) {
// This line is no good -- will create identical Domain structures instead of
// Need a way to define a coarse structure for the coarse domain (see above)
Dm[i].reset( new Domain(domain_db, comm) );
Dm[i].reset( new Domain(multidomain_db[i], comm) );
int N = (Nx[i]+2)*(Ny[i]+2)*(Nz[i]+2);
for (int n=0; n<N; n++){
Dm[i]->id[n] = 1;
@@ -135,6 +136,7 @@ int main(int argc, char **argv)
// array containing a distance mask
Array<float> MASK(Nx[0]+2,Ny[0]+2,Nz[0]+2);
MASK.fill(0);
// Create the level data
std::vector<Array<char>> ID(N_levels);
@@ -168,18 +170,17 @@ int main(int argc, char **argv)
PROFILE_START("ReadVolume");
int fid = netcdf::open(InputFile,netcdf::READ);
std::string varname("VOLUME");
netcdf::VariableType type = netcdf::getVarType( fid, varname );
std::vector<size_t> dim = netcdf::getVarDim( fid, varname );
auto type = netcdf::getVarType( fid, varname );
auto dim = netcdf::getVarDim( fid, varname );
if ( rank == 0 ) {
printf("Reading %s (%s)\n",varname.c_str(),netcdf::VariableTypeName(type).c_str());
printf(" dims = %i x %i x %i \n",int(dim[0]),int(dim[1]),int(dim[2]));
}
{
RankInfoStruct info( rank, nprocx, nprocy, nprocz );
int x = info.ix*nx;
int y = info.jy*ny;
int z = info.kz*nz;
// Read the local data
int x = Dm[0]->iproc()*nx + offset[0];
int y = Dm[0]->jproc()*ny + offset[1];
int z = Dm[0]->kproc()*nz + offset[2];
Array<short> VOLUME = netcdf::getVar<short>( fid, varname, {x,y,z}, {nx,ny,nz}, {1,1,1} );
// Copy the data and fill the halos
LOCVOL[0].fill(0);
@@ -191,30 +192,23 @@ int main(int argc, char **argv)
PROFILE_STOP("ReadVolume");
if (rank==0) printf("Read complete\n");
// Filter the original data
filter_src( *Dm[0], LOCVOL[0] );
// Set up the mask to be distance to cylinder (crop outside cylinder)
// float CylRad=900;
if (rank==0) printf("Cropping with cylinder: %i, %i, %i, radius=%f \n",Dm[0]->nprocx()*Nx[0],Dm[0]->nprocy()*Ny[0],Dm[0]->nprocz()*Nz[0],CylRad);
for (int k=0;k<Nz[0]+2;k++) {
for (int j=0;j<Ny[0]+2;j++) {
for (int i=0;i<Nx[0]+2;i++) {
int iproc = Dm[0]->iproc();
int jproc = Dm[0]->jproc();
int kproc = Dm[0]->kproc();
int x=iproc*Nx[0]+i-1;
int y=jproc*Ny[0]+j-1;
int z=kproc*Nz[0]+k-1;
//int cx = 0.5*nprocx*Nx[0];
//int cy = 0.5*nprocy*Ny[0];
//int cz = 0.5*nprocz*Nz[0];
float x= float(Dm[0]->iproc()*Nx[0]+i-1);
float y= float (Dm[0]->jproc()*Ny[0]+j-1);
float z= float(Dm[0]->kproc()*Nz[0]+k-1);
float cx = float(center[0] - offset[0]);
float cy = float(center[1] - offset[1]);
float cz = float(center[2] - offset[2]);
// distance from the center line
MASK(i,j,k) = CylRad - sqrt(float((z-cz)*(z-cz) + (y-cy)*(y-cy)) );
MASK(i,j,k) = sqrt((z-cz)*(z-cz) + (y-cy)*(y-cy));
//if (sqrt(((z-cz)*(z-cz) + (y-cy)*(y-cy)) ) > CylRad) LOCVOL[0](i,j,k)=background;
}
}
}
@@ -222,15 +216,20 @@ int main(int argc, char **argv)
// Compute the means for the high/low regions
// (should use automated mixture model to approximate histograms)
//float THRESHOLD = 0.05*maxReduce( Dm[0]->Comm, std::max( LOCVOL[0].max(), fabs(LOCVOL[0].min()) ) );
float THRESHOLD=0.5*float(target+background);
float THRESHOLD=0.5*(target+background);
float mean_plus=0;
float mean_minus=0;
float min_value = LOCVOL[0](0);
float max_value = LOCVOL[0](0);
int count_plus=0;
int count_minus=0;
for (int k=1;k<Nz[0]+1;k++) {
for (int j=1;j<Ny[0]+1;j++) {
for (int i=1;i<Nx[0]+1;i++) {
if (MASK(i,j,k) > 0.f ){
//LOCVOL[0](i,j,k) = MASK(i,j,k);
if (MASK(i,j,k) < CylRad ){
auto tmp = LOCVOL[0](i,j,k);
/* if ((tmp-background)*(tmp-target) > 0){
// direction to background / target is the same
@@ -241,35 +240,47 @@ int main(int argc, char **argv)
if ( tmp > THRESHOLD ) {
mean_plus += tmp;
count_plus++;
} else if ( tmp < -THRESHOLD ) {
}
else {
mean_minus += tmp;
count_minus++;
}
if (tmp < min_value) min_value = tmp;
if (tmp > max_value) max_value = tmp;
}
}
}
}
mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / sumReduce( Dm[0]->Comm, count_plus );
mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / sumReduce( Dm[0]->Comm, count_minus );
count_plus=sumReduce( Dm[0]->Comm, count_plus);
count_minus=sumReduce( Dm[0]->Comm, count_minus);
if (rank==0) printf("minimum value=%f, max value=%f \n",min_value,max_value);
if (rank==0) printf("plus=%i, minus=%i \n",count_plus,count_minus);
ASSERT( count_plus > 0 && count_minus > 0 );
MPI_Barrier(comm);
mean_plus = sumReduce( Dm[0]->Comm, mean_plus ) / count_plus;
mean_minus = sumReduce( Dm[0]->Comm, mean_minus ) / count_minus;
MPI_Barrier(comm);
if (rank==0) printf(" Region 1 mean (+): %f, Region 2 mean (-): %f \n",mean_plus, mean_minus);
MPI_Barrier(comm);
// Scale the source data to +-1.0
//if (rank==0) printf("Scale the input data (size = %i) \n",LOCVOL[0].length());
for (size_t i=0; i<LOCVOL[0].length(); i++) {
if (MASK(i) < 0.f){
LOCVOL[0](i) = 1.0;
if ( MASK(i) > CylRad ){
LOCVOL[0](i)=background;
}
else if ( LOCVOL[0](i) >= 0 ) {
LOCVOL[0](i) /= mean_plus;
LOCVOL[0](i) = std::min( LOCVOL[0](i), 1.0f );
} else {
LOCVOL[0](i) /= -mean_minus;
LOCVOL[0](i) = std::max( LOCVOL[0](i), -1.0f );
if ( LOCVOL[0](i) >= THRESHOLD ) {
auto tmp = LOCVOL[0](i)/ mean_plus;
LOCVOL[0](i) = std::min( tmp, 1.0f );
}
else {
auto tmp = -LOCVOL[0](i)/mean_minus;
LOCVOL[0](i) = std::max( tmp, -1.0f );
}
//LOCVOL[0](i) = MASK(i);
}
// Fill the source data for the coarse meshes
if (rank==0) printf("Coarsen the mesh for N_levels=%i \n",N_levels);
MPI_Barrier(comm);
PROFILE_START("CoarsenMesh");
for (int i=1; i<N_levels; i++) {
Array<float> filter(ratio[0],ratio[1],ratio[2]);
@@ -290,21 +301,12 @@ int main(int argc, char **argv)
PROFILE_STOP("CoarsenMesh");
// Initialize the coarse level
PROFILE_START("Solve full mesh");
if (rank==0)
printf("Initialize full mesh\n");
solve( LOCVOL[0], Mean[0], ID[0], Dist[0], MultiScaleSmooth[0],
NonLocalMean[0], *fillFloat[0], *Dm[0], nprocx,
rough_cutoff, lamda, nlm_sigsq, nlm_depth);
PROFILE_STOP("Solve full mesh");
MPI_Barrier(comm);
/* // Initialize the coarse level
PROFILE_START("Solve coarse mesh");
if (rank==0)
printf("Initialize coarse mesh\n");
printf("Initialize full mesh\n");
solve( LOCVOL.back(), Mean.back(), ID.back(), Dist.back(), MultiScaleSmooth.back(),
NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx );
NonLocalMean.back(), *fillFloat.back(), *Dm.back(), nprocx,
rough_cutoff, lamda, nlm_sigsq, nlm_depth);
PROFILE_STOP("Solve coarse mesh");
MPI_Barrier(comm);
@@ -312,11 +314,12 @@ int main(int argc, char **argv)
PROFILE_START("Refine distance");
if (rank==0)
printf("Refine mesh\n");
for (int i=int(Nx.size())-2; i>=0; i--) {
for (int i=N_levels-2; i>=0; i--) {
if (rank==0)
printf(" Refining to level %i\n",int(i));
printf(" Refining to level %i\n",i);
refine( Dist[i+1], LOCVOL[i], Mean[i], ID[i], Dist[i], MultiScaleSmooth[i],
NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i );
NonLocalMean[i], *fillFloat[i], *Dm[i], nprocx, i,
rough_cutoff, lamda, nlm_sigsq, nlm_depth);
}
PROFILE_STOP("Refine distance");
MPI_Barrier(comm);
@@ -328,12 +331,10 @@ int main(int argc, char **argv)
Array<float> filter_Mean, filter_Dist1, filter_Dist2;
filter_final( ID[0], Dist[0], *fillFloat[0], *Dm[0], filter_Mean, filter_Dist1, filter_Dist2 );
PROFILE_STOP("Filtering final domains");
*/
//removeDisconnected( ID[0], *Dm[0] );
/*
// Write the distance function to a netcdf file
const char* netcdf_filename = "Distance.nc";
/* const char* netcdf_filename = "Distance.nc";
{
RankInfoStruct info( rank, nprocx, nprocy, nprocz );
std::vector<int> dim = { Nx[0]*nprocx, Ny[0]*nprocy, Nz[0]*nprocz };
@@ -343,77 +344,50 @@ int main(int argc, char **argv)
fillFloat[0]->copy( Dist[0], data );
netcdf::write( fid, "Distance", dims, data, info );
netcdf::close( fid );
}
*/
} */
{
// Write the results
if (rank==0) printf("Setting up visualization structure \n");
// std::vector<IO::MeshDataStruct> meshData(N_levels);
std::vector<IO::MeshDataStruct> meshData(1);
// for (size_t i=0; i<Nx.size(); i++) {
std::vector<IO::MeshDataStruct> meshData(N_levels);
for (size_t i=0; i<Nx.size(); i++) {
// Mesh
meshData[0].meshName = "image";
meshData[0].mesh = std::shared_ptr<IO::DomainMesh>( new IO::DomainMesh(Dm[0]->rank_info,Nx[0],Ny[0],Nz[0],Lx,Ly,Lz) );
meshData[i].meshName = "image_" + std::to_string( i );
meshData[i].mesh = std::make_shared<IO::DomainMesh>(Dm[i]->rank_info,Nx[i],Ny[i],Nz[i],Lx,Ly,Lz);
// Source data
std::shared_ptr<IO::Variable> OrigData( new IO::Variable() );
OrigData->name = "Source Data";
auto OrigData = std::make_shared<IO::Variable>();
OrigData->name = "Source_Data_" + std::to_string( i );
OrigData->type = IO::VariableType::VolumeVariable;
OrigData->dim = 1;
OrigData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(OrigData);
fillDouble[0]->copy( LOCVOL[0], OrigData->data );
OrigData->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(OrigData);
fillDouble[i]->copy( LOCVOL[i], OrigData->data );
// Non-Local Mean
std::shared_ptr<IO::Variable> NonLocMean( new IO::Variable() );
NonLocMean->name = "Non-Local Mean";
auto NonLocMean = std::make_shared<IO::Variable>();
NonLocMean->name = "NonLocal_Mean_" + std::to_string( i );
NonLocMean->type = IO::VariableType::VolumeVariable;
NonLocMean->dim = 1;
NonLocMean->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(NonLocMean);
fillDouble[0]->copy( NonLocalMean[0], NonLocMean->data );
std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
SegData->name = "Segmented Data";
SegData->type = IO::VariableType::VolumeVariable;
SegData->dim = 1;
SegData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(SegData);
fillDouble[0]->copy( ID[0], SegData->data );
// Signed Distance
std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
DistData->name = "Signed Distance";
DistData->type = IO::VariableType::VolumeVariable;
DistData->dim = 1;
DistData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(DistData);
fillDouble[0]->copy( Dist[0], DistData->data );
// Smoothed Data
std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
SmoothData->name = "Smoothed Data";
SmoothData->type = IO::VariableType::VolumeVariable;
SmoothData->dim = 1;
SmoothData->data.resize(Nx[0],Ny[0],Nz[0]);
meshData[0].vars.push_back(SmoothData);
fillDouble[0]->copy( MultiScaleSmooth[0], SmoothData->data );
/*// Segmented Data
std::shared_ptr<IO::Variable> SegData( new IO::Variable() );
SegData->name = "Segmented Data";
NonLocMean->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(NonLocMean);
fillDouble[i]->copy( NonLocalMean[i], NonLocMean->data );
// Segmented Data
auto SegData = std::make_shared<IO::Variable>();
SegData->name = "Segmented_Data_" + std::to_string( i );
SegData->type = IO::VariableType::VolumeVariable;
SegData->dim = 1;
SegData->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(SegData);
fillDouble[i]->copy( ID[i], SegData->data );
// Signed Distance
std::shared_ptr<IO::Variable> DistData( new IO::Variable() );
DistData->name = "Signed Distance";
auto DistData = std::make_shared<IO::Variable>();
DistData->name = "Signed_Distance_" + std::to_string( i );
DistData->type = IO::VariableType::VolumeVariable;
DistData->dim = 1;
DistData->data.resize(Nx[i],Ny[i],Nz[i]);
meshData[i].vars.push_back(DistData);
fillDouble[i]->copy( Dist[i], DistData->data );
// Smoothed Data
std::shared_ptr<IO::Variable> SmoothData( new IO::Variable() );
SmoothData->name = "Smoothed Data";
auto SmoothData = std::make_shared<IO::Variable>();
SmoothData->name = "Smoothed_Data_" + std::to_string( i );
SmoothData->type = IO::VariableType::VolumeVariable;
SmoothData->dim = 1;
SmoothData->data.resize(Nx[i],Ny[i],Nz[i]);
@@ -444,20 +418,18 @@ int main(int argc, char **argv)
meshData[0].vars.push_back(filter_Dist2_var);
fillDouble[0]->copy( filter_Dist2, filter_Dist2_var->data );
#endif
*/
MPI_Barrier(comm);
if (rank==0) printf("Writing output \n");
// Write visulization data
IO::writeData( 0, meshData, comm );
if (rank==0) printf("Finished. \n");
}
// Compute the Minkowski functionals
MPI_Barrier(comm);
std::shared_ptr<Minkowski> Averages(new Minkowski(Dm[0]));
auto Averages = std::make_shared<Minkowski>(Dm[0]);
Array <char> phase_label(Nx[0],Ny[0],Nz[0]);
Array <double> phase_distance(Nx[0],Ny[0],Nz[0]);
Array <char> phase_label(Nx[0]+2,Ny[0]+2,Nz[0]+2);
Array <double> phase_distance(Nx[0]+2,Ny[0]+2,Nz[0]+2);
// Analyze the wetting fluid
for (int k=1;k<Nz[0]+1;k++) {
for (int j=1;j<Ny[0]+1;j++) {