From e6fa7d40656ceb490140496f39955a2ee628e4ca Mon Sep 17 00:00:00 2001 From: James McClure Date: Wed, 1 Dec 2021 08:08:16 -0500 Subject: [PATCH 1/8] add SCAL file --- docs/source/examples/color/steadyState.rst | 123 ++++++++++++++++++++- docs/source/userGuide/models/mrt/mrt.rst | 3 +- models/ColorModel.cpp | 28 ++++- 3 files changed, 148 insertions(+), 6 deletions(-) diff --git a/docs/source/examples/color/steadyState.rst b/docs/source/examples/color/steadyState.rst index 997f882e..35eaf46b 100644 --- a/docs/source/examples/color/steadyState.rst +++ b/docs/source/examples/color/steadyState.rst @@ -1,5 +1,120 @@ -******************* -Steady-state flow -******************* +******************************** +Steady-state flow (color model) +******************************** + +In this example we simulate a steady-state flow with a constant driving force. This will enforce a periodic boundary condition +in all directions. While the driving force may be set in any direction, we will set it in the z-direction to be consistent +with the convention for pressure and velocity boundary conditions. + + +For the case considered in ``example/DiscPack`` we specify the following information in the input file + +.. code:: c + + Domain { + Filename = "discs_3x128x128.raw.morphdrain.raw" + ReadType = "8bit" // data type + N = 3, 128, 128 // size of original image + nproc = 1, 2, 2 // process grid + n = 3, 64, 64 // sub-domain size + voxel_length = 1.0 // voxel length (in microns) + ReadValues = 0, 1, 2 // labels within the original image + WriteValues = 0, 1, 2 // associated labels to be used by LBPM + BC = 0 // fully periodic BC + Sw = 0.35 // target saturation for morphological tools + } + + Color { + capillary_number = -1e-5 // capillary number for the displacement, positive="oil injection" + timestepMax = 20000 // maximum timtestep + alpha = 0.01 // controls interfacial tension + rhoA = 1.0 // controls the density of fluid A + rhoB = 1.0 // controls the density of fluid B + tauA = 0.7 // controls the viscosity of fluid A + tauB = 0.7 // controls the viscosity of fluid B + F = 0, 0, 1e-5 // body force + WettingConvention = "SCAL" + ComponentLabels = 0 // image labels for solid voxels + ComponentAffinity = 0.9 // controls the wetting affinity for each label + Restart = false + } + Analysis { + analysis_interval = 1000 // logging interval for timelog.csv + subphase_analysis_interval = 500000 // loggging interval for subphase.csv + N_threads = 4 // number of analysis threads (GPU version only) + visualization_interval = 10000 // interval to write visualization files + restart_interval = 10000000 // interval to write restart file + restart_file = "Restart" // base name of restart file + } + Visualization { + format = "hdf5" + write_silo = true // write SILO databases with assigned variables + save_8bit_raw = true // write labeled 8-bit binary files with phase assignments + save_phase_field = true // save phase field within SILO database + save_pressure = true // save pressure field within SILO database + save_velocity = false // save velocity field within SILO database + } + FlowAdaptor { + } + + +Once this has been set, we launch ``lbpm_color_simulator`` in the same way as other parallel tools + +.. code:: bash + + mpirun -np 4 $LBPM_BIN/lbpm_color_simulator input.db + +Successful output looks like the following + + +.. code:: bash + + ******************************************************** + Running Color LBM + ******************************************************** + voxel length = 1.000000 micron + voxel length = 1.000000 micron + Input media: discs_3x128x128.raw.morphdrain.raw + Relabeling 3 values + oldvalue=0, newvalue =0 + oldvalue=1, newvalue =1 + oldvalue=2, newvalue =2 + Dimensions of segmented image: 3 x 128 x 128 + Reading 8-bit input data + Read segmented data from discs_3x128x128.raw.morphdrain.raw + Label=0, Count=11862 + Label=1, Count=26430 + Label=2, Count=10860 + Distributing subdomains across 4 processors + Process grid: 1 x 2 x 2 + Subdomain size: 3 x 64 x 64 + Size of transition region: 0 + Media porosity = 0.758667 + Initialized solid phase -- Converting to Signed Distance function + Domain set. + Create ScaLBL_Communicator + Set up memory efficient layout, 9090 | 9120 | 21780 + Allocating distributions + Setting up device map and neighbor list + Component labels: 1 + label=0, affinity=-0.900000, volume fraction==0.417582 + Initializing distributions + Initializing phase field + Affinities - rank 0: + Main: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 2: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 3: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 4: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Affinities - rank 0: + Main: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 1: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 2: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 3: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + Thread 4: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 + ******************************************************** + CPU time = 0.001501 + Lattice update rate (per core)= 6.074861 MLUPS + Lattice update rate (per MPI process)= 6.074861 MLUPS + (flatten density field) -In this example we simulate a steady-state flow with a constant driving force. diff --git a/docs/source/userGuide/models/mrt/mrt.rst b/docs/source/userGuide/models/mrt/mrt.rst index eddf10d4..55b4f158 100644 --- a/docs/source/userGuide/models/mrt/mrt.rst +++ b/docs/source/userGuide/models/mrt/mrt.rst @@ -225,4 +225,5 @@ Example Input File InletLayers = 0, 0, 10 // specify 10 layers along the z-inlet BC = 0 // boundary condition type (0 for periodic) } - + Visualization { + } diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c7caa292..c5eff406 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -184,6 +184,7 @@ void ScaLBL_ColorModel::ReadParams(string filename){ domain_db->putScalar( "BC", BoundaryCondition ); } else if (protocol == "fractional flow"){ + if (rank == 0) printf("Using fractional flow protocol \n"); if (BoundaryCondition != 0 && BoundaryCondition != 5){ BoundaryCondition = 0; if (rank==0) printf("WARNING: protocol (fractional flow) supports only full periodic boundary condition \n"); @@ -191,6 +192,7 @@ void ScaLBL_ColorModel::ReadParams(string filename){ domain_db->putScalar( "BC", BoundaryCondition ); } else if (protocol == "centrifuge"){ + if (rank == 0) printf("Using centrifuge protocol \n"); if (BoundaryCondition != 3 ){ BoundaryCondition = 3; if (rank==0) printf("WARNING: protocol (centrifuge) supports only constant pressure boundary condition \n"); @@ -617,6 +619,9 @@ double ScaLBL_ColorModel::Run(int returntime){ tolerance = analysis_db->getScalar( "tolerance" ); } + auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); + + runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); auto t1 = std::chrono::system_clock::now(); int CURRENT_TIMESTEP = 0; @@ -872,7 +877,28 @@ double ScaLBL_ColorModel::Run(int returntime){ 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); - + + /* SCAL file */ + if (WettingConvention == "SCAL"){ + WriteHeader=false; + FILE * scal_log_file = fopen("SCAL.csv","r"); + if (scal_log_file != NULL) + fclose(scal_log_file); + else + WriteHeader=true; + scal_log_file = fopen("SCAL.csv","a"); + if (WriteHeader){ + fprintf(scal_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(scal_log_file,"eff.perm.oil.upper.bound eff.perm.water.lower.bound eff.perm.oil.lower.bound eff.perm.water.upper.bound "); + fprintf(scal_log_file,"cap.pressure cap.pressure.connected pressure.drop Ca M\n"); + } + fprintf(scal_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(scal_log_file,"%.5g %.5g %.5g %.5g ",kAeff_filmA, kBeff_filmA, kAeff_filmB,kBeff_filmB); + fprintf(scal_log_file,"%.5g %.5g %.5g %.5g %.5g\n",pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility); + fclose(scal_log_file); + /****************/ + } + printf(" Measured capillary number %f \n ",Ca); } if (SET_CAPILLARY_NUMBER ){ From 7d525e999b2de7d779f83421bc117ad676d4708f Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Thu, 2 Dec 2021 23:32:57 +1100 Subject: [PATCH 2/8] enable different types of solid BC for Poisson solver; to be built and tested --- common/ScaLBL.cpp | 8 ++++ common/ScaLBL.h | 3 ++ cpu/D3Q7BC.cpp | 20 +++++++++ models/PoissonSolver.cpp | 89 +++++++++++++++++++++------------------- models/PoissonSolver.h | 3 +- 5 files changed, 80 insertions(+), 43 deletions(-) diff --git a/common/ScaLBL.cpp b/common/ScaLBL.cpp index 6a405c0a..420e9879 100644 --- a/common/ScaLBL.cpp +++ b/common/ScaLBL.cpp @@ -1356,6 +1356,14 @@ void ScaLBL_Communicator::SolidNeumannD3Q7(double *fq, double *BoundaryValue){ ScaLBL_Solid_Neumann_D3Q7(fq,BoundaryValue,bb_dist,bb_interactions,n_bb_d3q7); } +void ScaLBL_Communicator::SolidDirichletAndNeumannD3Q7(double *fq, double *BoundaryValue, int *BoundaryLabel){ + // fq is a D3Q7 distribution + // BoundaryValues is a list of values to assign at bounce-back solid sites + // BoundaryLabel: is a list of integer labels indicating the type of BCs + // 1-> Dirichlet BC; 2-> Neumann BC. + ScaLBL_Solid_DirichletAndNeumann_D3Q7(fq,BoundaryValue,BoundaryLabel,bb_dist,bb_interactions,n_bb_d3q7); +} + void ScaLBL_Communicator::SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad, double epsilon_LB, double tau, double rho0, double den_scale,double h, double time_conv){ // fq is a D3Q19 distribution diff --git a/common/ScaLBL.h b/common/ScaLBL.h index 816c1551..30a74113 100644 --- a/common/ScaLBL.h +++ b/common/ScaLBL.h @@ -593,6 +593,8 @@ extern "C" void ScaLBL_Solid_Dirichlet_D3Q7(double *dist,double *BoundaryValue,i extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int *BounceBackDist_list,int *BounceBackSolid_list,int N); +extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist,double *BoundaryValue,int *BoundaryLabel,int *BounceBackDist_list,int *BounceBackSolid_list,int N); + extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad, double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv, int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list, @@ -700,6 +702,7 @@ public: void SetupBounceBackList(IntArray &Map, signed char *id, int Np, bool SlippingVelBC=false); void SolidDirichletD3Q7(double *fq, double *BoundaryValue); void SolidNeumannD3Q7(double *fq, double *BoundaryValue); + void SolidDirichletAndNeumannD3Q7(double *fq, double *BoundaryValue, int *BoundaryLabel); void SolidSlippingVelocityBCD3Q19(double *fq, double *zeta_potential, double *ElectricField, double *SolidGrad, double epslion_LB, double tau, double rho0, double den_scale,double h, double time_conv); diff --git a/cpu/D3Q7BC.cpp b/cpu/D3Q7BC.cpp index 1c255495..b2fa75f7 100644 --- a/cpu/D3Q7BC.cpp +++ b/cpu/D3Q7BC.cpp @@ -30,6 +30,26 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist,double *BoundaryValue,int } } +extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist,double *BoundaryValue,int* BoundaryLabel,int *BounceBackDist_list,int *BounceBackSolid_list,int N){ + + int idx; + int iq,ib; + double value_b,value_b_label,value_q; + for (idx=0; idxkeyExists( "BC_Solid" )){ - BoundaryConditionSolid = electric_db->getScalar( "BC_Solid" ); + // BC_solid=1: Dirichlet-type surfacen potential + // BC_solid=2: Neumann-type surfacen charge density + BoundaryConditionSolidList.push_back(1); + if (electric_db->keyExists( "BC_SolidList" )){ + BoundaryConditionSolidList.clear(); + BoundaryConditionSolidList = electric_db->getVector( "BC_SolidList" ); } // Read boundary condition for electric potential // BC = 0: normal periodic BC @@ -133,19 +136,8 @@ void ScaLBL_Poisson::ReadParams(string filename){ else{ if (rank==0) printf("LB-Poisson Solver: tolerance_method=%s cannot be identified!\n",tolerance_method.c_str()); } - - switch (BoundaryConditionSolid){ - case 1: - if (rank==0) printf("LB-Poisson Solver: solid boundary: Dirichlet-type surfacen potential is assigned\n"); - break; - case 2: - if (rank==0) printf("LB-Poisson Solver: solid boundary: Neumann-type surfacen charge density is assigned\n"); - break; - default: - if (rank==0) printf("LB-Poisson Solver: solid boundary: Dirichlet-type surfacen potential is assigned\n"); - break; - } } + void ScaLBL_Poisson::SetDomain(){ Dm = std::shared_ptr(new Domain(domain_db,comm)); // full domain for analysis Mask = std::shared_ptr(new Domain(domain_db,comm)); // mask domain removes immobile phases @@ -243,17 +235,18 @@ void ScaLBL_Poisson::ReadInput(){ if (rank == 0) cout << " Domain set." << endl; } -void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid) +void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid, int *poisson_solid_BClabel) { signed char VALUE=0; double AFFINITY=0.f; + int BoundaryConditionSolid=0; auto LabelList = electric_db->getVector( "SolidLabels" ); auto AffinityList = electric_db->getVector( "SolidValues" ); size_t NLABELS = LabelList.size(); - if (NLABELS != AffinityList.size()){ - ERROR("Error: LB-Poisson Solver: SolidLabels and SolidValues must be the same length! \n"); + if (NLABELS != AffinityList.size() || NLABELS != BoundaryConditionSolidList.size()){ + ERROR("Error: LB-Poisson Solver: BC_SolidList, SolidLabels and SolidValues all must be of the same length! \n"); } std::vector label_count( NLABELS, 0.0 ); @@ -268,10 +261,15 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid) int n = k*Nx*Ny+j*Nx+i; VALUE=Mask->id[n]; AFFINITY=0.f; + BoundaryConditionSolid=0; // Assign the affinity from the paired list for (unsigned int idx=0; idx < NLABELS; idx++){ if (VALUE == LabelList[idx]){ AFFINITY=AffinityList[idx]; + BoundaryConditionSolid=BoundaryConditionSolidList[idx]; + if (BoundaryConditionSolid!=1 && BoundaryConditionSolid!=2){ + ERROR("Error: LB-Poisson Solver: Note only BC_SolidList of 1 or 2 is supported!\n"); + } //NOTE need to convert the user input phys unit to LB unit if (BoundaryConditionSolid==2){ //for BCS=1, i.e. Dirichlet-type, no need for unit conversion @@ -283,6 +281,7 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid) } } poisson_solid[n] = AFFINITY; + poisson_solid_BClabel[n] = BoundaryConditionSolid; } } } @@ -295,17 +294,16 @@ void ScaLBL_Poisson::AssignSolidBoundary(double *poisson_solid) for (unsigned int idx=0; idxBarrier(); ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_Poisson_Init(dvcMap, fq, Psi, 0, ScaLBL_Comm->LastExterior(), Np); delete [] psi_host; + delete [] psi_BCLabel_host; //extra treatment for halo layer //if (BoundaryCondition==1){ @@ -749,23 +752,25 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){ ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); - if (BoundaryConditionSolid==1){ - ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); - } - else if (BoundaryConditionSolid==2){ - ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi); - } + ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel); + //if (BoundaryConditionSolid==1){ + // ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); + //} + //else if (BoundaryConditionSolid==2){ + // ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi); + //} } void ScaLBL_Poisson::SolvePoissonAAeven(double *ChargeDensity){ ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_AAeven_Poisson(dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); - if (BoundaryConditionSolid==1){ - ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); - } - else if (BoundaryConditionSolid==2){ - ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi); - } + ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel); + //if (BoundaryConditionSolid==1){ + // ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); + //} + //else if (BoundaryConditionSolid==2){ + // ScaLBL_Comm->SolidNeumannD3Q7(fq, Psi); + //} } void ScaLBL_Poisson::DummyChargeDensity(){ diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 81312cab..07e03b05 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -46,7 +46,7 @@ public: int analysis_interval; int BoundaryConditionInlet; int BoundaryConditionOutlet; - int BoundaryConditionSolid; + vector BoundaryConditionSolidList; double tau; double tolerance; std::string tolerance_method; @@ -86,6 +86,7 @@ public: //signed char *dvcID; double *fq; double *Psi; + int *Psi_BCLabel; double *ElectricField; double *ChargeDensityDummy;// for debugging double *ResidualError; From 95d01c4e286d8ae364aa932f4e45c465354248e4 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Thu, 2 Dec 2021 07:47:49 -0500 Subject: [PATCH 3/8] fix bug --- models/PoissonSolver.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index 07e03b05..e36a1ea4 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -104,7 +104,7 @@ private: //int rank,nprocs; void LoadParams(std::shared_ptr db0); - void AssignSolidBoundary(double *poisson_solid); + void AssignSolidBoundary(double *poisson_solid, int *poisson_solid_label); void Potential_Init(double *psi_init); void ElectricField_LB_to_Phys(DoubleArray &Efield_reg); void SolveElectricPotentialAAodd(int timestep_from_Study); From 8fce93fc47600ad98ffb8b1f308f32b348133cc2 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Mon, 6 Dec 2021 13:50:17 +1100 Subject: [PATCH 4/8] update solid BCs of Poisson solver in GPU;to be built and tested --- cuda/D3Q7BC.cu | 31 +++++++++++++++++++++++++++++++ hip/D3Q7BC.cu | 31 +++++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+) diff --git a/cuda/D3Q7BC.cu b/cuda/D3Q7BC.cu index b2af849d..0651694d 100644 --- a/cuda/D3Q7BC.cu +++ b/cuda/D3Q7BC.cu @@ -38,6 +38,28 @@ __global__ void dvc_ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValu } } +__global__ void dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count) +{ + + int idx; + int iq,ib; + double value_b,value_b_label,value_q; + idx = blockIdx.x*blockDim.x + threadIdx.x; + if (idx < count){ + iq = BounceBackDist_list[idx]; + ib = BounceBackSolid_list[idx]; + value_b = BoundaryValue[ib];//get boundary value from a solid site + value_b_label = BoundaryLabel[ib];//get boundary label (i.e. type of BC) from a solid site + value_q = dist[iq]; + if (value_b_label==1){//Dirichlet BC + dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice + } + if (value_b_label==2){//Neumann BC + dist[iq] = value_q + value_b; + } + } +} + __global__ void dvc_ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad, double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv, int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list, @@ -733,6 +755,15 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue, i } } +extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count){ + int GRID = count / 512 + 1; + dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7<<>>(dist, BoundaryValue, BoundaryLabel, BounceBackDist_list, BounceBackSolid_list, count); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("CUDA error in ScaLBL_Solid_DirichletAndNeumann_D3Q7 (kernel): %s \n",cudaGetErrorString(err)); + } +} + extern "C" void ScaLBL_Solid_SlippingVelocityBC_D3Q19(double *dist, double *zeta_potential, double *ElectricField, double *SolidGrad, double epsilon_LB, double tau, double rho0,double den_scale, double h, double time_conv, int *BounceBackDist_list, int *BounceBackSolid_list, int *FluidBoundary_list, diff --git a/hip/D3Q7BC.cu b/hip/D3Q7BC.cu index 9413a68a..c594791a 100644 --- a/hip/D3Q7BC.cu +++ b/hip/D3Q7BC.cu @@ -37,6 +37,28 @@ __global__ void dvc_ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValu } } +__global__ void dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count) +{ + + int idx; + int iq,ib; + double value_b,value_b_label,value_q; + idx = blockIdx.x*blockDim.x + threadIdx.x; + if (idx < count){ + iq = BounceBackDist_list[idx]; + ib = BounceBackSolid_list[idx]; + value_b = BoundaryValue[ib];//get boundary value from a solid site + value_b_label = BoundaryLabel[ib];//get boundary label (i.e. type of BC) from a solid site + value_q = dist[iq]; + if (value_b_label==1){//Dirichlet BC + dist[iq] = -1.0*value_q + value_b*0.25;//NOTE 0.25 is the speed of sound for D3Q7 lattice + } + if (value_b_label==2){//Neumann BC + dist[iq] = value_q + value_b; + } + } +} + __global__ void dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np) { int idx,n; @@ -409,6 +431,15 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue, i } } +extern "C" void ScaLBL_Solid_DirichletAndNeumann_D3Q7(double *dist, double *BoundaryValue,int *BoundaryLabel, int *BounceBackDist_list, int *BounceBackSolid_list, int count){ + int GRID = count / 512 + 1; + dvc_ScaLBL_Solid_DirichletAndNeumann_D3Q7<<>>(dist, BoundaryValue, BoundaryLabel, BounceBackDist_list, BounceBackSolid_list, count); + cudaError_t err = cudaGetLastError(); + if (cudaSuccess != err){ + printf("hip error in ScaLBL_Solid_DirichletAndNeumann_D3Q7 (kernel): %s \n",cudaGetErrorString(err)); + } +} + extern "C" void ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z(int *list, double *dist, double Vin, int count, int Np){ int GRID = count / 512 + 1; dvc_ScaLBL_D3Q7_AAeven_Poisson_Potential_BC_z<<>>(list, dist, Vin, count, Np); From ba88a78afb0590009f73bc94c001f44fe15312b4 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Mon, 6 Dec 2021 16:28:08 +1100 Subject: [PATCH 5/8] update periodic potential BC for inlet and outlet;to be built and tested --- models/PoissonSolver.cpp | 65 ++++++++++++++-------------------------- models/PoissonSolver.h | 6 ++-- 2 files changed, 25 insertions(+), 46 deletions(-) diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index a60433a5..b055b753 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -1,3 +1,4 @@ + sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); /* * Multi-relaxation time LBM Model */ @@ -18,7 +19,7 @@ ScaLBL_Poisson::ScaLBL_Poisson(int RANK, int NP, const Utilities::MPI& COMM): epsilon0(0),epsilon0_LB(0),epsilonR(0),epsilon_LB(0),Vin(0),Vout(0),Nx(0),Ny(0),Nz(0),N(0),Np(0),analysis_interval(0), chargeDen_dummy(0),WriteLog(0),nprocx(0),nprocy(0),nprocz(0), BoundaryConditionInlet(0),BoundaryConditionOutlet(0),BoundaryConditionSolidList(0),Lx(0),Ly(0),Lz(0), - Vin0(0),freqIn(0),t0_In(0),Vin_Type(0),Vout0(0),freqOut(0),t0_Out(0),Vout_Type(0), + Vin0(0),freqIn(0),PhaseShift_In(0),Vout0(0),freqOut(0),PhaseShift_Out(0), TestPeriodic(0),TestPeriodicTime(0),TestPeriodicTimeConv(0),TestPeriodicSaveInterval(0), comm(COMM) { @@ -403,8 +404,7 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ //set up default boundary input parameters Vin0 = Vout0 = 1.0; //unit: [V] freqIn = freqOut = 50.0; //unit: [Hz] - t0_In = t0_Out = 0.0; //unit: [sec] - Vin_Type = Vout_Type = 1; //1->sin; 2->cos + PhaseShift_In = PhaseShift_Out = 0.0; //unit: [radian] Vin = 1.0; //Boundary-z (inlet) electric potential Vout = 1.0; //Boundary-Z (outlet) electric potential @@ -423,24 +423,12 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ if (electric_db->keyExists( "freqIn" )){//unit: Hz freqIn = electric_db->getScalar( "freqIn" ); } - if (electric_db->keyExists( "t0_In" )){//timestep shift, unit: lt - t0_In = electric_db->getScalar( "t0_In" ); - } - if (electric_db->keyExists( "Vin_Type" )){ - //type=1 -> sine - //tyep=2 -> cosine - Vin_Type = electric_db->getScalar( "Vin_Type" ); - if (Vin_Type>2 || Vin_Type<=0) ERROR("Error: user-input Vin_Type is currently not supported! \n"); + if (electric_db->keyExists( "PhaseShift_In" )){//phase shift, unit: radian + PhaseShift_In = electric_db->getScalar( "PhaseShift_In" ); } if (rank==0){ - if (Vin_Type==1){ - printf("LB-Poisson Solver: inlet boundary; periodic electric potential Vin = %.3g*Sin[2*pi*%.3g*(t+%.3g)] [V]\n",Vin0,freqIn,t0_In); - printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vin0,freqIn,t0_In); - } - else if (Vin_Type==2){ - printf("LB-Poisson Solver: inlet boundary; periodic electric potential Vin = %.3g*Cos[2*pi*%.3g*(t+%.3g)] [V] \n",Vin0,freqIn,t0_In); - printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vin0,freqIn,t0_In); - } + printf("LB-Poisson Solver: inlet boundary; periodic electric potential Vin = %.3g*Cos[2*pi*%.3g*t+%.3g] [V] \n",Vin0,freqIn,PhaseShift_In); + printf(" V0 = %.3g [V], frequency = %.3g [Hz], phase shift = %.3g [radian] \n",Vin0,freqIn,PhaseShift_In); } break; } @@ -460,31 +448,19 @@ void ScaLBL_Poisson::Potential_Init(double *psi_init){ if (electric_db->keyExists( "freqOut" )){//unit: Hz freqOut = electric_db->getScalar( "freqOut" ); } - if (electric_db->keyExists( "t0_Out" )){//timestep shift, unit: lt - t0_Out = electric_db->getScalar( "t0_Out" ); - } - if (electric_db->keyExists( "Vout_Type" )){ - //type=1 -> sine - //tyep=2 -> cosine - Vout_Type = electric_db->getScalar( "Vout_Type" ); - if (Vout_Type>2 || Vin_Type<=0) ERROR("Error: user-input Vout_Type is currently not supported! \n"); + if (electric_db->keyExists( "PhaseShift_Out" )){//timestep shift, unit: lt + PhaseShift_Out = electric_db->getScalar( "PhaseShift_Out" ); } if (rank==0){ - if (Vout_Type==1){ - printf("LB-Poisson Solver: outlet boundary; periodic electric potential Vout = %.3g*Sin[2*pi*%.3g*(t+%.3g)] [V]\n",Vout0,freqOut,t0_Out); - printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vout0,freqOut,t0_Out); - } - else if (Vout_Type==2){ - printf("LB-Poisson Solver: outlet boundary; periodic electric potential Vout = %.3g*Cos[2*pi*%.3g*(t+%.3g)] [V]\n",Vout0,freqOut,t0_Out); - printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [sec] \n",Vout0,freqOut,t0_Out); - } + printf("LB-Poisson Solver: outlet boundary; periodic electric potential Vout = %.3g*Cos[2*pi*%.3g*t+%.3g] [V]\n",Vout0,freqOut,PhaseShift_Out); + printf(" V0 = %.3g [V], frequency = %.3g [Hz], timestep shift = %.3g [radian] \n",Vout0,freqOut,PhaseShift_Out); } break; } } //By default only periodic BC is applied and Vin=Vout=1.0, i.e. there is no potential gradient along Z-axis - if (BoundaryConditionInlet==2) Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,t0_In,Vin_Type,0); - if (BoundaryConditionOutlet==2) Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,t0_Out,Vout_Type,0); + if (BoundaryConditionInlet==2) Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,0); + if (BoundaryConditionOutlet==2) Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,0); double slope = (Vout-Vin)/(Nz-2); double psi_linearized; for (int k=0;kD3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); break; case 2: - Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,t0_In,Vin_Type,timestep_from_Study); + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); break; } @@ -708,7 +684,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAodd(int timestep_from_Study){ ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); break; case 2: - Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,t0_Out,Vout_Type,timestep_from_Study); + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); break; } @@ -729,7 +705,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); break; case 2: - Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,t0_In,Vin_Type,timestep_from_Study); + Vin = getBoundaryVoltagefromPeriodicBC(Vin0,freqIn,PhaseShift_In,timestep_from_Study); ScaLBL_Comm->D3Q7_Poisson_Potential_BC_z(NeighborList, fq, Vin, timestep); break; } @@ -740,7 +716,7 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); break; case 2: - Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,t0_Out,Vout_Type,timestep_from_Study); + Vout = getBoundaryVoltagefromPeriodicBC(Vout0,freqOut,PhaseShift_Out,timestep_from_Study); ScaLBL_Comm->D3Q7_Poisson_Potential_BC_Z(NeighborList, fq, Vout, timestep); break; } @@ -752,6 +728,9 @@ void ScaLBL_Poisson::SolveElectricPotentialAAeven(int timestep_from_Study){ void ScaLBL_Poisson::SolvePoissonAAodd(double *ChargeDensity){ ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); ScaLBL_D3Q7_AAodd_Poisson(NeighborList, dvcMap, fq, ChargeDensity, Psi, ElectricField, tau, epsilon_LB, 0, ScaLBL_Comm->LastExterior(), Np); + //TODO: perhaps add another ScaLBL_Comm routine to update Psi values on solid boundary nodes. + //something like: + //ScaLBL_Comm->SolidDirichletBoundaryUpdates(Psi, Psi_BCLabel, timestep); ScaLBL_Comm->SolidDirichletAndNeumannD3Q7(fq, Psi, Psi_BCLabel); //if (BoundaryConditionSolid==1){ // ScaLBL_Comm->SolidDirichletD3Q7(fq, Psi); diff --git a/models/PoissonSolver.h b/models/PoissonSolver.h index e36a1ea4..94a6f36f 100644 --- a/models/PoissonSolver.h +++ b/models/PoissonSolver.h @@ -55,8 +55,8 @@ public: double Vin, Vout; double chargeDen_dummy;//for debugging bool WriteLog; - double Vin0,freqIn,t0_In,Vin_Type; - double Vout0,freqOut,t0_Out,Vout_Type; + double Vin0,freqIn,t0_In,PhaseShift_In; + double Vout0,freqOut,t0_Out,PhaseShift_Out; bool TestPeriodic; double TestPeriodicTime;//unit: [sec] double TestPeriodicTimeConv; //unit [sec/lt] @@ -113,7 +113,7 @@ private: void SolvePoissonAAodd(double *ChargeDensity); void SolvePoissonAAeven(double *ChargeDensity); void getConvergenceLog(int timestep,double error); - double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int V_type,int time_step); + double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step); }; #endif From d61cb8571fd46359124e2d6c693cc0294e21da68 Mon Sep 17 00:00:00 2001 From: Rex Zhe Li Date: Mon, 6 Dec 2021 00:46:44 -0500 Subject: [PATCH 6/8] fix dumb typo --- models/PoissonSolver.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/models/PoissonSolver.cpp b/models/PoissonSolver.cpp index b055b753..8658ac07 100644 --- a/models/PoissonSolver.cpp +++ b/models/PoissonSolver.cpp @@ -1,4 +1,3 @@ - sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); /* * Multi-relaxation time LBM Model */ From 44965c17f38cd8eac413d0d9f2554bb0af9edcf4 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 9 Dec 2021 10:11:05 -0500 Subject: [PATCH 7/8] fix merge --- models/ColorModel.cpp | 2235 ++++++++++++++++++++++------------------- 1 file changed, 1177 insertions(+), 1058 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index c5eff406..f63d7ee5 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -9,1123 +9,1242 @@ color lattice boltzmann model #include #include - -ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, const Utilities::MPI& COMM): - rank(RANK), nprocs(NP), Restart(0), timestep(0), timestepMax(0), - tauA(0), tauB(0), rhoA(0), rhoB(0), alpha(0), beta(0), - Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), - inletA(0), inletB(0), outletA(0), outletB(0), - Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), - BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), - NeighborList(nullptr), dvcMap(nullptr), fq(nullptr), Aq(nullptr), Bq(nullptr), - Den(nullptr), Phi(nullptr), ColorGrad(nullptr), Velocity(nullptr), Pressure(nullptr), - comm(COMM) -{ - REVERSE_FLOW_DIRECTION = false; +ScaLBL_ColorModel::ScaLBL_ColorModel(int RANK, int NP, + const Utilities::MPI &COMM) + : rank(RANK), nprocs(NP), Restart(0), timestep(0), timestepMax(0), tauA(0), + tauB(0), rhoA(0), rhoB(0), alpha(0), beta(0), Fx(0), Fy(0), Fz(0), + flux(0), din(0), dout(0), inletA(0), inletB(0), outletA(0), outletB(0), + Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), + BoundaryCondition(0), Lx(0), Ly(0), Lz(0), id(nullptr), + NeighborList(nullptr), dvcMap(nullptr), fq(nullptr), Aq(nullptr), + Bq(nullptr), Den(nullptr), Phi(nullptr), ColorGrad(nullptr), + Velocity(nullptr), Pressure(nullptr), comm(COMM) { + REVERSE_FLOW_DIRECTION = false; } -ScaLBL_ColorModel::~ScaLBL_ColorModel() -{ - delete [] id; - ScaLBL_FreeDeviceMemory( NeighborList ); - ScaLBL_FreeDeviceMemory( dvcMap ); - ScaLBL_FreeDeviceMemory( fq ); - ScaLBL_FreeDeviceMemory( Aq ); - ScaLBL_FreeDeviceMemory( Bq ); - ScaLBL_FreeDeviceMemory( Den ); - ScaLBL_FreeDeviceMemory( Phi ); - ScaLBL_FreeDeviceMemory( Pressure ); - ScaLBL_FreeDeviceMemory( Velocity ); - ScaLBL_FreeDeviceMemory( ColorGrad ); +ScaLBL_ColorModel::~ScaLBL_ColorModel() { + delete[] id; + ScaLBL_FreeDeviceMemory(NeighborList); + ScaLBL_FreeDeviceMemory(dvcMap); + ScaLBL_FreeDeviceMemory(fq); + ScaLBL_FreeDeviceMemory(Aq); + ScaLBL_FreeDeviceMemory(Bq); + ScaLBL_FreeDeviceMemory(Den); + ScaLBL_FreeDeviceMemory(Phi); + ScaLBL_FreeDeviceMemory(Pressure); + ScaLBL_FreeDeviceMemory(Velocity); + ScaLBL_FreeDeviceMemory(ColorGrad); } -/*void ScaLBL_ColorModel::WriteCheckpoint(const char *FILENAME, const double *cPhi, const double *cfq, int Np) -{ - int q,n; - double value; - ofstream File(FILENAME,ios::binary); - for (n=0; n(filename); + domain_db = db->getDatabase("Domain"); + color_db = db->getDatabase("Color"); + analysis_db = db->getDatabase("Analysis"); + vis_db = db->getDatabase("Visualization"); + + // set defaults + timestepMax = 100000; + tauA = tauB = 1.0; + rhoA = rhoB = 1.0; + Fx = Fy = Fz = 0.0; + alpha = 1e-3; + beta = 0.95; + Restart = false; + din = dout = 1.0; + flux = 0.0; + + // Color Model parameters + if (color_db->keyExists("timestepMax")) { + timestepMax = color_db->getScalar("timestepMax"); + } + if (color_db->keyExists("tauA")) { + tauA = color_db->getScalar("tauA"); + } + if (color_db->keyExists("tauB")) { + tauB = color_db->getScalar("tauB"); + } + if (color_db->keyExists("rhoA")) { + rhoA = color_db->getScalar("rhoA"); + } + if (color_db->keyExists("rhoB")) { + rhoB = color_db->getScalar("rhoB"); + } + if (color_db->keyExists("F")) { + Fx = color_db->getVector("F")[0]; + Fy = color_db->getVector("F")[1]; + Fz = color_db->getVector("F")[2]; + } + if (color_db->keyExists("alpha")) { + alpha = color_db->getScalar("alpha"); + } + if (color_db->keyExists("beta")) { + beta = color_db->getScalar("beta"); + } + if (color_db->keyExists("Restart")) { + Restart = color_db->getScalar("Restart"); + } + if (color_db->keyExists("din")) { + din = color_db->getScalar("din"); + } + if (color_db->keyExists("dout")) { + dout = color_db->getScalar("dout"); + } + if (color_db->keyExists("flux")) { + flux = color_db->getScalar("flux"); + } + inletA = 1.f; + inletB = 0.f; + outletA = 0.f; + outletB = 1.f; + + + BoundaryCondition = 0; + if (color_db->keyExists("BC")) { + BoundaryCondition = color_db->getScalar("BC"); + } 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; } } - File.close(); - -} - -void ScaLBL_ColorModel::ReadCheckpoint(char *FILENAME, double *cPhi, double *cfq, int Np) -{ - int q=0, n=0; - double value=0; - ifstream File(FILENAME,ios::binary); - for (n=0; nkeyExists("OutletLayersPhase")) { + int outlet_layers_phase = + domain_db->getScalar("OutletLayersPhase"); + if (outlet_layers_phase == 1) { + inletA = 1.0; + inletB = 0.0; } } - File.close(); -} - */ -void ScaLBL_ColorModel::ReadParams(string filename){ - // read the input database - db = std::make_shared( filename ); - domain_db = db->getDatabase( "Domain" ); - color_db = db->getDatabase( "Color" ); - analysis_db = db->getDatabase( "Analysis" ); - vis_db = db->getDatabase( "Visualization" ); - // set defaults - timestepMax = 100000; - tauA = tauB = 1.0; - rhoA = rhoB = 1.0; - Fx = Fy = Fz = 0.0; - alpha=1e-3; - beta=0.95; - Restart=false; - din=dout=1.0; - flux=0.0; - - // Color Model parameters - if (color_db->keyExists( "timestepMax" )){ - timestepMax = color_db->getScalar( "timestepMax" ); - } - if (color_db->keyExists( "tauA" )){ - tauA = color_db->getScalar( "tauA" ); - } - if (color_db->keyExists( "tauB" )){ - tauB = color_db->getScalar( "tauB" ); - } - if (color_db->keyExists( "rhoA" )){ - rhoA = color_db->getScalar( "rhoA" ); - } - if (color_db->keyExists( "rhoB" )){ - rhoB = color_db->getScalar( "rhoB" ); - } - if (color_db->keyExists( "F" )){ - Fx = color_db->getVector( "F" )[0]; - Fy = color_db->getVector( "F" )[1]; - Fz = color_db->getVector( "F" )[2]; - } - if (color_db->keyExists( "alpha" )){ - alpha = color_db->getScalar( "alpha" ); - } - if (color_db->keyExists( "beta" )){ - beta = color_db->getScalar( "beta" ); - } - if (color_db->keyExists( "Restart" )){ - Restart = color_db->getScalar( "Restart" ); - } - if (color_db->keyExists( "din" )){ - din = color_db->getScalar( "din" ); - } - if (color_db->keyExists( "dout" )){ - dout = color_db->getScalar( "dout" ); - } - if (color_db->keyExists( "flux" )){ - flux = color_db->getScalar( "flux" ); - } - inletA=1.f; - inletB=0.f; - outletA=0.f; - outletB=1.f; - //if (BoundaryCondition==4) flux *= rhoA; // mass flux must adjust for density (see formulation for details) - - BoundaryCondition = 0; - if (color_db->keyExists( "BC" )){ - BoundaryCondition = color_db->getScalar( "BC" ); - } - 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" ); - if (protocol == "seed water"){ - if (BoundaryCondition != 0 && BoundaryCondition != 5){ - BoundaryCondition = 0; - if (rank==0) printf("WARNING: protocol (seed water) supports only full periodic boundary condition \n"); - } - domain_db->putScalar( "BC", BoundaryCondition ); - } - else if (protocol == "open connected oil"){ - if (BoundaryCondition != 0 && BoundaryCondition != 5){ - BoundaryCondition = 0; - if (rank==0) printf("WARNING: protocol (open connected oil) supports only full periodic boundary condition \n"); - } - domain_db->putScalar( "BC", BoundaryCondition ); - } - else if (protocol == "shell aggregation"){ - if (BoundaryCondition != 0 && BoundaryCondition != 5){ - BoundaryCondition = 0; - if (rank==0) printf("WARNING: protocol (shell aggregation) supports only full periodic boundary condition \n"); - } - domain_db->putScalar( "BC", BoundaryCondition ); - } - else if (protocol == "fractional flow"){ - if (rank == 0) printf("Using fractional flow protocol \n"); - 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 (rank == 0) printf("Using centrifuge protocol \n"); - 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 (rank == 0) printf("Using core flooding protocol \n"); - 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 ); - } + // Override user-specified boundary condition for specific protocols + auto protocol = color_db->getWithDefault("protocol", "none"); + if (protocol == "seed water") { + if (BoundaryCondition != 0 && BoundaryCondition != 5) { + BoundaryCondition = 0; + if (rank == 0) + printf("WARNING: protocol (seed water) supports only full " + "periodic boundary condition \n"); + } + domain_db->putScalar("BC", BoundaryCondition); + } else if (protocol == "open connected oil") { + if (BoundaryCondition != 0 && BoundaryCondition != 5) { + BoundaryCondition = 0; + if (rank == 0) + printf("WARNING: protocol (open connected oil) supports only " + "full periodic boundary condition \n"); + } + domain_db->putScalar("BC", BoundaryCondition); + } else if (protocol == "shell aggregation") { + if (BoundaryCondition != 0 && BoundaryCondition != 5) { + BoundaryCondition = 0; + if (rank == 0) + printf("WARNING: protocol (shell aggregation) supports only " + "full periodic boundary condition \n"); + } + 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 (rank == 0) + printf("Using core flooding protocol \n"); + 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); + } } -void ScaLBL_ColorModel::SetDomain(){ - Dm = std::shared_ptr(new Domain(domain_db,comm)); // full domain for analysis - Mask = std::shared_ptr(new Domain(domain_db,comm)); // mask domain removes immobile phases - // domain parameters - Nx = Dm->Nx; - Ny = Dm->Ny; - Nz = Dm->Nz; - Lx = Dm->Lx; - Ly = Dm->Ly; - Lz = Dm->Lz; - N = Nx*Ny*Nz; - id = new signed char [N]; - for (int i=0; iid[i] = 1; // initialize this way - //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object - Averages = std::shared_ptr ( new SubPhase(Dm) ); // TwoPhase analysis object - comm.barrier(); - Dm->CommInit(); - comm.barrier(); - // Read domain parameters - rank = Dm->rank(); - nprocx = Dm->nprocx(); - nprocy = Dm->nprocy(); - nprocz = Dm->nprocz(); + +void ScaLBL_ColorModel::SetDomain() { + Dm = std::shared_ptr( + new Domain(domain_db, comm)); // full domain for analysis + Mask = std::shared_ptr( + new Domain(domain_db, comm)); // mask domain removes immobile phases + // domain parameters + Nx = Dm->Nx; + Ny = Dm->Ny; + Nz = Dm->Nz; + Lx = Dm->Lx; + Ly = Dm->Ly; + Lz = Dm->Lz; + N = Nx * Ny * Nz; + id = new signed char[N]; + for (int i = 0; i < Nx * Ny * Nz; i++) + Dm->id[i] = 1; // initialize this way + //Averages = std::shared_ptr ( new TwoPhase(Dm) ); // TwoPhase analysis object + Averages = + std::shared_ptr(new SubPhase(Dm)); // TwoPhase analysis object + comm.barrier(); + Dm->CommInit(); + comm.barrier(); + // Read domain parameters + rank = Dm->rank(); + nprocx = Dm->nprocx(); + nprocy = Dm->nprocy(); + nprocz = Dm->nprocz(); } -void ScaLBL_ColorModel::ReadInput(){ - - sprintf(LocalRankString,"%05d",rank); - sprintf(LocalRankFilename,"%s%s","ID.",LocalRankString); - sprintf(LocalRestartFile,"%s%s","Restart.",LocalRankString); - - if (color_db->keyExists( "image_sequence" )){ - auto ImageList = color_db->getVector( "image_sequence"); - int IMAGE_INDEX = color_db->getWithDefault( "image_index", 0 ); - std::string first_image = ImageList[IMAGE_INDEX]; - Mask->Decomp(first_image); - IMAGE_INDEX++; - } - else if (domain_db->keyExists( "GridFile" )){ +void ScaLBL_ColorModel::ReadInput() { + + sprintf(LocalRankString, "%05d", rank); + sprintf(LocalRankFilename, "%s%s", "ID.", LocalRankString); + sprintf(LocalRestartFile, "%s%s", "Restart.", LocalRankString); + + if (color_db->keyExists("image_sequence")) { + auto ImageList = color_db->getVector("image_sequence"); + int IMAGE_INDEX = color_db->getWithDefault("image_index", 0); + std::string first_image = ImageList[IMAGE_INDEX]; + Mask->Decomp(first_image); + IMAGE_INDEX++; + } else if (domain_db->keyExists("GridFile")) { // Read the local domain data - auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD ); + auto input_id = readMicroCT(*domain_db, MPI_COMM_WORLD); // Fill the halo (assuming GCW of 1) - array size0 = { (int) input_id.size(0), (int) input_id.size(1), (int) input_id.size(2) }; - ArraySize size1 = { (size_t) Mask->Nx, (size_t) Mask->Ny, (size_t) Mask->Nz }; - ASSERT( (int) size1[0] == size0[0]+2 && (int) size1[1] == size0[1]+2 && (int) size1[2] == size0[2]+2 ); - fillHalo fill( MPI_COMM_WORLD, Mask->rank_info, size0, { 1, 1, 1 }, 0, 1 ); + array size0 = {(int)input_id.size(0), (int)input_id.size(1), + (int)input_id.size(2)}; + ArraySize size1 = {(size_t)Mask->Nx, (size_t)Mask->Ny, + (size_t)Mask->Nz}; + ASSERT((int)size1[0] == size0[0] + 2 && (int)size1[1] == size0[1] + 2 && + (int)size1[2] == size0[2] + 2); + fillHalo fill(MPI_COMM_WORLD, Mask->rank_info, size0, + {1, 1, 1}, 0, 1); Array id_view; - id_view.viewRaw( size1, Mask->id.data() ); - fill.copy( input_id, id_view ); - fill.fill( id_view ); - } - else if (domain_db->keyExists( "Filename" )){ - auto Filename = domain_db->getScalar( "Filename" ); - Mask->Decomp(Filename); - } - else{ - Mask->ReadIDs(); - } - for (int i=0; iid[i]; // save what was read - - // Generate the signed distance map - // Initialize the domain and communication - Array id_solid(Nx,Ny,Nz); - // Solve for the position of the solid phase - for (int k=0;kid[n]; - if (label > 0) id_solid(i,j,k) = 1; - else id_solid(i,j,k) = 0; - } - } - } - // Initialize the signed distance function - for (int k=0;kSDs(i,j,k) = 2.0*double(id_solid(i,j,k))-1.0; - } - } - } -// MeanFilter(Averages->SDs); - Minkowski Solid(Dm); - if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n"); - CalcDist(Averages->SDs,id_solid,*Mask); - Solid.ComputeScalar(Averages->SDs,0.0); - /* save averages */ - Averages->solid.V = Solid.Vi; - Averages->solid.A = Solid.Ai; - Averages->solid.H = Solid.Ji; - Averages->solid.X = Solid.Xi; - Averages->gsolid.V = Solid.Vi_global; - Averages->gsolid.A = Solid.Ai_global; - Averages->gsolid.H = Solid.Ji_global; - Averages->gsolid.X = Solid.Xi_global; - /* write to file */ - if (rank == 0) { - FILE *SOLID = fopen("solid.csv","w"); - fprintf(SOLID,"Vs As Hs Xs\n"); - fprintf(SOLID,"%.8g %.8g %.8g %.8g\n",Solid.Vi_global,Solid.Ai_global,Solid.Ji_global,Solid.Xi_global); - fclose(SOLID); - } - if (rank == 0) cout << "Domain set." << endl; - - Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta); + id_view.viewRaw(size1, Mask->id.data()); + fill.copy(input_id, id_view); + fill.fill(id_view); + } else if (domain_db->keyExists("Filename")) { + auto Filename = domain_db->getScalar("Filename"); + Mask->Decomp(Filename); + } else { + Mask->ReadIDs(); + } + for (int i = 0; i < Nx * Ny * Nz; i++) + id[i] = Mask->id[i]; // save what was read + + // Generate the signed distance map + // Initialize the domain and communication + Array id_solid(Nx, Ny, Nz); + // Solve for the position of the solid phase + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + int n = k * Nx * Ny + j * Nx + i; + // Initialize the solid phase + signed char label = Mask->id[n]; + if (label > 0) + id_solid(i, j, k) = 1; + else + id_solid(i, j, k) = 0; + } + } + } + // Initialize the signed distance function + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + // Initialize distance to +/- 1 + Averages->SDs(i, j, k) = 2.0 * double(id_solid(i, j, k)) - 1.0; + } + } + } + // MeanFilter(Averages->SDs); + Minkowski Solid(Dm); + if (rank == 0) + printf("Initialized solid phase -- Converting to Signed Distance " + "function \n"); + CalcDist(Averages->SDs, id_solid, *Mask); + Solid.ComputeScalar(Averages->SDs, 0.0); + /* save averages */ + Averages->solid.V = Solid.Vi; + Averages->solid.A = Solid.Ai; + Averages->solid.H = Solid.Ji; + Averages->solid.X = Solid.Xi; + Averages->gsolid.V = Solid.Vi_global; + Averages->gsolid.A = Solid.Ai_global; + Averages->gsolid.H = Solid.Ji_global; + Averages->gsolid.X = Solid.Xi_global; + /* write to file */ + if (rank == 0) { + FILE *SOLID = fopen("solid.csv", "w"); + fprintf(SOLID, "Vs As Hs Xs\n"); + fprintf(SOLID, "%.8g %.8g %.8g %.8g\n", Solid.Vi_global, + Solid.Ai_global, Solid.Ji_global, Solid.Xi_global); + fclose(SOLID); + } + if (rank == 0) + cout << "Domain set." << endl; + + Averages->SetParams(rhoA, rhoB, tauA, tauB, Fx, Fy, Fz, alpha, beta); } -void ScaLBL_ColorModel::AssignComponentLabels(double *phase) -{ - size_t NLABELS=0; - signed char VALUE=0; - double AFFINITY=0.f; +void ScaLBL_ColorModel::AssignComponentLabels(double *phase) { + size_t NLABELS = 0; + signed char VALUE = 0; + double AFFINITY = 0.f; - auto LabelList = color_db->getVector( "ComponentLabels" ); - auto AffinityList = color_db->getVector( "ComponentAffinity" ); - auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); + auto LabelList = color_db->getVector("ComponentLabels"); + auto AffinityList = color_db->getVector("ComponentAffinity"); + auto WettingConvention = + color_db->getWithDefault("WettingConvention", "none"); - NLABELS=LabelList.size(); - if (NLABELS != AffinityList.size()){ - ERROR("Error: ComponentLabels and ComponentAffinity must be the same length! \n"); - } + NLABELS = LabelList.size(); + if (NLABELS != AffinityList.size()) { + ERROR("Error: ComponentLabels and ComponentAffinity must be the same " + "length! \n"); + } - if (WettingConvention == "SCAL"){ - for (size_t idx=0; idxid[n] = 0; // set mask to zero since this is an immobile component - } - } - // fluid labels are reserved - if (VALUE == 1) AFFINITY=1.0; - else if (VALUE == 2) AFFINITY=-1.0; - phase[n] = AFFINITY; - } - } - } + double *label_count; + double *label_count_global; + label_count = new double[NLABELS]; + label_count_global = new double[NLABELS]; + // Assign the labels + for (size_t idx = 0; idx < NLABELS; idx++) + label_count[idx] = 0; - // Set Dm to match Mask - for (int i=0; iid[i] = Mask->id[i]; - - for (size_t idx=0; idxComm.sumReduce( label_count[idx]); + for (int k = 0; k < Nz; k++) { + for (int j = 0; j < Ny; j++) { + for (int i = 0; i < Nx; i++) { + int n = k * Nx * Ny + j * Nx + i; + VALUE = id[n]; + // Assign the affinity from the paired list + for (unsigned int idx = 0; idx < NLABELS; idx++) { + if (VALUE == LabelList[idx]) { + AFFINITY = AffinityList[idx]; + label_count[idx] += 1.0; + idx = NLABELS; + } + } + // fluid labels are reserved + if (VALUE == 1) + AFFINITY = 1.0; + else if (VALUE == 2) + AFFINITY = -1.0; + phase[n] = AFFINITY; + } + } + } - if (rank==0){ - printf("Component labels: %lu \n",NLABELS); - for (unsigned int idx=0; idxid[i] = Mask->id[i]; + for (size_t idx = 0; idx < NLABELS; idx++) + label_count_global[idx] = Dm->Comm.sumReduce(label_count[idx]); + + if (rank == 0) { + printf("Component labels: %lu \n", NLABELS); + for (unsigned int idx = 0; idx < NLABELS; idx++) { + VALUE = LabelList[idx]; + AFFINITY = AffinityList[idx]; + double volume_fraction = + double(label_count_global[idx]) / + double((Nx - 2) * (Ny - 2) * (Nz - 2) * nprocs); + printf(" label=%d, affinity=%f, volume fraction==%f\n", VALUE, + AFFINITY, volume_fraction); + } + } } +void ScaLBL_ColorModel::Create() { + /* + * This function creates the variables needed to run a LBM + */ + //......................................................... + // don't perform computations at the eight corners + //id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0; + //id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0; -void ScaLBL_ColorModel::Create(){ - /* - * This function creates the variables needed to run a LBM - */ - //......................................................... - // don't perform computations at the eight corners - //id[0] = id[Nx-1] = id[(Ny-1)*Nx] = id[(Ny-1)*Nx + Nx-1] = 0; - //id[(Nz-1)*Nx*Ny] = id[(Nz-1)*Nx*Ny+Nx-1] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx] = id[(Nz-1)*Nx*Ny+(Ny-1)*Nx + Nx-1] = 0; + //......................................................... + // Initialize communication structures in averaging domain + for (int i = 0; i < Nx * Ny * Nz; i++) + Dm->id[i] = Mask->id[i]; + Mask->CommInit(); + Np = Mask->PoreCount(); + //........................................................................... + if (rank == 0) + printf("Create ScaLBL_Communicator \n"); + // Create a communicator for the device (will use optimized layout) + // ScaLBL_Communicator ScaLBL_Comm(Mask); // original + ScaLBL_Comm = + std::shared_ptr(new ScaLBL_Communicator(Mask)); + ScaLBL_Comm_Regular = + std::shared_ptr(new ScaLBL_Communicator(Mask)); - //......................................................... - // Initialize communication structures in averaging domain - for (int i=0; iid[i] = Mask->id[i]; - Mask->CommInit(); - Np=Mask->PoreCount(); - //........................................................................... - if (rank==0) printf ("Create ScaLBL_Communicator \n"); - // Create a communicator for the device (will use optimized layout) - // ScaLBL_Communicator ScaLBL_Comm(Mask); // original - ScaLBL_Comm = std::shared_ptr(new ScaLBL_Communicator(Mask)); - ScaLBL_Comm_Regular = std::shared_ptr(new ScaLBL_Communicator(Mask)); + int Npad = (Np / 16 + 2) * 16; + if (rank == 0) + printf("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); + Map.resize(Nx, Ny, Nz); + Map.fill(-2); + auto neighborList = new int[18 * Npad]; + Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map, neighborList, + Mask->id.data(), Np, 1); + comm.barrier(); - int Npad=(Np/16 + 2)*16; - if (rank==0) printf ("Set up memory efficient layout, %i | %i | %i \n", Np, Npad, N); - Map.resize(Nx,Ny,Nz); Map.fill(-2); - auto neighborList= new int[18*Npad]; - Np = ScaLBL_Comm->MemoryOptimizedLayoutAA(Map,neighborList,Mask->id.data(),Np,1); - comm.barrier(); + //........................................................................... + // MAIN VARIABLES ALLOCATED HERE + //........................................................................... + // LBM variables + if (rank == 0) + printf("Allocating distributions \n"); + //......................device distributions................................. + dist_mem_size = Np * sizeof(double); + neighborSize = 18 * (Np * sizeof(int)); + //........................................................................... + ScaLBL_AllocateDeviceMemory((void **)&NeighborList, neighborSize); + ScaLBL_AllocateDeviceMemory((void **)&dvcMap, sizeof(int) * Np); + ScaLBL_AllocateDeviceMemory((void **)&fq, 19 * dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **)&Aq, 7 * dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **)&Bq, 7 * dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **)&Den, 2 * dist_mem_size); + ScaLBL_AllocateDeviceMemory((void **)&Phi, sizeof(double) * Nx * Ny * Nz); + ScaLBL_AllocateDeviceMemory((void **)&Pressure, sizeof(double) * Np); + ScaLBL_AllocateDeviceMemory((void **)&Velocity, 3 * sizeof(double) * Np); + ScaLBL_AllocateDeviceMemory((void **)&ColorGrad, 3 * sizeof(double) * Np); + //........................................................................... + // Update GPU data structures + if (rank == 0) + printf("Setting up device map and neighbor list \n"); + fflush(stdout); + int *TmpMap; + TmpMap = new int[Np]; + for (int k = 1; k < Nz - 1; k++) { + for (int j = 1; j < Ny - 1; j++) { + for (int i = 1; i < Nx - 1; i++) { + int idx = Map(i, j, k); + if (!(idx < 0)) + TmpMap[idx] = k * Nx * Ny + j * Nx + i; + } + } + } + // check that TmpMap is valid + for (int idx = 0; idx < ScaLBL_Comm->LastExterior(); idx++) { + auto n = TmpMap[idx]; + if (n > Nx * Ny * Nz) { + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nx * Ny * Nz - 1; + } + } + for (int idx = ScaLBL_Comm->FirstInterior(); + idx < ScaLBL_Comm->LastInterior(); idx++) { + auto n = TmpMap[idx]; + if (n > Nx * Ny * Nz) { + printf("Bad value! idx=%i \n", n); + TmpMap[idx] = Nx * Ny * Nz - 1; + } + } + ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int) * Np); + ScaLBL_Comm->Barrier(); + delete[] TmpMap; - //........................................................................... - // MAIN VARIABLES ALLOCATED HERE - //........................................................................... - // LBM variables - if (rank==0) printf ("Allocating distributions \n"); - //......................device distributions................................. - dist_mem_size = Np*sizeof(double); - neighborSize=18*(Np*sizeof(int)); - //........................................................................... - ScaLBL_AllocateDeviceMemory((void **) &NeighborList, neighborSize); - ScaLBL_AllocateDeviceMemory((void **) &dvcMap, sizeof(int)*Np); - ScaLBL_AllocateDeviceMemory((void **) &fq, 19*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Aq, 7*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Bq, 7*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Den, 2*dist_mem_size); - ScaLBL_AllocateDeviceMemory((void **) &Phi, sizeof(double)*Nx*Ny*Nz); - ScaLBL_AllocateDeviceMemory((void **) &Pressure, sizeof(double)*Np); - ScaLBL_AllocateDeviceMemory((void **) &Velocity, 3*sizeof(double)*Np); - ScaLBL_AllocateDeviceMemory((void **) &ColorGrad, 3*sizeof(double)*Np); - //........................................................................... - // Update GPU data structures - if (rank==0) printf ("Setting up device map and neighbor list \n"); - fflush(stdout); - int *TmpMap; - TmpMap=new int[Np]; - for (int k=1; kLastExterior(); idx++){ - auto n = TmpMap[idx]; - if (n > Nx*Ny*Nz){ - printf("Bad value! idx=%i \n", n); - TmpMap[idx] = Nx*Ny*Nz-1; - } - } - for (int idx=ScaLBL_Comm->FirstInterior(); idxLastInterior(); idx++){ - auto n = TmpMap[idx]; - if ( n > Nx*Ny*Nz ){ - printf("Bad value! idx=%i \n",n); - TmpMap[idx] = Nx*Ny*Nz-1; - } - } - ScaLBL_CopyToDevice(dvcMap, TmpMap, sizeof(int)*Np); - ScaLBL_Comm->Barrier(); - delete [] TmpMap; - - // copy the neighbor list - ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); - delete [] neighborList; - // initialize phi based on PhaseLabel (include solid component labels) - double *PhaseLabel; - PhaseLabel = new double[N]; - AssignComponentLabels(PhaseLabel); - ScaLBL_CopyToDevice(Phi, PhaseLabel, N*sizeof(double)); - delete [] PhaseLabel; -} + // copy the neighbor list + ScaLBL_CopyToDevice(NeighborList, neighborList, neighborSize); + delete[] neighborList; + // initialize phi based on PhaseLabel (include solid component labels) + double *PhaseLabel; + PhaseLabel = new double[N]; + AssignComponentLabels(PhaseLabel); + ScaLBL_CopyToDevice(Phi, PhaseLabel, N * sizeof(double)); + delete[] PhaseLabel; +} /******************************************************** * AssignComponentLabels * ********************************************************/ -void ScaLBL_ColorModel::Initialize(){ - - /* if both capillary number and flux BC are specified */ - if (color_db->keyExists( "capillary_number" ) && BoundaryCondition == 4){ - double capillary_number = color_db->getScalar( "capillary_number" ); - if (rank==0) printf(" set flux to achieve Ca=%f \n", 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 = Mask->Porosity()*CrossSectionalArea*(Ny-2)*IFT*capillary_number/MuB; - if (rank==0) printf(" flux=%f \n",flux); - } - color_db->putScalar( "flux", flux ); - - if (rank==0) printf ("Initializing distributions \n"); - ScaLBL_D3Q19_Init(fq, Np); - /* - * This function initializes model - */ - if (Restart == true){ - if (rank==0){ - printf("Reading restart file! \n"); - } +void ScaLBL_ColorModel::Initialize() { - // Read in the restart file to CPU buffers - int *TmpMap; - TmpMap = new int[Np]; - - double *cPhi, *cDist, *cDen; - cPhi = new double[N]; - cDen = new double[2*Np]; - cDist = new double[19*Np]; - ScaLBL_CopyToHost(TmpMap, dvcMap, Np*sizeof(int)); - ScaLBL_CopyToHost(cPhi, Phi, N*sizeof(double)); - - ifstream File(LocalRestartFile,ios::binary); - int idx; - double value,va,vb; - for (int n=0; nLastExterior(); 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(); + /* if both capillary number and flux BC are specified */ + if (color_db->keyExists("capillary_number") && BoundaryCondition == 4) { + double capillary_number = + color_db->getScalar("capillary_number"); + if (rank == 0) + printf(" set flux to achieve Ca=%f \n", 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 = Mask->Porosity() * CrossSectionalArea * (Ny - 2) * IFT * + capillary_number / MuB; + if (rank == 0) + printf(" flux=%f \n", flux); + } + color_db->putScalar("flux", flux); - comm.barrier(); - } + if (rank == 0) + printf("Initializing distributions \n"); + ScaLBL_D3Q19_Init(fq, Np); + /* + * This function initializes model + */ + if (Restart == true) { + if (rank == 0) { + printf("Reading restart file! \n"); + } - if (rank==0) printf ("Initializing phase field \n"); - 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); + // Read in the restart file to CPU buffers + int *TmpMap; + TmpMap = new int[Np]; - // establish reservoirs for external bC - 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); - } - } - ScaLBL_CopyToHost(Averages->Phi.data(),Phi,N*sizeof(double)); + double *cPhi, *cDist, *cDen; + cPhi = new double[N]; + cDen = new double[2 * Np]; + cDist = new double[19 * Np]; + ScaLBL_CopyToHost(TmpMap, dvcMap, Np * sizeof(int)); + ScaLBL_CopyToHost(cPhi, Phi, N * sizeof(double)); + + ifstream File(LocalRestartFile, ios::binary); + int idx; + double value, va, vb; + for (int n = 0; n < Np; n++) { + File.read((char *)&va, sizeof(va)); + File.read((char *)&vb, sizeof(vb)); + cDen[n] = va; + cDen[Np + n] = vb; + } + for (int n = 0; n < Np; n++) { + // Read the distributions + for (int q = 0; q < 19; q++) { + File.read((char *)&value, sizeof(value)); + cDist[q * Np + n] = value; + } + } + File.close(); + + for (int n = 0; n < ScaLBL_Comm->LastExterior(); n++) { + va = cDen[n]; + vb = cDen[Np + n]; + value = (va - vb) / (va + vb); + idx = TmpMap[n]; + if (!(idx < 0) && idx < N) + cPhi[idx] = value; + } + for (int n = ScaLBL_Comm->FirstInterior(); + n < ScaLBL_Comm->LastInterior(); n++) { + va = cDen[n]; + vb = cDen[Np + n]; + value = (va - vb) / (va + vb); + idx = TmpMap[n]; + if (!(idx < 0) && idx < N) + cPhi[idx] = value; + } + + // Copy the restart data to the GPU + ScaLBL_CopyToDevice(Den, cDen, 2 * Np * sizeof(double)); + ScaLBL_CopyToDevice(fq, cDist, 19 * Np * sizeof(double)); + ScaLBL_CopyToDevice(Phi, cPhi, N * sizeof(double)); + ScaLBL_Comm->Barrier(); + + comm.barrier(); + } + + if (rank == 0) + printf("Initializing phase field \n"); + 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); + + // establish reservoirs for external bC + 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); + } + } + ScaLBL_CopyToHost(Averages->Phi.data(), Phi, N * sizeof(double)); } -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; - bool TRIGGER_FORCE_RESCALE = 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; - int INITIAL_TIMESTEP = timestep; +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"); + bool Regular = false; + bool RESCALE_FORCE = false; + bool SET_CAPILLARY_NUMBER = false; + bool TRIGGER_FORCE_RESCALE = false; + double tolerance = 0.01; + auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); + 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; + int INITIAL_TIMESTEP = timestep; - double capillary_number = 1.0e-5; - double Ca_previous = 0.0; - double minCa = 8.0e-6; - double maxCa = 1.0; - if (color_db->keyExists( "capillary_number" )){ - capillary_number = color_db->getScalar( "capillary_number" ); - SET_CAPILLARY_NUMBER=true; - maxCa = 2.0*capillary_number; - minCa = 0.8*capillary_number; - } - 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" ); - } - - auto WettingConvention = color_db->getWithDefault( "WettingConvention", "none" ); - - - runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); - auto t1 = std::chrono::system_clock::now(); - int CURRENT_TIMESTEP = 0; - int EXIT_TIMESTEP = min(timestepMax,returntime); - while (timestep < EXIT_TIMESTEP ) { - //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } - PROFILE_START("Update"); - // *************ODD TIMESTEP************* - timestep++; - // Compute the Phase indicator field - // Read for Aq, Bq happens in this routine (requires communication) - ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); - - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL - if (BoundaryCondition > 0 && BoundaryCondition < 5){ - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - // Halo exchange for phase field - ScaLBL_Comm_Regular->SendHalo(Phi); - - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set BCs - if (BoundaryCondition == 3){ - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - if (BoundaryCondition == 4){ - din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - else if (BoundaryCondition == 5){ - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); - - // *************EVEN TIMESTEP************* - timestep++; - // Compute the Phase indicator field - ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); - - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL - // Halo exchange for phase field - if (BoundaryCondition > 0 && BoundaryCondition < 5){ - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set boundary conditions - if (BoundaryCondition == 3){ - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - else if (BoundaryCondition == 4){ - din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - else if (BoundaryCondition == 5){ - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - 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 && BoundaryCondition == 0){ - 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 (isSteady && (Ca > maxCa || Ca < minCa) && SET_CAPILLARY_NUMBER ){ - /* re-run the point if the actual Ca is too far from the target Ca */ - isSteady = false; - RESCALE_FORCE = true; - t1 = std::chrono::system_clock::now(); - CURRENT_TIMESTEP = 0; - timestep = INITIAL_TIMESTEP; - TRIGGER_FORCE_RESCALE = true; - if (rank == 0) printf(" Capillary number missed target value = %f (measured value was Ca = %f) \n ",capillary_number, Ca); - } - - if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP){ - TRIGGER_FORCE_RESCALE = true; - } - - if (TRIGGER_FORCE_RESCALE){ - RESCALE_FORCE = false; - TRIGGER_FORCE_RESCALE = 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); - - /* SCAL file */ - if (WettingConvention == "SCAL"){ - WriteHeader=false; - FILE * scal_log_file = fopen("SCAL.csv","r"); - if (scal_log_file != NULL) - fclose(scal_log_file); - else - WriteHeader=true; - scal_log_file = fopen("SCAL.csv","a"); - if (WriteHeader){ - fprintf(scal_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(scal_log_file,"eff.perm.oil.upper.bound eff.perm.water.lower.bound eff.perm.oil.lower.bound eff.perm.water.upper.bound "); - fprintf(scal_log_file,"cap.pressure cap.pressure.connected pressure.drop Ca M\n"); - } - fprintf(scal_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(scal_log_file,"%.5g %.5g %.5g %.5g ",kAeff_filmA, kBeff_filmA, kAeff_filmB,kBeff_filmB); - fprintf(scal_log_file,"%.5g %.5g %.5g %.5g %.5g\n",pAB,pAB_connected,viscous_pressure_drop,Ca,Mobility); - fclose(scal_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"); - PROFILE_SAVE("lbpm_color_simulator",1); - //************************************************************************ - // Compute the walltime per timestep - auto t2 = std::chrono::system_clock::now(); - double cputime = std::chrono::duration( t2 - t1 ).count() / CURRENT_TIMESTEP; - // Performance obtained from each node - double MLUPS = double(Np)/cputime/1000000; - - if (rank==0) printf("********************************************************\n"); - if (rank==0) printf("CPU time = %f \n", cputime); - if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); - return(MLUPS); - MLUPS *= nprocs; - -} - -void ScaLBL_ColorModel::Run(){ - int nprocs=nprocx*nprocy*nprocz; - const RankInfoStruct rank_info(rank,nprocx,nprocy,nprocz); - int analysis_interval = 1000; // number of timesteps in between in situ analysis - if (analysis_db->keyExists( "analysis_interval" )){ - analysis_interval = analysis_db->getScalar( "analysis_interval" ); - } - - //************ MAIN ITERATION LOOP ***************************************/ - comm.barrier(); - PROFILE_START("Loop"); - //std::shared_ptr analysis_db; - bool Regular = false; - auto current_db = db->cloneDatabase(); - runAnalysis analysis( current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, Map ); - //analysis.createThreads( analysis_method, 4 ); + double capillary_number = 1.0e-5; + double Ca_previous = 0.0; + double minCa = 8.0e-6; + double maxCa = 1.0; + if (color_db->keyExists("capillary_number")) { + capillary_number = color_db->getScalar("capillary_number"); + SET_CAPILLARY_NUMBER = true; + maxCa = 2.0 * capillary_number; + minCa = 0.8 * capillary_number; + } + 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(); - while (timestep < timestepMax ) { - //if ( rank==0 ) { printf("Running timestep %i (%i MB)\n",timestep+1,(int)(Utilities::getMemoryUsage()/1048576)); } - PROFILE_START("Update"); - // *************ODD TIMESTEP************* - timestep++; - // Compute the Phase indicator field - // Read for Aq, Bq happens in this routine (requires communication) - ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + int CURRENT_TIMESTEP = 0; + int EXIT_TIMESTEP = min(timestepMax, returntime); + while (timestep < EXIT_TIMESTEP) { + PROFILE_START("Update"); + // *************ODD TIMESTEP************* + timestep++; + // Compute the Phase indicator field + // Read for Aq, Bq happens in this routine (requires communication) + ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, + ScaLBL_Comm->LastExterior(), Np); - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL - if (BoundaryCondition > 0 && BoundaryCondition < 5){ - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - // Halo exchange for phase field - ScaLBL_Comm_Regular->SendHalo(Phi); + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + if (BoundaryCondition > 0 && BoundaryCondition < 5) { + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + // Halo exchange for phase field + ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set BCs - if (BoundaryCondition == 3){ - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - if (BoundaryCondition == 4){ - din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - else if (BoundaryCondition == 5){ - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); + ScaLBL_D3Q19_AAodd_Color( + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set BCs + if (BoundaryCondition == 3) { + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + if (BoundaryCondition == 4) { + din = + ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 5) { + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, + Fx, Fy, Fz, Nx, Nx * Ny, 0, + ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); - // *************EVEN TIMESTEP************* - timestep++; - // Compute the Phase indicator field - ScaLBL_Comm->BiSendD3Q7AA(Aq,Bq); //READ FROM NORMAL - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm->BiRecvD3Q7AA(Aq,Bq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, ScaLBL_Comm->LastExterior(), Np); + // *************EVEN TIMESTEP************* + timestep++; + // Compute the Phase indicator field + ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, + ScaLBL_Comm->LastExterior(), Np); - // Perform the collision operation - ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL - // Halo exchange for phase field - if (BoundaryCondition > 0 && BoundaryCondition < 5){ - ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); - ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); - } - ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); - ScaLBL_Comm_Regular->RecvHalo(Phi); - ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE - ScaLBL_Comm->Barrier(); - // Set boundary conditions - if (BoundaryCondition == 3){ - ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - else if (BoundaryCondition == 4){ - din = ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); - ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); - } - else if (BoundaryCondition == 5){ - ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); - ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); - } - ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, tauA, tauB, - alpha, beta, Fx, Fy, Fz, Nx, Nx*Ny, 0, ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->Barrier(); - //************************************************************************ - PROFILE_STOP("Update"); + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL + // Halo exchange for phase field + if (BoundaryCondition > 0 && BoundaryCondition < 5) { + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + ScaLBL_Comm_Regular->SendHalo(Phi); + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, + Nx * Ny, ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set boundary conditions + if (BoundaryCondition == 3) { + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 4) { + din = + ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 5) { + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + rhoB, tauA, tauB, 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 - if (rank==0 && timestep%analysis_interval == 0 && BoundaryCondition == 4){ - printf("%i %f \n",timestep,din); - } - // Run the analysis - analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, fq, Den ); - } - analysis.finish(); - PROFILE_STOP("Loop"); - PROFILE_SAVE("lbpm_color_simulator",1); - //************************************************************************ - ScaLBL_Comm->Barrier(); - if (rank==0) printf("-------------------------------------------------------------------\n"); - // Compute the walltime per timestep + CURRENT_TIMESTEP += 2; + if (CURRENT_TIMESTEP > MIN_STEADY_TIMESTEPS && BoundaryCondition == 0) { + 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 (isSteady && (Ca > maxCa || Ca < minCa) && + SET_CAPILLARY_NUMBER) { + /* re-run the point if the actual Ca is too far from the target Ca */ + isSteady = false; + RESCALE_FORCE = true; + t1 = std::chrono::system_clock::now(); + CURRENT_TIMESTEP = 0; + timestep = INITIAL_TIMESTEP; + TRIGGER_FORCE_RESCALE = true; + if (rank == 0) + printf(" Capillary number missed target value = %f " + "(measured value was Ca = %f) \n ", + capillary_number, Ca); + } + + if (RESCALE_FORCE == true && SET_CAPILLARY_NUMBER == true && + CURRENT_TIMESTEP > RESCALE_FORCE_AFTER_TIMESTEP) { + TRIGGER_FORCE_RESCALE = true; + } + + if (TRIGGER_FORCE_RESCALE) { + RESCALE_FORCE = false; + TRIGGER_FORCE_RESCALE = 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); + + // Saturation normalized effective permeability to account for decoupled phases and + // effective porosity. + double kAeff_connected_low = + (1.0 - current_saturation) * h * h * muA * + flow_rate_A_connected / (force_mag); + double kBeff_connected_low = current_saturation * 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); + + // Saturation normalized effective permeability to account for decoupled phases and + // effective porosity. + double kAeff_low = (1.0 - current_saturation) * h * h * + muA * (flow_rate_A) / (force_mag); + double kBeff_low = current_saturation * h * h * muB * + (flow_rate_B) / (force_mag); + + double viscous_pressure_drop = + (rhoA * volA + rhoB * volB) * force_mag; + double Mobility = muA / muB; // visc contrast + double eff_pres = + 1.0 / (kAeff + kBeff); // effective pressure drop + + 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 "); + fprintf(kr_log_file, "eff.perm.oil.upper.bound " + "eff.perm.water.upper.bound "); + fprintf(kr_log_file, "eff.perm.oil.lower.bound " + "eff.perm.water.lower.bound "); + fprintf(kr_log_file, + "eff.perm.oil.connected.upper.bound " + "eff.perm.water.connected.upper.bound "); + fprintf(kr_log_file, + "eff.perm.oil.connected.lower.bound " + "eff.perm.water.connected.lower.bound "); + fprintf(kr_log_file, "eff.perm.oil.disconnected " + "eff.perm.water.disconnected "); + fprintf(kr_log_file, + "cap.pressure cap.pressure.connected " + "pressure.drop Ca M eff.pressure\n"); + } + fprintf(kr_log_file, "%i %.5g ", CURRENT_TIMESTEP, + current_saturation); + fprintf(kr_log_file, "%.5g %.5g ", kAeff, kBeff); + fprintf(kr_log_file, "%.5g %.5g ", kAeff_low, kBeff_low); + fprintf(kr_log_file, "%.5g %.5g ", kAeff_connected, + kBeff_connected); + fprintf(kr_log_file, "%.5g %.5g ", kAeff_connected_low, + kBeff_connected_low); + fprintf(kr_log_file, "%.5g %.5g ", kAeff_disconnected, + kBeff_disconnected); + fprintf(kr_log_file, "%.5g %.5g %.5g %.5g %.5g ", pAB, + pAB_connected, viscous_pressure_drop, Ca, Mobility); + fprintf(kr_log_file, "%.5g\n", eff_pres); + fclose(kr_log_file); + + if (WettingConvention == "SCAL"){ + WriteHeader = false; + FILE *scal_log_file = fopen("SCAL.csv", "r"); + if (scal_log_file != NULL) + fclose(scal_log_file); + else + WriteHeader = true; + scal_log_file = fopen("SCAL.csv", "a"); + if (WriteHeader) { + fprintf(scal_log_file, "timesteps sat.water "); + fprintf(scal_log_file, "eff.perm.oil.upper.bound " + "eff.perm.water.upper.bound "); + fprintf(scal_log_file, "eff.perm.oil.lower.bound " + "eff.perm.water.lower.bound "); + fprintf(scal_log_file, + "eff.perm.oil.connected.upper.bound " + "eff.perm.water.connected.upper.bound "); + fprintf(scal_log_file, + "eff.perm.oil.connected.lower.bound " + "eff.perm.water.connected.lower.bound "); + fprintf(scal_log_file, "eff.perm.oil.disconnected " + "eff.perm.water.disconnected "); + fprintf(scal_log_file, + "cap.pressure cap.pressure.connected " + "pressure.drop Ca M eff.pressure\n"); + } + fprintf(scal_log_file, "%i %.5g ", CURRENT_TIMESTEP, + current_saturation); + fprintf(scal_log_file, "%.5g %.5g ", kAeff, kBeff); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_low, kBeff_low); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected, + kBeff_connected); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected_low, + kBeff_connected_low); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_disconnected, + kBeff_disconnected); + fprintf(scal_log_file, "%.5g %.5g %.5g %.5g %.5g ", pAB, + pAB_connected, viscous_pressure_drop, Ca, Mobility); + fprintf(scal_log_file, "%.5g\n", eff_pres); + fclose(scal_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"); + PROFILE_SAVE("lbpm_color_simulator", 1); + //************************************************************************ + // Compute the walltime per timestep auto t2 = std::chrono::system_clock::now(); - double cputime = std::chrono::duration( t2 - t1 ).count() / timestep; - // Performance obtained from each node - double MLUPS = double(Np)/cputime/1000000; + double cputime = + std::chrono::duration(t2 - t1).count() / CURRENT_TIMESTEP; + // Performance obtained from each node + double MLUPS = double(Np) / cputime / 1000000; - if (rank==0) printf("********************************************************\n"); - if (rank==0) printf("CPU time = %f \n", cputime); - if (rank==0) printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); - MLUPS *= nprocs; - if (rank==0) printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); - if (rank==0) printf("********************************************************\n"); - - // ************************************************************************ + if (rank == 0) + printf("********************************************************\n"); + if (rank == 0) + printf("CPU time = %f \n", cputime); + if (rank == 0) + printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + return (MLUPS); + MLUPS *= nprocs; } -void ScaLBL_ColorModel::WriteDebug(){ - // Copy back final phase indicator field and convert to regular layout - DoubleArray PhaseField(Nx,Ny,Nz); - //ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); - ScaLBL_CopyToHost(PhaseField.data(), Phi, sizeof(double)*N); +void ScaLBL_ColorModel::Run() { + int nprocs = nprocx * nprocy * nprocz; + const RankInfoStruct rank_info(rank, nprocx, nprocy, nprocz); + int analysis_interval = + 1000; // number of timesteps in between in situ analysis + if (analysis_db->keyExists("analysis_interval")) { + analysis_interval = analysis_db->getScalar("analysis_interval"); + } - FILE *OUTFILE; - sprintf(LocalRankFilename,"Phase.%05i.raw",rank); - OUTFILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,OUTFILE); - fclose(OUTFILE); + //************ MAIN ITERATION LOOP ***************************************/ + comm.barrier(); + PROFILE_START("Loop"); + //std::shared_ptr analysis_db; + bool Regular = false; + auto current_db = db->cloneDatabase(); + runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular, + Map); + //analysis.createThreads( analysis_method, 4 ); + auto t1 = std::chrono::system_clock::now(); + while (timestep < timestepMax) { + PROFILE_START("Update"); - ScaLBL_Comm->RegularLayout(Map,&Den[0],PhaseField); - FILE *AFILE; - sprintf(LocalRankFilename,"A.%05i.raw",rank); - AFILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,AFILE); - fclose(AFILE); + // *************ODD TIMESTEP************* + timestep++; + // Compute the Phase indicator field + // Read for Aq, Bq happens in this routine (requires communication) + ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAodd_PhaseField(NeighborList, dvcMap, Aq, Bq, Den, Phi, 0, + ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->RegularLayout(Map,&Den[Np],PhaseField); - FILE *BFILE; - sprintf(LocalRankFilename,"B.%05i.raw",rank); - BFILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,BFILE); - fclose(BFILE); + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL + if (BoundaryCondition > 0 && BoundaryCondition < 5) { + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + // Halo exchange for phase field + ScaLBL_Comm_Regular->SendHalo(Phi); - ScaLBL_Comm->RegularLayout(Map,Pressure,PhaseField); - FILE *PFILE; - sprintf(LocalRankFilename,"Pressure.%05i.raw",rank); - PFILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,PFILE); - fclose(PFILE); + ScaLBL_D3Q19_AAodd_Color( + NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, rhoB, + tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, Nx * Ny, + ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set BCs + if (BoundaryCondition == 3) { + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } + if (BoundaryCondition == 4) { + din = + ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 5) { + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAodd_Color(NeighborList, dvcMap, fq, Aq, Bq, Den, Phi, + Velocity, rhoA, rhoB, tauA, tauB, alpha, beta, + Fx, Fy, Fz, Nx, Nx * Ny, 0, + ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); - ScaLBL_Comm->RegularLayout(Map,&Velocity[0],PhaseField); - FILE *VELX_FILE; - sprintf(LocalRankFilename,"Velocity_X.%05i.raw",rank); - VELX_FILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,VELX_FILE); - fclose(VELX_FILE); + // *************EVEN TIMESTEP************* + timestep++; + // Compute the Phase indicator field + ScaLBL_Comm->BiSendD3Q7AA(Aq, Bq); //READ FROM NORMAL + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, + ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm->BiRecvD3Q7AA(Aq, Bq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + ScaLBL_D3Q7_AAeven_PhaseField(dvcMap, Aq, Bq, Den, Phi, 0, + ScaLBL_Comm->LastExterior(), Np); - ScaLBL_Comm->RegularLayout(Map,&Velocity[Np],PhaseField); - FILE *VELY_FILE; - sprintf(LocalRankFilename,"Velocity_Y.%05i.raw",rank); - VELY_FILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,VELY_FILE); - fclose(VELY_FILE); + // Perform the collision operation + ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL + // Halo exchange for phase field + if (BoundaryCondition > 0 && BoundaryCondition < 5) { + ScaLBL_Comm->Color_BC_z(dvcMap, Phi, Den, inletA, inletB); + ScaLBL_Comm->Color_BC_Z(dvcMap, Phi, Den, outletA, outletB); + } + ScaLBL_Comm_Regular->SendHalo(Phi); + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, + Nx * Ny, ScaLBL_Comm->FirstInterior(), + ScaLBL_Comm->LastInterior(), Np); + ScaLBL_Comm_Regular->RecvHalo(Phi); + ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE + ScaLBL_Comm->Barrier(); + // Set boundary conditions + if (BoundaryCondition == 3) { + ScaLBL_Comm->D3Q19_Pressure_BC_z(NeighborList, fq, din, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 4) { + din = + ScaLBL_Comm->D3Q19_Flux_BC_z(NeighborList, fq, flux, timestep); + ScaLBL_Comm->D3Q19_Pressure_BC_Z(NeighborList, fq, dout, timestep); + } else if (BoundaryCondition == 5) { + ScaLBL_Comm->D3Q19_Reflection_BC_z(fq); + ScaLBL_Comm->D3Q19_Reflection_BC_Z(fq); + } + ScaLBL_D3Q19_AAeven_Color(dvcMap, fq, Aq, Bq, Den, Phi, Velocity, rhoA, + rhoB, tauA, tauB, alpha, beta, Fx, Fy, Fz, Nx, + Nx * Ny, 0, ScaLBL_Comm->LastExterior(), Np); + ScaLBL_Comm->Barrier(); + //************************************************************************ + PROFILE_STOP("Update"); - ScaLBL_Comm->RegularLayout(Map,&Velocity[2*Np],PhaseField); - FILE *VELZ_FILE; - sprintf(LocalRankFilename,"Velocity_Z.%05i.raw",rank); - VELZ_FILE = fopen(LocalRankFilename,"wb"); - fwrite(PhaseField.data(),8,N,VELZ_FILE); - fclose(VELZ_FILE); + if (rank == 0 && timestep % analysis_interval == 0 && + BoundaryCondition == 4) { + printf("%i %f \n", timestep, din); + } + // Run the analysis + analysis.basic(timestep, current_db, *Averages, Phi, Pressure, Velocity, + fq, Den); + } + analysis.finish(); + PROFILE_STOP("Loop"); + PROFILE_SAVE("lbpm_color_simulator", 1); + //************************************************************************ + ScaLBL_Comm->Barrier(); + if (rank == 0) + printf("---------------------------------------------------------------" + "----\n"); + // Compute the walltime per timestep + auto t2 = std::chrono::system_clock::now(); + double cputime = std::chrono::duration(t2 - t1).count() / timestep; + // Performance obtained from each node + double MLUPS = double(Np) / cputime / 1000000; + if (rank == 0) + printf("********************************************************\n"); + if (rank == 0) + printf("CPU time = %f \n", cputime); + if (rank == 0) + printf("Lattice update rate (per core)= %f MLUPS \n", MLUPS); + MLUPS *= nprocs; + if (rank == 0) + printf("Lattice update rate (total)= %f MLUPS \n", MLUPS); + if (rank == 0) + printf("********************************************************\n"); + + // ************************************************************************ +} + +void ScaLBL_ColorModel::WriteDebug() { + // Copy back final phase indicator field and convert to regular layout + DoubleArray PhaseField(Nx, Ny, Nz); + //ScaLBL_Comm->RegularLayout(Map,Phi,PhaseField); + ScaLBL_CopyToHost(PhaseField.data(), Phi, sizeof(double) * N); + + FILE *OUTFILE; + sprintf(LocalRankFilename, "Phase.%05i.raw", rank); + OUTFILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, OUTFILE); + fclose(OUTFILE); + + ScaLBL_Comm->RegularLayout(Map, &Den[0], PhaseField); + FILE *AFILE; + sprintf(LocalRankFilename, "A.%05i.raw", rank); + AFILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, AFILE); + fclose(AFILE); + + ScaLBL_Comm->RegularLayout(Map, &Den[Np], PhaseField); + FILE *BFILE; + sprintf(LocalRankFilename, "B.%05i.raw", rank); + BFILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, BFILE); + fclose(BFILE); + + ScaLBL_Comm->RegularLayout(Map, Pressure, PhaseField); + FILE *PFILE; + sprintf(LocalRankFilename, "Pressure.%05i.raw", rank); + PFILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, PFILE); + fclose(PFILE); + + ScaLBL_Comm->RegularLayout(Map, &Velocity[0], PhaseField); + FILE *VELX_FILE; + sprintf(LocalRankFilename, "Velocity_X.%05i.raw", rank); + VELX_FILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, VELX_FILE); + fclose(VELX_FILE); + + ScaLBL_Comm->RegularLayout(Map, &Velocity[Np], PhaseField); + FILE *VELY_FILE; + sprintf(LocalRankFilename, "Velocity_Y.%05i.raw", rank); + VELY_FILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, VELY_FILE); + fclose(VELY_FILE); + + ScaLBL_Comm->RegularLayout(Map, &Velocity[2 * Np], PhaseField); + FILE *VELZ_FILE; + sprintf(LocalRankFilename, "Velocity_Z.%05i.raw", rank); + VELZ_FILE = fopen(LocalRankFilename, "wb"); + fwrite(PhaseField.data(), 8, N, VELZ_FILE); + fclose(VELZ_FILE); } From a65ceef7d57dd70be6d59db57ece20df111a8ca3 Mon Sep 17 00:00:00 2001 From: James McClure Date: Thu, 9 Dec 2021 10:16:47 -0500 Subject: [PATCH 8/8] try 2 --- models/ColorModel.cpp | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index f63d7ee5..7a356fc0 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -960,32 +960,24 @@ double ScaLBL_ColorModel::Run(int returntime) { fprintf(scal_log_file, "timesteps sat.water "); fprintf(scal_log_file, "eff.perm.oil.upper.bound " "eff.perm.water.upper.bound "); - fprintf(scal_log_file, "eff.perm.oil.lower.bound " + fprintf(scal_log_file, + "eff.perm.oil.lower.bound " "eff.perm.water.lower.bound "); - fprintf(scal_log_file, - "eff.perm.oil.connected.upper.bound " - "eff.perm.water.connected.upper.bound "); - fprintf(scal_log_file, - "eff.perm.oil.connected.lower.bound " - "eff.perm.water.connected.lower.bound "); fprintf(scal_log_file, "eff.perm.oil.disconnected " "eff.perm.water.disconnected "); fprintf(scal_log_file, "cap.pressure cap.pressure.connected " - "pressure.drop Ca M eff.pressure\n"); + "Ca eff.pressure\n"); } fprintf(scal_log_file, "%i %.5g ", CURRENT_TIMESTEP, current_saturation); - fprintf(scal_log_file, "%.5g %.5g ", kAeff, kBeff); fprintf(scal_log_file, "%.5g %.5g ", kAeff_low, kBeff_low); - fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected, - kBeff_connected); fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected_low, kBeff_connected_low); fprintf(scal_log_file, "%.5g %.5g ", kAeff_disconnected, kBeff_disconnected); - fprintf(scal_log_file, "%.5g %.5g %.5g %.5g %.5g ", pAB, - pAB_connected, viscous_pressure_drop, Ca, Mobility); + fprintf(scal_log_file, "%.5g %.5g %.5g ", pAB, + pAB_connected, Ca); fprintf(scal_log_file, "%.5g\n", eff_pres); fclose(scal_log_file);