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

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

Browse Source
This commit is contained in:
James E McClure
2019-04-04 20:25:50 -04:00
parent 44b746acc0 ace0c78ff3
commit e169e1a089
18 changed files with 527 additions and 339 deletions
Download Patch File Download Diff File Expand all files Collapse all files

226
analysis/SubPhase.cpp
View File

@@ -28,23 +28,32 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
//.........................................
if (Dm->rank()==0){
TIMELOG = fopen("subphase.csv","a+");
if (fseek(TIMELOG,0,SEEK_SET) == fseek(TIMELOG,0,SEEK_CUR))
bool WriteHeader=false;
SUBPHASE = fopen("subphase.csv","r");
if (SUBPHASE != NULL)
fclose(SUBPHASE);
else
WriteHeader=true;
SUBPHASE = fopen("subphase.csv","a+");
if (WriteHeader)
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn ");
fprintf(TIMELOG,"pwc pwd pnc pnd "); // pressures
fprintf(TIMELOG,"Mwc Mwd Mwi Mnc Mnd Mni "); // mass
fprintf(TIMELOG,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x "); // momentum
fprintf(TIMELOG,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");
fprintf(TIMELOG,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");
fprintf(TIMELOG,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
fprintf(TIMELOG,"Vwc Awc Hwc Xwc "); // wc region
fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd "); // wd region
fprintf(TIMELOG,"Vnc Anc Hnc Xnc "); // nc region
fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd "); // nd region
fprintf(TIMELOG,"Vi Ai Hi Xi\n"); // interface region
//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn ");
fprintf(SUBPHASE,"pwc pwd pnc pnd "); // pressures
fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni "); // mass
fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x "); // momentum
fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");
fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");
fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
fprintf(SUBPHASE,"Vwc Awc Hwc Xwc "); // wc region
fprintf(SUBPHASE,"Vwd Awd Hwd Xwd Nwd "); // wd region
fprintf(SUBPHASE,"Vnc Anc Hnc Xnc "); // nc region
fprintf(SUBPHASE,"Vnd And Hnd Xnd Nnd "); // nd region
fprintf(SUBPHASE,"Vi Ai Hi Xi "); // interface region
fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n"); // interface region
// stress tensor?
}
@@ -54,20 +63,38 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
sprintf(LocalRankString,"%05d",Dm->rank());
char LocalRankFilename[40];
sprintf(LocalRankFilename,"%s%s","subphase.csv.",LocalRankString);
TIMELOG = fopen(LocalRankFilename,"a+");
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"time rn rw nun nuw Fx Fy Fz iftwn ");
fprintf(TIMELOG,"pwc pwd pnc pnd "); // pressures
fprintf(TIMELOG,"Mwc Mwd Mwi Mnc Mnd Mni "); // mass
fprintf(TIMELOG,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x "); // momentum
fprintf(TIMELOG,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");
fprintf(TIMELOG,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");
fprintf(TIMELOG,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
fprintf(TIMELOG,"Vwc Awc Hwc Xwc "); // wc region
fprintf(TIMELOG,"Vwd Awd Hwd Xwd Nwd "); // wd region
fprintf(TIMELOG,"Vnc Anc Hnc Xnc "); // nc region
fprintf(TIMELOG,"Vnd And Hnd Xnd Nnd "); // nd region
fprintf(TIMELOG,"Vi Ai Hi Xi\n"); // interface region
SUBPHASE = fopen(LocalRankFilename,"a+");
//fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n");
fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn ");
fprintf(SUBPHASE,"pwc pwd pnc pnd "); // pressures
fprintf(SUBPHASE,"Mwc Mwd Mwi Mnc Mnd Mni "); // mass
fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x "); // momentum
fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y ");
fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z ");
fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy
fprintf(SUBPHASE,"Vwc Awc Hwc Xwc "); // wc region
fprintf(SUBPHASE,"Vwd Awd Hwd Xwd Nwd "); // wd region
fprintf(SUBPHASE,"Vnc Anc Hnc Xnc "); // nc region
fprintf(SUBPHASE,"Vnd And Hnd Xnd Nnd "); // nd region
fprintf(SUBPHASE,"Vi Ai Hi Xi "); // interface region
fprintf(SUBPHASE,"Vic Aic Hic Xic Nic\n"); // interface region
}
if (Dm->rank()==0){
bool WriteHeader=false;
TIMELOG = fopen("timelog.csv","r");
if (TIMELOG != NULL)
fclose(TIMELOG);
else
WriteHeader=true;
TIMELOG = fopen("timelog.csv","a+");
if (WriteHeader)
{
// If timelog is empty, write a short header to list the averages
//fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n");
fprintf(TIMELOG,"sw krw krn qw qn pw pn\n");
}
}
}
@@ -75,41 +102,42 @@ SubPhase::SubPhase(std::shared_ptr <Domain> dm):
// Destructor
SubPhase::~SubPhase()
{
if ( TIMELOG!=NULL ) { fclose(TIMELOG); }
if ( SUBPHASE!=NULL ) { fclose(SUBPHASE); }
}
void SubPhase::Write(int timestep)
{
if (Dm->rank()==0){
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 %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Px, gwd.Px, giwn.Px, gnc.Px, gnd.Px, giwn.Px);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Py, gnc.Py, gnd.Py, giwn.Py);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pz, gnc.Pz, gnd.Pz, giwn.Pz);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",giwn.V, giwn.A, giwn.H, giwn.X);
fflush(TIMELOG);
fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.p, gwd.p, gnc.p, gnd.p);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gwc.V, gwc.A, gwc.H, gwc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",gnc.V, gnc.A, gnc.H, gnc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",gnd.V, gnd.A, gnd.H, gnd.X, gnd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",giwn.V, giwn.A, giwn.H, giwn.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i\n",giwnc.V, giwnc.A, giwnc.H, giwnc.X, giwnc.Nc);
fflush(SUBPHASE);
}
else{
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 %.5g ",wc.p, wd.p, nc.p, nd.p);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Px, nc.Px, nd.Px, iwn.Px);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Py, nc.Py, nd.Py, iwn.Py);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pz, nc.Pz, nd.Pz, iwn.Pz);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g\n",iwn.V, iwn.A, iwn.H, iwn.X);
fflush(TIMELOG);
fprintf(SUBPHASE,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.p, wd.p, nc.p, nd.p);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",wc.V, wc.A, wc.H, wc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",nc.V, nc.A, nc.H, nc.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %i ",nd.V, nd.A, nd.H, nd.X, nd.Nc);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g ",iwn.V, iwn.A, iwn.H, iwn.X);
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g\n",iwnc.V, iwnc.A, iwnc.H, iwnc.X);
}
}
@@ -136,28 +164,17 @@ void SubPhase::Basic(){
if (Dm->BoundaryCondition > 0 && Dm->kproc() == Dm->nprocz()-1) kmax=Nz-4;
imin=jmin=1;
// If inlet layers exist use these as default
// If inlet/outlet layers exist use these as default
if (Dm->inlet_layers_x > 0) imin = Dm->inlet_layers_x;
if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z;
nb.reset(); wb.reset();
/*
//Dm->CommunicateMeshHalo(Phi);
for (int k=1; k<Nz-1; k++){
for (int j=1; j<Ny-1; j++){
for (int i=1; i<Nx-1; i++){
// Compute all of the derivatives using finite differences
double fx = 0.5*(Phi(i+1,j,k) - Phi(i-1,j,k));
double fy = 0.5*(Phi(i,j+1,k) - Phi(i,j-1,k));
double fz = 0.5*(Phi(i,j,k+1) - Phi(i,j,k-1));
DelPhi(i,j,k) = sqrt(fx*fx+fy*fy+fz*fz);
}
}
}
*/
if (Dm->outlet_layers_z > 0) kmax = Dm->outlet_layers_z;
nb.reset(); wb.reset();
double nA,nB;
double count_w = 0.0;
double count_n = 0.0;
for (k=kmin; k<kmax; k++){
for (j=jmin; j<Ny-1; j++){
for (i=imin; i<Nx-1; i++){
@@ -187,6 +204,14 @@ void SubPhase::Basic(){
wb.Py += rho_w*nB*Vel_y(n);
wb.Pz += rho_w*nB*Vel_z(n);
}
if ( phi > 0.99 ){
nb.p += Pressure(n);
count_n += 1.0;
}
else if ( phi < -0.99 ){
wb.p += Pressure(n);
count_w += 1.0;
}
}
}
}
@@ -201,14 +226,35 @@ void SubPhase::Basic(){
gnb.Px=sumReduce( Dm->Comm, nb.Px);
gnb.Py=sumReduce( Dm->Comm, nb.Py);
gnb.Pz=sumReduce( Dm->Comm, nb.Pz);
count_w=sumReduce( Dm->Comm, count_w);
count_n=sumReduce( Dm->Comm, count_n);
gwb.p=sumReduce( Dm->Comm, wb.p) / count_w;
gnb.p=sumReduce( Dm->Comm, nb.p) / count_n;
if (Dm->rank() == 0){
double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = Fx/force_mag;
double dir_y = Fy/force_mag;
double dir_z = Fz/force_mag;
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
force_mag = 1.0;
}
double saturation=gwb.V/(gwb.V + gnb.V);
double water_flow_rate=gwb.V*sqrt(gwb.Px*gwb.Px + gwb.Py*gwb.Py + gwb.Pz*gwb.Pz)/gwb.M;
double not_water_flow_rate=gnb.V*sqrt(gnb.Px*gnb.Px + gnb.Py*gnb.Py + gnb.Pz*gnb.Pz)/gnb.M;
double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M;
double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M;
double total_flow_rate = water_flow_rate + not_water_flow_rate;
double fractional_flow= water_flow_rate / total_flow_rate;
printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
double krn = nu_n*not_water_flow_rate / force_mag;
double krw = nu_w*water_flow_rate / force_mag;
//printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,water_flow_rate,not_water_flow_rate, gwb.p, gnb.p);
fflush(TIMELOG);
}
}
@@ -285,7 +331,7 @@ void SubPhase::Full(){
if (Dm->inlet_layers_y > 0) jmin = Dm->inlet_layers_y;
if (Dm->inlet_layers_z > 0) kmin = Dm->inlet_layers_z;
nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset();
nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); iwnc.reset();
Dm->CommunicateMeshHalo(Phi);
for (int k=1; k<Nz-1; k++){
@@ -429,7 +475,18 @@ void SubPhase::Full(){
giwn.A=sumReduce( Dm->Comm, iwn.A);
giwn.H=sumReduce( Dm->Comm, iwn.H);
giwn.X=sumReduce( Dm->Comm, iwn.X);
// measure only the connected part
iwnc.Nc = morph_i->MeasureConnectedPathway();
iwnc.V = morph_i->V();
iwnc.A = morph_i->A();
iwnc.H = morph_i->H();
iwnc.X = morph_i->X();
giwnc.V=sumReduce( Dm->Comm, iwnc.V);
giwnc.A=sumReduce( Dm->Comm, iwnc.A);
giwnc.H=sumReduce( Dm->Comm, iwnc.H);
giwnc.X=sumReduce( Dm->Comm, iwnc.X);
giwnc.Nc = iwnc.Nc;
double vol_nc_bulk = 0.0;
double vol_wc_bulk = 0.0;
double vol_nd_bulk = 0.0;
@@ -528,7 +585,18 @@ void SubPhase::Full(){
gwc.Pz=sumReduce( Dm->Comm, wc.Pz);
gwc.K=sumReduce( Dm->Comm, wc.K);
gwc.p=sumReduce( Dm->Comm, wc.p);
giwn.Mn=sumReduce( Dm->Comm, iwn.Mn);
giwn.Pnx=sumReduce( Dm->Comm, iwn.Pnx);
giwn.Pny=sumReduce( Dm->Comm, iwn.Pny);
giwn.Pnz=sumReduce( Dm->Comm, iwn.Pnz);
giwn.Kn=sumReduce( Dm->Comm, iwn.Kn);
giwn.Mw=sumReduce( Dm->Comm, iwn.Mw);
giwn.Pwx=sumReduce( Dm->Comm, iwn.Pwx);
giwn.Pwy=sumReduce( Dm->Comm, iwn.Pwy);
giwn.Pwz=sumReduce( Dm->Comm, iwn.Pwz);
giwn.Kw=sumReduce( Dm->Comm, iwn.Kw);
// pressure averaging
if (vol_wc_bulk > 0.0)
wc.p = wc.p /vol_wc_bulk;

10
analysis/SubPhase.h
View File

@@ -11,8 +11,6 @@
#include "analysis/analysis.h"
#include "analysis/distance.h"
#include "analysis/Minkowski.h"
#include "shared_ptr.h"
#include "common/Utilities.h"
#include "common/MPI_Helpers.h"
#include "IO/MeshDatabase.h"
@@ -37,10 +35,12 @@ private:
class interface{
public:
int Nc;
double M,Px,Py,Pz,K;
double Mw,Mn,Pnx,Pny,Pnz,Pwx,Pwy,Pwz,Kw,Kn;
double V,A,H,X;
void reset(){
Nc = 0;
M=Px=Py=Pz=K=0.0;
V=A=H=X=0.0;
Mw=Mn=Pnx=Pny=Pnz=Pwx=Pwy=Pwz=Kw=Kn=0.0;
@@ -69,11 +69,11 @@ public:
*/
// local entities
phase wc,wd,wb,nc,nd,nb;
interface iwn;
interface iwn,iwnc;
// global entities
phase gwc,gwd,gwb,gnc,gnd,gnb;
interface giwn;
interface giwn,giwnc;
//...........................................................................
int Nx,Ny,Nz;
IntArray PhaseID; // Phase ID array (solid=0, non-wetting=1, wetting=2)
@@ -104,7 +104,7 @@ public:
private:
FILE *TIMELOG;
FILE *SUBPHASE;
};
#endif

196
analysis/morphology.cpp
View File

@@ -1,7 +1,7 @@
#include <analysis/morphology.h>
// Implementation of morphological opening routine
inline void PackID(int *list, int count, char *sendbuf, char *ID){
inline void PackID(int *list, int count, signed char *sendbuf, signed char *ID){
// Fill in the phase ID values from neighboring processors
// This packs up the values that need to be sent from one processor to another
int idx,n;
@@ -13,7 +13,7 @@ inline void PackID(int *list, int count, char *sendbuf, char *ID){
}
//***************************************************************************************
inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
inline void UnpackID(int *list, int count, signed char *recvbuf, signed char *ID){
// Fill in the phase ID values from neighboring processors
// This unpacks the values once they have been recieved from neighbors
int idx,n;
@@ -25,7 +25,7 @@ inline void UnpackID(int *list, int count, char *recvbuf, char *ID){
}
//***************************************************************************************
double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
// SignDist is the distance to the object that you want to constaing the morphological opening
// VoidFraction is the the empty space where the object inst
// id is a labeled map
@@ -73,51 +73,51 @@ double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, do
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
// Communication buffers
char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
// send buffers
sendID_x = new char [Dm->sendCount_x];
sendID_y = new char [Dm->sendCount_y];
sendID_z = new char [Dm->sendCount_z];
sendID_X = new char [Dm->sendCount_X];
sendID_Y = new char [Dm->sendCount_Y];
sendID_Z = new char [Dm->sendCount_Z];
sendID_xy = new char [Dm->sendCount_xy];
sendID_yz = new char [Dm->sendCount_yz];
sendID_xz = new char [Dm->sendCount_xz];
sendID_Xy = new char [Dm->sendCount_Xy];
sendID_Yz = new char [Dm->sendCount_Yz];
sendID_xZ = new char [Dm->sendCount_xZ];
sendID_xY = new char [Dm->sendCount_xY];
sendID_yZ = new char [Dm->sendCount_yZ];
sendID_Xz = new char [Dm->sendCount_Xz];
sendID_XY = new char [Dm->sendCount_XY];
sendID_YZ = new char [Dm->sendCount_YZ];
sendID_XZ = new char [Dm->sendCount_XZ];
sendID_x = new signed char [Dm->sendCount_x];
sendID_y = new signed char [Dm->sendCount_y];
sendID_z = new signed char [Dm->sendCount_z];
sendID_X = new signed char [Dm->sendCount_X];
sendID_Y = new signed char [Dm->sendCount_Y];
sendID_Z = new signed char [Dm->sendCount_Z];
sendID_xy = new signed char [Dm->sendCount_xy];
sendID_yz = new signed char [Dm->sendCount_yz];
sendID_xz = new signed char [Dm->sendCount_xz];
sendID_Xy = new signed char [Dm->sendCount_Xy];
sendID_Yz = new signed char [Dm->sendCount_Yz];
sendID_xZ = new signed char [Dm->sendCount_xZ];
sendID_xY = new signed char [Dm->sendCount_xY];
sendID_yZ = new signed char [Dm->sendCount_yZ];
sendID_Xz = new signed char [Dm->sendCount_Xz];
sendID_XY = new signed char [Dm->sendCount_XY];
sendID_YZ = new signed char [Dm->sendCount_YZ];
sendID_XZ = new signed char [Dm->sendCount_XZ];
//......................................................................................
// recv buffers
recvID_x = new char [Dm->recvCount_x];
recvID_y = new char [Dm->recvCount_y];
recvID_z = new char [Dm->recvCount_z];
recvID_X = new char [Dm->recvCount_X];
recvID_Y = new char [Dm->recvCount_Y];
recvID_Z = new char [Dm->recvCount_Z];
recvID_xy = new char [Dm->recvCount_xy];
recvID_yz = new char [Dm->recvCount_yz];
recvID_xz = new char [Dm->recvCount_xz];
recvID_Xy = new char [Dm->recvCount_Xy];
recvID_xZ = new char [Dm->recvCount_xZ];
recvID_xY = new char [Dm->recvCount_xY];
recvID_yZ = new char [Dm->recvCount_yZ];
recvID_Yz = new char [Dm->recvCount_Yz];
recvID_Xz = new char [Dm->recvCount_Xz];
recvID_XY = new char [Dm->recvCount_XY];
recvID_YZ = new char [Dm->recvCount_YZ];
recvID_XZ = new char [Dm->recvCount_XZ];
recvID_x = new signed char [Dm->recvCount_x];
recvID_y = new signed char [Dm->recvCount_y];
recvID_z = new signed char [Dm->recvCount_z];
recvID_X = new signed char [Dm->recvCount_X];
recvID_Y = new signed char [Dm->recvCount_Y];
recvID_Z = new signed char [Dm->recvCount_Z];
recvID_xy = new signed char [Dm->recvCount_xy];
recvID_yz = new signed char [Dm->recvCount_yz];
recvID_xz = new signed char [Dm->recvCount_xz];
recvID_Xy = new signed char [Dm->recvCount_Xy];
recvID_xZ = new signed char [Dm->recvCount_xZ];
recvID_xY = new signed char [Dm->recvCount_xY];
recvID_yZ = new signed char [Dm->recvCount_yZ];
recvID_Yz = new signed char [Dm->recvCount_Yz];
recvID_Xz = new signed char [Dm->recvCount_Xz];
recvID_XY = new signed char [Dm->recvCount_XY];
recvID_YZ = new signed char [Dm->recvCount_YZ];
recvID_XZ = new signed char [Dm->recvCount_XZ];
//......................................................................................
int sendtag,recvtag;
sendtag = recvtag = 7;
@@ -327,7 +327,7 @@ double morph_open()
*/
//***************************************************************************************
double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction){
// SignDist is the distance to the object that you want to constaing the morphological opening
// VoidFraction is the the empty space where the object inst
// id is a labeled map
@@ -378,51 +378,51 @@ double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, d
if (rank==0) printf("Maximum pore size: %f \n",maxdistGlobal);
// Communication buffers
char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
signed char *sendID_x, *sendID_y, *sendID_z, *sendID_X, *sendID_Y, *sendID_Z;
signed char *sendID_xy, *sendID_yz, *sendID_xz, *sendID_Xy, *sendID_Yz, *sendID_xZ;
signed char *sendID_xY, *sendID_yZ, *sendID_Xz, *sendID_XY, *sendID_YZ, *sendID_XZ;
signed char *recvID_x, *recvID_y, *recvID_z, *recvID_X, *recvID_Y, *recvID_Z;
signed char *recvID_xy, *recvID_yz, *recvID_xz, *recvID_Xy, *recvID_Yz, *recvID_xZ;
signed char *recvID_xY, *recvID_yZ, *recvID_Xz, *recvID_XY, *recvID_YZ, *recvID_XZ;
// send buffers
sendID_x = new char [Dm->sendCount_x];
sendID_y = new char [Dm->sendCount_y];
sendID_z = new char [Dm->sendCount_z];
sendID_X = new char [Dm->sendCount_X];
sendID_Y = new char [Dm->sendCount_Y];
sendID_Z = new char [Dm->sendCount_Z];
sendID_xy = new char [Dm->sendCount_xy];
sendID_yz = new char [Dm->sendCount_yz];
sendID_xz = new char [Dm->sendCount_xz];
sendID_Xy = new char [Dm->sendCount_Xy];
sendID_Yz = new char [Dm->sendCount_Yz];
sendID_xZ = new char [Dm->sendCount_xZ];
sendID_xY = new char [Dm->sendCount_xY];
sendID_yZ = new char [Dm->sendCount_yZ];
sendID_Xz = new char [Dm->sendCount_Xz];
sendID_XY = new char [Dm->sendCount_XY];
sendID_YZ = new char [Dm->sendCount_YZ];
sendID_XZ = new char [Dm->sendCount_XZ];
sendID_x = new signed char [Dm->sendCount_x];
sendID_y = new signed char [Dm->sendCount_y];
sendID_z = new signed char [Dm->sendCount_z];
sendID_X = new signed char [Dm->sendCount_X];
sendID_Y = new signed char [Dm->sendCount_Y];
sendID_Z = new signed char [Dm->sendCount_Z];
sendID_xy = new signed char [Dm->sendCount_xy];
sendID_yz = new signed char [Dm->sendCount_yz];
sendID_xz = new signed char [Dm->sendCount_xz];
sendID_Xy = new signed char [Dm->sendCount_Xy];
sendID_Yz = new signed char [Dm->sendCount_Yz];
sendID_xZ = new signed char [Dm->sendCount_xZ];
sendID_xY = new signed char [Dm->sendCount_xY];
sendID_yZ = new signed char [Dm->sendCount_yZ];
sendID_Xz = new signed char [Dm->sendCount_Xz];
sendID_XY = new signed char [Dm->sendCount_XY];
sendID_YZ = new signed char [Dm->sendCount_YZ];
sendID_XZ = new signed char [Dm->sendCount_XZ];
//......................................................................................
// recv buffers
recvID_x = new char [Dm->recvCount_x];
recvID_y = new char [Dm->recvCount_y];
recvID_z = new char [Dm->recvCount_z];
recvID_X = new char [Dm->recvCount_X];
recvID_Y = new char [Dm->recvCount_Y];
recvID_Z = new char [Dm->recvCount_Z];
recvID_xy = new char [Dm->recvCount_xy];
recvID_yz = new char [Dm->recvCount_yz];
recvID_xz = new char [Dm->recvCount_xz];
recvID_Xy = new char [Dm->recvCount_Xy];
recvID_xZ = new char [Dm->recvCount_xZ];
recvID_xY = new char [Dm->recvCount_xY];
recvID_yZ = new char [Dm->recvCount_yZ];
recvID_Yz = new char [Dm->recvCount_Yz];
recvID_Xz = new char [Dm->recvCount_Xz];
recvID_XY = new char [Dm->recvCount_XY];
recvID_YZ = new char [Dm->recvCount_YZ];
recvID_XZ = new char [Dm->recvCount_XZ];
recvID_x = new signed char [Dm->recvCount_x];
recvID_y = new signed char [Dm->recvCount_y];
recvID_z = new signed char [Dm->recvCount_z];
recvID_X = new signed char [Dm->recvCount_X];
recvID_Y = new signed char [Dm->recvCount_Y];
recvID_Z = new signed char [Dm->recvCount_Z];
recvID_xy = new signed char [Dm->recvCount_xy];
recvID_yz = new signed char [Dm->recvCount_yz];
recvID_xz = new signed char [Dm->recvCount_xz];
recvID_Xy = new signed char [Dm->recvCount_Xy];
recvID_xZ = new signed char [Dm->recvCount_xZ];
recvID_xY = new signed char [Dm->recvCount_xY];
recvID_yZ = new signed char [Dm->recvCount_yZ];
recvID_Yz = new signed char [Dm->recvCount_Yz];
recvID_Xz = new signed char [Dm->recvCount_Xz];
recvID_XY = new signed char [Dm->recvCount_XY];
recvID_YZ = new signed char [Dm->recvCount_YZ];
recvID_XZ = new signed char [Dm->recvCount_XZ];
//......................................................................................
int sendtag,recvtag;
sendtag = recvtag = 7;
@@ -754,12 +754,20 @@ double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
if (morph_delta / morph_delta_previous > 2.0 ) morph_delta = morph_delta_previous*2.0;
//MAX_DISPLACEMENT *= max(TargetGrowth/GrowthEstimate,1.25);
if (MAX_DISPLACEMENT > 1.0 ){
if (morph_delta > 0.0 ) morph_delta = 1.0;
else morph_delta = -1.0;
//if (COUNT_FOR_LOOP > 2) COUNT_FOR_LOOP = 100;
COUNT_FOR_LOOP = 100; // exit loop if displacement is too large
if (morph_delta > 0.0 ){
// object is growing
if (MAX_DISPLACEMENT > 3.0 ){
morph_delta = 3.0;
COUNT_FOR_LOOP = 100; // exit loop if displacement is too large
}
}
else{
// object is shrinking
if (MAX_DISPLACEMENT > 1.0 ){
morph_delta = -1.0;
COUNT_FOR_LOOP = 100; // exit loop if displacement is too large
}
}
}
if (rank == 0) printf("Final delta=%f \n",morph_delta);

6
analysis/morphology.h
View File

@@ -3,6 +3,6 @@
#include "common/Domain.h"
#include "analysis/runAnalysis.h"
double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol);
double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr<Domain> Dm, double VoidFraction);
double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id, std::shared_ptr<Domain> Dm, double TargetVol);

68
analysis/runAnalysis.cpp
View File

@@ -404,9 +404,9 @@ public:
// Averages.ColorToSignedDistance(beta,Averages.Phase,Averages.Phase_tplus);
}
if ( matches(type,AnalysisType::ComputeAverages) ) {
PROFILE_START("Compute subphase",1);
PROFILE_START("Compute basic averages",1);
Averages.Basic();
PROFILE_STOP("Compute subphase",1);
PROFILE_STOP("Compute basic averages",1);
}
}
private:
@@ -417,6 +417,28 @@ private:
double beta;
};
class SubphaseWorkItem: public ThreadPool::WorkItemRet<void>
{
public:
SubphaseWorkItem( AnalysisType type_, int timestep_, SubPhase& Averages_ ):
type(type_), timestep(timestep_), Averages(Averages_){ }
~SubphaseWorkItem() { }
virtual void run() {
PROFILE_START("Compute subphase",1);
Averages.Full();
Averages.Write(timestep);
PROFILE_STOP("Compute subphase",1);
}
private:
SubphaseWorkItem();
AnalysisType type;
int timestep;
SubPhase& Averages;
double beta;
};
/******************************************************************
* MPI comm wrapper for use with analysis *
@@ -483,16 +505,29 @@ runAnalysis::runAnalysis( std::shared_ptr<Database> db,
ThreadPool::thread_id_t d_wait_analysis;
ThreadPool::thread_id_t d_wait_vis;
ThreadPool::thread_id_t d_wait_restart;
ThreadPool::thread_id_t d_wait_subphase;
char rankString[20];
sprintf(rankString,"%05d",Dm->rank());
d_N[0] = Dm->Nx;
d_N[1] = Dm->Ny;
d_N[2] = Dm->Nz;
d_restart_interval = db->getScalar<int>( "restart_interval" );
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
d_blobid_interval = db->getScalar<int>( "blobid_interval" );
d_visualization_interval = db->getScalar<int>( "visualization_interval" );
d_analysis_interval = db->getScalar<int>( "analysis_interval" );
d_subphase_analysis_interval = INT_MAX;
d_visualization_interval = INT_MAX;
d_blobid_interval = INT_MAX;
if (db->keyExists( "blobid_interval" )){
d_blobid_interval = db->getScalar<int>( "blobid_interval" );
}
if (db->keyExists( "visualization_interval" )){
d_visualization_interval = db->getScalar<int>( "visualization_interval" );
}
if (db->keyExists( "subphase_analysis_interval" )){
d_subphase_analysis_interval = db->getScalar<int>( "subphase_analysis_interval" );
}
auto restart_file = db->getScalar<std::string>( "restart_file" );
d_restartFile = restart_file + "." + rankString;
d_rank = MPI_WORLD_RANK();
@@ -579,6 +614,7 @@ void runAnalysis::finish( )
d_wait_blobID.reset();
d_wait_analysis.reset();
d_wait_vis.reset();
d_wait_subphase.reset();
d_wait_restart.reset();
// Syncronize
MPI_Barrier( d_comm );
@@ -881,6 +917,7 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
//if ( matches(type,AnalysisType::CopySimState) ) {
if ( timestep%d_analysis_interval == 0 ) {
finish(); // can't copy if threads are still working on data
// Copy the members of Averages to the cpu (phase was copied above)
PROFILE_START("Copy-Pressure",1);
ScaLBL_D3Q19_Pressure(fq,Pressure,d_Np);
@@ -907,9 +944,19 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
// if ( matches(type,AnalysisType::ComputeAverages) ) {
if ( timestep%d_analysis_interval == 0 ) {
auto work = new BasicWorkItem(type,timestep,Averages);
work->add_dependency(d_wait_analysis); // Make sure we are done using analysis before modifying
work->add_dependency(d_wait_subphase); // Make sure we are done using analysis before modifying
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_vis);
d_wait_analysis = d_tpool.add_work(work);
}
if ( timestep%d_subphase_analysis_interval == 0 ) {
auto work = new SubphaseWorkItem(type,timestep,Averages);
work->add_dependency(d_wait_subphase); // Make sure we are done using analysis before modifying
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_vis);
d_wait_subphase = d_tpool.add_work(work);
}
if (timestep%d_restart_interval==0){
std::shared_ptr<double> cfq,cDen;
@@ -930,6 +977,15 @@ void runAnalysis::basic( int timestep, SubPhase &Averages, const double *Phi, do
d_wait_restart = d_tpool.add_work(work1);
}
if (timestep%d_visualization_interval==0){
// Write the vis files
auto work = new IOWorkItem( timestep, d_meshData, Averages, d_fillData, getComm() );
work->add_dependency(d_wait_analysis);
work->add_dependency(d_wait_subphase);
work->add_dependency(d_wait_vis);
d_wait_vis = d_tpool.add_work(work);
}
PROFILE_STOP("run");
}

6
analysis/runAnalysis.h
View File

@@ -22,7 +22,7 @@
#include "common/Communication.h"
#include "common/ScaLBL.h"
#include "threadpool/thread_pool.h"
#include <limits.h>
typedef std::shared_ptr<std::pair<int,IntArray>> BlobIDstruct;
typedef std::shared_ptr<std::vector<BlobIDType>> BlobIDList;
@@ -30,7 +30,7 @@ typedef std::shared_ptr<std::vector<BlobIDType>> BlobIDList;
// Types of analysis
enum class AnalysisType : uint64_t { AnalyzeNone=0, IdentifyBlobs=0x01, CopyPhaseIndicator=0x02,
CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20 };
CopySimState=0x04, ComputeAverages=0x08, CreateRestart=0x10, WriteVis=0x20, ComputeSubphase=0x40 };
//! Class to run the analysis in multiple threads
@@ -103,6 +103,7 @@ private:
int d_Np;
int d_rank;
int d_restart_interval, d_analysis_interval, d_blobid_interval, d_visualization_interval;
int d_subphase_analysis_interval;
double d_beta;
bool d_regular;
ThreadPool d_tpool;
@@ -122,6 +123,7 @@ private:
// 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_subphase;
ThreadPool::thread_id_t d_wait_vis;
ThreadPool::thread_id_t d_wait_restart;

91
common/Domain.cpp
View File

@@ -92,6 +92,7 @@ Domain::Domain( std::shared_ptr<Database> db, MPI_Comm Communicator):
Lx(0), Ly(0), Lz(0), Volume(0), BoundaryCondition(0),
Comm(MPI_COMM_NULL),
inlet_layers_x(0), inlet_layers_y(0), inlet_layers_z(0),
outlet_layers_x(0), outlet_layers_y(0), outlet_layers_z(0),
sendCount_x(0), sendCount_y(0), sendCount_z(0), sendCount_X(0), sendCount_Y(0), sendCount_Z(0),
sendCount_xy(0), sendCount_yz(0), sendCount_xz(0), sendCount_Xy(0), sendCount_Yz(0), sendCount_xZ(0),
sendCount_xY(0), sendCount_yZ(0), sendCount_Xz(0), sendCount_XY(0), sendCount_YZ(0), sendCount_XZ(0),
@@ -141,6 +142,12 @@ void Domain::initialize( std::shared_ptr<Database> db )
inlet_layers_y = InletCount[1];
inlet_layers_z = InletCount[2];
}
if (d_db->keyExists( "OutletLayers" )){
auto OutletCount = d_db->getVector<int>( "OutletLayers" );
outlet_layers_x = OutletCount[0];
outlet_layers_y = OutletCount[1];
outlet_layers_z = OutletCount[2];
}
ASSERT( n.size() == 3u );
ASSERT( L.size() == 3u );
@@ -159,13 +166,17 @@ void Domain::initialize( std::shared_ptr<Database> db )
MPI_Comm_rank( Comm, &myrank );
rank_info = RankInfoStruct(myrank,nproc[0],nproc[1],nproc[2]);
// inlet layers only apply to lower part of domain
if (nproc[0] > 0) inlet_layers_x = 0;
if (nproc[1] > 0) inlet_layers_y = 0;
if (nproc[2] > 0) inlet_layers_z = 0;
if (rank_info.ix > 0) inlet_layers_x = 0;
if (rank_info.jy > 0) inlet_layers_y = 0;
if (rank_info.kz > 0) inlet_layers_z = 0;
// outlet layers only apply to top part of domain
if (rank_info.ix < nproc[0]-1 ) outlet_layers_x = 0;
if (rank_info.jy < nproc[1]-1) outlet_layers_y = 0;
if (rank_info.kz < nproc[2]-1) outlet_layers_z = 0;
// Fill remaining variables
N = Nx*Ny*Nz;
Volume = nx*ny*nx*nproc[0]*nproc[1]*nproc[2]*1.0;
id = new char[N];
id = new signed char[N];
memset(id,0,N);
BoundaryCondition = d_db->getScalar<int>("BC");
int nprocs;
@@ -529,43 +540,45 @@ void Domain::CommInit()
void Domain::ReadIDs(){
// Read the IDs from input file
int nprocs=nprocx()*nprocy()*nprocz();
size_t readID;
char LocalRankString[8];
char LocalRankFilename[40];
//.......................................................................
if (rank() == 0) printf("Read input media... \n");
//.......................................................................
sprintf(LocalRankString,"%05d",rank());
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// .......... READ THE INPUT FILE .......................................
if (rank()==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank());
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Domain::ReadIDs -- Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
fclose(IDFILE);
// Compute the porosity
double sum;
double porosity;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
//.........................................................
// If external boundary conditions are applied remove solid
if (BoundaryCondition > 0 && kproc() == 0){
for (int k=0; k<3; k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
}
}
int nprocs=nprocx()*nprocy()*nprocz();
size_t readID;
char LocalRankString[8];
char LocalRankFilename[40];
//.......................................................................
if (rank() == 0) printf("Read input media... \n");
//.......................................................................
sprintf(LocalRankString,"%05d",rank());
sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString);
// .......... READ THE INPUT FILE .......................................
if (rank()==0) printf("Initialize from segmented data: solid=0, NWP=1, WP=2 \n");
sprintf(LocalRankFilename,"ID.%05i",rank());
FILE *IDFILE = fopen(LocalRankFilename,"rb");
if (IDFILE==NULL) ERROR("Domain::ReadIDs -- Error opening file: ID.xxxxx");
readID=fread(id,1,N,IDFILE);
if (readID != size_t(N)) printf("Domain::ReadIDs -- Error reading ID (rank=%i) \n",rank());
fclose(IDFILE);
// Compute the porosity
double sum;
double sum_local=0.0;
double iVol_global = 1.0/(1.0*(Nx-2)*(Ny-2)*(Nz-2)*nprocs);
if (BoundaryCondition > 0) iVol_global = 1.0/(1.0*(Nx-2)*nprocx()*(Ny-2)*nprocy()*((Nz-2)*nprocz()-6));
//.........................................................
// If external boundary conditions are applied remove solid
if (BoundaryCondition > 0 && kproc() == 0){
if (inlet_layers_z < 4) inlet_layers_z=4;
for (int k=0; k<inlet_layers_z; k++){
for (int j=0;j<Ny;j++){
for (int i=0;i<Nx;i++){
int n = k*Nx*Ny+j*Nx+i;
id[n] = 1;
}
}
}
}
if (BoundaryCondition > 0 && kproc() == nprocz()-1){
for (int k=Nz-3; k<Nz; k++){
if (outlet_layers_z < 4) outlet_layers_z=4;
for (int k=Nz-outlet_layers_z; 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;
@@ -574,7 +587,7 @@ void Domain::ReadIDs(){
}
}
}
for (int k=1;k<Nz-1;k++){
for (int k=inlet_layers_z+1; k<Nz-outlet_layers_z-1;k++){
for (int j=1;j<Ny-1;j++){
for (int i=1;i<Nx-1;i++){
int n = k*Nx*Ny+j*Nx+i;

9
common/Domain.h
View File

@@ -125,6 +125,8 @@ public: // Public variables (need to create accessors instead)
double Lx,Ly,Lz,Volume;
int Nx,Ny,Nz,N;
int inlet_layers_x, inlet_layers_y, inlet_layers_z;
int outlet_layers_x, outlet_layers_y, outlet_layers_z;
double porosity;
RankInfoStruct rank_info;
MPI_Comm Comm; // MPI Communicator for this domain
@@ -136,6 +138,7 @@ public: // Public variables (need to create accessors instead)
//**********************************
// MPI ranks for all 18 neighbors
//**********************************
inline double Porosity() const { return porosity; }
inline int iproc() const { return rank_info.ix; }
inline int jproc() const { return rank_info.jy; }
inline int kproc() const { return rank_info.kz; }
@@ -184,7 +187,7 @@ public: // Public variables (need to create accessors instead)
int *recvList_xY, *recvList_yZ, *recvList_Xz, *recvList_XY, *recvList_YZ, *recvList_XZ;
//......................................................................................
// Solid indicator function
char *id;
signed char *id;
void ReadIDs();
void CommunicateMeshHalo(DoubleArray &Mesh);
@@ -193,8 +196,8 @@ public: // Public variables (need to create accessors instead)
private:
void PackID(int *list, int count, char *sendbuf, char *ID);
void UnpackID(int *list, int count, char *recvbuf, char *ID);
void PackID(int *list, int count, signed char *sendbuf, signed char *ID);
void UnpackID(int *list, int count, signed char *recvbuf, signed char *ID);
void CommHaloIDs();
//......................................................................................

2
common/ScaLBL.cpp
View File

@@ -367,7 +367,7 @@ void ScaLBL_Communicator::D3Q19_MapRecv(int Cqx, int Cqy, int Cqz, int *list, i
delete [] ReturnDist;
}
int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, char *id, int Np){
int ScaLBL_Communicator::MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np){
/*
* Generate a memory optimized layout
* id[n] == 0 implies that site n should be ignored (treat as a mask)

2
common/ScaLBL.h
View File

@@ -173,7 +173,7 @@ public:
int FirstInterior();
int LastInterior();
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, char *id, int Np);
int MemoryOptimizedLayoutAA(IntArray &Map, int *neighborList, signed char *id, int Np);
// void MemoryOptimizedLayout(IntArray &Map, int *neighborList, char *id, int Np);
// void MemoryOptimizedLayoutFull(IntArray &Map, int *neighborList, char *id, int Np);
// void MemoryDenseLayout(IntArray &Map, int *neighborList, char *id, int Np);

121
models/ColorModel.cpp
View File

@@ -25,7 +25,7 @@ rank(RANK), nprocs(NP), Restart(0),timestep(0),timestepMax(0),tauA(0),tauB(0),rh
Fx(0),Fy(0),Fz(0),flux(0),din(0),dout(0),inletA(0),inletB(0),outletA(0),outletB(0),
Nx(0),Ny(0),Nz(0),N(0),Np(0),nprocx(0),nprocy(0),nprocz(0),BoundaryCondition(0),Lx(0),Ly(0),Lz(0),comm(COMM)
{
REVERSE_FLOW_DIRECTION = false;
}
ScaLBL_ColorModel::~ScaLBL_ColorModel(){
@@ -119,7 +119,7 @@ void ScaLBL_ColorModel::SetDomain(){
Mask = std::shared_ptr<Domain>(new Domain(domain_db,comm)); // mask domain removes immobile phases
Nx+=2; Ny+=2; Nz += 2;
N = Nx*Ny*Nz;
id = new char [N];
id = new signed char [N];
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
Averages = std::shared_ptr<SubPhase> ( new SubPhase(Dm) ); // TwoPhase analysis object
@@ -175,10 +175,10 @@ void ScaLBL_ColorModel::ReadInput(){
void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
{
size_t NLABELS=0;
char VALUE=0;
signed char VALUE=0;
double AFFINITY=0.f;
auto LabelList = color_db->getVector<char>( "ComponentLabels" );
auto LabelList = color_db->getVector<int>( "ComponentLabels" );
auto AffinityList = color_db->getVector<double>( "ComponentAffinity" );
NLABELS=LabelList.size();
@@ -220,12 +220,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase)
for (int idx=0; idx<NLABELS; idx++) label_count_global[idx]=sumReduce( Dm->Comm, label_count[idx]);
if (rank==0){
printf("Components labels: %lu \n",NLABELS);
printf("Component labels: %lu \n",NLABELS);
for (unsigned int idx=0; idx<NLABELS; idx++){
VALUE=LabelList[idx];
AFFINITY=AffinityList[idx];
double volume_fraction = double(label_count_global[idx])/double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
printf(" label=%hhd, affinity=%f, volume fraction==%f\n",VALUE,AFFINITY,volume_fraction);
printf(" label=%d, affinity=%f, volume fraction==%f\n",VALUE,AFFINITY,volume_fraction);
}
}
@@ -418,7 +418,7 @@ void ScaLBL_ColorModel::Initialize(){
ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3);
}
}
ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double));
}
void ScaLBL_ColorModel::Run(){
@@ -430,6 +430,7 @@ void ScaLBL_ColorModel::Run(){
bool USE_MORPH = false;
int analysis_interval = 1000; // number of timesteps in between in situ analysis
int MAX_MORPH_TIMESTEPS = 50000; // maximum number of LBM timesteps to spend in morphological adaptation routine
int MAX_STEADY_TIMESTEPS = 200000;
int RAMP_TIMESTEPS = 0;//50000; // number of timesteps to run initially (to get a reasonable velocity field before other pieces kick in)
int morph_interval = 1000000;
int CURRENT_MORPH_TIMESTEPS=0; // counter for number of timesteps spent in morphological adaptation routine (reset each time)
@@ -467,7 +468,12 @@ void ScaLBL_ColorModel::Run(){
if (analysis_db->keyExists( "analysis_interval" )){
analysis_interval = analysis_db->getScalar<int>( "analysis_interval" );
}
if (analysis_db->keyExists( "max_steady_timesteps" )){
MAX_STEADY_TIMESTEPS = analysis_db->getScalar<int>( "max_steady_timesteps" );
}
if (analysis_db->keyExists( "max_morph_timesteps" )){
MAX_MORPH_TIMESTEPS = analysis_db->getScalar<int>( "max_morph_timesteps" );
}
if (rank==0){
printf("********************************************************\n");
printf("No. of timesteps: %i \n", timestepMax);
@@ -569,6 +575,10 @@ void ScaLBL_ColorModel::Run(){
//analysis.run( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
analysis.basic( timestep, *Averages, Phi, Pressure, Velocity, fq, Den );
if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition > 0){
printf("....inlet pressure=%f \n",din);
}
// allow initial ramp-up to get closer to steady state
if (timestep > RAMP_TIMESTEPS && timestep%analysis_interval == 0 && USE_MORPH){
analysis.finish();
@@ -586,14 +596,18 @@ void ScaLBL_ColorModel::Run(){
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 force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz);
double dir_x = Fx/force_mag;
double dir_y = Fy/force_mag;
double dir_z = Fz/force_mag;
double flow_rate_A = (vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
double flow_rate_B = (vB_x*dir_x + vB_y*dir_y + vB_z*dir_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-2)*(Ny-2)*(Nz-2)*nprocs));
double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz);
if (fabs((Ca - Ca_previous)/Ca) < tolerance ){
if (fabs((Ca - Ca_previous)/Ca) < tolerance || CURRENT_STEADY_TIMESTEPS > MAX_STEADY_TIMESTEPS){
MORPH_ADAPT = true;
CURRENT_MORPH_TIMESTEPS=0;
delta_volume_target = (volA + volB)*morph_delta; // set target volume change
@@ -636,11 +650,15 @@ void ScaLBL_ColorModel::Run(){
// flow reversal criteria based on fractional flow rate
if (delta_volume_target < 0.0 &&
volA*flow_rate_A/(volA*flow_rate_A+volB*flow_rate_B) < RESIDUAL_ENDPOINT_THRESHOLD){
delta_volume_target *= (-1.0);
REVERSE_FLOW_DIRECTION = true;
}
else if (delta_volume_target > 0.0 &&
volB*flow_rate_B/(volA*flow_rate_A+volB*flow_rate_B) < RESIDUAL_ENDPOINT_THRESHOLD){
REVERSE_FLOW_DIRECTION = true;
}
if ( REVERSE_FLOW_DIRECTION ){
delta_volume_target *= (-1.0);
REVERSE_FLOW_DIRECTION = false;
}
}
else{
@@ -663,9 +681,18 @@ void ScaLBL_ColorModel::Run(){
CURRENT_STEADY_TIMESTEPS=0;
}
else if (CURRENT_MORPH_TIMESTEPS > MAX_MORPH_TIMESTEPS) {
delta_volume = 0.0;
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
}
if ( REVERSE_FLOW_DIRECTION ){
if (rank==0) printf("*****REVERSE FLOW DIRECTION***** \n");
delta_volume = 0.0;
// flow direction will reverse after next steady point
MORPH_ADAPT = false;
CURRENT_STEADY_TIMESTEPS=0;
}
MPI_Barrier(comm);
}
morph_timesteps += analysis_interval;
@@ -774,47 +801,47 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
}
}
if (volume_connected < 0.025*volume_initial){
// if connected volume is less than 2.5% just delete the whole thing
if (rank==0) printf("Connected region has shrunk to less than 2.5% of total fluid volume (remove the whole thing) \n");
if (volume_connected < 0.02*volume_initial && target_delta_volume < 0.0){
// if connected volume is less than 2% just delete the whole thing
if (rank==0) printf("Connected region has shrunk to less than 2 %% of total fluid volume \n");
REVERSE_FLOW_DIRECTION = true;
}
else {
if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
double target_delta_volume_incremental = target_delta_volume;
if (fabs(target_delta_volume) > 0.01*volume_initial)
target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental);
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) < 0.0 ) phase_id(i,j,k) = 0;
else phase_id(i,j,k) = 1;
//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
}
}
}
if (rank==0) printf("MorphGrow with target volume fraction change %f \n", target_delta_volume/volume_initial);
double target_delta_volume_incremental = target_delta_volume;
if (fabs(target_delta_volume) > 0.01*volume_initial)
target_delta_volume_incremental = 0.01*volume_initial*target_delta_volume/fabs(target_delta_volume);
delta_volume = MorphGrow(Averages->SDs,phase_distance,phase_id,Averages->Dm, target_delta_volume_incremental);
CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
// 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){
if (d < 3.f){
//phase(i,j,k) = -1.0;
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
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) < 0.0 ) phase_id(i,j,k) = 0;
else phase_id(i,j,k) = 1;
//if (phase_distance(i,j,k) < 0.0 ) phase(i,j,k) = 1.0;
}
}
fillDouble.fill(phase);
}
CalcDist(phase_distance,phase_id,*Dm); // re-calculate distance
// 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){
if (d < 3.f){
//phase(i,j,k) = -1.0;
phase(i,j,k) = (2.f*(exp(-2.f*beta*d))/(1.f+exp(-2.f*beta*d))-1.f);
}
}
}
}
}
fillDouble.fill(phase);
count = 0.f;
for (int k=1; k<Nz-1; k++){

3
models/ColorModel.h
View File

@@ -47,6 +47,7 @@ public:
void WriteDebug();
bool Restart,pBC;
bool REVERSE_FLOW_DIRECTION;
int timestep,timestepMax;
int BoundaryCondition;
double tauA,tauB,rhoA,rhoB,alpha,beta;
@@ -71,7 +72,7 @@ public:
std::shared_ptr<Database> analysis_db;
IntArray Map;
char *id;
signed char *id;
int *NeighborList;
int *dvcMap;
double *fq, *Aq, *Bq;

4
models/MRTModel.cpp
View File

@@ -256,8 +256,8 @@ void ScaLBL_MRTModel::Run(){
Hs=sumReduce( Dm->Comm, Hs);
Xs=sumReduce( Dm->Comm, Xs);
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);
double h = Lz/double((Nz-2)*nprocz);
double absperm = h*h*mu*Mask->Porosity()*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,

4
tests/CMakeLists.txt
View File

@@ -12,9 +12,9 @@ ADD_LBPM_EXECUTABLE( lbpm_refine_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphdrain_pp )
ADD_LBPM_EXECUTABLE( lbpm_morphopen_pp )
#ADD_LBPM_EXECUTABLE( lbpm_morph_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
#ADD_LBPM_EXECUTABLE( lbpm_segmented_pp )
#ADD_LBPM_EXECUTABLE( lbpm_block_pp )
ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )
#ADD_LBPM_EXECUTABLE( lbpm_segmented_decomp )
ADD_LBPM_EXECUTABLE( lbpm_serial_decomp )
#ADD_LBPM_EXECUTABLE( lbpm_disc_pp )
#ADD_LBPM_EXECUTABLE( lbpm_juanes_bench_disc_pp )

35
tests/TestColorBubble.cpp
View File

@@ -14,13 +14,13 @@ using namespace std;
inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius){
// initialize a bubble
int i,j,k,n;
int rank = ColorModel.Mask->rank();
int nprocx = ColorModel.Mask->nprocx();
int nprocy = ColorModel.Mask->nprocy();
int nprocz = ColorModel.Mask->nprocz();
int Nx = ColorModel.Mask->Nx;
int Ny = ColorModel.Mask->Ny;
int Nz = ColorModel.Mask->Nz;
int rank = ColorModel.Dm->rank();
int nprocx = ColorModel.Dm->nprocx();
int nprocy = ColorModel.Dm->nprocy();
int nprocz = ColorModel.Dm->nprocz();
int Nx = ColorModel.Dm->Nx;
int Ny = ColorModel.Dm->Ny;
int Nz = ColorModel.Dm->Nz;
if (rank == 0) cout << "Setting up bubble..." << endl;
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
@@ -34,14 +34,12 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
for (k=0;k<Nz;k++){
for (j=0;j<Ny;j++){
for (i=0;i<Nx;i++){
n = k*Nx*Ny + j*Nz + i;
int iglobal= i+(Nx-2)*ColorModel.Mask->iproc();
int jglobal= j+(Ny-2)*ColorModel.Mask->jproc();
int kglobal= k+(Nz-2)*ColorModel.Mask->kproc();
n = k*Nx*Ny + j*Nx + i;
double iglobal= double(i+(Nx-2)*ColorModel.Dm->iproc())-double((Nx-2)*nprocx)*0.5;
double jglobal= double(j+(Ny-2)*ColorModel.Dm->jproc())-double((Ny-2)*nprocy)*0.5;
double kglobal= double(k+(Nz-2)*ColorModel.Dm->kproc())-double((Nz-2)*nprocz)*0.5;
// Initialize phase position field for parallel bubble test
if ((iglobal-0.5*(Nx-2)*nprocx)*(iglobal-0.5*(Nx-2)*nprocx)
+(jglobal-0.5*(Ny-2)*nprocy)*(jglobal-0.5*(Ny-2)*nprocy)
+(kglobal-0.5*(Nz-2)*nprocz)*(kglobal-0.5*(Nz-2)*nprocz) < BubbleRadius*BubbleRadius){
if ((iglobal*iglobal)+(jglobal*jglobal)+(kglobal*kglobal) < BubbleRadius*BubbleRadius){
ColorModel.Mask->id[n] = 2;
ColorModel.Mask->id[n] = 2;
}
@@ -49,9 +47,18 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
ColorModel.Mask->id[n]=1;
ColorModel.Mask->id[n]=1;
}
ColorModel.id[n] = ColorModel.Mask->id[n];
ColorModel.Dm->id[n] = ColorModel.Mask->id[n];
}
}
}
FILE *OUTFILE;
char LocalRankFilename[40];
sprintf(LocalRankFilename,"Bubble.%05i.raw",rank);
OUTFILE = fopen(LocalRankFilename,"wb");
fwrite(ColorModel.id,1,Nx*Ny*Nz,OUTFILE);
fclose(OUTFILE);
// initialize the phase indicator field
}
//***************************************************************************************

22
tests/lbpm_morphdrain_pp.cpp
View File

@@ -59,8 +59,8 @@ int main(int argc, char **argv)
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
SW = domain_db->getScalar<double>("Sw");
// Generate the NWP configuration
@@ -82,8 +82,8 @@ int main(int argc, char **argv)
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit();
char *id;
id = new char [N];
signed char *id;
id = new signed char [N];
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
@@ -131,14 +131,14 @@ int main(int argc, char **argv)
// calculate distance to non-wetting fluid
if (domain_db->keyExists( "HistoryLabels" )){
if (rank==0) printf("Relabel solid components that touch fluid 1 \n");
auto LabelList = domain_db->getVector<char>( "ComponentLabels" );
auto HistoryLabels = domain_db->getVector<char>( "HistoryLabels" );
auto LabelList = domain_db->getVector<int>( "ComponentLabels" );
auto HistoryLabels = domain_db->getVector<int>( "HistoryLabels" );
size_t NLABELS=LabelList.size();
if (rank==0){
for (unsigned int idx=0; idx < NLABELS; idx++){
char VALUE = LabelList[idx];
char NEWVAL = HistoryLabels[idx];
printf(" Relabel component %d as %d \n", VALUE, NEWVAL);
signed char VALUE = LabelList[idx];
signed char NEWVAL = HistoryLabels[idx];
printf(" Relabel component %hhd as %hhd \n", VALUE, NEWVAL);
}
}
for (int k=0;k<nz;k++){
@@ -169,8 +169,8 @@ int main(int argc, char **argv)
int n = k*nx*ny+j*nx+i;
signed char LOCVAL = id[n];
for (unsigned int idx=0; idx < NLABELS; idx++){
char VALUE=LabelList[idx];
char NEWVALUE=HistoryLabels[idx];
signed char VALUE=LabelList[idx];
signed char NEWVALUE=HistoryLabels[idx];
if (LOCVAL == VALUE){
idx = NLABELS;
if (SignDist(i,j,k) < 1.0){

20
tests/lbpm_morphopen_pp.cpp
View File

@@ -59,8 +59,8 @@ int main(int argc, char **argv)
auto L = domain_db->getVector<double>( "L" );
auto size = domain_db->getVector<int>( "n" );
auto nproc = domain_db->getVector<int>( "nproc" );
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
SW = domain_db->getScalar<double>("Sw");
// Generate the NWP configuration
@@ -82,8 +82,8 @@ int main(int argc, char **argv)
for (n=0; n<N; n++) Dm->id[n]=1;
Dm->CommInit();
char *id;
id = new char [N];
signed char *id;
id = new signed char [N];
sprintf(LocalRankFilename,"ID.%05i",rank);
size_t readID;
FILE *IDFILE = fopen(LocalRankFilename,"rb");
@@ -131,13 +131,13 @@ int main(int argc, char **argv)
// calculate distance to non-wetting fluid
if (domain_db->keyExists( "HistoryLabels" )){
if (rank==0) printf("Relabel solid components that touch fluid 1 \n");
auto LabelList = domain_db->getVector<char>( "ComponentLabels" );
auto HistoryLabels = domain_db->getVector<char>( "HistoryLabels" );
auto LabelList = domain_db->getVector<int>( "ComponentLabels" );
auto HistoryLabels = domain_db->getVector<int>( "HistoryLabels" );
size_t NLABELS=LabelList.size();
if (rank==0){
for (unsigned int idx=0; idx < NLABELS; idx++){
char VALUE = LabelList[idx];
char NEWVAL = HistoryLabels[idx];
signed char VALUE = LabelList[idx];
signed char NEWVAL = HistoryLabels[idx];
printf(" Relabel component %d as %d \n", VALUE, NEWVAL);
}
}
@@ -169,8 +169,8 @@ int main(int argc, char **argv)
int n = k*nx*ny+j*nx+i;
signed char LOCVAL = id[n];
for (unsigned int idx=0; idx < NLABELS; idx++){
char VALUE=LabelList[idx];
char NEWVALUE=HistoryLabels[idx];
signed char VALUE=LabelList[idx];
signed char NEWVALUE=HistoryLabels[idx];
if (LOCVAL == VALUE){
idx = NLABELS;
if (SignDist(i,j,k) < 1.0){

41
tests/lbpm_serial_decomp.cpp
View File

@@ -87,8 +87,11 @@ int main(int argc, char **argv)
if (domain_db->keyExists( "checkerSize" )){
checkerSize = domain_db->getScalar<int>( "checkerSize" );
}
auto ReadValues = domain_db->getVector<char>( "ReadValues" );
auto WriteValues = domain_db->getVector<char>( "WriteValues" );
else {
checkerSize = SIZE[0];
}
auto ReadValues = domain_db->getVector<int>( "ReadValues" );
auto WriteValues = domain_db->getVector<int>( "WriteValues" );
auto ReadType = domain_db->getScalar<std::string>( "ReadType" );
if (ReadType == "8bit"){
}
@@ -112,8 +115,8 @@ int main(int argc, char **argv)
printf("Input media: %s\n",Filename.c_str());
printf("Relabeling %lu values\n",ReadValues.size());
for (int idx=0; idx<ReadValues.size(); idx++){
char oldvalue=ReadValues[idx];
char newvalue=WriteValues[idx];
int oldvalue=ReadValues[idx];
int newvalue=WriteValues[idx];
printf("oldvalue=%d, newvalue =%d \n",oldvalue,newvalue);
}
@@ -158,13 +161,13 @@ int main(int argc, char **argv)
for (int j = 0; j<Ny; j++){
for (int i = xStart; i < xStart+inlet_count_x; i++){
if ( (j/checkerSize + k/checkerSize)%2 == 0){
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
else{
// void checkers
SegData[k*Nx*Ny+j*Nx+i] = 2;
}
else{
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
}
}
}
@@ -177,13 +180,13 @@ int main(int argc, char **argv)
for (int j = yStart; i < yStart+inlet_count_y; j++){
for (int i = 0; i<Nx; i++){
if ( (i/checkerSize + k/checkerSize)%2 == 0){
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
else{
// void checkers
SegData[k*Nx*Ny+j*Nx+i] = 2;
}
else{
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
}
}
}
@@ -196,13 +199,13 @@ int main(int argc, char **argv)
for (int j = 0; j<Ny; j++){
for (int i = 0; i<Nx; i++){
if ( (i/checkerSize+j/checkerSize)%2 == 0){
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
else{
// void checkers
SegData[k*Nx*Ny+j*Nx+i] = 2;
}
else{
// solid checkers
SegData[k*Nx*Ny+j*Nx+i] = 0;
}
}
}
}
@@ -259,8 +262,8 @@ int main(int argc, char **argv)
n = k*(nx+2)*(ny+2) + j*(nx+2) + i;;
char locval = loc_id[n];
for (int idx=0; idx<ReadValues.size(); idx++){
char oldvalue=ReadValues[idx];
char newvalue=WriteValues[idx];
signed char oldvalue=ReadValues[idx];
signed char newvalue=WriteValues[idx];
if (locval == oldvalue){
loc_id[n] = newvalue;
LabelCount[idx]++;
@@ -285,7 +288,7 @@ int main(int argc, char **argv)
}
}
for (int idx=0; idx<ReadValues.size(); idx++){
char label=ReadValues[idx];
int label=ReadValues[idx];
int count=LabelCount[idx];
printf("Label=%d, Count=%d \n",label,count);
}