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

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

Browse Source
This commit is contained in:
JamesEMcclure 2021-05-19 14:24:54 -04:00
commit 3742568ac9
3 changed files with 375 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

352
models/ColorModel.cpp
View File

@ -185,6 +185,13 @@ void ScaLBL_ColorModel::ReadParams(string filename){
}
domain_db->putScalar<int>( "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<int>( "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<Database> 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<int>( "min_steady_timesteps", 1000000 );
int MAX_STEADY_TIMESTEPS = flow_db->getWithDefault<int>( "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<double>( "capillary_number" );
SET_CAPILLARY_NUMBER=true;
}
if (color_db->keyExists( "rescale_force_after_timestep" )){
RESCALE_FORCE_AFTER_TIMESTEP = color_db->getScalar<int>( "rescale_force_after_timestep" );
RESCALE_FORCE = true;
}
if (analysis_db->keyExists( "tolerance" )){
tolerance = analysis_db->getScalar<double>( "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<double>("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<double>("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<double>( "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<double>( "move_interface_cutoff", 0.975 );

1
models/ColorModel.h
View File

@ -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:

29
tests/lbpm_color_simulator.cpp
View File

@ -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<int>( "analysis_interval" );
/* flow adaptor keys to control */
auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
int MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
int SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 100000 );
int ANALYSIS_INTERVAL = ColorModel.timestepMax;
if (ColorModel.analysis_db->keyExists( "analysis_interval" )){
ANALYSIS_INTERVAL = ColorModel.analysis_db->getScalar<int>( "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();