diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 9cc977e5..4890a8d3 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -28,23 +28,32 @@ SubPhase::SubPhase(std::shared_ptr 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 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 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; koutlet_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 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; kComm, 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; diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h index ac5cbf76..738d2188 100644 --- a/analysis/SubPhase.h +++ b/analysis/SubPhase.h @@ -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 diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index f86e2de8..57ead57e 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -1,7 +1,7 @@ #include // 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 Dm, double VoidFraction){ +double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr 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 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 Dm, double VoidFraction){ +double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr 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 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 &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); diff --git a/analysis/morphology.h b/analysis/morphology.h index 4692475d..f69af5f2 100644 --- a/analysis/morphology.h +++ b/analysis/morphology.h @@ -3,6 +3,6 @@ #include "common/Domain.h" #include "analysis/runAnalysis.h" -double MorphOpen(DoubleArray &SignDist, char *id, std::shared_ptr Dm, double VoidFraction); -double MorphDrain(DoubleArray &SignDist, char *id, std::shared_ptr Dm, double VoidFraction); -double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array &id, std::shared_ptr Dm, double TargetVol); \ No newline at end of file +double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr Dm, double VoidFraction); +double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr Dm, double VoidFraction); +double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array &id, std::shared_ptr Dm, double TargetVol); diff --git a/analysis/runAnalysis.cpp b/analysis/runAnalysis.cpp index 89c15e4b..4697aa20 100644 --- a/analysis/runAnalysis.cpp +++ b/analysis/runAnalysis.cpp @@ -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 +{ +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 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( "restart_interval" ); - d_analysis_interval = db->getScalar( "analysis_interval" ); - d_blobid_interval = db->getScalar( "blobid_interval" ); - d_visualization_interval = db->getScalar( "visualization_interval" ); + d_analysis_interval = db->getScalar( "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( "blobid_interval" ); + } + if (db->keyExists( "visualization_interval" )){ + d_visualization_interval = db->getScalar( "visualization_interval" ); + } + if (db->keyExists( "subphase_analysis_interval" )){ + d_subphase_analysis_interval = db->getScalar( "subphase_analysis_interval" ); + } + auto restart_file = db->getScalar( "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 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"); } diff --git a/analysis/runAnalysis.h b/analysis/runAnalysis.h index 83892d68..e837e4c5 100644 --- a/analysis/runAnalysis.h +++ b/analysis/runAnalysis.h @@ -22,7 +22,7 @@ #include "common/Communication.h" #include "common/ScaLBL.h" #include "threadpool/thread_pool.h" - +#include typedef std::shared_ptr> BlobIDstruct; typedef std::shared_ptr> BlobIDList; @@ -30,7 +30,7 @@ typedef std::shared_ptr> 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; diff --git a/common/Domain.cpp b/common/Domain.cpp index f16dbb8a..d65770de 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -92,6 +92,7 @@ Domain::Domain( std::shared_ptr 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 db ) inlet_layers_y = InletCount[1]; inlet_layers_z = InletCount[2]; } + if (d_db->keyExists( "OutletLayers" )){ + auto OutletCount = d_db->getVector( "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 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("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 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 0 && kproc() == nprocz()-1){ - for (int k=Nz-3; k(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; iid[i] = 1; // initialize this way //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object Averages = std::shared_ptr ( 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( "ComponentLabels" ); + auto LabelList = color_db->getVector( "ComponentLabels" ); auto AffinityList = color_db->getVector( "ComponentAffinity" ); NLABELS=LabelList.size(); @@ -220,12 +220,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) for (int idx=0; idxComm, label_count[idx]); if (rank==0){ - printf("Components labels: %lu \n",NLABELS); + printf("Component labels: %lu \n",NLABELS); for (unsigned int idx=0; idxPhi.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( "analysis_interval" ); } - + if (analysis_db->keyExists( "max_steady_timesteps" )){ + MAX_STEADY_TIMESTEPS = analysis_db->getScalar( "max_steady_timesteps" ); + } + if (analysis_db->keyExists( "max_morph_timesteps" )){ + MAX_MORPH_TIMESTEPS = analysis_db->getScalar( "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 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; kSDs(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; kSDs(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 analysis_db; IntArray Map; - char *id; + signed char *id; int *NeighborList; int *dvcMap; double *fq, *Aq, *Bq; diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index c57cddee..ff6b0843 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -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, diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c6775e68..37561fbe 100755 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -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 ) diff --git a/tests/TestColorBubble.cpp b/tests/TestColorBubble.cpp index 8025176b..0e6ea25a 100644 --- a/tests/TestColorBubble.cpp +++ b/tests/TestColorBubble.cpp @@ -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;kiproc(); - 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 } //*************************************************************************************** diff --git a/tests/lbpm_morphdrain_pp.cpp b/tests/lbpm_morphdrain_pp.cpp index c9a58222..675a76de 100644 --- a/tests/lbpm_morphdrain_pp.cpp +++ b/tests/lbpm_morphdrain_pp.cpp @@ -59,8 +59,8 @@ int main(int argc, char **argv) auto L = domain_db->getVector( "L" ); auto size = domain_db->getVector( "n" ); auto nproc = domain_db->getVector( "nproc" ); - auto ReadValues = domain_db->getVector( "ReadValues" ); - auto WriteValues = domain_db->getVector( "WriteValues" ); + auto ReadValues = domain_db->getVector( "ReadValues" ); + auto WriteValues = domain_db->getVector( "WriteValues" ); SW = domain_db->getScalar("Sw"); // Generate the NWP configuration @@ -82,8 +82,8 @@ int main(int argc, char **argv) for (n=0; nid[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( "ComponentLabels" ); - auto HistoryLabels = domain_db->getVector( "HistoryLabels" ); + auto LabelList = domain_db->getVector( "ComponentLabels" ); + auto HistoryLabels = domain_db->getVector( "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;kgetVector( "L" ); auto size = domain_db->getVector( "n" ); auto nproc = domain_db->getVector( "nproc" ); - auto ReadValues = domain_db->getVector( "ReadValues" ); - auto WriteValues = domain_db->getVector( "WriteValues" ); + auto ReadValues = domain_db->getVector( "ReadValues" ); + auto WriteValues = domain_db->getVector( "WriteValues" ); SW = domain_db->getScalar("Sw"); // Generate the NWP configuration @@ -82,8 +82,8 @@ int main(int argc, char **argv) for (n=0; nid[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( "ComponentLabels" ); - auto HistoryLabels = domain_db->getVector( "HistoryLabels" ); + auto LabelList = domain_db->getVector( "ComponentLabels" ); + auto HistoryLabels = domain_db->getVector( "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){ diff --git a/tests/lbpm_serial_decomp.cpp b/tests/lbpm_serial_decomp.cpp index 738e387e..41cb9686 100644 --- a/tests/lbpm_serial_decomp.cpp +++ b/tests/lbpm_serial_decomp.cpp @@ -87,8 +87,11 @@ int main(int argc, char **argv) if (domain_db->keyExists( "checkerSize" )){ checkerSize = domain_db->getScalar( "checkerSize" ); } - auto ReadValues = domain_db->getVector( "ReadValues" ); - auto WriteValues = domain_db->getVector( "WriteValues" ); + else { + checkerSize = SIZE[0]; + } + auto ReadValues = domain_db->getVector( "ReadValues" ); + auto WriteValues = domain_db->getVector( "WriteValues" ); auto ReadType = domain_db->getScalar( "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