Merge branch 'morphLBM' of github.com:JamesEMcClure/LBPM-WIA into morphLBM
This commit is contained in:
James E McClure
commit
82afdbc744
@ -28,23 +28,25 @@ 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))
|
||||
SUBPHASE = fopen("subphase.csv","a+");
|
||||
if (ftell(SUBPHASE) == 0)
|
||||
{
|
||||
// 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 +56,21 @@ 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
|
||||
}
|
||||
}
|
||||
|
||||
@ -75,41 +78,43 @@ 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 ",gwc.Px, gwd.Px, giwn.Px, gnc.Px, gnd.Px, giwn.Px);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Py, gwd.Py, giwn.Py, gnc.Py, gnd.Py, giwn.Py);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",gwc.Pz, gwd.Pz, giwn.Pz, gnc.Pz, gnd.Pz, giwn.Pz);
|
||||
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.Px, nc.Px, nd.Px, iwn.Px);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Py, wd.Py, iwn.Py, nc.Py, nd.Py, iwn.Py);
|
||||
fprintf(SUBPHASE,"%.5g %.5g %.5g %.5g %.5g %.5g ",wc.Pz, wd.Pz, iwn.Pz, nc.Pz, nd.Pz, iwn.Pz);
|
||||
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);
|
||||
fflush(SUBPHASE);
|
||||
}
|
||||
|
||||
}
|
||||
@ -203,9 +208,14 @@ void SubPhase::Basic(){
|
||||
gnb.Pz=sumReduce( Dm->Comm, nb.Pz);
|
||||
|
||||
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;
|
||||
|
||||
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*sqrt(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);
|
||||
@ -285,7 +295,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,6 +439,16 @@ 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);
|
||||
|
||||
double vol_nc_bulk = 0.0;
|
||||
double vol_wc_bulk = 0.0;
|
||||
@ -528,7 +548,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;
|
||||
|
||||
@ -35,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;
|
||||
@ -67,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)
|
||||
@ -102,7 +104,7 @@ public:
|
||||
|
||||
private:
|
||||
FILE *TIMELOG;
|
||||
|
||||
FILE *SUBPHASE;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -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);
|
||||
|
||||
@ -415,6 +415,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)
|
||||
@ -452,7 +453,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);
|
||||
@ -571,14 +577,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
|
||||
@ -772,10 +782,9 @@ double ScaLBL_ColorModel::MorphInit(const double beta, const double target_delta
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (volume_connected < 0.02*volume_initial){
|
||||
// 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");
|
||||
if (rank==0) printf("Connected region has shrunk to less than 2 %% of total fluid volume \n");
|
||||
REVERSE_FLOW_DIRECTION = true;
|
||||
}
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
@ -50,9 +48,17 @@ inline void InitializeBubble(ScaLBL_ColorModel &ColorModel, double BubbleRadius)
|
||||
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
|
||||
}
|
||||
//***************************************************************************************
|
||||
|
||||
Loading…
Reference in New Issue
Block a user