diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 484a5656..8e32d5ac 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -185,6 +185,13 @@ 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 ); + } } void ScaLBL_ColorModel::SetDomain(){ @@ -549,14 +556,39 @@ void ScaLBL_ColorModel::Initialize(){ double ScaLBL_ColorModel::Run(int returntime){ int nprocs=nprocx*nprocy*nprocz; - + const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); //************ MAIN ITERATION LOOP ***************************************/ comm.barrier(); PROFILE_START("Loop"); //std::shared_ptr 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 ) { @@ -642,7 +674,158 @@ 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; + + 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); + + 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 cap.pressure cap.pressure.connected pressure.drop Ca M\n"); + + fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_TIMESTEP,current_saturation,kAeff,kBeff,kAeff_connected,kBeff_connected,kAeff_disconnected,kBeff_disconnected,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"); @@ -818,6 +1001,13 @@ 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); + printf(" morph_delta = %f \n",morph_delta); + } printf("No. of timesteps: %i \n", timestepMax); fflush(stdout); } @@ -1752,6 +1942,166 @@ FlowAdaptor::~FlowAdaptor(){ } +double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){ + + double MASS_FRACTION_CHANGE = 0.05; + if (M.db->keyExists( "FlowAdaptor" )){ + auto flow_db = M.db->getDatabase( "FlowAdaptor" ); + MASS_FRACTION_CHANGE = flow_db->getWithDefault( "fractional_flow_increment", 0.05); + } + + 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)); + + /* compute the total momentum */ + vax = vay = vaz = 0.0; + vbx = vby = vbz = 0.0; + mass_a = mass_b = 0.0; + 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]; + 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]; + } + } + 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]; + 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]; + } + } + 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 n=0; n < M.ScaLBL_Comm->LastExterior(); n++){ + phi = Phase[n]; + vx = Vel_x[n]; + vy = Vel_y[n]; + vz = Vel_z[n]; + double local_momentum = sqrt(vx*vx+vy*vy+vz*vz); + 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; + } + 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; + } + } + + for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){ + phi = Phase[n]; + vx = Vel_x[n]; + vy = Vel_y[n]; + vz = Vel_z[n]; + double local_momentum = sqrt(vx*vx+vy*vy+vz*vz); + 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; + } + 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; + } + } + + 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); +} + double FlowAdaptor::MoveInterface(ScaLBL_ColorModel &M){ double INTERFACE_CUTOFF = M.color_db->getWithDefault( "move_interface_cutoff", 0.975 ); diff --git a/models/ColorModel.h b/models/ColorModel.h index c76a4568..1ec72b73 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -113,6 +113,7 @@ public: FlowAdaptor(ScaLBL_ColorModel &M); ~FlowAdaptor(); double MoveInterface(ScaLBL_ColorModel &M); + double UpdateFractionalFlow(ScaLBL_ColorModel &M); DoubleArray phi; DoubleArray phi_t; private: diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index ac76c3d6..507156a1 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -71,21 +71,38 @@ int main( int argc, char **argv ) 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" ); + /* flow adaptor keys to control */ + auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" ); + int MAX_STEADY_TIME = flow_db->getWithDefault( "max_steady_timesteps", 1000000 ); + int SKIP_TIMESTEPS = flow_db->getWithDefault( "skip_timesteps", 100000 ); + + 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 (ColorModel.timestep < ColorModel.timestepMax){ - timestep += analysis_interval; + /* 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); + Adapt.UpdateFractionalFlow(ColorModel); - Adapt.MoveInterface(ColorModel); + /* apply timestep skipping algorithm to accelerate steady-state */ + int skip_time = 0; + timestep = ColorModel.timestep; + while (skip_time < SKIP_TIMESTEPS){ + timestep += ANALYSIS_INTERVAL; + MLUPS = ColorModel.Run(timestep); + Adapt.MoveInterface(ColorModel); + skip_time += ANALYSIS_INTERVAL; + } } ColorModel.WriteDebug(); - } //Analysis.WriteVis(LeeModel,LeeModel.db, timestep); + } else ColorModel.Run();