diff --git a/CMakeLists.txt b/CMakeLists.txt index 31518d86..c9b7e66a 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,8 +18,7 @@ SET( TEST_MAX_PROCS 16 ) # Initialize the project -PROJECT( ${PROJ} LANGUAGES CXX ) - +PROJECT( ${PROJ} LANGUAGES CXX) # Prevent users from building in place IF ("${CMAKE_CURRENT_SOURCE_DIR}" STREQUAL "${CMAKE_CURRENT_BINARY_DIR}" ) diff --git a/analysis/ElectroChemistry.cpp b/analysis/ElectroChemistry.cpp index 2a0b7169..dc932f91 100644 --- a/analysis/ElectroChemistry.cpp +++ b/analysis/ElectroChemistry.cpp @@ -53,19 +53,32 @@ void ElectroChemistryAnalyzer::Basic(ScaLBL_IonModel &Ion, ScaLBL_Poisson &Poiss Poisson.getElectricPotential(ElectricalPotential); /* local sub-domain averages */ - double rho_avg_local[Ion.number_ion_species]; - double rho_mu_avg_local[Ion.number_ion_species]; - double rho_mu_fluctuation_local[Ion.number_ion_species]; - double rho_psi_avg_local[Ion.number_ion_species]; - double rho_psi_fluctuation_local[Ion.number_ion_species]; + double *rho_avg_local; + double *rho_mu_avg_local; + double *rho_mu_fluctuation_local; + double *rho_psi_avg_local; + double *rho_psi_fluctuation_local; /* global averages */ - double rho_avg_global[Ion.number_ion_species]; - double rho_mu_avg_global[Ion.number_ion_species]; - double rho_mu_fluctuation_global[Ion.number_ion_species]; - double rho_psi_avg_global[Ion.number_ion_species]; - double rho_psi_fluctuation_global[Ion.number_ion_species]; + double *rho_avg_global; + double *rho_mu_avg_global; + double *rho_mu_fluctuation_global; + double *rho_psi_avg_global; + double *rho_psi_fluctuation_global; - for (int ion=0; ionrank()==0){ fprintf(TIMELOG,"%i ",timestep); - for (int ion=0; ion( Dm->rank_info,Dm->Nx-2,Dm->Ny-2,Dm->Nz-2,Dm->Lx,Dm->Ly,Dm->Lz ); auto ElectricPotential = std::make_shared(); std::vector> IonConcentration; - for (int ion=0; ion()); } auto VxVar = std::make_shared(); @@ -163,8 +176,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P } if (vis_db->getWithDefault( "save_concentration", true )){ - for (int ion=0; ionname = VisName; IonConcentration[ion]->type = IO::VariableType::VolumeVariable; IonConcentration[ion]->dim = 1; @@ -199,8 +212,8 @@ void ElectroChemistryAnalyzer::WriteVis( ScaLBL_IonModel &Ion, ScaLBL_Poisson &P } if (vis_db->getWithDefault( "save_concentration", true )){ - for (int ion=0; ionname = VisName; ASSERT(visData[0].vars[1+ion]->name==VisName); Array& IonConcentrationData = visData[0].vars[1+ion]->data; diff --git a/analysis/FreeEnergy.cpp b/analysis/FreeEnergy.cpp index 6a641a95..b49f2ec5 100644 --- a/analysis/FreeEnergy.cpp +++ b/analysis/FreeEnergy.cpp @@ -47,8 +47,6 @@ void FreeEnergyAnalyzer::SetParams(){ void FreeEnergyAnalyzer::Basic(ScaLBL_FreeLeeModel &LeeModel, int timestep){ - int i,j,k; - if (Dm->rank()==0){ fprintf(TIMELOG,"%i ",timestep); /*for (int ion=0; ion input_db, int timestep){ auto vis_db = input_db->getDatabase( "Visualization" ); - char VisName[40]; std::vector visData; fillHalo fillData(Dm->Comm,Dm->rank_info,{Dm->Nx-2,Dm->Ny-2,Dm->Nz-2},{1,1,1},0,1); diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 5349ae38..42ea8d7c 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -23,6 +23,7 @@ SubPhase::SubPhase(std::shared_ptr dm): Vel_x.resize(Nx,Ny,Nz); Vel_x.fill(0); // Gradient of the phase indicator field Vel_y.resize(Nx,Ny,Nz); Vel_y.fill(0); Vel_z.resize(Nx,Ny,Nz); Vel_z.fill(0); + Dissipation.resize(Nx,Ny,Nz); Dissipation.fill(0); SDs.resize(Nx,Ny,Nz); SDs.fill(0); //......................................... @@ -42,11 +43,12 @@ SubPhase::SubPhase(std::shared_ptr dm): //fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n"); fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn wet "); 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,"Mwc Mwd Mwi Mnc Mnd Mni Msw Msn "); // mass + fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x "); // momentum + fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y "); + fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z "); fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy + fprintf(SUBPHASE,"Dwc Dwd Dnc Dnd "); // viscous dissipation 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 @@ -67,11 +69,12 @@ SubPhase::SubPhase(std::shared_ptr dm): //fprintf(SUBPHASE,"--------------------------------------------------------------------------------------\n"); fprintf(SUBPHASE,"time rn rw nun nuw Fx Fy Fz iftwn wet "); 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,"Mwc Mwd Mwi Mnc Mnd Mni Msw Msn "); // mass + fprintf(SUBPHASE,"Pwc_x Pwd_x Pwi_x Pnc_x Pnd_x Pni_x Psw_x Psn_x "); // momentum + fprintf(SUBPHASE,"Pwc_y Pwd_y Pwi_y Pnc_y Pnd_y Pni_y Psw_y Psn_y "); + fprintf(SUBPHASE,"Pwc_z Pwd_z Pwi_z Pnc_z Pnd_z Pni_z Psw_z Psn_z "); fprintf(SUBPHASE,"Kwc Kwd Kwi Knc Knd Kni "); // kinetic energy + fprintf(SUBPHASE,"Dwc Dwd Dnc Dnd "); // viscous dissipation 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 @@ -93,7 +96,7 @@ SubPhase::SubPhase(std::shared_ptr dm): { // If timelog is empty, write a short header to list the averages //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"sw krw krn vw vn force pw pn wet\n"); + fprintf(TIMELOG,"sw krw krn krwf krnf vw vn force pw pn wet\n"); } } } @@ -111,11 +114,12 @@ void SubPhase::Write(int timestep) if (Dm->rank()==0){ fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn,total_wetting_interaction_global); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.p, gwd.p, gnc.p, gnd.p); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.M, gwd.M, giwn.Mw, gnc.M, gnd.M, giwn.Mn, gifs.Mw, gifs.Mn); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Px, gwd.Px, giwn.Pwx, gnc.Px, gnd.Px, giwn.Pnx, gifs.Pwx, gifs.Pnx); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Py, gwd.Py, giwn.Pwy, gnc.Py, gnd.Py, giwn.Pny, gifs.Pwy, gifs.Pny); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",gwc.Pz, gwd.Pz, giwn.Pwz, gnc.Pz, gnd.Pz, giwn.Pnz, gifs.Pwz, gifs.Pnz); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",gwc.K, gwd.K, giwn.Kw, gnc.K, gnd.K, giwn.Kn); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.visc, gwd.visc, gnc.visc, gnd.visc); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gwc.V, gwc.A, gwc.H, gwc.X); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",gwd.V, gwd.A, gwd.H, gwd.X, gwd.Nc); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",gnc.V, gnc.A, gnc.H, gnc.X); @@ -127,11 +131,12 @@ void SubPhase::Write(int timestep) else{ fprintf(SUBPHASE,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",timestep,rho_n,rho_w,nu_n,nu_w,Fx,Fy,Fz,gamma_wn,total_wetting_interaction); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.p, wd.p, nc.p, nd.p); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny); - fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.M, wd.M, iwn.Mw, nc.M, nd.M, iwn.Mn, ifs.Mw, ifs.Mn); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Px, wd.Px, iwn.Pwx, nc.Px, nd.Px, iwn.Pnx, ifs.Pwx, ifs.Pnx); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Py, wd.Py, iwn.Pwy, nc.Py, nd.Py, iwn.Pny, ifs.Pwy, ifs.Pny); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g ",wc.Pz, wd.Pz, iwn.Pwz, nc.Pz, nd.Pz, iwn.Pnz, ifs.Pwz, ifs.Pnz); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %.8g %.8g ",wc.K, wd.K, iwn.Kw, nc.K, nd.K, iwn.Kn); + fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.visc, wd.visc, nc.visc, nd.visc); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",wc.V, wc.A, wc.H, wc.X); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g %i ",wd.V, wd.A, wd.H, wd.X, wd.Nc); fprintf(SUBPHASE,"%.8g %.8g %.8g %.8g ",nc.V, nc.A, nc.H, nc.X); @@ -167,7 +172,7 @@ void SubPhase::Basic(){ if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z; if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z; */ - nb.reset(); wb.reset(); + nb.reset(); wb.reset(); iwn.reset(); double count_w = 0.0; double count_n = 0.0; @@ -245,6 +250,27 @@ void SubPhase::Basic(){ wb.p += Pressure(n); count_w += 1.0; } + /* compute the film contribution */ + else if (SDs(i,j,k) < 2.0){ + if ( phi > 0.0 ){ + nA = 1.0; + iwn.V += 1.0; + iwn.Mn += nA*rho_n; + // velocity + iwn.Pnx += rho_n*nA*Vel_x(n); + iwn.Pny += rho_n*nA*Vel_y(n); + iwn.Pnz += rho_n*nA*Vel_z(n); + } + else{ + nB = 1.0; + iwn.Mw += nB*rho_w; + iwn.V += 1.0; + + iwn.Pwx += rho_w*nB*Vel_x(n); + iwn.Pwy += rho_w*nB*Vel_y(n); + iwn.Pwz += rho_w*nB*Vel_z(n); + } + } } } } @@ -267,12 +293,6 @@ void SubPhase::Basic(){ //printf("wetting interaction = %f, count = %f\n",total_wetting_interaction,count_wetting_interaction); total_wetting_interaction_global=Dm->Comm.sumReduce( total_wetting_interaction); count_wetting_interaction_global=Dm->Comm.sumReduce( count_wetting_interaction); - /* normalize wetting interactions <-- Don't do this if normalizing laplacian (use solid surface area) - if (count_wetting_interaction > 0.0) - total_wetting_interaction /= count_wetting_interaction; - if (count_wetting_interaction_global > 0.0) - total_wetting_interaction_global /= count_wetting_interaction_global; - */ gwb.V=Dm->Comm.sumReduce( wb.V); gnb.V=Dm->Comm.sumReduce( nb.V); @@ -285,6 +305,16 @@ void SubPhase::Basic(){ gnb.Py=Dm->Comm.sumReduce( nb.Py); gnb.Pz=Dm->Comm.sumReduce( nb.Pz); + giwn.Mw=Dm->Comm.sumReduce( iwn.Mw); + giwn.Pwx=Dm->Comm.sumReduce( iwn.Pwx); + giwn.Pwy=Dm->Comm.sumReduce( iwn.Pwy); + giwn.Pwz=Dm->Comm.sumReduce( iwn.Pwz); + + giwn.Mn=Dm->Comm.sumReduce( iwn.Mn); + giwn.Pnx=Dm->Comm.sumReduce( iwn.Pnx); + giwn.Pny=Dm->Comm.sumReduce( iwn.Pny); + giwn.Pnz=Dm->Comm.sumReduce( iwn.Pnz); + count_w=Dm->Comm.sumReduce( count_w); count_n=Dm->Comm.sumReduce( count_n); if (count_w > 0.0) @@ -310,20 +340,21 @@ void SubPhase::Basic(){ if (gnb.Pz != gnb.Pz) err=true; if (Dm->rank() == 0){ + /* align flow direction based on total mass flux */ + double dir_x = gwb.Px + gnb.Px; + double dir_y = gwb.Py + gnb.Py; + double dir_z = gwb.Pz + gnb.Pz; + double flow_magnitude = dir_x*dir_x + dir_y*dir_y + dir_z*dir_z; double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz); - double dir_x = 0.0; - double dir_y = 0.0; - double dir_z = 0.0; if (force_mag > 0.0){ dir_x = Fx/force_mag; dir_y = Fy/force_mag; dir_z = Fz/force_mag; } else { - // default to z direction - dir_x = 0.0; - dir_y = 0.0; - dir_z = 1.0; + dir_x /= flow_magnitude; + dir_y /= flow_magnitude; + dir_z /= flow_magnitude; } if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 || Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4 ){ // compute the pressure drop @@ -331,7 +362,7 @@ void SubPhase::Basic(){ double length = ((Nz-2)*Dm->nprocz()); force_mag -= pressure_drop/length; } - if (force_mag == 0.0){ + if (force_mag == 0.0 && flow_magnitude == 0.0){ // default to z direction dir_x = 0.0; dir_y = 0.0; @@ -341,14 +372,21 @@ void SubPhase::Basic(){ double saturation=gwb.V/(gwb.V + gnb.V); double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume; double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume; + + /* contribution from water films */ + double water_film_flow_rate=gwb.V*(giwn.Pwx*dir_x + giwn.Pwy*dir_y + giwn.Pwz*dir_z)/gwb.M / Dm->Volume; + double not_water_film_flow_rate=gnb.V*(giwn.Pnx*dir_x + giwn.Pny*dir_y + giwn.Pnz*dir_z)/gnb.M / Dm->Volume; //double total_flow_rate = water_flow_rate + not_water_flow_rate; //double fractional_flow = water_flow_rate / total_flow_rate; double h = Dm->voxel_length; double krn = h*h*nu_n*not_water_flow_rate / force_mag ; double krw = h*h*nu_w*water_flow_rate / force_mag; + /* not counting films */ + double krnf = krn - h*h*nu_n*not_water_film_flow_rate / force_mag ; + double krwf = krw - h*h*nu_w*water_film_flow_rate / force_mag; //printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow); - fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global); + fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,krwf,krnf,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global); fflush(TIMELOG); } if (err==true){ @@ -364,6 +402,7 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl double A1,A2,A3,A4,A5,A6; double B1,B2,B3,B4,B5,B6; double nAB,delta; + double phi = (nA-nB)/(nA+nB); // Instantiate mass transport distributions // Stationary value - distribution 0 nAB = 1.0/(nA+nB); @@ -402,7 +441,7 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl double uwx = (B1-B2); double uwy = (B3-B4); double uwz = (B5-B6); - + /* I.Mn += rA*nA; I.Mw += rB*nB; I.Pnx += rA*nA*unx; @@ -413,6 +452,21 @@ inline void InterfaceTransportMeasures( double beta, double rA, double rB, doubl I.Pwz += rB*nB*uwz; I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz); I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz); + */ + if (phi > 0.0){ + I.Mn += rA; + I.Pnx += rA*ux; + I.Pny += rA*uy; + I.Pnz += rA*uz; + } + else { + I.Mw += rB; + I.Pwx += rB*ux; + I.Pwy += rB*uy; + I.Pwz += rB*uz; + } + I.Kn += rA*nA*(unx*unx + uny*uny + unz*unz); + I.Kw += rB*nB*(uwx*uwx + uwy*uwy + uwz*uwz); } @@ -432,7 +486,7 @@ void SubPhase::Full(){ if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z; if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z; */ - nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); iwnc.reset(); + nd.reset(); nc.reset(); wd.reset(); wc.reset(); iwn.reset(); iwnc.reset(); ifs.reset(); Dm->CommunicateMeshHalo(Phi); for (int k=1; kCommunicateMeshHalo(DelPhi); + + + Dm->CommunicateMeshHalo(Vel_x); + Dm->CommunicateMeshHalo(Vel_y); + Dm->CommunicateMeshHalo(Vel_z); + for (int k=1; k 2.0){ + Dissipation(i,j,k) = 2*rho*nu*( ux*ux + vy*vy + wz*wz + 0.5*(vx + uy)*(vx + uy)+ 0.5*(vz + wy)*(vz + wy)+ 0.5*(uz + wx)*(uz + wx)); + } + } + } + } /* Set up geometric analysis of each region */ @@ -605,13 +686,33 @@ void SubPhase::Full(){ double ux = Vel_x(n); double uy = Vel_y(n); double uz = Vel_z(n); + double visc = Dissipation(n); - if (DelPhi(n) > 1e-3){ - // interface region + if (DelPhi(n) > 1e-3 ){ + // get the normal vector double nx = 0.5*(Phi(i+1,j,k)-Phi(i-1,j,k)); double ny = 0.5*(Phi(i,j+1,k)-Phi(i,j-1,k)); double nz = 0.5*(Phi(i,j,k+1)-Phi(i,j,k-1)); - InterfaceTransportMeasures( beta, rho_w, rho_n, nA, nB, nx, ny, nz, ux, uy, uz, iwn); + if (SDs(n) > 2.5){ + // not a film region + InterfaceTransportMeasures( beta, rho_w, rho_n, nA, nB, nx, ny, nz, ux, uy, uz, iwn); + } + else{ + // films that are close to the wetting fluid + if ( morph_w->distance(i,j,k) < 2.5 && phi > 0.0){ + ifs.Mw += rho_w; + ifs.Pwx += rho_w*ux; + ifs.Pwy += rho_w*uy; + ifs.Pwz += rho_w*uz; + } + // films that are close to the NWP + if ( morph_n->distance(i,j,k) < 2.5 && phi < 0.0){ + ifs.Mn += rho_n; + ifs.Pnx += rho_n*ux; + ifs.Pny += rho_n*uy; + ifs.Pnz += rho_n*uz; + } + } } else if ( phi > 0.0){ if (morph_n->label(i,j,k) > 0 ){ @@ -642,6 +743,7 @@ void SubPhase::Full(){ nd.Py += nA*rho_n*uy; nd.Pz += nA*rho_n*uz; nd.K += nA*rho_n*(ux*ux + uy*uy + uz*uz); + nd.visc += visc; } else{ nA = 1.0; @@ -650,6 +752,7 @@ void SubPhase::Full(){ nc.Py += nA*rho_n*uy; nc.Pz += nA*rho_n*uz; nc.K += nA*rho_n*(ux*ux + uy*uy + uz*uz); + nc.visc += visc; } } else{ @@ -661,6 +764,7 @@ void SubPhase::Full(){ wd.Py += nB*rho_w*uy; wd.Pz += nB*rho_w*uz; wd.K += nB*rho_w*(ux*ux + uy*uy + uz*uz); + wd.visc += visc; } else{ nB = 1.0; @@ -669,6 +773,7 @@ void SubPhase::Full(){ wc.Py += nB*rho_w*uy; wc.Pz += nB*rho_w*uz; wc.K += nB*rho_w*(ux*ux + uy*uy + uz*uz); + wc.visc += visc; } } } @@ -681,25 +786,29 @@ void SubPhase::Full(){ gnd.Py=Dm->Comm.sumReduce( nd.Py); gnd.Pz=Dm->Comm.sumReduce( nd.Pz); gnd.K=Dm->Comm.sumReduce( nd.K); + gnd.visc=Dm->Comm.sumReduce( nd.visc); gwd.M=Dm->Comm.sumReduce( wd.M); gwd.Px=Dm->Comm.sumReduce( wd.Px); gwd.Py=Dm->Comm.sumReduce( wd.Py); gwd.Pz=Dm->Comm.sumReduce( wd.Pz); gwd.K=Dm->Comm.sumReduce( wd.K); + gwd.visc=Dm->Comm.sumReduce( wd.visc); gnc.M=Dm->Comm.sumReduce( nc.M); gnc.Px=Dm->Comm.sumReduce( nc.Px); gnc.Py=Dm->Comm.sumReduce( nc.Py); gnc.Pz=Dm->Comm.sumReduce( nc.Pz); gnc.K=Dm->Comm.sumReduce( nc.K); + gnc.visc=Dm->Comm.sumReduce( nc.visc); gwc.M=Dm->Comm.sumReduce( wc.M); gwc.Px=Dm->Comm.sumReduce( wc.Px); gwc.Py=Dm->Comm.sumReduce( wc.Py); gwc.Pz=Dm->Comm.sumReduce( wc.Pz); gwc.K=Dm->Comm.sumReduce( wc.K); - + gwc.visc=Dm->Comm.sumReduce( wc.visc); + giwn.Mn=Dm->Comm.sumReduce( iwn.Mn); giwn.Pnx=Dm->Comm.sumReduce( iwn.Pnx); giwn.Pny=Dm->Comm.sumReduce( iwn.Pny); @@ -711,6 +820,15 @@ void SubPhase::Full(){ giwn.Pwz=Dm->Comm.sumReduce( iwn.Pwz); giwn.Kw=Dm->Comm.sumReduce( iwn.Kw); + gifs.Mn= Dm->Comm.sumReduce( ifs.Mn); + gifs.Pnx=Dm->Comm.sumReduce( ifs.Pnx); + gifs.Pny=Dm->Comm.sumReduce( ifs.Pny); + gifs.Pnz=Dm->Comm.sumReduce( ifs.Pnz); + gifs.Mw= Dm->Comm.sumReduce( ifs.Mw); + gifs.Pwx=Dm->Comm.sumReduce( ifs.Pwx); + gifs.Pwy=Dm->Comm.sumReduce( ifs.Pwy); + gifs.Pwz=Dm->Comm.sumReduce( ifs.Pwz); + // pressure averaging gnc.p=Dm->Comm.sumReduce( nc.p); gnd.p=Dm->Comm.sumReduce( nd.p); diff --git a/analysis/SubPhase.h b/analysis/SubPhase.h index a6d35edd..c00f00de 100644 --- a/analysis/SubPhase.h +++ b/analysis/SubPhase.h @@ -22,10 +22,11 @@ class phase{ public: int Nc; double p; - double M,Px,Py,Pz,K; + double M,Px,Py,Pz,K,visc; double V,A,H,X; void reset(){ p=M=Px=Py=Pz=K=0.0; + visc=0.0; V=A=H=X=0.0; Nc=1; } @@ -70,10 +71,12 @@ public: // local entities phase wc,wd,wb,nc,nd,nb,solid; interface iwn,iwnc; + interface ifs; // global entities phase gwc,gwd,gwb,gnc,gnd,gnb,gsolid; interface giwn,giwnc; + interface gifs; /* fluid-solid wetting interaction */ double total_wetting_interaction, count_wetting_interaction; double total_wetting_interaction_global, count_wetting_interaction_global; @@ -92,6 +95,7 @@ public: DoubleArray Vel_x; // velocity field DoubleArray Vel_y; DoubleArray Vel_z; + DoubleArray Dissipation; DoubleArray SDs; std::shared_ptr morph_w; diff --git a/analysis/morphology.cpp b/analysis/morphology.cpp index 324c45cb..4b023db6 100644 --- a/analysis/morphology.cpp +++ b/analysis/morphology.cpp @@ -120,6 +120,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr sendtag = recvtag = 7; int ii,jj,kk; + int imin,jmin,kmin,imax,jmax,kmax; int Nx = nx; int Ny = ny; int Nz = nz; @@ -128,17 +129,13 @@ double MorphOpen(DoubleArray &SignDist, signed char *id, std::shared_ptr double void_fraction_new=1.0; double void_fraction_diff_old = 1.0; double void_fraction_diff_new = 1.0; - - // Increase the critical radius until the target saturation is met - double deltaR=0.05; // amount to change the radius in voxel units - double Rcrit_old; - - int imin,jmin,kmin,imax,jmax,kmax; - if (ErodeLabel == 1){ VoidFraction = 1.0 - VoidFraction; } + // Increase the critical radius until the target saturation is met + double deltaR=0.05; // amount to change the radius in voxel units + double Rcrit_old = maxdistGlobal; double Rcrit_new = maxdistGlobal; while (void_fraction_new > VoidFraction) @@ -406,6 +403,7 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptr2){ // Rcrit_new = strtod(argv[2],NULL); @@ -457,7 +452,6 @@ double MorphDrain(DoubleArray &SignDist, signed char *id, std::shared_ptrgetWithDefault( "save_dissipation", false ) ) { + ASSERT( visData[0].vars[5]->name == "ViscousDissipation" ); + Array &ViscousDissipation = visData[0].vars[5]->data; + fillData.copy( Averages.Dissipation, ViscousDissipation ); + } + if ( vis_db->getWithDefault( "save_distance", false ) ) { - ASSERT( visData[0].vars[5]->name == "SignDist" ); - Array &SignData = visData[0].vars[5]->data; + ASSERT( visData[0].vars[6]->name == "SignDist" ); + Array &SignData = visData[0].vars[6]->data; fillData.copy( Averages.SDs, SignData ); } if ( vis_db->getWithDefault( "save_connected_components", false ) ) { - ASSERT( visData[0].vars[6]->name == "BlobID" ); - Array &BlobData = visData[0].vars[6]->data; + ASSERT( visData[0].vars[7]->name == "BlobID" ); + Array &BlobData = visData[0].vars[7]->data; fillData.copy( Averages.morph_n->label, BlobData ); } @@ -652,6 +658,7 @@ runAnalysis::runAnalysis( std::shared_ptr input_db, const RankInfoStru auto VxVar = std::make_shared(); auto VyVar = std::make_shared(); auto VzVar = std::make_shared(); + auto ViscousDissipationVar = std::make_shared(); auto SignDistVar = std::make_shared(); auto BlobIDVar = std::make_shared(); @@ -688,6 +695,14 @@ runAnalysis::runAnalysis( std::shared_ptr input_db, const RankInfoStru VzVar->data.resize( d_n[0], d_n[1], d_n[2] ); d_meshData[0].vars.push_back( VzVar ); } + + if ( vis_db->getWithDefault( "save_dissipation", false ) ) { + ViscousDissipationVar->name = "ViscousDissipation"; + ViscousDissipationVar->type = IO::VariableType::VolumeVariable; + ViscousDissipationVar->dim = 1; + ViscousDissipationVar->data.resize( d_n[0], d_n[1], d_n[2] ); + d_meshData[0].vars.push_back( ViscousDissipationVar ); + } if ( vis_db->getWithDefault( "save_distance", false ) ) { SignDistVar->name = "SignDist"; @@ -729,7 +744,6 @@ runAnalysis::runAnalysis( ScaLBL_ColorModel &ColorModel) d_comm = ColorModel.Dm->Comm.dup(); d_Np = ColorModel.Np; - bool Regular = false; auto input_db = ColorModel.db; auto db = input_db->getDatabase( "Analysis" ); diff --git a/common/Domain.cpp b/common/Domain.cpp index b38add0c..6a4c61fb 100644 --- a/common/Domain.cpp +++ b/common/Domain.cpp @@ -383,7 +383,7 @@ void Domain::Decomp( const std::string& Filename ) for (int i = 0; i0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2068,6 +2082,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2088,6 +2106,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2109,9 +2131,9 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map //CP: capillary penalty // also turn off recoloring for grey nodes extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi,double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Velocity,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, - double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){ + double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; double fq; @@ -2127,18 +2149,16 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do double a1,b1,a2,b2,nAB,delta; double C,nx,ny,nz; //color gradient magnitude and direction double phi,tau,rho0,rlx_setA,rlx_setB; - + double W;//greyscale wetting strength + double Sn_grey,Sw_grey; + //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) double porosity; double perm;//voxel permeability //double c0, c1; //Guo's model parameters double tau_eff; double mu_eff;//kinematic viscosity - double nx_gs,ny_gs,nz_gs;//grey-solid color gradient - double nx_phase,ny_phase,nz_phase,C_phase; double Fx,Fy,Fz; - double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18; - double gp3,gp5,gp7; double Fcpx,Fcpy,Fcpz;//capillary penalty force const double mrt_V1=0.05263157894736842; @@ -2155,15 +2175,15 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do const double mrt_V12=0.04166666666666666; for (n=start; n0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2764,6 +2798,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2779,6 +2817,10 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2790,6 +2832,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do Aq[6*Np+n] = a2; Bq[6*Np+n] = b2; //............................................... + } } diff --git a/cuda/GreyscaleColor.cu b/cuda/GreyscaleColor.cu index d5ab9460..6132a30c 100644 --- a/cuda/GreyscaleColor.cu +++ b/cuda/GreyscaleColor.cu @@ -1450,7 +1450,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ @@ -1478,6 +1478,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int double Fx,Fy,Fz; double Fcpx,Fcpy,Fcpz;//capillary penalty force double W;//greyscale wetting strength + double Sn_grey,Sw_grey; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1504,6 +1505,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int porosity = Poros[n]; perm = Perm[n]; W = GreySolidW[n]; + Sn_grey = GreySn[n]; + Sw_grey = GreySw[n]; // compute phase indicator field phi=(nA-nB)/(nA+nB); @@ -2165,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2184,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2204,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2226,7 +2241,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; @@ -2249,6 +2264,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis double Fx,Fy,Fz; double Fcpx,Fcpy,Fcpz;//capillary penalty force double W;//greyscale wetting strength + double Sn_grey,Sw_grey; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -2276,6 +2292,8 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis porosity = Poros[n]; perm = Perm[n]; W = GreySolidW[n]; + Sn_grey = GreySn[n]; + Sw_grey = GreySw[n]; // compute phase indicator field phi=(nA-nB)/(nA+nB); @@ -2870,6 +2888,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2886,6 +2908,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2901,6 +2927,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -4500,12 +4530,13 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm, Vel, Pressure, + dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure, rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); + cudaError_t err = cudaGetLastError(); if (cudaSuccess != err){ printf("CUDA error in ScaLBL_D3Q19_AAeven_GreyscaleColor_CP: %s \n",cudaGetErrorString(err)); @@ -4515,11 +4546,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *GreySolidW, double *Poros,double *Perm,double *Vel,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, Poros, Perm,Vel,Pressure, + dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure, rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); cudaError_t err = cudaGetLastError(); diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 00000000..7fad3c1f --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = $(HOME)/local/doc/build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/README.md b/docs/README.md new file mode 100644 index 00000000..bd31ba05 --- /dev/null +++ b/docs/README.md @@ -0,0 +1,19 @@ +Dependencies for LBPM documentation + +# install sphinx +pip install Sphinx + +# foamatting requires sphinx read-the-docs-theme +pip install sphinx-rtd-theme + +# equation rendering requires latex and dvipng command +sudo apt-get install dvipng +sudo apt-get install texlive texstudio +sudo apt-get install texlive-latex-recommended texlive-pictures texlive-latex-extra + + +# To build the docs +Step 1) install dependencies listed above +Step 2) type 'make html' from the docs/ directory +Step 3) point your browser at ~/local/doc/build/html/index.html +# \ No newline at end of file diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 00000000..fa5df96e --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,70 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +# sys.path.insert(0, os.path.abspath('.')) + + +# -- Project information ----------------------------------------------------- + +project = 'LBPM' +copyright = '2021, James E McClure' +author = 'James E McClure' + +# The full version, including alpha/beta/rc tags +release = '1.0' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.imgmath' +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'alabaster' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] + + +## Read the docs style: +if os.environ.get('READTHEDOCS') != 'True': + try: + import sphinx_rtd_theme + except ImportError: + pass # assume we have sphinx >= 1.3 + else: + html_theme_path = [sphinx_rtd_theme.get_html_theme_path()] + html_theme = 'sphinx_rtd_theme' + +#def setup(app): +# app.add_stylesheet("fix_rtd.css") diff --git a/docs/source/developerGuide/buildingModels/overview.rst b/docs/source/developerGuide/buildingModels/overview.rst new file mode 100644 index 00000000..223b7c99 --- /dev/null +++ b/docs/source/developerGuide/buildingModels/overview.rst @@ -0,0 +1,20 @@ +=========================== +Implementing a new LB model +=========================== + +While LBPM includes a range of fully-functioning lattice Boltzmann models, the commonly used +Bhatnager-Gross-Krook (BGK) model has been deliberately excluded. While the physical limitations +of this model are well-known, implementing the BGK model is an excellent way to understand +how to implement new LB models within the more general framework of LBPM. In this excercise +you will + +* learn "what goes where" + + + +* don't modify core data structures (unless you have a really good reason) + +* re-use existing components whenever possible + + + diff --git a/docs/source/developerGuide/designOverview/dataStructures.rst b/docs/source/developerGuide/designOverview/dataStructures.rst new file mode 100644 index 00000000..bea4f91b --- /dev/null +++ b/docs/source/developerGuide/designOverview/dataStructures.rst @@ -0,0 +1,6 @@ +=============== +Data Structures +=============== + +LBPM includes a variety of generalized data structures to facilitate the implementation +of different lattice Boltzmann models. diff --git a/docs/source/developerGuide/index.rst b/docs/source/developerGuide/index.rst new file mode 100644 index 00000000..d8448fa1 --- /dev/null +++ b/docs/source/developerGuide/index.rst @@ -0,0 +1,18 @@ +############################################################################### +Developer guide +############################################################################### + +The LBPM developer guide provides essential information on how to add new physics +into the framework. + +.. toctree:: + :glob: + :maxdepth: 2 + + designOverview/* + + buildingModels/* + + testingModels/* + + diff --git a/docs/source/developerGuide/testingModels/unitTests.rst b/docs/source/developerGuide/testingModels/unitTests.rst new file mode 100644 index 00000000..4a4e84bf --- /dev/null +++ b/docs/source/developerGuide/testingModels/unitTests.rst @@ -0,0 +1,9 @@ +================= +Adding unit tests +================= + +Unit tests in LBPM are implemented using ctest + +* general overview + +* launching unit tests for GPU (MPI flags etc.) diff --git a/docs/source/examples/running.rst b/docs/source/examples/running.rst new file mode 100644 index 00000000..83993f59 --- /dev/null +++ b/docs/source/examples/running.rst @@ -0,0 +1,9 @@ +============ +Running LBPM +============ + +There are two main components to running LBPM simulators. +First is understanding how to launch MPI tasks on your system, +which depends on the particular implementation of MPI that you are using, +as well as other details of the local configuration. The second component is +understanding the LBPM input file structure. diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 00000000..fc5225d1 --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,25 @@ +.. LBPM documentation master file, created by + sphinx-quickstart on Thu May 20 12:19:14 2021. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +LBPM -- Documentation +=================================================== + +.. toctree:: + :glob: + :maxdepth: 2 + :caption: Contents: + + install + examples/* + userGuide/* + developerGuide/* + publications/* + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/source/install.rst b/docs/source/install.rst new file mode 100644 index 00000000..635c92d0 --- /dev/null +++ b/docs/source/install.rst @@ -0,0 +1,164 @@ +============ +Installation +============ + +Sample scripts to configure and build LBPM are included in the `sample_scripts` directory. To build this package, MPI and cmake are required, along with a compiler that supports the C++-14 standard. Building in source is not supported. + +The essential dependencies needed to build LBPM are: + +1. cmake (version 3.9 or higher) +2. MPI +3. HDF5 +4. silo + +************************* +Building Dependencies +************************* +(skip if they are already installed on your system) + +1. Download third-party library source files + +.. code-block:: bash + + zlib-1.2.11.tar.gz + hdf5-1.8.12.tar.gz + silo-4.10.2.tar.gz + + +2. Set up path where you want to install + +.. code-block:: bash + + export MPI_DIR=/path/to/mpi/ + export LBPM_ZLIB_DIR=/path/to/zlib + export LBPM_HDF5_DIR=/path/to/hdf5 + export LBPM_SILO_DIR=/path/to/silo + +3. Build third-party library dependencies + +.. code-block:: bash + + tar -xzvf zlib-1.2.11.tar.gz + tar -xzvf hdf5-1.8.12.tar.gz + tar -xzvf silo-4.10.2.tar.gz + cd zlib-1.2.11 + + ./configure --prefix=$LBPM_ZLIB_DIR && make && make install + + cd ../hdf5-1.8.12 + + CC=$MPI_DIR/bin/mpicc CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \ + ./configure --prefix=$LBPM_HDF5_DIR --enable-parallel --enable-shared --with-zlib=$LBPM_ZLIB_DIR + make && make install + + + cd ../silo-4.10.2 + + CC=$MPI_DIR/bin/mpicc CXX=$MPI_DIR/bin/mpicxx CXXFLAGS="-fPIC -O3 -std=c++14" \ + ./configure --prefix=$LBPM_SILO_DIR -with-hdf5=$LBPM_HDF5_DIR/include,$LBPM_HDF5_DIR/lib --enable-static + make && make install + +************************* +Building LBPM +************************* + +Many HPC systems will include all of these dependencies, and LBPM can be built simply by setting the paths to the associated libraries in the cmake configure line. + +The steps typically used to build LBPM are as follows: + +1. Set environment variables for the source directory `$LBPM_WIA_SOURCE` and build directory `$LBPM_WIA_DIR` + +.. code-block:: bash + + export LBPM_SOURCE=/path/to/source/LBPM + export LBPM_DIR=/path/to/build/LBPM + +2. Set environment variables for the path to HDF5 and SILO (required), and optionally for TimerUtility and NetCDF (optional) + +.. code-block:: bash + + export LBPM_HDF5_DIR=/path/to/hdf5 + export LBPM_SILO_DIR=/path/to/silo + export LBPM_TIMER_DIR=/path/to/timer + export LBPM_NETCDF_DIR=/path/to/netcdf + +3. Create the build directory and navigate to it + +.. code-block:: bash + + mkdir $LBPM_WIA_DIR + cd $LBPM_WIA_DIR + +4. Configure the project. Numerous scripts exist to build LBPM on different HPC clusters, which are available in the `$LBPM_SOURCE/sample_scripts/` directory. It is often possible to run these scripts directly if one exists for a system similar to the one you are building on. For a standard CPU build: + +.. code-block:: bash + + cmake \ + -D CMAKE_BUILD_TYPE:STRING=Release \ + -D CMAKE_C_COMPILER:PATH=mpicc \ + -D CMAKE_CXX_COMPILER:PATH=mpicxx \ + -D CMAKE_C_FLAGS="-fPIC" \ + -D CMAKE_CXX_FLAGS="-fPIC" \ + -D CMAKE_CXX_STD=14 \ + -D USE_TIMER=0 \ + -D TIMER_DIRECTORY=$LBPM_TIMER_DIR \ + -D USE_NETCDF=0 \ + -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR \ + -D USE_SILO=1 \ + -D HDF5_DIRECTORY=$LBPM_HDF5_DIR \ + -D SILO_DIRECTORY=$LBPM_SILO_DIR \ + -D USE_CUDA=0 \ + $LBPM_SOURCE + +For GPU support, it is necessary to have CUDA along with a GPU-aware MPI implementation. Otherwise, the LBPM routines should behave identically irrespective of the underlying hardware. + +.. code-block:: bash + + cmake \ + -D CMAKE_BUILD_TYPE:STRING=Release \ + -D CMAKE_C_COMPILER:PATH=mpicc \ + -D CMAKE_CXX_COMPILER:PATH=mpicxx \ + -D CMAKE_C_FLAGS="-fPIC" \ + -D CMAKE_CXX_FLAGS="-fPIC" \ + -D CMAKE_CXX_STD=14 \ + -D USE_TIMER=0 \ + -D TIMER_DIRECTORY=$LBPM_TIMER_DIR \ + -D USE_NETCDF=0 \ + -D NETCDF_DIRECTORY=$LBPM_NETCDF_DIR \ + -D USE_SILO=1 \ + -D HDF5_DIRECTORY=$LBPM_HDF5_DIR \ + -D SILO_DIRECTORY=$LBPM_SILO_DIR \ + -D USE_CUDA=1 \ + -D CMAKE_CUDA_FLAGS="-arch sm_70" \ + $LBPM_SOURCE + +5. Build the project (using four cores to build) + +.. code-block:: bash + + make -j4 + +6. Install the project + +.. code-block:: bash + + make install + +7. Run the tests to make sure they execute correctly (on a cluster, it is recommended to run these using the batch system rather than on the head node) + +.. code-block:: bash + + ctest + + +************************* +Sample Scripts +************************* + +The LBPM repository contains sample scripts showing successful CMake configuration, build and +install steps for a range of systems. Refer to the project sub-directory below for these examples. + +.. code-block:: bash + + ls $LBPM_SOURCE/sample_scripts + diff --git a/docs/source/publications/publications.rst b/docs/source/publications/publications.rst new file mode 100644 index 00000000..aea8e5b6 --- /dev/null +++ b/docs/source/publications/publications.rst @@ -0,0 +1,13 @@ +============ +Publications +============ + +* James E McClure, Zhe Li, Mark Berrill, Thomas Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9 + + +* James E. McClure, Zhe Li, Adrian P. Sheppard, Cass T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670 + + +* Y.D. Wang, T. Chung, R.T. Armstrong, J. McClure, T. Ramstad, P. Mostaghimi, "Accelerated Computation of Relative Permeability by Coupled Morphological and Direct Multiphase Flow Simulation" Journal of Computational Physics (401) (2020) https://doi.org/10.1016/j.jcp.2019.108966 + + diff --git a/docs/source/userGuide/IO/fileformat.rst b/docs/source/userGuide/IO/fileformat.rst new file mode 100644 index 00000000..faec04a6 --- /dev/null +++ b/docs/source/userGuide/IO/fileformat.rst @@ -0,0 +1,13 @@ +======================== +I/O conventions for LBPM +======================== + +There are three main kinds of output file that are supported by LBPM. + + +* CSV files -- + +* formatted binary files -- + +* unformatted binary files -- + diff --git a/docs/source/userGuide/analysis/analysis.rst b/docs/source/userGuide/analysis/analysis.rst new file mode 100644 index 00000000..10d5282a --- /dev/null +++ b/docs/source/userGuide/analysis/analysis.rst @@ -0,0 +1,5 @@ +=========================== +Internal Analysis Framework +=========================== + +placeholder for analysis diff --git a/docs/source/userGuide/index.rst b/docs/source/userGuide/index.rst new file mode 100644 index 00000000..49025202 --- /dev/null +++ b/docs/source/userGuide/index.rst @@ -0,0 +1,17 @@ +############################################################################### +User Guide +############################################################################### + +Welcome to the LBPM user guide. + +.. toctree:: + :glob: + :maxdepth: 2 + + models/* + + analysis/* + + visualization/* + + IO/* diff --git a/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst b/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst new file mode 100644 index 00000000..93dbaff0 --- /dev/null +++ b/docs/source/userGuide/models/PoissonBoltzmann/PoissonBoltzmann.rst @@ -0,0 +1,6 @@ +============================================= +Poisson-Boltzmann model +============================================= + +The LBPM Poisson-Boltzmann solver is designed to solve the Poisson-Boltzmann equation +to solve for the electric field in an ionic fluid. diff --git a/docs/source/userGuide/models/color/index.rst b/docs/source/userGuide/models/color/index.rst new file mode 100644 index 00000000..23bb5949 --- /dev/null +++ b/docs/source/userGuide/models/color/index.rst @@ -0,0 +1,91 @@ +############################################################################### +Color model +############################################################################### + +The LBPM color model is implemented by combining a multi-relaxation time D3Q19 +lattice Boltzmann equation (LBE) to solve for the momentum transport with two D3Q7 +LBEs for the mass transport. The color model will obey strict mass and momentum +conservation while minimizing diffusive fluxes across the interface between fluids. +The color model is a good choice for modeling dense fluids that are strongly immiscible +(e.g. water-oil systems). Due to the strong anti-diffusion in the interface region, +the color model is not suitable for modeling processes such as Ostwald ripening that +depend on diffusive fluxes between fluid phases. + +A typical command to launch the LBPM color simulator is as follows + +``` +mpirun -np $NUMPROCS lbpm_color_simulator input.db +``` + +where ``$NUMPROCS`` is the number of MPI processors to be used and ``input.db`` is +the name of the input database that provides the simulation parameters. +Note that the specific syntax to launch MPI tasks may vary depending on your system. +For additional details please refer to your local system documentation. + +**************************** +Simulation protocols +**************************** + +Simulation protocols are designed to make it simpler to design and execute common +computational experiments. Protocols will automatically determine boundary conditions +needed to perform a particular simulation. LBPM will internall set default simulation paramaters +that can be over-ridden to develop customized simulations. + +.. toctree:: + :glob: + :maxdepth: 2 + + protocols/* + + +*************************** +Model parameters +*************************** + +The essential model parameters for the color model are + +- :math:`\alpha` -- control the interfacial tension between fluids with key ``alpha`` +- :math:`\beta` -- control the width of the interface with key ``beta`` +- :math:`\tau_A` -- control the viscosity of fluid A with key ``tauA`` +- :math:`\tau_B` -- control the viscosity of fluid B with key ``tauB`` + +**************************** +Model Formulation +**************************** + + + +**************************** +Boundary Conditions +**************************** + +The following external boundary conditions are supported by ``lbpm_color_simulator`` +and can be set by setting the ``BC`` key values in the ``Domain`` section of the +input file database + +- ``BC = 0`` -- fully periodic boundary conditions +- ``BC = 3`` -- constant pressure boundary condition +- ``BC = 4`` -- constant volumetric flux boundary condition + +For ``BC = 0`` any mass that exits on one side of the domain will re-enter at the other +side. If the pore-structure for the image is tight, the mismatch between the inlet and +outlet can artificially reduce the permeability of the sample due to the blockage of +flow pathways at the boundary. LBPM includes an internal utility that will reduce the impact +of the boundary mismatch by eroding the solid labels within the inlet and outlet layers +(https://doi.org/10.1007/s10596-020-10028-9) to create a mixing layer. +The number mixing layers to use can be set using the key values in the ``Domain`` section +of the input database + +- ``InletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the inlet +- ``OUtletLayers = 5`` -- set the number of mixing layers to ``5`` voxels at the outlet + +For the other boundary conditions a thin reservoir of fluid (default ``3`` voxels) +is established at either side of the domain. The inlet is defined as the boundary face +where ``z = 0`` and the outlet is the boundary face where ``z = nprocz*nz``. By default a +reservoir of fluid A is established at the inlet and a reservoir of fluid B is established at +the outlet, each with a default thickness of three voxels. To over-ride the default label at +the inlet or outlet, the ``Domain`` section of the database may specify the following key values + +- ``InletLayerPhase = 2`` -- establish a reservoir of component B at the inlet +- ``OutletLayerPhase = 1`` -- establish a reservoir of component A at the outlet + diff --git a/docs/source/userGuide/models/color/protocols/centrifuge.rst b/docs/source/userGuide/models/color/protocols/centrifuge.rst new file mode 100644 index 00000000..29ebd186 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/centrifuge.rst @@ -0,0 +1,43 @@ +====================================== +Color model -- Centrifuge Protocol +====================================== + +The centrifuge protocol is designed to mimic SCAL centrifuge experiments that +are used to infer the capillary pressure. The LBPM centrifuge protocol is +constructed as an unsteady simulation with constant pressure boundary conditions +and zero pressure drop across the sample. This will enforce the following key values + +* ``BC = 3`` -- constant pressure boundary condition +* ``din = 1.0`` -- inlet pressure value +* ``dout = 1.0`` -- outlet pressure value + +By default LBPM will populate the inlet reservoir with fluid A (usually the non-wetting fluid) +and the outlet reservoir with fluid B (usually water). Flow is induced by setting an external +body force to generate displacement in the ``z`` direction. If the body force is set to +zero, e.g. ``F = 0, 0, 0``, the simulation will produce spontaneous imbibition, with the +balancing point being determined based on zero pressure drop across the sample. Setting +an external body force will shift the capillary pressure. Setting a positive force will +cause fluid A to be forced into the sample. Once steady conditions are achieved, +the pressure of fluid A will be larger than fluid B. Alternatively, if the driving force is +negative then fluid B will be forced into the sample, and the steady-state configuration +will stabilize to a configuration where fluid B has a larger pressure compared to fluid A. +The capillary pressure is thereby inferred based on the body force. + +In a conventional SCAL experiment the centrifugal forces are proportional to the density +difference between fluids. While this is still true for LBPM simulation, the body force will +still be effective even if there is no difference in density between the fluids. +This is because a positive body force will favor a larger saturation of fluid A +(positive capillary pressure ) whereas a negative body force will favor a lower +saturation of fluid A (negative capillary pressure). + + +To enable the ``centrifuge`` protocol such that the pressure of fluid B is higher than +fluid A, the following keys can be set. Increasing the body force will lead to a larger +capillary pressure + +.. code-block:: bash + + protocol = "centrifuge" + F = 0, 0, -1.0e-5 + + diff --git a/docs/source/userGuide/models/color/protocols/coreFlooding.rst b/docs/source/userGuide/models/color/protocols/coreFlooding.rst new file mode 100644 index 00000000..12d1e799 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/coreFlooding.rst @@ -0,0 +1,63 @@ +====================================== +Color model -- Core Flooding +====================================== + +The core flooding protocol is designed to mimic SCAL experiments where one +immiscible fluid is injected into the sample at a constant rate, displacing the +other fluid. The core flooding protocol relies on a flux boundary condition +to ensure that fluid is injected into the sample at a constant rate. The flux +boundary condition implements a time-varying pressure boundary condition that +adapts to ensure a constant volumetric flux. Details for the flux boundary +condition are available +(see: https://doi.org/10.1016/j.compfluid.2020.104670) + +.. code-block:: bash + + protocol = "core flooding" + + +To match experimental conditions, it is usually important to match the capillary +number, which is + +.. math:: + \mbox{Ca} = \frac{\mu u_z}{\gamma} + + +where :math:`\mu` is the dynamic viscosity, :math:`u_z` is the fluid +(usually water) velocity and :math:`\gamma` is the interfacial tension. +The volumetric flow rate is related to the fluid velocity based on + +.. math:: + Q_z = \epsilon C_{xy} u_z + +where :math:`C_{xy}` is the cross-sectional area and :math:`\epsilon` +is the porosity. Given a particular experimental system +self-similar conditions can be determined for the lattice Boltzmann +simulation by matching the non-dimensional :math:`mbox{Ca}`. It is nearly +awlays advantageous for the timestep to be as large as possible so +that time-to-solution will be more favorable. This is accomplished by + +* use a high value for the numerical surface tension (e.g. ``alpha=1.0e-2``) +* use a small value for the fluid viscosity (e.g. ``tau_w = 0.7`` and ``tau_n = 0.7`` ) +* determine the volumetric flow rate needed to match :math:`\mbox{Ca}` + +For the color LBM the interfacial tension is +:math:`\gamma = 6 \alpha` and the dynamic viscosity is :math:`\mu = \rho(\tau-1/2)/3`, +where the units are relative to the lattice spacing, timestep and mass +density. Agreemetn between the experimental and simulated values for +:math:`\mbox{Ca}` is ensured by setting the volumetric flux + +.. math:: + Q_z = \frac{\epsilon C_{xy} \gamma }{\mu} \mbox{Ca} + +where the LB units of the volumetric flux will be voxels per timestep. + +In some situations it may also be important to match other non-dimensional numbers, +such as the viscosity ratio, density ratio, and/or Ohnesorge/Laplace number. This +can be accomplished with an analogous procedure. Enforcing additional constraints +will necessarily restrict the LB parameters that can be used, which are ultimately +manifested as a constraint on the size of a timestep. + + + + diff --git a/docs/source/userGuide/models/color/protocols/fractionalFlow.rst b/docs/source/userGuide/models/color/protocols/fractionalFlow.rst new file mode 100644 index 00000000..7f5a26f0 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/fractionalFlow.rst @@ -0,0 +1,17 @@ +========================================== +Color model -- Fractional Flow Protocol +========================================== + +The fractional flow protocol is designed to perform steady-state relative +permeability simulations by using an internal routine to change the fluid +saturation by adding and subtracting mass to the fluid phases. The +mass density is updated for each fluid based on the locations where +the local rate of flow is high. + +.. code-block:: bash + + protocol = "fractional flow" + + + + diff --git a/docs/source/userGuide/models/color/protocols/imageSequence.rst b/docs/source/userGuide/models/color/protocols/imageSequence.rst new file mode 100644 index 00000000..25c91e6a --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/imageSequence.rst @@ -0,0 +1,27 @@ +====================================== +Color model -- Image Sequence Protocol +====================================== + +The image sequence protocol is designed to perform a set steady-state +simulations based on a sequence of 3D (8-bit) images provided by the user. +The images might be the output of a previous LBPM simulation, a sequence of +(segmented) experimental data, or data generated from a custom routine. +The image sequence protocol will apply the same set of flow conditions +to all images in the sequence. This means + +* the image labels and any associated properties are the same +* the external boundary conditions are the same +* the physical simulation parameters are the same + +The image sequence protocol does not set boundary conditions by default. +It is up to the user to determine the flow condition, with the understanding +that the same set of will be applied to each image in the sequence. + +To enable the image sequence protocol, the following keys should be set +within the ``Color`` section of the input database + +.. code-block:: bash + + protocol = "image sequence" + image_sequence = "image1.raw", "image2.raw" + diff --git a/docs/source/userGuide/models/color/protocols/shellAggregation.rst b/docs/source/userGuide/models/color/protocols/shellAggregation.rst new file mode 100644 index 00000000..1ab54131 --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/shellAggregation.rst @@ -0,0 +1,18 @@ +========================================== +Color model -- Shell Aggregation Protocol +========================================== + +The shell aggregation protocol is designed to perform steady-state relative +permeability simulations by using an internal routine to change the fluid +saturation by moving the interface. The basic design for the shell aggregation +protocol is described by Wang et al. ( https://doi.org/10.1016/j.jcp.2019.108966 ). + + +.. code-block:: bash + + + protocol = "shell aggregation" + + + + diff --git a/docs/source/userGuide/models/freeEnergy/freeEnergy.rst b/docs/source/userGuide/models/freeEnergy/freeEnergy.rst new file mode 100644 index 00000000..f53cfaad --- /dev/null +++ b/docs/source/userGuide/models/freeEnergy/freeEnergy.rst @@ -0,0 +1,6 @@ +============================================= +Free energy model +============================================= + +The LBPM free energy model is constructed to solve the Allen-Cahn equations, +which are typically used to describe liquid-gas systems. diff --git a/docs/source/userGuide/models/greyscale/greyscale.rst b/docs/source/userGuide/models/greyscale/greyscale.rst new file mode 100644 index 00000000..47c0621f --- /dev/null +++ b/docs/source/userGuide/models/greyscale/greyscale.rst @@ -0,0 +1,7 @@ +============================================= +Greyscale model model +============================================= + +The LBPM greyscale lattice Boltzmann model is constructed to approximate the +solution of the Darcy-Brinkman equations in grey regions, coupled to a Navier-Stokes +solution in open regions. diff --git a/docs/source/userGuide/models/index.rst b/docs/source/userGuide/models/index.rst new file mode 100644 index 00000000..6a158aa6 --- /dev/null +++ b/docs/source/userGuide/models/index.rst @@ -0,0 +1,23 @@ +############################################################################### +LBPM model summary +############################################################################### + +Currently supported lattice Boltzmann models + +.. toctree:: + :glob: + :maxdepth: 2 + + color/* + + mrt/* + + nernstPlanck/* + + PoissonBoltzmann/* + + greyscale/* + + freeEnergy/* + + diff --git a/docs/source/userGuide/models/mrt/mrt.rst b/docs/source/userGuide/models/mrt/mrt.rst new file mode 100644 index 00000000..1e232ea6 --- /dev/null +++ b/docs/source/userGuide/models/mrt/mrt.rst @@ -0,0 +1,6 @@ +============================================= +MRT model +============================================= + +The multi-relaxation time (MRT) lattice Boltzmann model is constructed to approximate the +solution of the Navier-Stokes equations. diff --git a/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst b/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst new file mode 100644 index 00000000..3842e2e6 --- /dev/null +++ b/docs/source/userGuide/models/nernstPlanck/nerstPlanck.rst @@ -0,0 +1,6 @@ +============================================= +Nernst-Planck model +============================================= + +The Nernst-Planck model is designed to model ion transport based on the +Nernst-Planck equation. diff --git a/docs/source/userGuide/visualization/visit.rst b/docs/source/userGuide/visualization/visit.rst new file mode 100644 index 00000000..f9ff2c8a --- /dev/null +++ b/docs/source/userGuide/visualization/visit.rst @@ -0,0 +1,5 @@ +====================================== +Visualizing simulation data with visit +====================================== + +placeholder for visit diff --git a/example/DiscPack/DiscPack.ipynb b/example/DiscPack/DiscPack.ipynb new file mode 100644 index 00000000..1a6c152d --- /dev/null +++ b/example/DiscPack/DiscPack.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 37, + "metadata": {}, + "outputs": [], + "source": [ + "import random\n", + "import numpy as np\n", + "import matplotlib.pylab as plt\n" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "473.53100883248396" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Radius = 36\n", + "random.seed(23)\n", + "random.uniform(0,512)" + ] + }, + { + "cell_type": "code", + "execution_count": 203, + "metadata": {}, + "outputs": [], + "source": [ + "cx = [312, 70, 173, 250, 100, 30, 70, 210, 200, 112, 332, 346, 452, 431, 500, 480]\n", + "cy = [327, 64, 112, 20, 350, 450, 210, 438, 261, 500, 172, 455, 349, 60, 210, 460]\n", + "radius = [46, 42, 49, 48, 43, 48, 64, 56, 46, 30, 86, 62, 52, 38, 59, 41]\n", + "#for i in range(1,30):\n", + "# cx.append(random.uniform(0,512))\n", + "# cy.append(random.uniform(0,512))\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 244, + "metadata": {}, + "outputs": [], + "source": [ + "Lid = np.ones((512,512),dtype=np.uint8)\n", + "ID = np.ones((512,512),dtype=np.uint8)\n", + "for i in range(len(radius)):\n", + " Cx = cx[i]\n", + " Cy = cy[i]\n", + " Radius = radius[i]\n", + " for i in range(np.int(Cx-Radius),np.int(Cx+Radius)):\n", + " for j in range(np.int(Cy-Radius),np.int(Cy+Radius)):\n", + " if i>1 and i<512 and j>1 and j<512:\n", + " if ((i-Cx)*(i-Cx) + (j-Cy)*(j-Cy) < Radius*Radius):\n", + " ID[i,j] = 0\n", + " \n", + "for i in range(512):\n", + " ID[1,i] = 0\n", + " ID[2,i] = 0\n", + " ID[510,i] = 0\n", + " ID[511,i] = 0" + ] + }, + { + "cell_type": "code", + "execution_count": 245, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dfbAldX3n8ffH4UHChTAIjuMMcdBMiGCtTu6IZtlsZgIKGgzsGrNjLTqpYE3VFkatSEWI5VNcdjU1umsKqJVyNGxQbyaoC2GLIMG5a9yVIHdHIo9hFAKXGRkfZkquazDgd//ovkzP4T6ch374dZ/Pq+rUOafPOd3f/vXpz+nz69N9FBGYmVm3PKfpAszMrHwOdzOzDnK4m5l1kMPdzKyDHO5mZh3kcDcz6yCHu1lJJP2upK8t8tgvSJqTtKLuumw8OdzNcpKmJb2tinFHxCMRMRERT1cxfrNeDnczsw5yuFtpJL1H0mOSnpD0gKSz8+HPkXSZpG9L+oGknZJOLLzurZL+MX/sfZIelnRO/tgHJf2lpOvy8X5L0i9JulzSfkmPSnptYVw/L2mHpH15Lf9xvitkvttE0nZJByQ9JOl1+WNXAL8GXJl3n1yZD/9lSbdK+mE+T79TmNbzJN0o6UeS7gBeskTbrJMUko7I70/ntf2ffHp/lY/vs/n4viFpXeH1n8jn9UeSZiT9WuGxYyRdm8/TfZL+UNJs4fEXSvqCpO/l8/yOIRextYjD3Uoh6TTg7cArI+I44Fzg4fzhdwAXAr8OvBA4AFyVv+504Grg3wOrgZ8H1vSM/g3AnwMrgd3ALWTv3TXAHwOfLDz3WuAp4BeBDcBrgWJXy6uAB4CTgD8BdkhSRLwX+Fvg7Xn3ydslHQvcCnwOeD7wZuBqSWfk47oK+Ke87t/LL4PYArwln4+XAF8HPgOcCNwHfKDw3G8Ar8gf+xzwl5Kemz/2AWAd8GLgNcBF8y+S9Bzgr4C78umcDbxL0rkD1mptExG++DLyhSxM9wPnAEf2PHYfcHbh/mrgn4EjgPcDny889nPAT4Fz8vsfBG4tPP4GYA5Ykd8/DgjgBGAV8CRwTOH5bwZ25bd/F9jTM60AXpDfnwbeVnj83wF/2zMvnyQL0xX5PPxy4bH/BHxtkfZZl0/riMK03lt4/GPAzT3z+c0l2vsA8PL89neAcwuPvQ2YzW+/Cnik57WXA59p+j3jS7WXI/r9EDBbSkTskfQusjA+Q9ItwB9ExF7gRcCXJP2s8JKnycL4hcCjhfH8P0k/6Bn944XbPwG+H4d2TP4kv57Ix3UksE/S/POfUxw/8N2eac2/diEvAl4l6WBh2BFk3yJOzm8Xx/2Pi4xnMb3z1Xv/mbokvZsstF9I9iFxPNm3D+hpw57bLwJe2DMPK8i+pViHOdytNBHxOeBzko4n28L9KFm3w6PA70XE/+59jaR9wGmF+8cAzxuyhEfJttxPioinhnh97ylSHwX+V0S8pveJeT/+U8ApwP354F8YYprLyvvX30PWpXJPRPxM0gFg/hNsH7AWuDe/f0rh5Y8CD0XE+ipqs3S5z91KIek0Sb8h6WiyfuifkG2dA/w34ApJL8qfe7KkC/LHrgfeIOlfSjoK+BCHQmsgEbEP+DLwMUnH5ztyXyLp1/scxeNk/dbzbgJ+SdJbJB2ZX14p6aX5N4cvAh+U9HP5voOtw9Tdh+PIPki+Bxwh6f1kW+7zdgKXS1opaQ3Zvo95dwA/ynd2HyNphaSXSXplRbVaIhzuVpajgY8A3yfr+ng+8Ef5Y58AbgS+LOkJ4HayvmAi4h7g94Epsi3QJ8j67p8cso63AkeRbcUeIPvwWN3naz8B/Hb+q5M/jYgnyHbIbgH25vP10XxeIQvRiXz4n5HtDK3CLcDNwD+Qdf38E4d3vfwxMAs8BPwN2Tw/CZB/CL2BbGfsQ2TL51NkO66twxThP+uwdEiaAA4C6yPioabraSNJ/wHYEhH9fmOxDvKWuzVO0hvyro1jge3Atzj0M0pbhqTVks7Ku6FOA94NfKnpuqxZDndLwQVk3R57gfVkW53+Stm/o8h2YD8BfAW4gezYARtj7pYxM+sgb7mbmXVQEr9zP+mkk2LdunUjjePHP/4xxx57bDkFVWjf3r18d9++pstY1gtWr25FndCeWttSJ7Sn1sXq3DA52UA1S6sio2ZmZr4fEScv9FgS4b5u3TruvPPOkcYxPT3Npk2byimoQv/lYx/jfZde2nQZy7r03e9uRZ3QnlrbUie0p9bF6hw1T6pQRUZJWvSo6L7CXdLDZDtrngaeioiNys7q9xdk58x4GPidiDiQP/9y4OL8+e+IiFtGqL9TNkxOMteC/RzT09OV1DmhoY5PSs4wbVNVm1ZhmFrrWLa9NbWpTes2SJ/75oh4RURszO9fBtyWH9Z8W35//ix/W4AzgPPIzqLnf58x66gJqbYP7a5sHNRhlB2qF5CdXpX8+sLC8KmIeDI/CGUPcOYI0zEzDoXoUhezeX39FFLSQ2SHcgfwyYi4RtLBiDih8JwDEbEy/5OD2yPiunz4DrJTmV7fM85twDaAVatWTU5NTY00I3Nzc0xMLHZyv3SMe527Z2ZKH+eatWt5bHZ2+SeWaJgddsO26bBtNspOxX5rrWJ59mN+3tqyPkE1tW7evHmm0JtymH53qJ4VEXslPR+4VdL9Szx3oc2HZ32CRMQ1wDUAGzdujFF3NLRlh+q413n+5s2lj/PD27fXvvNv2D73ftu07K3wQevtt9Yqlme/5iJasz5B/et+X90y+Tm5iYj9ZIc1nwk8Lmk1ZIc/k53sCbITGBVPObqW7MhDM1tGVd0rVYzX3UBpWzbcJR0r6bj522Rnybub7Cx/86c43Up2yDP58C2SjpZ0Ktnh5HeUXbhZ19QRlg7k8dHPlvsq4GuS7iIL6f8ZEX9NdnrX10h6kOx/Gz8Cz5zCdSfZKVf/Grik8K85Ztaj7p2hZUwvhQ+JFGpI2bJ97hHxHeDlCwz/Adk/wyz0miuAK0auzixBZf6u2gFlVfG5Zcwa0nSw++eT3eZwt1r5aEKzejjczRqQ0hZzSrVYeRzuZgMo45tHimGaYk02Goe71c5dM2bVc7ib1SjlLeSUa7PBOdytEW3cem9jzTa+HO5mfSgj2Js6ydYgvPXeHUn8E5ONbqmVMtUtzrkIh4lZRbzl3lKDnMc75fN9p/rBU9TVX8iMIoXllkINKXO4t8yoIZ1iyKe8kqZcm9lSHO4tUmYopxbwlo5+3xv+4Eubw70lqjrHdypSDIoUa7KMl83yHO4tUGUIO+AXllItKWuinbxs+uNwT9y4/YHDXESjK2/T02+jOtvLy6Z/DveE1f0HDinxFmG71NF2Xj6Dcbhbsupamb21Xo4q29DLZ3A+iMmeMSEltxIV6yn720Vq89oF821a5rLychqOwz1RqXWTpKCs4HBYVK+Mo4+9nEbjcLfDpLj13mvQrfnU56erRgl4L7PROdwT5K32/s1FMD097TAo0XybljUua4Z3qJrVwCFndXO4m5l1kMPdrCYbJiebLmFZ/obRHQ53M7MOcribmXWQw92sRil3e6Rcmw3O4W5mDvYOcrgnyCtat3n5Wh0c7nYYB089UmrnlGqx8vgI1USVcW6OVPmcI5kUlnFX2tKezVvu9oyqVvT5P+Uu68+5yx5fk5r+YxLrLoe7Vaau8N09M9PqkK/7fPI+f/14cLdMwur82l7Wyt5kyM5P28G1OLfN+Oh7y13SCkm7Jd2U3z9R0q2SHsyvVxaee7mkPZIekHRuFYWPi7b8fVlKXSQp1TKI+S3qKpZ53VvrvV1n/V6sPIN0y7wTuK9w/zLgtohYD9yW30fS6cAW4AzgPOBqSSvKKXc8pf73ZamulG0OjLLCuKlQb+r1dkhf4S5pLfCbwKcKgy8Ars1vXwtcWBg+FRFPRsRDwB7gzHLKHV+pfp1uw4rYhhoXM2w4ty3Uqx7fOFL08QaQdD3wn4HjgEsj4nxJByPihMJzDkTESklXArdHxHX58B3AzRFxfc84twHbAFatWjU5NTU10ozMzc0xMTEx0jjqMGqdu2dmSqljuTMULldnWXWUYc3atTw2O7vs85o+K2Nb3qPQX611vwcWWn5da9NBbd68eSYiNi702LI7VCWdD+yPiBlJm/qY3kIft8/6BImIa4BrADZu3BibNvUz6sVNT08z6jjqMGqd86+t+u/LlqoztS2qD2/fzvsuvbTv5zf1Lagt71FYvtam3gO9y65LbVq2fn4tcxbwW5JeDzwXOF7SdcDjklZHxD5Jq4H9+fNngVMKr18L7C2zaBvuz6K73L8+iDb8T2zKmv5FlJddf5btc4+IyyNibUSsI9tR+pWIuAi4EdiaP20rcEN++0Zgi6SjJZ0KrAfuKL1yAw7/hcVyl1F1IdjndWle6pJKP3gqdaRulN+5fwTYKeli4BHgTQARcY+kncC9wFPAJRHx9MiVWqO6uDJ5K7B/XVz+XTfQEaoRMR0R5+e3fxARZ0fE+vz6h4XnXRERL4mI0yLi5rKLtnp1ecXu8rx1nZfd0nz6AVuSV6DxlnoXSNtPPVElh7stalxWmnGZTxsvDndbUEq/Y6+DA/7Z2tQmbaq1Lg53s5wD4pA2frh7+R3O4W5m1kEOd3uWcd4CGud5n9fmNmhz7WVzuNthvHJY2/k9nHG4m/VwOFgXONztGQ4183ugOxzuZgtwyFnbOdzNzDrI4W6At1TN74GuGeWskGadVvdZI5cLV5/B0gbhcDdr0CBby8XnOuhtOQ5389fxmpXR3vPjcMjbYtznbraEsj/4Uh+fdYfD3awGE1JlJ+NK/Zzr1gyH+5hzKFSvrjb2srQih7tZheoOXAe8zXO4m1WkqaB1wBs43M06yQFvDnezCjhcrWkOd7OSpRLsqdRhzXC4m5l1kMPdrMMG2Xr30a7d4nA3K5G7QiwVDncze4a33rvD4W5m1kEOd7OSpNolk2pdVfG3j4zDfcx5RbBebX5PtLn2sjnczcw6yOFutoRx3RLcMDnZdAkDG9dltRiHu5ktqE1h2aZa67JsuEt6rqQ7JN0l6R5JH8qHnyjpVkkP5tcrC6+5XNIeSQ9IOrfKGbDRecUw655+ttyfBH4jIl4OvAI4T9KrgcuA2yJiPXBbfh9JpwNbgDOA84CrJa2oonizKvlDrx1t0IYam7BsuEdmLr97ZH4J4ALg2nz4tcCF+e0LgKmIeDIiHgL2AGeWWrWVziuILSbl90Yb9w3URdHHgsu3vGeAXwSuioj3SDoYEScUnnMgIlZKuhK4PSKuy4fvAG6OiOt7xrkN2AawatWqyampqZFmZG5ujomJiZHGUYeU6yz+x+eatWt5bHa2wWr6V1WtwwTHUv+T2lSbDjMfC71Pq/oP2GFtmJxMen3qVUWtmzdvnomIjQs+GBF9X4ATgF3Ay4CDPY8dyK+vAi4qDN8BvHGp8U5OTsaodu3aNfI46pByncfCM5ePb99+2P2UL1XUWkYb1lFnVfOy2Pu06WXdO08pr0+9qqgVuDMWydUjBvmUiIiDkqbJ+tIfl7Q6IvZJWg3sz582C5xSeNlaYO8g0xkHgx41mPJXYxsfcxGNHvHq9aB//fxa5mRJJ+S3jwHOAe4HbgS25k/bCtyQ374R2CLpaEmnAuuBO8ouvK12z8wMtXJMSM9cquIVJzNKO6TWhlXUMxfRyHym1rap62fLfTVwbd7v/hxgZ0TcJOnrwE5JFwOPAG8CiIh7JO0E7gWeAi6JiKerKX88zQd8VSvuuJ2LxIYz//6r+v3iUB/OsuEeEX8PbFhg+A+Asxd5zRXAFSNXZ0uakPzGL1kZ7ZnKB2Rd742qQt7v7dH4CNWWq6KrxivV6MaxDee7a0bt1mqq26drBtqhatZlXQqUpuel6embt9w7o+yt93E7OKSq/RdNcLAaONw7pewumnEJiSrns+42HJdlZstzuNuSHBajq3vHphk43G3MdSV4HezWy+HeQVX8eqaL4VH3PG2YnCx9ml1dNjY6h3tHpfA765Q1GYhlBLJD3Zbjn0Ja3+o6IrFqqYRibx3LtWsqdVs7ONw7rKojWFM5AnNQqYdj6vUtpp/3Qlvnrc0c7jaUtga8lWPQZd/7fId99RzuNrS2dNM4SMpR5nIujsvLpxreoWojS3Xl9E7H8lT5AZ76xkFbOdytFKkFaUq1tFnV/yFQnI6Vy90yVqpiqNa9wjrQy1X38qvyfwrGkbfcrTJ1bc1XcXDQOJuQGv0zbG/Fl8PhbpWrIuR93u9uc8CPzt0yVhsHcfpSClX/09hoHO5mPQYJuC6FT0rBPs8BPzyHu1lumHDryk7AFIPdRuNwt7FWVqj5oJzqeOt9ON6hamOrqq3Vtm0Ft6HeNtSYGoe7jaWqw6Kug39G1YYabTgO9w7zV9lnqzt0HZ7lcVsOxuFuNqYclt3mcO8ob7U/W1Nh5hC1JjjcO8jB/mxNB2zT07fx43A3s9bwh2T/HO4dU9ZW++6ZmWd2Pi53SV0qNTZ5Mq5eqbSJVcfh3iGjBvuwgd2moG+a28jq4nDviFGCvcxgTi28UqvHrC4+/UAHDBvsVR+h6R27Zs3xlnvLpRbsdU/DzBbmcG+xlIO9OC3/vvxwqdZl3bJsuEs6RdIuSfdJukfSO/PhJ0q6VdKD+fXKwmsul7RH0gOSzq1yBsbNqP9A5KA1Gw/9bLk/Bbw7Il4KvBq4RNLpwGXAbRGxHrgtv0/+2BbgDOA84GpJK6oovo2G/b/PMv5SzgFrNj6W3aEaEfuAffntJyTdB6wBLgA25U+7FpgG3pMPn4qIJ4GHJO0BzgS+XnbxbTaOOxt9Xm6z+igGWNkkrQO+CrwMeCQiTig8diAiVkq6Erg9Iq7Lh+8Abo6I63vGtQ3YBrBq1arJqampkWZkbm6OiYmJkcZRh6bqHPQAmjVr1/LY7GwltWyYnCx1fIu1aUoHDcHhbVp2GwxjqfapcvmPqth2bVnvoZpaN2/ePBMRGxd8MCL6ugATwAzwb/P7B3seP5BfXwVcVBi+A3jjUuOenJyMUe3atWvkcdShiTqPhYEvH9++fajX9XMp22JtWlX9w16KbZqCppZ/me+ftqz3EdXUCtwZi+RqX7+WkXQk8AXgsxHxxXzw45JW54+vBvbnw2eBUwovXwvs7Wc6Vr4U+9lTrMnawd16/evn1zIi2/q+LyI+XnjoRmBrfnsrcENh+BZJR0s6FVgP3FFeydYFDnizavVzhOpZwFuAb0n6Zj7sj4CPADslXQw8ArwJICLukbQTuJfslzaXRMTTpVdutoy5iCQ/RLz1aXXo59cyXwMWW0POXuQ1VwBXjFCXmVUs1Q8/K4ePUO0wr7jWJf7GMxiHu3WaA2Fpbp/ucrhbY8bxm0WKYZpiTb3aUGNqHO5mljQH+3Ac7tZ5qYRDCkelLiaVNrLyONxtLDQdXk1Pvx8p1phiTW3hcDezZ6QUpinV0kYOdxsbTYWFQ2pwbrPROdxtrJRxXvxBp9c2cxGN7R+oe/l0mcPdrCJtD6m66297e6Wmn3PLmFWiyZV5ftpV/Na+SyFVZTv1TsPK5S33DvNKs7wyuwG63KVQ1bx1tb1S4C13M0bfQh2XkCprS35c2qtJDnezAodOf9xO6XO3jDXC4WBWLYd7xzlEzcaTw30MpBbwqdVj1kUOd6uVg92sHg73MeFQNRsvDvcx0nTANz39NpuQBr50Ue887p6ZGYv5HobDfcz45FntMkpgdSnsBpmXrszzqBzuVjkH++DKDOa2h13b62+Kw30M1XWYfB3TWeyreZu/qldRb1vboW01p8RHqI6xuYjKVp4qQ32YmouvSfmbRNVhNiElPf/zHOqjc7iPubLP+ldVcJS5ss+PK7WQqyvQ2hLwNhp3yxhwqAtlmK6UYV/Xjyq/mqfUbVN3DSnM82JSrq1NvOVuC9owOdno1l1TYdfEPDvMDnFblMdb7pacJlfwcdqCdpB2m8PdkpFSF0kKddRhXOZzHDncLQkOGbNyOdytcakGe9Vb8KnMdyp1WLkc7taoNgRLG2o06+Vfy3REvwHk3zebjYdlt9wlfVrSfkl3F4adKOlWSQ/m1ysLj10uaY+kBySdW1XhNtzvtMf5t92jKLvWNs27tVM/3TJ/BpzXM+wy4LaIWA/clt9H0unAFuCM/DVXS1pRWrUGlNcX3GTItzHc2lizja9lwz0ivgr8sGfwBcC1+e1rgQsLw6ci4smIeAjYA5xZUq1GdSeVMrNuUfTRBytpHXBTRLwsv38wIk4oPH4gIlZKuhK4PSKuy4fvAG6OiOsXGOc2YBvAqlWrJqempkaakbm5OSYmJkYaRx2GrXP3zEwF1Rxuw+TkM7erbM+y52XN2rU8Njtb6jiXUmynQRTbtI7lOYjeeWpqfRq0XRZb9sMuoypV0aabN2+eiYiNCz4YEctegHXA3YX7B3seP5BfXwVcVBi+A3jjcuOfnJyMUe3atWvkcdRh2DqPhVouo9bZxHx8fPv22tqn2EaDKrZpnfUOM09Nrk9lLPsUVdGmwJ2xSK4O+1PIxyWtBsiv9+fDZ4FTCs9bC+wdchqWq7PbxF00y3MbWRsMG+43Alvz21uBGwrDt0g6WtKpwHrgjtFKHG9NBInDy5rkn+uWY9nfuUv6PLAJOEnSLPAB4CPATkkXA48AbwKIiHsk7QTuBZ4CLomIpyuqvfO6FrJdmh+fE71aVf6RzLhYNtwj4s2LPHT2Is+/ArhilKKsebtnZti0aVPTZXSWw2t5bqPR+PQDZpasuv7vt4sc7olKYYulir+2s0NSCa1U6ljKICHfhvmpg88tYzYk97vXr7e9p6envQwW4S33BKW0lZtSLVY+B2N3ecvdrEFd3mm43Hz5g6VaDnezhjUV8GWH66Dz0Pt8h3253C1jloC6g63M6XXhLKVd5HA3GzNlB3vZHPDlcLeMWSLmQ7fKcCsr2KsO4Pnxu6tmeN5yN0tMVYHWlmBvalpd43A3S1CZAV/mUZ4O2/Zwt0yCUvp5nL8WN6fY9sO8H5r+NUzZ0/V7cTAOd7Mh1Rk2/fbHV1VTKhsb1j+Hu1mLjPPWq0/3MBj3uScqhTdx2f2+1k4pbbWnVEvqHO5mtiiHaXs53BPW5NZuiv8ebwb+wOmXw91q06WumS7Ni3WTwz1xTYSIg8us/RzuLdDET+5scePSRil3f6RcWyoc7i1RR6B0ZRpVanv9Nj4c7i1SZbD424FZtzjcW6bsf4Nv6t/l2xjwbazZxpfDvaVGDZqmQt3M6uHTD7TYMCeWSinQUzpB2nJSajezfjjcO6Kt4dOGgG9r29p4c7eMNS7l8Ey5NrOlONzNFuFgtzZzuFsSUtvBm1ItZsNwuFtSUgjVFGowG5XD3ZLTZLg62A9JuS1Sri0V/rWMJanfv5Urc1pmXeJwt6RVGfIOdesyh7u1Qm8QDxP2DvPBpXgcgpdjfyrrc5d0nqQHJO2RdFlV07HxNP/rmrkINkxOHnZ/sYu1n5dj/yoJd0krgKuA1wGnA2+WdHoV0zKzajlQ26mqbpkzgT0R8R0ASVPABcC9FU2vNXbPzHD+5s1Nl7GsD2/f3oo6oT21tqVOSLfW3i6ixer0BxIoKmgESb8NnBcRb8vvvwV4VUS8vfCcbcC2/O5pwAMjTvYk4PsjjqMOrrN8bam1LXVCe2ptS51QTa0vioiTF3qgqi33hfbAHPYpEhHXANeUNkHpzojYWNb4quI6y9eWWttSJ7Sn1rbUCfXXWtUO1VnglML9tcDeiqZlZmY9qgr3bwDrJZ0q6ShgC3BjRdMyM7MelXTLRMRTkt4O3AKsAD4dEfdUMa2C0rp4KuY6y9eWWttSJ7Sn1rbUCTXXWskOVTMza5ZPHGZm1kEOdzOzDmpduEt6k6R7JP1M0saexy7PT3fwgKRzC8MnJX0rf+xPpfpPlpHa6RgkfVrSfkl3F4adKOlWSQ/m1ysLjy3YtjXUeYqkXZLuy5f7OxOu9bmS7pB0V17rh1KtNZ/2Ckm7Jd2UeJ0P5+vvNyXdmWqtkk6QdL2k+/P36682WmdEtOoCvJTsoKdpYGNh+OnAXcDRwKnAt4EV+WN3AL9K9vv7m4HX1VzziryeFwNH5XWe3nA7/mvgV4C7C8P+BLgsv30Z8NHl2raGOlcDv5LfPg74h7yeFGsVMJHfPhL4O+DVKdaaT/8PgM8BN6W6/PPpPwyc1DMsuVqBa4G35bePAk5oss7WbblHxH0RsdDRrBcAUxHxZEQ8BOwBzpS0Gjg+Ir4eWav+d+DCGkuGwukYIuKnwPzpGBoTEV8Fftgz+AKyNyj59YWF4c9q25rq3BcR/ze//QRwH7Am0VojIubyu0fml0ixVklrgd8EPlUYnFydS0iqVknHk20w7QCIiJ9GxMEm62xduC9hDfBo4f5sPmxNfrt3eJ0Wqy01qyJiH2ShCjw/H55E/ZLWARvItoiTrDXv6vgmsB+4NSJSrfW/An8I/KwwLMU6IfuA/LKkmfy0JZBerS8Gvgd8Ju/q+pSkY5usM8nzuUv6G+AFCzz03oi4YbGXLTAslhhepxRqGEXj9UuaAL4AvCsifrTEbpNGa42Ip4FXSDoB+JKkly3x9EZqlXQ+sD8iZiRt6uclCwyrc/mfFRF7JT0fuFXS/Us8t6lajyDr5vz9iPg7SZ8g64ZZTOV1JhnuEXHOEC9b7JQHs/nt3uF1asvpGB6XtDoi9uXdWfvz4Y3WL+lIsmD/bER8MeVa50XEQUnTwHmkV+tZwG9Jej3wXOB4SdclWCcAEbE3v94v6Utk3Rep1ToLzObf1ACuJwv3xursUrfMjcAWSUdLOhVYD9yRfxV6QtKr81/JvBVYbOu/Km05HcONwNb89lYOtdOCbVtHQfky2wHcFxEfT7zWk/MtdiQdA5wD3J9arRFxeUSsjYh1ZO/Fr0TERanVCSDpWEnHzd8GXgvcnVqtEfFd4FFJp+WDziY7xXlzddaxF7nMC/BvyD71ngQeB24pPPZesr3OD1D4RQywkewN8W3gSvIjc2uu+/Vkv/T4Nln3UtPt+HlgH/DPeXteDDwPuA14ML8+cbm2raHOf0X2dfXvgW/ml9cnWuu/AHbnteNEedIAAABdSURBVN4NvD8fnlythelv4tCvZZKrk6wv+678cs/8upNora8A7syX//8AVjZZp08/YGbWQV3qljEzs5zD3cysgxzuZmYd5HA3M+sgh7uZWQc53M3MOsjhbmbWQf8fBEplud2yjAwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(1)\n", + "plt.title('segmented image')\n", + "plt.pcolormesh(ID[:,:],cmap='hot')\n", + "plt.grid(True)\n", + "plt.axis('equal')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 246, + "metadata": {}, + "outputs": [], + "source": [ + "Data = np.transpose(np.array([Lid,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,ID,Lid]))" + ] + }, + { + "cell_type": "code", + "execution_count": 248, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df/BddX3n8efL8EOaLzRB5Gsg1KBNqeCs0i+iXbetKShotWHXdTfOoulUJzM72OqoU6HOWu1sdrWTZmoHmNUx2uyifidFXVJ2UCkma92VRr4FKhAp0Vj4kpj4Ixn5ujYWfO8f53zlcPn+uD/Oj8859/WYuXPvPefec97n1+ue+zn3nKuIwMzMuuUZTRdgZmblc7ibmXWQw93MrIMc7mZmHeRwNzPrIIe7mVkHOdzNSiLpdyR9ZZF+vyBpTtKKuuuy8eRwN8tJ2ivprVUMOyIejoiJiHiiiuGb9XK4m5l1kMPdSiPpPZIelfSYpAclXZZ3f4akayV9U9L3Je2SdGbhfW+W9I95v/8k6duSLs/7vV/SX0q6KR/u1yX9kqTrJB2V9IikVxWG9fOSdkg6nNfyn+ebQuabTSRtk3RM0kFJr877bQV+Dbg+bz65Pu/+y5Jul/SDfJr+XWFcz5K0W9IPJe0Dnr/EvFknKSSdlD/fm9f2f/Px/VU+vE/mw/uapHWF9384n9YfSpqR9GuFfqdJ2plP035JfyBpttD/HEmfkfTdfJp/f8hFbC3icLdSSLoAeBvwkog4HbgC+Hbe+/eBq4DfAM4BjgE35O+7ELgR+A/AGuDngXN7Bv864H8Aq4G7gS+QrbvnAn8MfKTw2p3A48AvAhcDrwKKTS0vBR4EzgL+BNghSRHxXuBvgLflzSdvk7QSuB34FHA28EbgRkkX5cO6AfinvO7fzW+D2AS8KZ+O5wNfBT4BnAnsB/6o8NqvAS/O+30K+EtJz8z7/RGwDnge8Erg6vk3SXoG8FfAvfl4LgPeIemKAWu1tokI33wb+UYWpkeBy4GTe/rtBy4rPF8D/DNwEvA+4NOFfj8H/AS4PH/+fuD2Qv/XAXPAivz56UAAq4BJ4ARwWuH1bwT25I9/BzjQM64AnpM/3wu8tdD/3wN/0zMtHyEL0xX5NPxyod9/Ab6yyPxZl4/rpMK43lvo/6fAbT3Tec8S8/sY8KL88beAKwr93grM5o9fCjzc897rgE80vc74Vu3tpH4/BMyWEhEHJL2DLIwvkvQF4J0RcQh4LvA5ST8tvOUJsjA+B3ikMJz/J+n7PYM/Unj8Y+B78eSByR/n9xP5sE4GDkuaf/0zisMHvtMzrvn3LuS5wEslHS90O4nsW8Sz88fFYf/jIsNZTO909T7/WV2S3kUW2ueQfUicQfbtA3rmYc/j5wLn9EzDCrJvKdZhDncrTUR8CviUpDPI9nA/RNbs8AjwuxHxf3rfI+kwcEHh+WnAs4Ys4RGyPfezIuLxId7fe4nUR4D/HRGv7H1h3o7/OHAe8I288y8MMc5l5e3r7yFrUrk/In4q6Rgw/wl2GFgLPJA/P6/w9keAgxGxvoraLF1uc7dSSLpA0m9KOpWsHfrHZHvnAP8N2Crpuflrny1pY97vZuB1kv6lpFOAD/BkaA0kIg4DXwT+VNIZ+YHc50v6jT4HcYSs3XrercAvSXqTpJPz20skvSD/5vBZ4P2Sfi4/drB5mLr7cDrZB8l3gZMkvY9sz33eLuA6SaslnUt27GPePuCH+cHu0yStkPRCSS+pqFZLhMPdynIq8EHge2RNH2cDf5j3+zCwG/iipMeAO8nagomI+4HfA6bJ9kAfI2u7PzFkHW8GTiHbiz1G9uGxps/3fhj4t/mvTv48Ih4jOyC7CTiUT9eH8mmFLEQn8u5/QXYwtApfAG4D/oGs6eefeGrTyx8Ds8BB4K/JpvkEQP4h9Dqyg7EHyZbPx8gOXFuHKcJ/1mHpkDQBHAfWR8TBputpI0n/EdgUEf1+Y7EO8p67NU7S6/KmjZXANuDrPPkzSluGpDWSXp43Q10AvAv4XNN1WbMc7paCjWTNHoeA9WR7nf5K2b9TyA5gPwZ8CbiF7NwBG2NuljEz6yDvuZuZdVASv3M/66yzYt26dSMN40c/+hErV64sp6AKHT50iDXnnNN0GT9z98zMgt2fs2YN3zl8eKRhXzw1NdL7+9WWZb9cnYstizIMuizaOk/n52Fd694gqpinMzMz34uIZy/Ys+lTZCOCqampGNWePXtGHkYdtm/b1nQJP7MSFr1t37Ztyf793OrSlmW/XJ2jzu8yl0db52nd694gqpinwF2xSK721Syj7Cp9X5d0j6S78m5n5lfLeyi/X114/XWSDuRX0fMFisyWMCExoaHO2xp4PDY+Bmlz3xARL46IS/Ln1wJ3RHZa8x358/mr/G0CLgKuJLuKnv99Zgw5TJZWV6g3PU5rxigHVDeSXV6V/P6qQvfpiDgR2UkoB4BLRxiPWec0HbBNj9+q19dPISUdJDuVO4CPRMRHJR2PiFWF1xyLiNX5nxzcGRE35d13kF3K9OaeYW4BtgBMTk5OTU9PjzQhc3NzTEwsdnG/dBw9coSzJyebLgNY+gDeuWvX8ujs7KL9+1XHga22LPv5Oqs8cDqoxZZP2+bpvJQPqFYxTzds2DBTaE15qsUa44s34Jz8/myyi/7/OnC85zXH8vsbgKsL3XcAr19q+D6g2oyqD6jWdWCrLct+z549lR80LWs51TlPR6nXB1RHPKAa2TW5iYijZKc1XwockbQGstOfyS72BNkFjIqXHF1LduahLWO+PbSfm7VPSnvsKRhlXfZ2sLxlw13SSkmnzz8mu0refWRX+Zu/xOlmslOeybtvknSqpPPJTiffV3bhXTFsYDvo2yX15dTEgd0y+ANzcf3suU8CX5F0L1lI/6+I+DzZ5V1fKekhsv9t/CAwfwnXXWSXXP08cE08+a85liszmFMN+bmKLm3R+y3m7pmZpL/1pLhsFlLXzzHLHk+q63/Tlj1DNSK+Bbxoge7fJ/tnmIXesxXYOnJ1HVXVijghDRSocxGt2ijK/DCE6j582mzQdWjQYVepTetyHZK4/ICVp4vBVeWHIVQ7rxw4Gc+H+vnCYR3V9MZUVmD6zM36VdFsYvVzuHdYPxtVFXutZQyz7nbUqtqCzZricO+4ugO+rGBvStnt+m1V1q9Q2j4f2szhbkA6bfQphEEKNZiNyuE+BvoNq1ECfi5i5A+IlEI1pVqaMuo88DxslsN9TAwS8HMRfV+bo4xQhzSDYJSzJ82a5p9C2qJSaaoxs8F5z90al/Kebsq1pczzrXkO9zGS4gaXYk292lBjVcZ52tvO4W6NcXCYVcfhbtYHfxBZ2zjcx4xDymw8ONytEW38kGljzTa+HO5mZh3kcDez0vkcieY53M3MOsjhbrVz27VZ9RzuZlYJN800y+FuZotyQLeXw33MeGM1Gw8OdzOrjHcmmuNwNyuZA+2pPD+a4XAfI97IbBCp/beuDcbhblYBh9nTVT1PPM+fyuFuZk/T798sDqqqAHawP53DfUyktPKnVEuVxmU6B1X2fPF8XpjDfQx45bdB1LG+lPHH6oP8kfs4cribVcgfrEsbJuTL+GAYByc1XYBVK9WNYC6iddeYGXZetmlam1pfUl1P28x77h3mDcYG4fWlWxzuHeSvreUro304ZanXZ4NzuHdMmzbSNtVahlQP/o3bchgXDvcG+Le+T2pDzV0+UzO1eqw8fYe7pBWS7pZ0a/78TEm3S3oov19deO11kg5IelDSFVUU3nZlNp20vRkm5dqrqC2V5ZVCDVadQfbc3w7sLzy/FrgjItYDd+TPkXQhsAm4CLgSuFHSinLK7Z75DX3Yn4N5A22vJped15vu6yvcJa0Ffgv4WKHzRmBn/ngncFWh+3REnIiIg8AB4NJyyu22YmAvd+uaFKeprpN56tTV9ceeTtHHgpZ0M/BfgdOBd0fEayUdj4hVhdcci4jVkq4H7oyIm/LuO4DbIuLmnmFuAbYATE5OTk1PT480IXNzc0xMTIw0jDocPXKEsycnmy5jWU3Nz7tnZgZ+z7lr1/Lo7GypdVRx8HO5eTrMtPdr0OkZdfkPOi3Dzu/eOufHm+LB6yq2qQ0bNsxExCUL9oyIJW/Aa4Eb88evAG7NHx/ved2x/P4G4OpC9x3A65cax9TUVIxqz549Iw+jDtu3bWu6hL40OT9XwkC37du2DfyepW5VGWSeNj0twyz/JurtrbPqZTiKKrYp4K5YJFf7OUP15cBvS3oN8EzgDEk3AUckrYmIw5LWAEfz188C5xXevxY41OcHkVmjZ3Sm0mQxX8cw86HuaShzWc0PK5Xl0GbLtrlHxHURsTYi1pEdKP1SRFwN7AY25y/bDNySP94NbJJ0qqTzgfXAvtIrt1JNSE+53T0z85Tndau7bTjVtuhBjsPUPQ1VrhtNrXddMsq1ZT4I7JL0FuBh4A0AEXG/pF3AA8DjwDUR8cTIlVol+t2Aiq+rO3Sr3shTDPXU1RW8E5KXz5AGCveI2AvszR9/H7hskddtBbaOWJtVZNQNs+6vzlUGvINjcHXvUTvgh+OrQo6ZtraPFscx6jQ4KIbXVFOJ2+IH53AfI1W2j9bdVFMc9yCvt+G5DbxdHO5joI6NsqmvznMR7N271wFesVSC3U00/fOFwzquzo0ylQCwcqW2XFOrJ1UOdyuVNzyrg9ez5TncO6zpg1/Wfl6W7eVw7yhvlGlZ6MSwpW62PM+npTncrRLe8DLDhrVD3kblcO8gh0LzygrnJvfm27AeVXklzbZzuJuVrMrzCcz65XC3yoxjGFU9zXXNU+8Rt5/DvWPGMVBTUGfTiZfxU3l+LMzhbtZCDjRbjsPdbEQ+n8BS5HA3G0HTAdv0+C1dDnezITlYLWUOdzOzDnK4mw0hpb32lGqxdDjczTrAAW+9HO4d4z8yqJ6D1NrA4W5mreYdmoU53K0y3ujq5W8UVuRwN7OnuXhqqukSbEQO9w7yHrONC38ILc7hbpXwB0z7eRm2m8O9o5rcMB0KVgevZ0tzuHeYV36z8eVwt1L5A6VbvDzby+HecXMRtW2gDoJuSnG5plhTahzuY6LqjcEbW/PGZRmMy3SOyuE+RqraKLyxdV8qyziVOtrA4T5mytw46mzySck4TjM0P91Nj79tTmq6AKtfcSMZ5pR1b2TpqfO4ShOXOfA6N7hl99wlPVPSPkn3Srpf0gfy7mdKul3SQ/n96sJ7rpN0QNKDkq6ocgJsNPN73723i6emFu1n4x02da8H4zyvR9FPs8wJ4Dcj4kXAi4ErJb0MuBa4IyLWA3fkz5F0IbAJuAi4ErhR0ooqirfxMiEteLt7ZuZp3cZJU+FXx0F6B/vwlg33yMzlT0/ObwFsBHbm3XcCV+WPNwLTEXEiIg4CB4BLS63axsowgT2uQV83H6RPl6KPmZjvec8AvwjcEBHvkXQ8IlYVXnMsIlZLuh64MyJuyrvvAG6LiJt7hrkF2AIwOTk5NT09PdKEzM3NMTExMdIw6nD0yBHOnpxsuoxlNT0/756Z6fu1565dy6Ozs8u+roqLTFVR5yCqunDWsMt/kPnRa5hp6a1zfvwpXlCsim1qw4YNMxFxyYI9I6LvG7AK2AO8EDje0+9Yfn8DcHWh+w7g9UsNd2pqKka1Z8+ekYdRh+3btjVdQl+anJ8rYaDb9m3b+n5tk/UOUmdT0zKvjOVfxzT01lnHvBlWFdsUcFcskqsD/VomIo5L2kvWln5E0pqIOCxpDXA0f9kscF7hbWuBQ4OMx8ZX1c0oE1LpX/mb+AVJG5ot2lBjl/Xza5lnS1qVPz4NuBz4BrAb2Jy/bDNwS/54N7BJ0qmSzgfWA/vKLty6pc728SrG5V+PWGr62XNfA+zM292fAeyKiFslfRXYJektwMPAGwAi4n5Ju4AHgMeBayLiiWrKN0tHHXvwDnbr17LhHhF/D1y8QPfvA5ct8p6twNaRq7Ox0NQvWtrWRONgt0H4DFVrVNM/Vawq4IvDL3N4Zv3ytWXMKjRKMPskHhuFw90a0/Re+7w62snnL+nQ72sd6jYqN8uYUU3zzEIc2lYX77lbI1LZazfrKoe7mVkHuVmmQ3xt9tHU1TRjVgeHe8uN2rxRfH9dweYmGbPquVmmhaq6nK0vk2spWuw6/gtdy9+e5HBvmbpW4FEu3WpWBgf2aNws0xJNrOTz43Q7tNXJgV4O77m3QNMre9Pjt/HgPfVyec89YSmt6P4liVUlpfW8Sxzu1jcHfDVGDTcvE1uIwz1BKe/JuB1+dFX8yqmoLcsm5fW8C9zmnphxWOHbEj5lKv5sr45xpb4epV5fFzjcbSjeOPvTZNCmGvIp1tRFDveEtG2lb1u9yynzG0VKwZpaLVYPh3si2rrSD1t3l5tmUl2WTYd8qvOlqxzuZpT3YdOGAGtDjTY6h3sC2r6xtb3+MjS9Vzyoumtt07zpCoe7NSaVpplR62hrcLXtA8kG43BvWFc2rnFse+9KOFY9DV2YR23kcLfGNRnwbf5wKZMDuHsc7ja2Rgn2LoZhFdPUxfnUFg73BnVtxR9leuYiat2LHnZcXWmKWUyXp23cONwtKXUEvJtibBw43C05VYXvqN8OxmWvdlyms+sc7g3xBrS0MptpyhjWuC2vMv5mcdzmWWoc7pa0UYPZTTDDczi3m6/nbq1QDOnlQseBbuZwtxbqDe+9e/dWGujjvAfrf99qLzfLmC1hnIPd2s3hbmZL8gdcOy0b7pLOk7RH0n5J90t6e979TEm3S3oov19deM91kg5IelDSFVVOgFlVHGrWZv3suT8OvCsiXgC8DLhG0oXAtcAdEbEeuCN/Tt5vE3ARcCVwo6QVVRRvZmYLWzbcI+JwRPxd/vgxYD9wLrAR2Jm/bCdwVf54IzAdESci4iBwALi07MLNrD7+FtM+igGOhEtaB3wZeCHwcESsKvQ7FhGrJV0P3BkRN+XddwC3RcTNPcPaAmwBmJycnJqenh5pQubm5piYmBhpGHU4euQIZ09OlnKSSJXOXbuWR2dnB37fxVNTFVSztCqWfRXLZ9h52oSFah1m2Va9ni82T5tYD5dTxXq6YcOGmYi4ZMGeEdHXDZgAZoB/kz8/3tP/WH5/A3B1ofsO4PVLDXtqaipGtWfPnpGHUYft27b97PFKSPa2fdu2od7XhCqWfUrzNKXln8J8rKLOOlSxngJ3xSK52tevZSSdDHwG+GREfDbvfETSmrz/GuBo3n0WOK/w9rXAoX7GY5YCN0FYF/TzaxmR7X3vj4jthV67gc35483ALYXumySdKul8YD2wr7ySzcxsOf2cofpy4E3A1yXdk3f7Q+CDwC5JbwEeBt4AEBH3S9oFPED2S5trIuKJ0ivvgLmITu0l+kzG7vMZq+2xbLhHxFeAxRLoskXesxXYOkJdZtYBXduBaROfodqwruwFdWU6HETWFQ53M6tUVz7428bhbiPzxmvL8TpSP4d7Arzim1nZHO5mZh3kcE9EW/fe21q31a/M/8W15TncE9K2Fb9t9ZqNE4e7DcXBbsPyulMPh3ti2vDVNfX6LH1eh6rncE9Uqit/qnVZ+7RhR6bNHO4JS23FT60e6wavV9Xo58Jh1qBUrs1RxQY46HQ5BLprftmmsK53hffcW6DpUCt7/BPSUBvxsO8bRNPzetx5/pfHe+4tUfeezcVTU6VuaGXWXRyWw6B7isvUe/LD8557y9RxEKqKPfWqeOOvV90fpvPruw++Ds7h3lJVrOxVDLOO8E39z8atPL1hP/8N0+H/dG6WabmFVuh+ArWODaHOver5cZUxXakcxDYbhcO9g1LYg2kqHP03cNXxfG0XN8tY6bqw1+sgs7ZzuFupUgj2FGroGn/YtY/D3TrJAW/jzuFupXGgmqXD4W6lSDHYR63JTREZz4d2crhbpzngbVw53M1sUf5way+Hu40sxSaZMjngrI0c7ma2IH+otZvD3TqvjG8W4xZ0F09NNV2CjcjhbtancQn4cZnOrnO4m5l1kMPdbABdv7Rsl6dt3DjczYbQxRDs4jSNM4e7mTnYO2jZcJf0cUlHJd1X6HampNslPZTfry70u07SAUkPSrqiqsLNmtaVJpouTIM9XT977n8BXNnT7VrgjohYD9yRP0fShcAm4KL8PTdKWlFatdaXCWnJm5WrreHYlQ8nW9iy4R4RXwZ+0NN5I7Azf7wTuKrQfToiTkTEQeAAcGlJtdoy+g1vB335HJKWGkUfK6WkdcCtEfHC/PnxiFhV6H8sIlZLuh64MyJuyrvvAG6LiJsXGOYWYAvA5OTk1PT09EgTMjc3x8TExEjDqMPRI0c4e3KytOGV+efQxRNXBpmfTf9B9blr1/Lo7Oyi/es+IWex+bFcnXXpZ360ZXvqrXN+3qd4ElYV83TDhg0zEXHJgj0jYtkbsA64r/D8eE//Y/n9DcDVhe47gNcvN/ypqakY1Z49e0YeRh22b9tW2rBWQum3eYPMzyrqGOS2fdu2vqapTsPUWfVtEG3ZnnrrbHKZL6eKeQrcFYvk6rB/kH1E0pqIOCxpDXA07z4LnFd43Vrg0JDjsCVU1aTiP5gux/w8TKXpq8xlOsg0eV1qzrA/hdwNbM4fbwZuKXTfJOlUSecD64F9o5VoRXW0lQ86fG/Ai2vyoOX8uMsa/zDrno/vNKefn0J+GvgqcIGkWUlvAT4IvFLSQ8Ar8+dExP3ALuAB4PPANRHxRFXFW3WabkcvSyofPHMRXDw1VUs9VXyglBHODvh6LdssExFvXKTXZYu8fiuwdZSibGF1bxyDNNHMRXjj7VMVTTZVfWiUvUzd7FefYdvcrWZNBWebAz71EEm9Ph/XaTdffsA6yeExmqqb5VLaCegqh3sLNL0htO3XESnU0GZ1rW9Nr9dd53C30jlcrV8O+Oo43K0STf78z8wc7slLZc9mmDrqDFpfBKscTaxvqazjXeNwt0rVEbgpXkekjRyy3eJwt8pVGfDeWzdbmH/nbrUohvCoe4hVBPpyNflDpFr+7Xv5HO5Wu2HP0GzylPriax1C1gYOd2tMUyE56jeH+fc75C1lDncbG1VcJwUc8pYmH1C1zqv6krP+lYmlyOFuVgIHvKXGzTLWWU1cIhncTGNp8J574lIJilTq6Jf3pNulbetXGzjczUrmDxZLgcO9BZreq2l6/INKIVxTqGFQbVvOtjSHu3VKSqGaUi0p84dKNRzu1hkO09E5aLvD4d4Svj56O/kDZ2lev6rjcG+RujcEb3jjycu9GxzuLVPXhudrpI+3qpe//1yleg73Fqp6o2jjRpdy80fKtS2ljeuBPcnh3lJVbXjeoK2oivXB61g9fPmBFhv2uuhLDcusV1nrmdexejncO2CUjc8bnPUrlT9Zsf443Duk343PG5uNop/1zOtY8xzuHeQNy+rg9SxtPqBqZtZBDnczsw5yuJuZdZDD3cysgxzuZmYdVFm4S7pS0oOSDki6tqrxmEHav9xIuTbrrkrCXdIK4Abg1cCFwBslXVjFuMzM7Omq+p37pcCBiPgWgKRpYCPwQEXjM/MesllBVeF+LvBI4fks8NLiCyRtAbbkT+ckPTjiOM8CvjfiMOpw1jvf/e5W1Ek75ie0p9a21AntqXXBOpXmlTirmKfPXaxHVeG+0Jx9ym5VRHwU+GhpI5TuiohLyhpeVVxn+dpSa1vqhPbU2pY6of5aqzqgOgucV3i+FjhU0bjMzKxHVeH+NWC9pPMlnQJsAnZXNC4zM+tRSbNMRDwu6W3AF4AVwMcj4v4qxlVQWhNPxVxn+dpSa1vqhPbU2pY6oeZaFf6FgZlZ5/gMVTOzDnK4m5l1UOvCXdIbJN0v6aeSLunpd11+uYMHJV1R6D4l6et5vz9XAz+CTe1yDJI+LumopPsK3c6UdLukh/L71YV+C87bGuo8T9IeSfvz5f72hGt9pqR9ku7Na/1AqrXm414h6W5JtyZe57fz7fceSXelWqukVZJulvSNfH391UbrjIhW3YAXABcAe4FLCt0vBO4FTgXOB74JrMj77QN+lez397cBr6655hV5Pc8DTsnrvLDh+fjrwK8A9xW6/Qlwbf74WuBDy83bGupcA/xK/vh04B/yelKsVcBE/vhk4G+Bl6VYaz7+dwKfAm5Ndfnn4/82cFZPt+RqBXYCb80fnwKsarLO1u25R8T+iFjobNaNwHREnIiIg8AB4FJJa4AzIuKrkc3V/w5cVWPJULgcQ0T8BJi/HENjIuLLwA96Om8kW0HJ768qdH/avK2pzsMR8Xf548eA/WRnQKdYa0TEXP705PwWKdYqaS3wW8DHCp2Tq3MJSdUq6QyyHaYdABHxk4g43mSdrQv3JSx0yYNz89vsAt3rtFhtqZmMiMOQhSpwdt49ifolrQMuJtsjTrLWvKnjHuAocHtEpFrrnwF/APy00C3FOiH7gPyipJn8siWQXq3PA74LfCJv6vqYpJVN1pnkH2RL+mvgOQv0em9E3LLY2xboFkt0r1MKNYyi8folTQCfAd4RET9c4rBJo7VGxBPAiyWtAj4n6YVLvLyRWiW9FjgaETOSXtHPWxboVufyf3lEHJJ0NnC7pG8s8dqmaj2JrJnz9yLibyV9mKwZZjGV15lkuEfE5UO8bbFLHszmj3u716ktl2M4ImlNRBzOm7OO5t0brV/SyWTB/smI+GzKtc6LiOOS9gJXkl6tLwd+W9JrgGcCZ0i6KcE6AYiIQ/n9UUmfI2u+SK3WWWA2/6YGcDNZuDdWZ5eaZXYDmySdKul8YD2wL/8q9Jikl+W/knkzsNjef1XacjmG3cDm/PFmnpxPC87bOgrKl9kOYH9EbE+81mfne+xIOg24HPhGarVGxHURsTYi1pGti1+KiKtTqxNA0kpJp88/Bl4F3JdarRHxHeARSRfknS4ju8R5c3XWcRS5zBvwr8k+9U4AR4AvFPq9l+yo84MUfhEDXEK2QnwTuJ78zNya634N2S89vknWvNT0fPw0cBj453x+vgV4FnAH8FB+f+Zy87aGOv8V2dfVvwfuyW+vSbTWfwHcndd6H/C+vHtytRbG/wqe/LVMcnWStWXfm9/un992Eq31xcBd+cQ41IsAAAA7SURBVPL/n8DqJuv05QfMzDqoS80yZmaWc7ibmXWQw93MrIMc7mZmHeRwNzPrIIe7mVkHOdzNzDro/wPD8AlMiG6OYAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(1)\n", + "plt.title('segmented image')\n", + "plt.pcolormesh(Data[:,:,9],cmap='hot')\n", + "plt.grid(True)\n", + "plt.axis('equal')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 249, + "metadata": {}, + "outputs": [], + "source": [ + "Data.tofile(\"discs512x512x18.raw\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/hip/GreyscaleColor.cu b/hip/GreyscaleColor.cu index c9dafb37..aff7b39a 100644 --- a/hip/GreyscaleColor.cu +++ b/hip/GreyscaleColor.cu @@ -1450,9 +1450,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff,double alpha, double beta, - double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){ + double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ int n,nn,ijk,nread; int nr1,nr2,nr3,nr4,nr5,nr6; @@ -1462,8 +1462,6 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int double fq; // conserved momemnts double rho,jx,jy,jz; - //double vx,vy,vz,v_mag; - //double ux,uy,uz,u_mag; double ux,uy,uz; // non-conserved moments double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; @@ -1473,18 +1471,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int double C,nx,ny,nz; //color gradient magnitude and direction double phi,tau,rho0,rlx_setA,rlx_setB; - //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) double porosity; double perm;//voxel permeability - //double c0, c1; //Guo's model parameters double tau_eff; double mu_eff;//kinematic viscosity - double nx_gs,ny_gs,nz_gs;//grey-solid color gradient - double nx_phase,ny_phase,nz_phase,C_phase; double Fx,Fy,Fz; - double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18; - double gp3,gp5,gp7; double Fcpx,Fcpy,Fcpz;//capillary penalty force + double W;//greyscale wetting strength + double Sn_grey,Sw_grey; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -1510,9 +1504,9 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int porosity = Poros[n]; perm = Perm[n]; - nx_gs = GreySolidGrad[n+0*Np]; - ny_gs = GreySolidGrad[n+1*Np]; - nz_gs = GreySolidGrad[n+2*Np]; + W = GreySolidW[n]; + Sn_grey = GreySn[n]; + Sw_grey = GreySw[n]; // compute phase indicator field phi=(nA-nB)/(nA+nB); @@ -1534,98 +1528,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int //........................................................................ nn = ijk-1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 - gp1 = Psi[nn]; //........................................................................ nn = ijk+1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 - gp2 = Psi[nn]; //........................................................................ nn = ijk-strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 - gp3 = Psi[nn]; //........................................................................ nn = ijk+strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 - gp4 = Psi[nn]; //........................................................................ nn = ijk-strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 - gp5 = Psi[nn]; //........................................................................ nn = ijk+strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 - gp6 = Psi[nn]; //........................................................................ nn = ijk-strideY-1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 - gp7 = Psi[nn]; //........................................................................ nn = ijk+strideY+1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 - gp8 = Psi[nn]; //........................................................................ nn = ijk+strideY-1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 - gp9 = Psi[nn]; //........................................................................ nn = ijk-strideY+1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 - gp10 = Psi[nn]; //........................................................................ nn = ijk-strideZ-1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 - gp11 = Psi[nn]; //........................................................................ nn = ijk+strideZ+1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 - gp12 = Psi[nn]; //........................................................................ nn = ijk+strideZ-1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 - gp13 = Psi[nn]; //........................................................................ nn = ijk-strideZ+1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 - gp14 = Psi[nn]; //........................................................................ nn = ijk-strideZ-strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 - gp15 = Psi[nn]; //........................................................................ nn = ijk+strideZ+strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 - gp16 = Psi[nn]; //........................................................................ nn = ijk+strideZ-strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 - gp17 = Psi[nn]; //........................................................................ nn = ijk-strideZ+strideY; // neighbor index (get convention) m18 = Phi[nn]; // get neighbor for phi - 18 - gp18 = Psi[nn]; //............Compute the Color Gradient................................... - nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); - ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); - nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); - C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase); - - //correct the normal color gradient by considering the effect of grey solid - nx = nx_phase + (1.0-porosity)*nx_gs; - ny = ny_phase + (1.0-porosity)*ny_gs; - nz = nz_phase + (1.0-porosity)*nz_gs; - if (C_phase==0.0){ - nx = nx_phase; - ny = ny_phase; - nz = nz_phase; - } - - //...........Normalize the Color Gradient................................. - C = sqrt(nx*nx+ny*ny+nz*nz); - double ColorMag = C; - if (C==0.0) ColorMag=1.0; - nx = nx/ColorMag; - ny = ny/ColorMag; - nz = nz/ColorMag; + nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); + nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); //............Compute the Greyscale Potential Gradient..................... // Fcpx = 0.0; @@ -1648,14 +1605,24 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // //ny = Fcpy/Fcp_mag; // //nz = Fcpz/Fcp_mag; // } - Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); - Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); - Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); + Fcpx = nx; + Fcpy = ny; + Fcpz = nz; + double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz); + if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0; //NOTE for open node (porosity=1.0),Fcp=0.0 Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm); Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm); Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm); + //...........Normalize the Color Gradient................................. + C = sqrt(nx*nx+ny*ny+nz*nz); + double ColorMag = C; + if (C==0.0) ColorMag=1.0; + nx = nx/ColorMag; + ny = ny/ColorMag; + nz = nz/ColorMag; + // q=0 fq = dist[n]; rho = fq; @@ -2201,6 +2168,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2220,6 +2191,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2240,6 +2215,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2262,15 +2241,13 @@ __global__ void dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *neighborList, int //CP: capillary penalty // also turn off recoloring for grey nodes __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm, double *Velocity, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm, double *Velocity, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, - double Gx, double Gy, double Gz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){ + double Gx, double Gy, double Gz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ int ijk,nn,n; double fq; // conserved momemnts double rho,jx,jy,jz; - //double vx,vy,vz,v_mag; - //double ux,uy,uz,u_mag; double ux,uy,uz; // non-conserved moments double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18; @@ -2280,18 +2257,14 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis double C,nx,ny,nz; //color gradient magnitude and direction double phi,tau,rho0,rlx_setA,rlx_setB; - //double GeoFun=0.0;//geometric function from Guo's PRE 66, 036304 (2002) double porosity; double perm;//voxel permeability - //double c0, c1; //Guo's model parameters double tau_eff; double mu_eff;//kinematic viscosity - double nx_gs,ny_gs,nz_gs;//grey-solid color gradient - double nx_phase,ny_phase,nz_phase,C_phase; double Fx,Fy,Fz; - double gp1,gp2,gp4,gp6,gp8,gp9,gp10,gp11,gp12,gp13,gp14,gp15,gp16,gp17,gp18; - double gp3,gp5,gp7; double Fcpx,Fcpy,Fcpz;//capillary penalty force + double W;//greyscale wetting strength + double Sn_grey,Sw_grey; const double mrt_V1=0.05263157894736842; const double mrt_V2=0.012531328320802; @@ -2315,11 +2288,12 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // read the component number densities nA = Den[n]; nB = Den[Np + n]; + porosity = Poros[n]; perm = Perm[n]; - nx_gs = GreySolidGrad[n+0*Np]; - ny_gs = GreySolidGrad[n+1*Np]; - nz_gs = GreySolidGrad[n+2*Np]; + W = GreySolidW[n]; + Sn_grey = GreySn[n]; + Sw_grey = GreySw[n]; // compute phase indicator field phi=(nA-nB)/(nA+nB); @@ -2341,98 +2315,61 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis //........................................................................ nn = ijk-1; // neighbor index (get convention) m1 = Phi[nn]; // get neighbor for phi - 1 - gp1 = Psi[nn]; //........................................................................ nn = ijk+1; // neighbor index (get convention) m2 = Phi[nn]; // get neighbor for phi - 2 - gp2 = Psi[nn]; //........................................................................ nn = ijk-strideY; // neighbor index (get convention) m3 = Phi[nn]; // get neighbor for phi - 3 - gp3 = Psi[nn]; //........................................................................ nn = ijk+strideY; // neighbor index (get convention) m4 = Phi[nn]; // get neighbor for phi - 4 - gp4 = Psi[nn]; //........................................................................ nn = ijk-strideZ; // neighbor index (get convention) m5 = Phi[nn]; // get neighbor for phi - 5 - gp5 = Psi[nn]; //........................................................................ nn = ijk+strideZ; // neighbor index (get convention) m6 = Phi[nn]; // get neighbor for phi - 6 - gp6 = Psi[nn]; //........................................................................ nn = ijk-strideY-1; // neighbor index (get convention) m7 = Phi[nn]; // get neighbor for phi - 7 - gp7 = Psi[nn]; //........................................................................ nn = ijk+strideY+1; // neighbor index (get convention) m8 = Phi[nn]; // get neighbor for phi - 8 - gp8 = Psi[nn]; //........................................................................ nn = ijk+strideY-1; // neighbor index (get convention) m9 = Phi[nn]; // get neighbor for phi - 9 - gp9 = Psi[nn]; //........................................................................ nn = ijk-strideY+1; // neighbor index (get convention) m10 = Phi[nn]; // get neighbor for phi - 10 - gp10 = Psi[nn]; //........................................................................ nn = ijk-strideZ-1; // neighbor index (get convention) m11 = Phi[nn]; // get neighbor for phi - 11 - gp11 = Psi[nn]; //........................................................................ nn = ijk+strideZ+1; // neighbor index (get convention) m12 = Phi[nn]; // get neighbor for phi - 12 - gp12 = Psi[nn]; //........................................................................ nn = ijk+strideZ-1; // neighbor index (get convention) m13 = Phi[nn]; // get neighbor for phi - 13 - gp13 = Psi[nn]; //........................................................................ nn = ijk-strideZ+1; // neighbor index (get convention) m14 = Phi[nn]; // get neighbor for phi - 14 - gp14 = Psi[nn]; //........................................................................ nn = ijk-strideZ-strideY; // neighbor index (get convention) m15 = Phi[nn]; // get neighbor for phi - 15 - gp15 = Psi[nn]; //........................................................................ nn = ijk+strideZ+strideY; // neighbor index (get convention) m16 = Phi[nn]; // get neighbor for phi - 16 - gp16 = Psi[nn]; //........................................................................ nn = ijk+strideZ-strideY; // neighbor index (get convention) m17 = Phi[nn]; // get neighbor for phi - 17 - gp17 = Psi[nn]; //........................................................................ nn = ijk-strideZ+strideY; // neighbor index (get convention) m18 = Phi[nn]; // get neighbor for phi - 18 - gp18 = Psi[nn]; //............Compute the Color Gradient................................... - nx_phase = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); - ny_phase = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); - nz_phase = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); - C_phase = sqrt(nx_phase*nx_phase+ny_phase*ny_phase+nz_phase*nz_phase); - - //correct the normal color gradient by considering the effect of grey solid - nx = nx_phase + (1.0-porosity)*nx_gs; - ny = ny_phase + (1.0-porosity)*ny_gs; - nz = nz_phase + (1.0-porosity)*nz_gs; - if (C_phase==0.0){ - nx = nx_phase; - ny = ny_phase; - nz = nz_phase; - } - - //...........Normalize the Color Gradient................................. - C = sqrt(nx*nx+ny*ny+nz*nz); - double ColorMag = C; - if (C==0.0) ColorMag=1.0; - nx = nx/ColorMag; - ny = ny/ColorMag; - nz = nz/ColorMag; + nx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); + ny = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); + nz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); //............Compute the Greyscale Potential Gradient..................... // Fcpx = 0.0; @@ -2455,14 +2392,23 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // ny = Fcpy/Fcp_mag; // nz = Fcpz/Fcp_mag; // } - Fcpx = -3.0/18.0*(m1-m2+0.5*(m7-m8+m9-m10+m11-m12+m13-m14)); - Fcpy = -3.0/18.0*(m3-m4+0.5*(m7-m8-m9+m10+m15-m16+m17-m18)); - Fcpz = -3.0/18.0*(m5-m6+0.5*(m11-m12-m13+m14+m15-m16-m17+m18)); + Fcpx = nx; + Fcpy = ny; + Fcpz = nz; + double Fcp_mag=sqrt(Fcpx*Fcpx+Fcpy*Fcpy+Fcpz*Fcpz); + if (Fcp_mag==0.0); Fcpx=Fcpy=Fcpz=0.0; //NOTE for open node (porosity=1.0),Fcp=0.0 Fcpx *= alpha*W*(1.0-porosity)/sqrt(perm); Fcpy *= alpha*W*(1.0-porosity)/sqrt(perm); Fcpz *= alpha*W*(1.0-porosity)/sqrt(perm); + //...........Normalize the Color Gradient................................. + C = sqrt(nx*nx+ny*ny+nz*nz); + double ColorMag = C; + if (C==0.0) ColorMag=1.0; + nx = nx/ColorMag; + ny = ny/ColorMag; + nz = nz/ColorMag; // q=0 fq = dist[n]; @@ -2942,6 +2888,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {1,0,0}, {0,1,0}, {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nx; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*ux))+delta; b1 = nB*(0.1111111111111111*(1+4.5*ux))-delta; @@ -2958,6 +2908,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {0,1,0} delta = beta*nA*nB*nAB*0.1111111111111111*ny; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uy))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uy))-delta; @@ -2973,6 +2927,10 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis // Cq = {0,0,1} delta = beta*nA*nB*nAB*0.1111111111111111*nz; if (!(nA*nB*nAB>0)) delta=0; + //----------------newly added for better control of recoloring---------------// + if (nA/(nA+nB)>=Sn_grey && nA/(nA+nB) <= Sw_grey && porosity !=1.0) delta = 0.0; + if (nA/(nA+nB)>Sw_grey && porosity !=1.0) delta = -1.0*delta; + //---------------------------------------------------------------------------// if (RecoloringOff==true && porosity !=1.0) delta=0; a1 = nA*(0.1111111111111111*(1+4.5*uz))+delta; b1 = nB*(0.1111111111111111*(1+4.5*uz))-delta; @@ -2989,6 +2947,7 @@ __global__ void dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dis } } + __global__ void dvc_ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, double *Bq, int start, int finish, int Np){ int idx; double nA,nB; @@ -4572,12 +4531,12 @@ extern "C" void ScaLBL_PhaseField_InitFromRestart(double *Den, double *Aq, doubl //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi,double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel, double *Pressure, double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, - double Fx, double Fy, double Fz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){ + double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, Psi, GreySolidGrad, Poros, Perm, Vel, Pressure, - rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, W, strideY, strideZ, start, finish, Np); + dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor_CP<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm, Vel, Pressure, + rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ printf("hip error in ScaLBL_D3Q19_AAeven_GreyscaleColor_CP: %s \n",hipGetErrorString(err)); @@ -4587,12 +4546,12 @@ extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(int *Map, double *dist, do //Model-1 & 4 with capillary pressure penalty extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, - double *Phi, double *Psi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel,double *Pressure, + double *Phi, double *GreySolidW, double *GreySn, double *GreySw, double *Poros,double *Perm,double *Vel,double *Pressure, double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, - double Fx, double Fy, double Fz, bool RecoloringOff, double W, int strideY, int strideZ, int start, int finish, int Np){ + double Fx, double Fy, double Fz, bool RecoloringOff, int strideY, int strideZ, int start, int finish, int Np){ - dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, Psi, GreySolidGrad, Poros, Perm,Vel,Pressure, - rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, W, strideY, strideZ, start, finish, Np); + dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor_CP<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidW, GreySn, GreySw, Poros, Perm,Vel,Pressure, + rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, RecoloringOff, strideY, strideZ, start, finish, Np); hipError_t err = hipGetLastError(); if (hipSuccess != err){ @@ -4600,51 +4559,3 @@ extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(int *d_neighborList, int *M } } -//extern "C" void ScaLBL_Update_GreyscalePotential(int *Map, double *Phi, double *Psi, double *Poro, double *Perm, double alpha, double W, -// int start, int finish, int Np){ -// -// dvc_ScaLBL_Update_GreyscalePotential<<>>(Map, Phi, Psi, Poro, Perm, alpha, W, start, finish, Np); -// -// hipError_t err = hipGetLastError(); -// if (hipSuccess != err){ -// printf("hip error in ScaLBL_Update_GreyscalePotential: %s \n",hipGetErrorString(err)); -// } -//} - -////Model-2&3 -//extern "C" void ScaLBL_D3Q19_AAeven_GreyscaleColor(int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *Phi,double *GreySolidGrad, double *Poros,double *Perm,double *Vel, -// double rhoA, double rhoB, double tauA, double tauB,double tauA_eff,double tauB_eff, double alpha, double beta, -// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ -// -// //cudaProfilerStart(); -// //cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor, cudaFuncCachePreferL1); -// -// dvc_ScaLBL_D3Q19_AAeven_GreyscaleColor<<>>(Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm, Vel, -// rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); -// hipError_t err = hipGetLastError(); -// if (hipSuccess != err){ -// printf("hip error in ScaLBL_D3Q19_AAeven_GreyscaleColor: %s \n",hipGetErrorString(err)); -// } -// //cudaProfilerStop(); -// -//} -// -////Model-2&3 -//extern "C" void ScaLBL_D3Q19_AAodd_GreyscaleColor(int *d_neighborList, int *Map, double *dist, double *Aq, double *Bq, double *Den, -// double *Phi, double *GreySolidGrad, double *Poros,double *Perm,double *Vel, -// double rhoA, double rhoB, double tauA, double tauB, double tauA_eff,double tauB_eff, double alpha, double beta, -// double Fx, double Fy, double Fz, int strideY, int strideZ, int start, int finish, int Np){ -// -// //cudaProfilerStart(); -// //cudaFuncSetCacheConfig(dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor, cudaFuncCachePreferL1); -// -// dvc_ScaLBL_D3Q19_AAodd_GreyscaleColor<<>>(d_neighborList, Map, dist, Aq, Bq, Den, Phi, GreySolidGrad, Poros, Perm,Vel, -// rhoA, rhoB, tauA, tauB, tauA_eff, tauB_eff,alpha, beta, Fx, Fy, Fz, strideY, strideZ, start, finish, Np); -// -// hipError_t err = hipGetLastError(); -// if (hipSuccess != err){ -// printf("hip error in ScaLBL_D3Q19_AAodd_GreyscaleColor: %s \n",hipGetErrorString(err)); -// } -// //cudaProfilerStop(); -//} diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 66497226..f473c4d4 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -145,6 +145,20 @@ void ScaLBL_ColorModel::ReadParams(string filename){ else if (domain_db->keyExists( "BC" )){ BoundaryCondition = domain_db->getScalar( "BC" ); } + if (domain_db->keyExists( "InletLayersPhase" )){ + int inlet_layers_phase = domain_db->getScalar( "InletLayersPhase" ); + if (inlet_layers_phase == 2 ) { + inletA = 0.0; + inletB = 1.0; + } + } + if (domain_db->keyExists( "OutletLayersPhase" )){ + int outlet_layers_phase = domain_db->getScalar( "OutletLayersPhase" ); + if (outlet_layers_phase == 1 ) { + inletA = 1.0; + inletB = 0.0; + } + } // Override user-specified boundary condition for specific protocols auto protocol = color_db->getWithDefault( "protocol", "none" ); @@ -169,6 +183,36 @@ void ScaLBL_ColorModel::ReadParams(string filename){ } domain_db->putScalar( "BC", BoundaryCondition ); } + else if (protocol == "fractional flow"){ + if (BoundaryCondition != 0 && BoundaryCondition != 5){ + BoundaryCondition = 0; + if (rank==0) printf("WARNING: protocol (fractional flow) supports only full periodic boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + } + else if (protocol == "centrifuge"){ + if (BoundaryCondition != 3 ){ + BoundaryCondition = 3; + if (rank==0) printf("WARNING: protocol (centrifuge) supports only constant pressure boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + } + else if (protocol == "core flooding"){ + if (BoundaryCondition != 4){ + BoundaryCondition = 4; + if (rank==0) printf("WARNING: protocol (core flooding) supports only volumetric flux boundary condition \n"); + } + domain_db->putScalar( "BC", BoundaryCondition ); + if (color_db->keyExists( "capillary_number" )){ + double capillary_number = color_db->getScalar( "capillary_number" ); + double MuB = rhoB*(tauB - 0.5)/3.0; + double IFT = 6.0*alpha; + double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2)); + flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB; + if (rank==0) printf(" protocol (core flooding): set flux=%f to achieve Ca=%f \n",flux, capillary_number); + } + color_db->putScalar( "flux", flux ); + } } void ScaLBL_ColorModel::SetDomain(){ @@ -299,11 +343,12 @@ void ScaLBL_ColorModel::AssignComponentLabels(double *phase) if (WettingConvention == "SCAL"){ for (size_t idx=0; idx analysis_db; bool Regular = false; + bool RESCALE_FORCE = false; + bool SET_CAPILLARY_NUMBER = false; + double tolerance = 0.01; auto current_db = db->cloneDatabase(); + auto flow_db = db->getDatabase( "FlowAdaptor" ); + int MIN_STEADY_TIMESTEPS = flow_db->getWithDefault( "min_steady_timesteps", 1000000 ); + int MAX_STEADY_TIMESTEPS = flow_db->getWithDefault( "max_steady_timesteps", 1000000 ); + + int RESCALE_FORCE_AFTER_TIMESTEP = MAX_STEADY_TIMESTEPS*2; + + double capillary_number = 1.0e-5; + double Ca_previous = 0.0; + if (color_db->keyExists( "capillary_number" )){ + capillary_number = color_db->getScalar( "capillary_number" ); + SET_CAPILLARY_NUMBER=true; + } + if (color_db->keyExists( "rescale_force_after_timestep" )){ + RESCALE_FORCE_AFTER_TIMESTEP = color_db->getScalar( "rescale_force_after_timestep" ); + RESCALE_FORCE = true; + } + if (analysis_db->keyExists( "tolerance" )){ + tolerance = analysis_db->getScalar( "tolerance" ); + } + + runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); auto t1 = std::chrono::system_clock::now(); + int CURRENT_TIMESTEP = 0; int START_TIMESTEP = timestep; int EXIT_TIMESTEP = min(timestepMax,returntime); while (timestep < EXIT_TIMESTEP ) { @@ -626,7 +696,183 @@ double ScaLBL_ColorModel::Run(int returntime){ alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); ScaLBL_Comm->Barrier(); //************************************************************************ + analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den ); // allow initial ramp-up to get closer to steady state + + CURRENT_TIMESTEP += 2; + if (CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS){ + analysis.finish(); + + double volB = Averages->gwb.V; + double volA = Averages->gnb.V; + volA /= Dm->Volume; + volB /= Dm->Volume;; + //initial_volume = volA*Dm->Volume; + double vA_x = Averages->gnb.Px/Averages->gnb.M; + double vA_y = Averages->gnb.Py/Averages->gnb.M; + double vA_z = Averages->gnb.Pz/Averages->gnb.M; + double vB_x = Averages->gwb.Px/Averages->gwb.M; + double vB_y = Averages->gwb.Py/Averages->gwb.M; + double vB_z = Averages->gwb.Pz/Averages->gwb.M; + double muA = rhoA*(tauA-0.5)/3.f; + double muB = rhoB*(tauB-0.5)/3.f; + 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 current_saturation = volB/(volA+volB); + double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z); + double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z); + double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha); + + bool isSteady = false; + if ( (fabs((Ca - Ca_previous)/Ca) < tolerance && CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS)) + isSteady = true; + if (CURRENT_TIMESTEP >= MAX_STEADY_TIMESTEPS) + isSteady = true; + if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP){ + RESCALE_FORCE = false; + double RESCALE_FORCE_FACTOR = capillary_number / Ca; + if (RESCALE_FORCE_FACTOR > 2.0) RESCALE_FORCE_FACTOR = 2.0; + if (RESCALE_FORCE_FACTOR < 0.5) RESCALE_FORCE_FACTOR = 0.5; + Fx *= RESCALE_FORCE_FACTOR; + Fy *= RESCALE_FORCE_FACTOR; + Fz *= RESCALE_FORCE_FACTOR; + force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz); + if (force_mag > 1e-3){ + Fx *= 1e-3/force_mag; // impose ceiling for stability + Fy *= 1e-3/force_mag; + Fz *= 1e-3/force_mag; + } + if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca); + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); + color_db->putVector("F",{Fx,Fy,Fz}); + } + if ( isSteady ){ + Averages->Full(); + Averages->Write(timestep); + analysis.WriteVisData(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den ); + analysis.finish(); + + if (rank==0){ + printf("** WRITE STEADY POINT *** "); + printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); + double h = Dm->voxel_length; + // pressures + double pA = Averages->gnb.p; + double pB = Averages->gwb.p; + double pAc = Averages->gnc.p; + double pBc = Averages->gwc.p; + double pAB = (pA-pB)/(h*6.0*alpha); + double pAB_connected = (pAc-pBc)/(h*6.0*alpha); + // connected contribution + double Vol_nc = Averages->gnc.V/Dm->Volume; + double Vol_wc = Averages->gwc.V/Dm->Volume; + double Vol_nd = Averages->gnd.V/Dm->Volume; + double Vol_wd = Averages->gwd.V/Dm->Volume; + double Mass_n = Averages->gnc.M + Averages->gnd.M; + double Mass_w = Averages->gwc.M + Averages->gwd.M; + double vAc_x = Averages->gnc.Px/Mass_n; + double vAc_y = Averages->gnc.Py/Mass_n; + double vAc_z = Averages->gnc.Pz/Mass_n; + double vBc_x = Averages->gwc.Px/Mass_w; + double vBc_y = Averages->gwc.Py/Mass_w; + double vBc_z = Averages->gwc.Pz/Mass_w; + // disconnected contribution + double vAd_x = Averages->gnd.Px/Mass_n; + double vAd_y = Averages->gnd.Py/Mass_n; + double vAd_z = Averages->gnd.Pz/Mass_n; + double vBd_x = Averages->gwd.Px/Mass_w; + double vBd_y = Averages->gwd.Py/Mass_w; + double vBd_z = Averages->gwd.Pz/Mass_w; + // film contribution + double Mfn = Averages->gifs.Mn; + double Mfw = Averages->gifs.Mw; + double vfn_x = Averages->gifs.Pnx; + double vfn_y = Averages->gifs.Pny; + double vfn_z = Averages->gifs.Pnz; + double vfw_x = Averages->gifs.Pwx; + double vfw_y = Averages->gifs.Pwy; + double vfw_z = Averages->gifs.Pwz; + + double flow_rate_A_connected = Vol_nc*(vAc_x*dir_x + vAc_y*dir_y + vAc_z*dir_z); + double flow_rate_B_connected = Vol_wc*(vBc_x*dir_x + vBc_y*dir_y + vBc_z*dir_z); + double flow_rate_A_disconnected = (Vol_nd)*(vAd_x*dir_x + vAd_y*dir_y + vAd_z*dir_z); + double flow_rate_B_disconnected = (Vol_wd)*(vBd_x*dir_x + vBd_y*dir_y + vBd_z*dir_z); + + double kAeff_connected = h*h*muA*flow_rate_A_connected/(force_mag); + double kBeff_connected = h*h*muB*flow_rate_B_connected/(force_mag); + + double kAeff_disconnected = h*h*muA*flow_rate_A_disconnected/(force_mag); + double kBeff_disconnected = h*h*muB*flow_rate_B_disconnected/(force_mag); + + double kAeff = h*h*muA*(flow_rate_A)/(force_mag); + double kBeff = h*h*muB*(flow_rate_B)/(force_mag); + + /* flow rate = [volume of fluid] x [momentum of fluid] / [mass of fluid] */ + /* fluid A eats the films */ + double flow_rate_A_filmA = (flow_rate_A*Averages->gnb.M + volA*(vfn_x*dir_x + vfn_y*dir_y + vfn_z*dir_z))/(Averages->gnb.M + Mfn); + double flow_rate_B_filmA = (flow_rate_B*Averages->gwb.M - volB*(vfn_x*dir_x + vfn_y*dir_y + vfn_z*dir_z))/(Averages->gwb.M - (rhoB/rhoA)*Mfn); + /* fluid B eats the films */ + double flow_rate_A_filmB = (flow_rate_A*Averages->gnb.M - volA*(vfw_x*dir_x + vfw_y*dir_y + vfw_z*dir_z))/(Averages->gnb.M - (rhoA/rhoB)*Mfw); + double flow_rate_B_filmB = (flow_rate_B*Averages->gwb.M + volB*(vfw_x*dir_x + vfw_y*dir_y + vfw_z*dir_z))/(Averages->gwb.M + Mfw); + /* effective permeability uncertainty limits */ + double kAeff_filmA = h*h*muA*(flow_rate_A_filmA)/(force_mag); + double kBeff_filmA = h*h*muB*(flow_rate_B_filmA)/(force_mag); + double kAeff_filmB = h*h*muA*(flow_rate_A_filmB)/(force_mag); + double kBeff_filmB = h*h*muB*(flow_rate_B_filmB)/(force_mag); + + double viscous_pressure_drop = (rhoA*volA + rhoB*volB)*force_mag; + double Mobility = muA/muB; + + bool WriteHeader=false; + FILE * kr_log_file = fopen("relperm.csv","r"); + if (kr_log_file != NULL) + fclose(kr_log_file); + else + WriteHeader=true; + kr_log_file = fopen("relperm.csv","a"); + if (WriteHeader){ + fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water eff.perm.oil.connected eff.perm.water.connected eff.perm.oil.disconnected eff.perm.water.disconnected "); + fprintf(kr_log_file,"eff.perm.oil.upper.bound eff.perm.water.lower.bound eff.perm.oil.lower.bound eff.perm.water.upper.bound "); + fprintf(kr_log_file,"cap.pressure cap.pressure.connected pressure.drop Ca M\n"); + } + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",CURRENT_TIMESTEP,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected); + fprintf(kr_log_file,"%.5g %.5g %.5g %.5g ",kAeff_filmA, kBeff_filmA, kAeff_filmB,kBeff_filmB); + fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g\n",pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility); + fclose(kr_log_file); + + printf(" Measured capillary number %f \n ",Ca); + } + if (SET_CAPILLARY_NUMBER ){ + Fx *= capillary_number / Ca; + Fy *= capillary_number / Ca; + Fz *= capillary_number / Ca; + if (force_mag > 1e-3){ + Fx *= 1e-3/force_mag; // impose ceiling for stability + Fy *= 1e-3/force_mag; + Fz *= 1e-3/force_mag; + } + if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca); + Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); + color_db->putVector("F",{Fx,Fy,Fz}); + } + else{ + if (rank==0){ + printf("** Continue to simulate steady *** \n "); + printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous); + } + } + } + } } + analysis.finish(); PROFILE_STOP("Update"); PROFILE_STOP("Loop"); @@ -695,32 +941,7 @@ void ScaLBL_ColorModel::Run(){ if (color_db->keyExists( "krA_morph_factor" )){ KRA_MORPH_FACTOR = color_db->getScalar( "krA_morph_factor" ); } - /* defaults for simulation protocols */ - auto protocol = color_db->getWithDefault( "protocol", "none" ); - if (protocol == "image sequence"){ - // Get the list of images - USE_DIRECT = true; - ImageList = color_db->getVector( "image_sequence"); - IMAGE_INDEX = color_db->getWithDefault( "image_index", 0 ); - IMAGE_COUNT = ImageList.size(); - morph_interval = 10000; - USE_MORPH = true; - } - else if (protocol == "seed water"){ - morph_delta = -0.05; - seed_water = 0.01; - USE_SEED = true; - USE_MORPH = true; - } - else if (protocol == "open connected oil"){ - morph_delta = -0.05; - USE_MORPH = true; - USE_MORPHOPEN_OIL = true; - } - else if (protocol == "shell aggregation"){ - morph_delta = -0.05; - USE_MORPH = true; - } + if (color_db->keyExists( "capillary_number" )){ capillary_number = color_db->getScalar( "capillary_number" ); SET_CAPILLARY_NUMBER=true; @@ -739,7 +960,6 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "seed_water" )){ seed_water = analysis_db->getScalar( "seed_water" ); if (rank == 0) printf("Seed water in oil %f (seed_water) \n",seed_water); - ASSERT(protocol == "seed water"); } if (analysis_db->keyExists( "morph_delta" )){ morph_delta = analysis_db->getScalar( "morph_delta" ); @@ -769,7 +989,54 @@ void ScaLBL_ColorModel::Run(){ if (analysis_db->keyExists( "max_morph_timesteps" )){ MAX_MORPH_TIMESTEPS = analysis_db->getScalar( "max_morph_timesteps" ); } - + + /* defaults for simulation protocols */ + auto protocol = color_db->getWithDefault( "protocol", "none" ); + if (protocol == "image sequence"){ + // Get the list of images + USE_DIRECT = true; + ImageList = color_db->getVector( "image_sequence"); + IMAGE_INDEX = color_db->getWithDefault( "image_index", 0 ); + IMAGE_COUNT = ImageList.size(); + morph_interval = 10000; + USE_MORPH = true; + USE_SEED = false; + } + else if (protocol == "seed water"){ + morph_delta = -0.05; + seed_water = 0.01; + USE_SEED = true; + USE_MORPH = true; + } + else if (protocol == "open connected oil"){ + morph_delta = -0.05; + USE_SEED = false; + USE_MORPH = true; + USE_MORPHOPEN_OIL = true; + } + else if (protocol == "shell aggregation"){ + morph_delta = -0.05; + USE_MORPH = true; + USE_SEED = false; + } + else if (protocol == "fractional flow"){ + USE_MORPH = false; + USE_SEED = false; + } + else if (protocol == "centrifuge"){ + USE_MORPH = false; + USE_SEED = false; + } + else if (protocol == "core flooding"){ + USE_MORPH = false; + USE_SEED = false; + if (SET_CAPILLARY_NUMBER){ + double MuB = rhoB*(tauB - 0.5)/3.0; + double IFT = 6.0*alpha; + double CrossSectionalArea = (double) (nprocx*(Nx-2)*nprocy*(Ny-2)); + flux = Dm->Porosity()*CrossSectionalArea*IFT*capillary_number/MuB; + } + } if (rank==0){ printf("********************************************************\n"); if (protocol == "image sequence"){ @@ -802,6 +1069,26 @@ void ScaLBL_ColorModel::Run(){ printf(" tolerance = %f \n",tolerance); printf(" morph_delta = %f \n",morph_delta); } + else if (protocol == "fractional flow"){ + printf(" using protocol = fractional flow \n"); + printf(" min_steady_timesteps = %i \n",MIN_STEADY_TIMESTEPS); + printf(" max_steady_timesteps = %i \n",MAX_STEADY_TIMESTEPS); + printf(" tolerance = %f \n",tolerance); + } + else if (protocol == "centrifuge"){ + printf(" using protocol = centrifuge \n"); + printf(" driving force = %f \n",Fz); + if (Fz < 0){ + printf(" Component B displacing component A \n"); + } + else if (Fz > 0){ + printf(" Component A displacing component B \n"); + } + } + else if (protocol == "core flooding"){ + printf(" using protocol = core flooding \n"); + printf(" capillary number = %f \n", capillary_number); + } printf("No. of timesteps: %i \n", timestepMax); fflush(stdout); } @@ -914,7 +1201,7 @@ void ScaLBL_ColorModel::Run(){ double volB = Averages->gwb.V; double volA = Averages->gnb.V; volA /= Dm->Volume; - volB /= Dm->Volume;; + volB /= Dm->Volume; //initial_volume = volA*Dm->Volume; double vA_x = Averages->gnb.Px/Averages->gnb.M; double vA_y = Averages->gnb.Py/Averages->gnb.M; @@ -1727,61 +2014,565 @@ FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M){ Nz = M.Dm->Nz; timestep=-1; timestep_previous=-1; - + phi.resize(Nx,Ny,Nz); phi.fill(0); // phase indicator field phi_t.resize(Nx,Ny,Nz); phi_t.fill(0); // time derivative for the phase indicator field } FlowAdaptor::~FlowAdaptor(){ + +} + +double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename){ + int rank = M.rank; + int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz; + if (rank==0) printf("Re-initializing fluids from file: %s \n", Filename.c_str()); + M.Mask->Decomp(Filename); + for (int i=0; iid[i]; // save what was read + for (int i=0; iid[i] = M.Mask->id[i]; // save what was read + + double *PhaseLabel; + PhaseLabel = new double[Nx*Ny*Nz]; + M.AssignComponentLabels(PhaseLabel); + + double Count = 0.0; + double PoreCount = 0.0; + for (int k=1; kComm.sumReduce( Count); + PoreCount=M.Dm->Comm.sumReduce( PoreCount); + if (rank==0) printf(" new saturation: %f (%f / %f) \n", Count / PoreCount, Count, PoreCount); + ScaLBL_CopyToDevice(M.Phi, PhaseLabel, Nx*Ny*Nz*sizeof(double)); + M.Dm->Comm.barrier(); + + ScaLBL_D3Q19_Init(M.fq, M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np); + M.Dm->Comm.barrier(); + + ScaLBL_CopyToHost(M.Averages->Phi.data(),M.Phi,Nx*Ny*Nz*sizeof(double)); + + double saturation = Count/PoreCount; + return saturation; +} + + +double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ + + double MASS_FRACTION_CHANGE = 0.0005; + if (M.db->keyExists( "FlowAdaptor" )){ + auto flow_db = M.db->getDatabase( "FlowAdaptor" ); + MASS_FRACTION_CHANGE = flow_db->getWithDefault( "mass_fraction_factor", 0.0005); + } + int Np = M.Np; + double dA, dB, phi; + double vx,vy,vz; + double vax,vay,vaz; + double vbx,vby,vbz; + double vax_global,vay_global,vaz_global; + double vbx_global,vby_global,vbz_global; + + double mass_a, mass_b, mass_a_global, mass_b_global; + + double *Aq_tmp, *Bq_tmp; + double *Vel_x, *Vel_y, *Vel_z, *Phase; + + Aq_tmp = new double [7*Np]; + Bq_tmp = new double [7*Np]; + Phase = new double [Np]; + Vel_x = new double [Np]; + Vel_y = new double [Np]; + Vel_z = new double [Np]; + + ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Vel_x, &M.Velocity[0], Np*sizeof(double)); + ScaLBL_CopyToHost(Vel_y, &M.Velocity[Np], Np*sizeof(double)); + ScaLBL_CopyToHost(Vel_z, &M.Velocity[2*Np], Np*sizeof(double)); + + int Nx = M.Nx; int Ny = M.Ny; int Nz = M.Nz; + + /* compute the total momentum */ + vax = vay = vaz = 0.0; + vbx = vby = vbz = 0.0; + mass_a = mass_b = 0.0; + double maxSpeed = 0.0; + double localMaxSpeed = 0.0; + for (int k=1; kSDs(i,j,k); + if (!(n<0) ){ + dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + phi = (dA - dB) / (dA + dB); + Phase[n] = phi; + mass_a += dA; + mass_b += dB; + if (phi > 0.0){ + vax += Vel_x[n]; + vay += Vel_y[n]; + vaz += Vel_z[n]; + } + else { + vbx += Vel_x[n]; + vby += Vel_y[n]; + vbz += Vel_z[n]; + } + double speed = sqrt(vax*vax + vay*vay + vaz*vaz + vbx*vbx + vby*vby + vbz*vbz); + if (distance > 3.0 && speed > localMaxSpeed){ + localMaxSpeed = speed; + } + } + } + } + } + maxSpeed = M.Dm->Comm.maxReduce(localMaxSpeed); + mass_a_global = M.Dm->Comm.sumReduce(mass_a); + mass_b_global = M.Dm->Comm.sumReduce(mass_b); + vax_global = M.Dm->Comm.sumReduce(vax); + vay_global = M.Dm->Comm.sumReduce(vay); + vaz_global = M.Dm->Comm.sumReduce(vaz); + vbx_global = M.Dm->Comm.sumReduce(vbx); + vby_global = M.Dm->Comm.sumReduce(vby); + vbz_global = M.Dm->Comm.sumReduce(vbz); + + double total_momentum_A = sqrt(vax_global*vax_global+vay_global*vay_global+vaz_global*vaz_global); + double total_momentum_B = sqrt(vbx_global*vbx_global+vby_global*vby_global+vbz_global*vbz_global); + /* compute the total mass change */ + double TOTAL_MASS_CHANGE = MASS_FRACTION_CHANGE*(mass_a_global + mass_b_global); + if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_a_global ) + TOTAL_MASS_CHANGE = 0.1*mass_a_global; + if (fabs(TOTAL_MASS_CHANGE) > 0.1*mass_b_global ) + TOTAL_MASS_CHANGE = 0.1*mass_b_global; + + double LOCAL_MASS_CHANGE = 0.0; + for (int k=1; k maxSpeed) local_momentum = maxSpeed; + if (phi > 0.0){ + LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_A; + Aq_tmp[n] -= 0.3333333333333333*LOCAL_MASS_CHANGE; + Aq_tmp[n+Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+2*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+3*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+4*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+5*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + Aq_tmp[n+6*Np] -= 0.1111111111111111*LOCAL_MASS_CHANGE; + //DebugMassA[n] = (-1.0)*LOCAL_MASS_CHANGE; + } + else{ + LOCAL_MASS_CHANGE = TOTAL_MASS_CHANGE*local_momentum/total_momentum_B; + Bq_tmp[n] += 0.3333333333333333*LOCAL_MASS_CHANGE; + Bq_tmp[n+Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+2*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+3*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+4*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+5*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + Bq_tmp[n+6*Np] += 0.1111111111111111*LOCAL_MASS_CHANGE; + //DebugMassB[n] = LOCAL_MASS_CHANGE; + } + } + } + } + } + + if (M.rank == 0) printf("Update Fractional Flow: change mass of fluid B by %f \n",TOTAL_MASS_CHANGE/mass_b_global); + + // Need to initialize Aq, Bq, Den, Phi directly + //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + + return(TOTAL_MASS_CHANGE); +} + +void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){ + + int Np = M.Np; + double dA, dB; + + double *Aq_tmp, *Bq_tmp; + + Aq_tmp = new double [7*Np]; + Bq_tmp = new double [7*Np]; + + ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + + for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ + dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + if (dA > 1.0){ + double mass_change = dA - 1.0; + Aq_tmp[n] -= 0.333333333333333*mass_change; + Aq_tmp[n+Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + if (dB > 1.0){ + double mass_change = dB - 1.0; + Bq_tmp[n] -= 0.333333333333333*mass_change; + Bq_tmp[n+Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + } + for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ + dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + if (dA > 1.0){ + double mass_change = dA - 1.0; + Aq_tmp[n] -= 0.333333333333333*mass_change; + Aq_tmp[n+Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Aq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + if (dB > 1.0){ + double mass_change = dB - 1.0; + Bq_tmp[n] -= 0.333333333333333*mass_change; + Bq_tmp[n+Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+2*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+3*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+4*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+5*Np] -= 0.111111111111111*mass_change; + Bq_tmp[n+6*Np] -= 0.111111111111111*mass_change; + } + } + + ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); } double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){ - - double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.975 ); - double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault( "move_interface_factor", 10.0 ); - ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) ); - /* compute the local derivative of phase indicator field */ - double beta = M.beta; - double factor = 0.5/beta; - double total_interface_displacement = 0.0; - double total_interface_sites = 0.0; - for (int n=0; nPhi(n); - double dist1 = factor*log((1.0+value1)/(1.0-value1)); - double value2 = phi(n); - double dist2 = factor*log((1.0+value2)/(1.0-value2)); - phi_t(n) = value2; - if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){ - /* time derivative of distance */ - double dxdt = 0.125*(dist2-dist1); - /* extrapolate to move the distance further */ - double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt; - /* compute the new phase interface */ - phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f); - total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt); - total_interface_sites += 1.0; - } - } - ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) ); - - -/* ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_PhaseField_Init(dvcMap, Phi, Den, Aq, Bq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){ - if (Dm->kproc()==0){ - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,0); - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,1); - ScaLBL_SetSlice_z(Phi,1.0,Nx,Ny,Nz,2); - } - if (Dm->kproc() == nprocz-1){ - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-1); - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-2); - ScaLBL_SetSlice_z(Phi,-1.0,Nx,Ny,Nz,Nz-3); + double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.1 ); + double MOVE_INTERFACE_FACTOR = M.color_db->getWithDefault( "move_interface_factor", 10.0 ); + + ScaLBL_CopyToHost( phi.data(), M.Phi, Nx*Ny*Nz* sizeof( double ) ); + /* compute the local derivative of phase indicator field */ + double beta = M.beta; + double factor = 0.5/beta; + double total_interface_displacement = 0.0; + double total_interface_sites = 0.0; + for (int n=0; nPhi(n); + double dist1 = factor*log((1.0+value1)/(1.0-value1)); + double value2 = phi(n); + double dist2 = factor*log((1.0+value2)/(1.0-value2)); + phi_t(n) = value2; + if (value1 < INTERFACE_CUTOFF && value1 > -1*INTERFACE_CUTOFF && value2 < INTERFACE_CUTOFF && value2 > -1*INTERFACE_CUTOFF ){ + /* time derivative of distance */ + double dxdt = 0.125*(dist2-dist1); + /* extrapolate to move the distance further */ + double dist3 = dist2 + MOVE_INTERFACE_FACTOR*dxdt; + /* compute the new phase interface */ + phi_t(n) = (2.f*(exp(-2.f*beta*(dist3)))/(1.f+exp(-2.f*beta*(dist3))) - 1.f); + total_interface_displacement += fabs(MOVE_INTERFACE_FACTOR*dxdt); + total_interface_sites += 1.0; } } - */ + ScaLBL_CopyToDevice( M.Phi, phi_t.data(), Nx*Ny*Nz* sizeof( double ) ); + return total_interface_sites; } +double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_delta_volume){ + + const RankInfoStruct rank_info(M.rank,M.nprocx,M.nprocy,M.nprocz); + auto rank = M.rank; + auto Nx = M.Nx; auto Ny = M.Ny; auto Nz = M.Nz; + auto N = Nx*Ny*Nz; + double vF = 0.f; + double vS = 0.f; + double delta_volume; + double WallFactor = 1.0; + bool USE_CONNECTED_NWP = false; + + DoubleArray phase(Nx,Ny,Nz); + IntArray phase_label(Nx,Ny,Nz);; + DoubleArray phase_distance(Nx,Ny,Nz); + Array phase_id(Nx,Ny,Nz); + fillHalo fillDouble(M.Dm->Comm,M.Dm->rank_info,{Nx-2,Ny-2,Nz-2},{1,1,1},0,1); + + // Basic algorithm to + // 1. Copy phase field to CPU + ScaLBL_CopyToHost(phase.data(), M.Phi, N*sizeof(double)); + + double count = 0.f; + for (int k=1; k 0.f && M.Averages->SDs(i,j,k) > 0.f) count+=1.f; + } + } + } + double volume_initial = M.Dm->Comm.sumReduce( count); + double PoreVolume = M.Dm->Volume*M.Dm->Porosity(); + /*ensure target isn't an absurdly small fraction of pore volume */ + if (volume_initial < target_delta_volume*PoreVolume){ + volume_initial = target_delta_volume*PoreVolume; + } + + // 2. Identify connected components of phase field -> phase_label + + double volume_connected = 0.0; + double second_biggest = 0.0; + if (USE_CONNECTED_NWP){ + ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm); + M.Dm->Comm.barrier(); + + // only operate on component "0"ScaLBL_ColorModel &M, + count = 0.0; + + for (int k=0; kComm.sumReduce( count); + second_biggest = M.Dm->Comm.sumReduce( second_biggest); + } + else { + // use the whole NWP + for (int k=0; kSDs(i,j,k) > 0.f){ + if (phase(i,j,k) > 0.f ){ + phase_id(i,j,k) = 0; + } + else { + phase_id(i,j,k) = 1; + } + } + else { + phase_id(i,j,k) = 1; + } + } + } + } + } + + // 3. Generate a distance map to the largest object -> phase_distance + CalcDist(phase_distance,phase_id,*M.Dm); + + double temp,value; + double factor=0.5/M.beta; + for (int k=0; k 1.f) value=1.f; + if (value < -1.f) value=-1.f; + // temp -- distance based on analytical form McClure, Prins et al, Comp. Phys. Comm. + temp = -factor*log((1.0+value)/(1.0-value)); + /// use this approximation close to the object + if (fabs(value) < 0.8 && M.Averages->SDs(i,j,k) > 1.f ){ + phase_distance(i,j,k) = temp; + } + // erase the original object + phase(i,j,k) = -1.0; + } + } + } + } + + + if (rank==0) printf("Pathway volume / next largest ganglion %f \n",volume_connected/second_biggest ); + + 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(M.Averages->SDs,phase_distance,phase_id,M.Averages->Dm, target_delta_volume_incremental, WallFactor); + + 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*M.beta*d))/(1.f+exp(-2.f*M.beta*d))-1.f); + } + } + } + } + } + fillDouble.fill(phase); + + count = 0.f; + for (int k=1; k 0.f && M.Averages->SDs(i,j,k) > 0.f){ + count+=1.f; + } + } + } + } + double volume_final= M.Dm->Comm.sumReduce( count); + + delta_volume = (volume_final-volume_initial); + if (rank == 0) printf("Shell Aggregation: change fluid volume fraction by %f \n", delta_volume/volume_initial); + if (rank == 0) printf(" new saturation = %f \n", volume_final/(M.Mask->Porosity()*double((Nx-2)*(Ny-2)*(Nz-2)*M.nprocs))); + + // 6. copy back to the device + //if (rank==0) printf("MorphInit: copy data back to device\n"); + ScaLBL_CopyToDevice(M.Phi,phase.data(),N*sizeof(double)); + + // 7. Re-initialize phase field and density + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np); + ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np); + auto BoundaryCondition = M.BoundaryCondition; + if (BoundaryCondition == 1 || BoundaryCondition == 2 || BoundaryCondition == 3 || BoundaryCondition == 4){ + if (M.Dm->kproc()==0){ + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,0); + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,1); + ScaLBL_SetSlice_z(M.Phi,1.0,Nx,Ny,Nz,2); + } + if (M.Dm->kproc() == M.nprocz-1){ + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-1); + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-2); + ScaLBL_SetSlice_z(M.Phi,-1.0,Nx,Ny,Nz,Nz-3); + } + } + return delta_volume; +} + + +double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil){ + srand(time(NULL)); + auto rank = M.rank; + auto Np = M.Np; + double mass_loss =0.f; + double count =0.f; + double *Aq_tmp, *Bq_tmp; + + Aq_tmp = new double [7*Np]; + Bq_tmp = new double [7*Np]; + + ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double)); + ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double)); + + + for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ + double random_value = seed_water_in_oil*double(rand())/ RAND_MAX; + double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + double phase_id = (dA - dB) / (dA + dB); + if (phase_id > 0.0){ + Aq_tmp[n] -= 0.3333333333333333*random_value; + Aq_tmp[n+Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value; + + Bq_tmp[n] += 0.3333333333333333*random_value; + Bq_tmp[n+Np] += 0.1111111111111111*random_value; + Bq_tmp[n+2*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+3*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+4*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+5*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+6*Np] += 0.1111111111111111*random_value; + } + mass_loss += random_value*seed_water_in_oil; + } + + for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ + double random_value = seed_water_in_oil*double(rand())/ RAND_MAX; + double dA = Aq_tmp[n] + Aq_tmp[n+Np] + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np]; + double dB = Bq_tmp[n] + Bq_tmp[n+Np] + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np]; + double phase_id = (dA - dB) / (dA + dB); + if (phase_id > 0.0){ + Aq_tmp[n] -= 0.3333333333333333*random_value; + Aq_tmp[n+Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value; + Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value; + + Bq_tmp[n] += 0.3333333333333333*random_value; + Bq_tmp[n+Np] += 0.1111111111111111*random_value; + Bq_tmp[n+2*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+3*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+4*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+5*Np] += 0.1111111111111111*random_value; + Bq_tmp[n+6*Np] += 0.1111111111111111*random_value; + } + mass_loss += random_value*seed_water_in_oil; + } + + count= M.Dm->Comm.sumReduce( count); + mass_loss= M.Dm->Comm.sumReduce( mass_loss); + if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count); + + // Need to initialize Aq, Bq, Den, Phi directly + //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double)); + ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double)); + + return(mass_loss); +} diff --git a/models/ColorModel.h b/models/ColorModel.h index 7d3c858a..d69e3980 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -72,6 +72,8 @@ public: double *ColorGrad; double *Velocity; double *Pressure; + + void AssignComponentLabels(double *phase); private: Utilities::MPI comm; @@ -85,7 +87,6 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); - void AssignComponentLabels(double *phase); double ImageInit(std::string filename); double MorphInit(const double beta, const double morph_delta); double SeedPhaseField(const double seed_water_in_oil); @@ -97,6 +98,11 @@ public: FlowAdaptor(ScaLBL_ColorModel &M); ~FlowAdaptor(); double MoveInterface(ScaLBL_ColorModel &M); + double ImageInit(ScaLBL_ColorModel &M, std::string Filename); + double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume); + double UpdateFractionalFlow(ScaLBL_ColorModel &M); + double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil); + void Flatten(ScaLBL_ColorModel &M); DoubleArray phi; DoubleArray phi_t; private: diff --git a/models/FreeLeeModel.cpp b/models/FreeLeeModel.cpp index b1c0e99f..19a2afa8 100644 --- a/models/FreeLeeModel.cpp +++ b/models/FreeLeeModel.cpp @@ -63,7 +63,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray // Gets data (in optimized layout) from the HOST and stores in regular layout // Primarly for debugging int i,j,k,idx; - int n; // initialize the array regdata.fill(0.f); @@ -72,7 +71,6 @@ void ScaLBL_FreeLeeModel::getData_RegularLayout(const double *data, DoubleArray for (k=0; kLastExterior(); n++){ -// va = cDen[n]; -// vb = cDen[Np + n]; -// value = (va-vb)/(va+vb); -// idx = TmpMap[n]; -// if (!(idx < 0) && idxFirstInterior(); nLastInterior(); n++){ -// va = cDen[n]; -// vb = cDen[Np + n]; -// value = (va-vb)/(va+vb); -// idx = TmpMap[n]; -// if (!(idx < 0) && idxBarrier(); -// comm.barrier(); -// -// if (rank==0) printf ("Initializing phase and density fields on device from Restart\n"); -// //TODO the following function is to be updated. -// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, 0, ScaLBL_Comm->LastExterior(), Np); -// //ScaLBL_FreeLeeModel_PhaseField_InitFromRestart(Den, hq, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); } } double ScaLBL_FreeLeeModel::Run_TwoFluid(int returntime){ - int nprocs=nprocx*nprocy*nprocz; int START_TIME = timestep; int EXIT_TIME = min(returntime, timestepMax); diff --git a/models/GreyscaleColorModel.cpp b/models/GreyscaleColorModel.cpp index 73c1f878..9e3a9bd5 100644 --- a/models/GreyscaleColorModel.cpp +++ b/models/GreyscaleColorModel.cpp @@ -309,18 +309,26 @@ void ScaLBL_GreyscaleColorModel::AssignGreySolidLabels()//apply capillary penalt // oil-wet < 0 // neutral = 0 (i.e. no penalty) double *GreySolidW_host = new double [Np]; + double *GreySn_host = new double [Np]; + double *GreySw_host = new double [Np]; size_t NLABELS=0; signed char VALUE=0; double AFFINITY=0.f; + double Sn,Sw;//end-point saturation of greynodes set by users auto LabelList = greyscaleColor_db->getVector( "GreySolidLabels" ); auto AffinityList = greyscaleColor_db->getVector( "GreySolidAffinity" ); + auto SnList = greyscaleColor_db->getVector( "GreySnList" ); + auto SwList = greyscaleColor_db->getVector( "GreySwList" ); NLABELS=LabelList.size(); if (NLABELS != AffinityList.size()){ ERROR("Error: GreySolidLabels and GreySolidAffinity must be the same length! \n"); } + if (NLABELS != SnList.size() || NLABELS != SwList.size()){ + ERROR("Error: GreySolidLabels, GreySnList, and GreySwList must be the same length! \n"); + } for (int k=0;k0: water-wet || grey-solid affinity<0: oil-wet \n"); } ScaLBL_CopyToDevice(GreySolidW, GreySolidW_host, Np*sizeof(double)); + ScaLBL_CopyToDevice(GreySn, GreySn_host, Np*sizeof(double)); + ScaLBL_CopyToDevice(GreySw, GreySw_host, Np*sizeof(double)); ScaLBL_Comm->Barrier(); delete [] GreySolidW_host; + delete [] GreySn_host; + delete [] GreySw_host; } ////----------------------------------------------------------------------------------------------------------// @@ -598,6 +619,8 @@ void ScaLBL_GreyscaleColorModel::Create(){ //ScaLBL_AllocateDeviceMemory((void **) &GreySolidPhi, sizeof(double)*Nx*Ny*Nz); //ScaLBL_AllocateDeviceMemory((void **) &GreySolidGrad, 3*sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &GreySolidW, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &GreySn, sizeof(double)*Np); + ScaLBL_AllocateDeviceMemory((void **) &GreySw, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Porosity_dvc, sizeof(double)*Np); ScaLBL_AllocateDeviceMemory((void **) &Permeability_dvc, sizeof(double)*Np); //........................................................................... @@ -921,7 +944,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ // Halo exchange for phase field ScaLBL_Comm_Regular->SendHalo(Phi); //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //Model-1&4 @@ -950,7 +973,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAodd_GreyscaleColor_CP(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); //Model-1&4 @@ -983,7 +1006,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } ScaLBL_Comm_Regular->SendHalo(Phi); //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); //Model-1&4 @@ -1012,7 +1035,7 @@ void ScaLBL_GreyscaleColorModel::Run(){ } //Model-1&4 with capillary pressure penalty for grey nodes - ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,Porosity_dvc,Permeability_dvc,Velocity,Pressure, + ScaLBL_D3Q19_AAeven_GreyscaleColor_CP(dvcMap, fq, Aq, Bq, Den, Phi, GreySolidW,GreySn,GreySw,Porosity_dvc,Permeability_dvc,Velocity,Pressure, rhoA, rhoB, tauA, tauB,tauA_eff, tauB_eff, alpha, beta, Fx, Fy, Fz, RecoloringOff, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); //Model-1&4 diff --git a/models/GreyscaleColorModel.h b/models/GreyscaleColorModel.h index a9a0bfc9..bff186e1 100644 --- a/models/GreyscaleColorModel.h +++ b/models/GreyscaleColorModel.h @@ -68,6 +68,8 @@ public: //double *GreySolidPhi; //Model 2 & 3 //double *GreySolidGrad;//Model 1 & 4 double *GreySolidW; + double *GreySn; + double *GreySw; //double *ColorGrad; double *Velocity; double *Pressure; diff --git a/models/GreyscaleModel.cpp b/models/GreyscaleModel.cpp index 308cc1e6..11f8ed63 100644 --- a/models/GreyscaleModel.cpp +++ b/models/GreyscaleModel.cpp @@ -150,7 +150,6 @@ void ScaLBL_GreyscaleModel::ReadInput(){ // Generate the signed distance map // Initialize the domain and communication Array id_solid(Nx,Ny,Nz); - int count = 0; // Solve for the position of the solid phase for (int k=0;kid[i] = Mask->id[i]; - for (int idx=0; idxComm.sumReduce( label_count[idx]); + for (size_t idx=0; idxComm.sumReduce( label_count[idx]); //Initialize a weighted porosity after considering grey voxels GreyPorosity=0.0; for (unsigned int idx=0; idx0 @@ -691,23 +691,22 @@ void ScaLBL_GreyscaleModel::Run(){ //double absperm = h*h*mu*Mask->Porosity()*flow_rate / force_mag; double absperm = h*h*mu*GreyPorosity*flow_rate / force_mag; - if (rank==0){ + if (rank==0){ printf(" AbsPerm = %.5g [micron^2]\n",absperm); - bool WriteHeader=false; - FILE * log_file = fopen("Permeability.csv","r"); - if (log_file != NULL) - fclose(log_file); - else - WriteHeader=true; - log_file = fopen("Permeability.csv","a"); - if (WriteHeader) - fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n", - timestep,Fx,Fy,Fz,mu,h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz,absperm); + bool WriteHeader=false; + FILE * log_file = fopen("Permeability.csv","r"); + if (log_file != NULL) + fclose(log_file); + else + WriteHeader=true; + log_file = fopen("Permeability.csv","a"); + if (WriteHeader) + fprintf(log_file,"timestep Fx Fy Fz mu Vs As Hs Xs vax vay vaz AbsPerm \n"); - fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, - h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm); - fclose(log_file); - } + fprintf(log_file,"%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",timestep, Fx, Fy, Fz, mu, + h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm); + fclose(log_file); + } } if (timestep%visualization_interval==0){ diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 8b6a9ff6..45c6a1ea 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -97,7 +97,7 @@ void ScaLBL_IonModel::ReadParams(string filename,vector &num_iter){ } else{ time_conv.clear(); - for (int i=0; i &num_iter){ ERROR("Error: number_ion_species and IonDiffusivityList must be the same length! \n"); } else{ - for (int i=0; i &num_iter){ ERROR("Error: number_ion_species and IonConcentrationList must be the same length! \n"); } else{ - for (int i=0; i &num_iter){ else { ERROR("Error: Non-periodic BCs are specified but InletValueList cannot be found! \n"); } - for (unsigned int i=0;i &num_iter){ else { ERROR("Error: Non-periodic BCs are specified but OutletValueList cannot be found! \n"); } - for (unsigned int i=0;iid[n]; AFFINITY=0.f; // Assign the affinity from the paired list - for (unsigned int idx=0; idx < NLABELS; idx++){ + for (size_t idx=0; idx < NLABELS; idx++){ //printf("idx=%i, value=%i, %i, \n",idx, VALUE,LabelList[idx]); if (VALUE == LabelList[idx]){ AFFINITY=AffinityList[idx]; @@ -725,23 +727,23 @@ void ScaLBL_IonModel::Initialize(){ auto File_ion = ion_db->getVector( "IonConcentrationFile" ); double *Ci_host; Ci_host = new double[number_ion_species*Np]; - for (int ic=0; icFirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np); } @@ -758,13 +760,13 @@ void ScaLBL_IonModel::Initialize(){ break; } - for (int i=0; i rlx; - for (unsigned int ic=0;icBarrier(); comm.barrier(); //auto t1 = std::chrono::system_clock::now(); - for (int ic=0; icFirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_Ion_ChargeDensity(Ci, ChargeDensity, IonValence[ic], ic, 0, ScaLBL_Comm->LastExterior(), Np); } @@ -961,7 +963,7 @@ void ScaLBL_IonModel::Run(double *Velocity, double *ElectricField){ //if (rank==0) printf("********************************************************\n"); } -void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const int ic){ +void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const size_t ic){ //This function wirte out the data in a normal layout (by aggregating all decomposed domains) ScaLBL_Comm->RegularLayout(Map,&Ci[ic*Np],IonConcentration); @@ -973,13 +975,13 @@ void ScaLBL_IonModel::getIonConcentration(DoubleArray &IonConcentration, const i void ScaLBL_IonModel::getIonConcentration_debug(int timestep){ //This function write out decomposed data DoubleArray PhaseField(Nx,Ny,Nz); - for (int ic=0; icRegularLayout(Map,&Ci[ic*Np],PhaseField); ScaLBL_Comm->Barrier(); comm.barrier(); IonConcentration_LB_to_Phys(PhaseField); FILE *OUTFILE; - sprintf(LocalRankFilename,"Ion%02i_Time_%i.%05i.raw",ic+1,timestep,rank); + sprintf(LocalRankFilename,"Ion%02zu_Time_%i.%05i.raw",ic+1,timestep,rank); OUTFILE = fopen(LocalRankFilename,"wb"); fwrite(PhaseField.data(),8,N,OUTFILE); fclose(OUTFILE); @@ -1046,7 +1048,7 @@ double ScaLBL_IonModel::CalIonDenConvergence(vector &ci_avg_previous){ Ci_host = new double[Np]; vector error(number_ion_species,0.0); - for (int ic=0; ic BoundaryConditionInlet; vector BoundaryConditionOutlet; vector IonDiffusivity;//User input unit [m^2/sec] diff --git a/models/StokesModel.cpp b/models/StokesModel.cpp index bbd02e85..8385cbff 100644 --- a/models/StokesModel.cpp +++ b/models/StokesModel.cpp @@ -59,6 +59,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){ BoundaryCondition = 0; if (stokes_db->keyExists( "BC" )){ BoundaryCondition = stokes_db->getScalar( "BC" ); + } if (stokes_db->keyExists( "tolerance" )){ tolerance = stokes_db->getScalar( "tolerance" ); @@ -110,6 +111,7 @@ void ScaLBL_StokesModel::ReadParams(string filename,int num_iter){ void ScaLBL_StokesModel::ReadParams(string filename){ //NOTE the max time step is left unspecified + // read the input database db = std::make_shared( filename ); domain_db = db->getDatabase( "Domain" ); @@ -298,8 +300,10 @@ void ScaLBL_StokesModel::AssignZetaPotentialSolid(double *zeta_potential_solid) ERROR("Error: LB Single-Fluid Solver: SolidLabels and ZetaPotentialSolidList must be the same length! \n"); } - double label_count[NLABELS]; - double label_count_global[NLABELS]; + double *label_count; + double *label_count_global; + label_count = new double [NLABELS]; + label_count_global = new double [NLABELS]; for (size_t idx=0; idx9){ - printf("Input arguments accepted \n"); - Nx = atoi(argv[1]); - Ny = atoi(argv[2]); - Nz = atoi(argv[3]); - x0 = atoi(argv[4]); - y0 = atoi(argv[5]); - z0 = atoi(argv[6]); - NX = atol(argv[7]); - NY = atol(argv[8]); - NZ = atol(argv[9]); - printf("Size %i X %i X %i \n",NX,NY,NZ); - fflush(stdout); + printf("Input arguments accepted \n"); + Nx = atoi(argv[1]); + Ny = atoi(argv[2]); + Nz = atoi(argv[3]); + x0 = atoi(argv[4]); + y0 = atoi(argv[5]); + z0 = atoi(argv[6]); + NX = atol(argv[7]); + NY = atol(argv[8]); + NZ = atol(argv[9]); + printf("Size %llu X %llu X %llu \n",(unsigned long long) NX, (unsigned long long) NY, (unsigned long long) NZ); + fflush(stdout); } else{ - printf("setting defaults \n"); - Nx=Ny=Nz=1024; - x0=y0=z0=0; - NX=NY=8640; - NZ=6480; + printf("setting defaults \n"); + Nx=Ny=Nz=1024; + x0=y0=z0=0; + NX=NY=8640; + NZ=6480; } //Bx = By = Bz = 9; //Nx = Ny = Nz = 1024; @@ -53,29 +53,28 @@ int main(int argc, char **argv){ if (By>8) By=8; if (Bz>8) Bz=8; - printf("System size (output) is: %i x %i x %i \n",NX,NY,NZ); - printf("Block size (read) is: %i x %i x %i \n",Nx,Ny,Nz); - printf("Starting location (read) is: %i, %i, %i \n", x0,y0,z0); - printf("Block number (read): %i x %i x %i \n",Bx,By,Bz); + printf("System size (output) is: %llu x %llu x %llu \n",(unsigned long long) NX,(unsigned long long) NY, (unsigned long long) NZ); + printf("Block size (read) is: %llu x %llu x %llu \n",(unsigned long long) Nx,(unsigned long long) Ny,(unsigned long long) Nz); + printf("Starting location (read) is: %llu, %llu, %llu \n", (unsigned long long) x0,(unsigned long long) y0,(unsigned long long) z0); + printf("Block number (read): %llu x %llu x %llu \n",(unsigned long long) Bx,(unsigned long long) By,(unsigned long long) Bz); fflush(stdout); - + // Filenames used //char LocalRankString[8]; char LocalRankFilename[40]; char sx[2]; char sy[2]; char sz[2]; - char tmpstr[10]; //sprintf(LocalRankString,"%05d",rank); N = Nx*Ny*Nz; N_full=NX*NY*NZ; - + char *id; id = new char [N]; char *ID; ID = new char [N_full]; - + for (unsigned int bz=0; bz SW) { sw_diff_old = sw_diff_new; diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index ac76c3d6..3dade4d5 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -9,8 +9,6 @@ #include "common/Utilities.h" #include "models/ColorModel.h" -//#define WRE_SURFACES - /* * Simulator for two-phase flow in porous media * James E. McClure 2013-2014 @@ -24,81 +22,170 @@ int main( int argc, char **argv ) { - // Initialize - Utilities::startup( argc, argv ); + // Initialize + Utilities::startup( argc, argv ); - { // Limit scope so variables that contain communicators will free before MPI_Finialize + { // Limit scope so variables that contain communicators will free before MPI_Finialize - Utilities::MPI comm( MPI_COMM_WORLD ); - int rank = comm.getRank(); - int nprocs = comm.getSize(); - std::string SimulationMode = "production"; - // Load the input database - auto db = std::make_shared( argv[1] ); - if (argc > 2) { - SimulationMode = "development"; - } + Utilities::MPI comm( MPI_COMM_WORLD ); + int rank = comm.getRank(); + int nprocs = comm.getSize(); + std::string SimulationMode = "production"; + // Load the input database + auto db = std::make_shared( argv[1] ); + if (argc > 2) { + SimulationMode = "development"; + } - if ( rank == 0 ) { - printf( "********************************************************\n" ); - printf( "Running Color LBM \n" ); - printf( "********************************************************\n" ); - if (SimulationMode == "development") - printf("**** DEVELOPMENT MODE ENABLED *************\n"); - } - // Initialize compute device - int device = ScaLBL_SetDevice( rank ); - NULL_USE( device ); - ScaLBL_DeviceBarrier(); - comm.barrier(); + if ( rank == 0 ) { + printf( "********************************************************\n" ); + printf( "Running Color LBM \n" ); + printf( "********************************************************\n" ); + if (SimulationMode == "development") + printf("**** DEVELOPMENT MODE ENABLED *************\n"); + } + // Initialize compute device + int device = ScaLBL_SetDevice( rank ); + NULL_USE( device ); + ScaLBL_DeviceBarrier(); + comm.barrier(); - PROFILE_ENABLE( 1 ); - // PROFILE_ENABLE_TRACE(); - // PROFILE_ENABLE_MEMORY(); - PROFILE_SYNCHRONIZE(); - PROFILE_START( "Main" ); - Utilities::setErrorHandlers(); + PROFILE_ENABLE( 1 ); + // PROFILE_ENABLE_TRACE(); + // PROFILE_ENABLE_MEMORY(); + PROFILE_SYNCHRONIZE(); + PROFILE_START( "Main" ); + Utilities::setErrorHandlers(); - auto filename = argv[1]; - ScaLBL_ColorModel ColorModel( rank, nprocs, comm ); - ColorModel.ReadParams( filename ); - ColorModel.SetDomain(); - ColorModel.ReadInput(); - ColorModel.Create(); // creating the model will create data structure to match the pore - // structure and allocate variables - ColorModel.Initialize(); // initializing the model will set initial conditions for variables - - if (SimulationMode == "development"){ - double MLUPS=0.0; - int timestep = 0; - int analysis_interval = ColorModel.timestepMax; - if (ColorModel.analysis_db->keyExists( "" )){ - analysis_interval = ColorModel.analysis_db->getScalar( "analysis_interval" ); - } - FlowAdaptor Adapt(ColorModel); - runAnalysis analysis(ColorModel); - while (ColorModel.timestep < ColorModel.timestepMax){ - timestep += analysis_interval; - MLUPS = ColorModel.Run(timestep); - if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS); - - Adapt.MoveInterface(ColorModel); - } - ColorModel.WriteDebug(); - } //Analysis.WriteVis(LeeModel,LeeModel.db, timestep); + auto filename = argv[1]; + ScaLBL_ColorModel ColorModel( rank, nprocs, comm ); + ColorModel.ReadParams( filename ); + ColorModel.SetDomain(); + ColorModel.ReadInput(); + ColorModel.Create(); // creating the model will create data structure to match the pore + // structure and allocate variables + ColorModel.Initialize(); // initializing the model will set initial conditions for variables - else - ColorModel.Run(); + if (SimulationMode == "development"){ + double MLUPS=0.0; + int timestep = 0; + bool ContinueSimulation = true; + + /* Variables for simulation protocols */ + auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "none" ); + /* image sequence protocol */ + int IMAGE_INDEX = 0; + int IMAGE_COUNT = 0; + std::vector ImageList; + /* flow adaptor keys to control behavior */ + int SKIP_TIMESTEPS = 0; + int MAX_STEADY_TIME = 1000000; + double ENDPOINT_THRESHOLD = 0.1; + double FRACTIONAL_FLOW_INCREMENT = 0.0; // this will skip the flow adaptor if not enabled + double SEED_WATER = 0.0; + if (ColorModel.db->keyExists( "FlowAdaptor" )){ + auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" ); + MAX_STEADY_TIME = flow_db->getWithDefault( "max_steady_timesteps", 1000000 ); + SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 50000 ); + ENDPOINT_THRESHOLD = flow_db->getWithDefault( "endpoint_threshold", 0.1); + /* protocol specific key values */ + if (PROTOCOL == "fractional flow") + FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault( "fractional_flow_increment", 0.05); + if (PROTOCOL == "seed water") + SEED_WATER = flow_db->getWithDefault( "seed_water", 0.01); + } + /* analysis keys*/ + int ANALYSIS_INTERVAL = ColorModel.timestepMax; + if (ColorModel.analysis_db->keyExists( "analysis_interval" )){ + ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar( "analysis_interval" ); + } + /* Launch the simulation */ + FlowAdaptor Adapt(ColorModel); + runAnalysis analysis(ColorModel); + while (ContinueSimulation){ + /* this will run steady points */ + timestep += MAX_STEADY_TIME; + MLUPS = ColorModel.Run(timestep); + if (rank==0) printf("Lattice update rate (per MPI process)= %f MLUPS \n", MLUPS); + if (ColorModel.timestep > ColorModel.timestepMax){ + ContinueSimulation = false; + } + + /* Load a new image if image sequence is specified */ + if (PROTOCOL == "image sequence"){ + IMAGE_INDEX++; + if (IMAGE_INDEX < IMAGE_COUNT){ + std::string next_image = ImageList[IMAGE_INDEX]; + if (rank==0) printf("***Loading next image in sequence (%i) ***\n",IMAGE_INDEX); + ColorModel.color_db->putScalar("image_index",IMAGE_INDEX); + Adapt.ImageInit(ColorModel, next_image); + } + else{ + if (rank==0) printf("Finished simulating image sequence \n"); + ColorModel.timestep = ColorModel.timestepMax; + } + } + /*********************************************************/ + /* update the fluid configuration with the flow adapter */ + int skip_time = 0; + timestep = ColorModel.timestep; + /* get the averaged flow measures computed internally for the last simulation point*/ + double SaturationChange = 0.0; + double volB = ColorModel.Averages->gwb.V; + double volA = ColorModel.Averages->gnb.V; + double initialSaturation = volB/(volA + volB); + double vA_x = ColorModel.Averages->gnb.Px/ColorModel.Averages->gnb.M; + double vA_y = ColorModel.Averages->gnb.Py/ColorModel.Averages->gnb.M; + double vA_z = ColorModel.Averages->gnb.Pz/ColorModel.Averages->gnb.M; + double vB_x = ColorModel.Averages->gwb.Px/ColorModel.Averages->gwb.M; + double vB_y = ColorModel.Averages->gwb.Py/ColorModel.Averages->gwb.M; + double vB_z = ColorModel.Averages->gwb.Pz/ColorModel.Averages->gwb.M; + double speedA = sqrt(vA_x*vA_x + vA_y*vA_y + vA_z*vA_z); + double speedB = sqrt(vB_x*vB_x + vB_y*vB_y + vB_z*vB_z); + /* stop simulation if previous point was sufficiently close to the endpoint*/ + if (volA*speedA < ENDPOINT_THRESHOLD*volB*speedB) ContinueSimulation = false; + if (ContinueSimulation){ + while (skip_time < SKIP_TIMESTEPS && fabs(SaturationChange) < fabs(FRACTIONAL_FLOW_INCREMENT) ){ + timestep += ANALYSIS_INTERVAL; + if (PROTOCOL == "fractional flow") { + Adapt.UpdateFractionalFlow(ColorModel); + } + else if (PROTOCOL == "shell aggregation"){ + double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange; + Adapt.ShellAggregation(ColorModel,target_volume_change); + } + else if (PROTOCOL == "seed water"){ + Adapt.SeedPhaseField(ColorModel,SEED_WATER); + } + /* Run some LBM timesteps to let the system relax a bit */ + MLUPS = ColorModel.Run(timestep); + /* Recompute the volume fraction now that the system has adjusted */ + double volB = ColorModel.Averages->gwb.V; + double volA = ColorModel.Averages->gnb.V; + SaturationChange = volB/(volA + volB) - initialSaturation; + skip_time += ANALYSIS_INTERVAL; + } + if (rank==0) printf(" ********************************************************************* \n"); + if (rank==0) printf(" Updated fractional flow with saturation change = %f \n", SaturationChange); + if (rank==0) printf(" Used protocol = %s \n", PROTOCOL.c_str()); + if (rank==0) printf(" ********************************************************************* \n"); + } + /*********************************************************/ + } + } + else + ColorModel.Run(); - PROFILE_STOP( "Main" ); - auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); - auto level = db->getWithDefault( "TimerLevel", 1 ); - PROFILE_SAVE( file, level ); - // **************************************************** + PROFILE_STOP( "Main" ); + auto file = db->getWithDefault( "TimerFile", "lbpm_color_simulator" ); + auto level = db->getWithDefault( "TimerLevel", 1 ); + NULL_USE(level); + PROFILE_SAVE( file, level ); + // **************************************************** - } // Limit scope so variables that contain communicators will free before MPI_Finialize + } // Limit scope so variables that contain communicators will free before MPI_Finialize - Utilities::shutdown(); - return 0; + Utilities::shutdown(); + return 0; } diff --git a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp index 19d99b9c..30ec4fec 100644 --- a/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp +++ b/tests/lbpm_freelee_SingleFluidBGK_simulator.cpp @@ -59,6 +59,7 @@ int main( int argc, char **argv ) PROFILE_STOP("Main"); auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_SingleFluidBGK_simulator" ); auto level = db->getWithDefault( "TimerLevel", 1 ); + NULL_USE(level); PROFILE_SAVE( file,level ); // **************************************************** diff --git a/tests/lbpm_freelee_simulator.cpp b/tests/lbpm_freelee_simulator.cpp index cdba5082..e2ca2d64 100644 --- a/tests/lbpm_freelee_simulator.cpp +++ b/tests/lbpm_freelee_simulator.cpp @@ -83,6 +83,7 @@ int main( int argc, char **argv ) PROFILE_STOP("Main"); auto file = db->getWithDefault( "TimerFile", "lbpm_freelee_simulator" ); auto level = db->getWithDefault( "TimerLevel", 1 ); + NULL_USE(level); PROFILE_SAVE( file,level ); // ****************************************************