diff --git a/docs/source/userGuide/models/color/protocols/sw_steady.rst b/docs/source/userGuide/models/color/protocols/sw_steady.rst new file mode 100644 index 00000000..d63d066a --- /dev/null +++ b/docs/source/userGuide/models/color/protocols/sw_steady.rst @@ -0,0 +1,72 @@ +====================================== +Color model -- Sw Steady +====================================== + +The water saturation steady state protocol is identical to the centrifuge protocol, with +the exception that the simulation explicity converges on the saturation state of fluid +A. + +That is, the simulation exits when + +.. math:: + :nowrap: + + $$ + \frac{\left | S_{w, i+1} - S_{w, i} \right |}{S_{w, i}} \le \epsilon + $$ + +averaged over the ``analysis_interval`` where :math:`S_{w,i}` is the saturation of fluid +A at step :math:`i` and :math:`\epsilon` is an allowed convergence threshold, or when +the ``timestepMax`` is reached, whichever occurs first. + +By default, :math:`\epsilon` is set to zero (i.e., the simulation continues through the +``timestepMax``), but can be set via ``tolerance`` within the ``Analysis`` section of +the input file database as shown below. + + +.. code-block:: c + + Color { + protocol = "sw_steady" + timestepMax = 1000000 // maximum timtestep + alpha = 0.005 // 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, -1.0e-5 // body force + din = 1.0 // inlet density (controls pressure) + dout = 1.0 // outlet density (controls pressure) + WettingConvention = "SCAL" // convention for sign of wetting affinity + ComponentLabels = 0, -1, -2 // image labels for solid voxels + ComponentAffinity = 1.0, 1.0, 0.6 // controls the wetting affinity for each label + Restart = false + } + Domain { + Filename = "Bentheimer_LB_sim_intermediate_oil_wet_Sw_0p37.raw" + ReadType = "8bit" // data type + N = 900, 900, 1600 // size of original image + nproc = 2, 2, 2 // process grid + n = 200, 200, 200 // sub-domain size + offset = 300, 300, 300 // offset to read sub-domain + voxel_length = 1.66 // voxel length (in microns) + ReadValues = -2, -1, 0, 1, 2 // labels within the original image + WriteValues = -2, -1, 0, 1, 2 // associated labels to be used by LBPM + BC = 3 // boundary condition type (0 for periodic) + } + Analysis { + analysis_interval = 1000 // logging interval for timelog.csv + subphase_analysis_interval = 5000 // loggging interval for subphase.csv + visualization_interval = 100000 // interval to write visualization files + N_threads = 4 // number of analysis threads (GPU version only) + restart_interval = 1000000 // interval to write restart file + restart_file = "Restart" // base name of restart file + tolerance = 1e-9 // Sw convergence tolerance + } + Visualization { + 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 = false // save pressure field within SILO database + save_velocity = false // save velocity field within SILO database + } diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 88caa2ea..4c1a0ad1 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -172,7 +172,7 @@ void ScaLBL_ColorModel::ReadParams(string filename) { "periodic boundary condition \n"); } domain_db->putScalar("BC", BoundaryCondition); - } else if (protocol == "centrifuge") { + } else if (protocol == "centrifuge" || protocol == "sw_steady") { if (BoundaryCondition != 3) { BoundaryCondition = 3; if (rank == 0) @@ -1116,7 +1116,11 @@ void ScaLBL_ColorModel::Run() { Map); //analysis.createThreads( analysis_method, 4 ); auto t1 = std::chrono::system_clock::now(); - while (timestep < timestepMax) { + + double delta_sw = 1.0; + double sw_prev = -1.0; + + while (timestep < timestepMax && delta_sw > tolerance) { PROFILE_START("Update"); // *************ODD TIMESTEP************* @@ -1213,13 +1217,23 @@ void ScaLBL_ColorModel::Run() { //************************************************************************ PROFILE_STOP("Update"); - 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); + + if (timestep % analysis_interval == 0){ + analysis.finish(); + + double volA = Averages->gnb.V / Dm->Volume; + double volB = Averages->gwb.V / Dm->Volume; + double sw = volB / (volA + volB); + + delta_sw = fabs(sw - sw_prev) / analysis_interval / sw; + if (rank == 0) + printf("t: %d sw: %0.5e dSw/dt: %.5e\n", timestep, sw, delta_sw); + + sw_prev = sw; + } } analysis.finish(); PROFILE_STOP("Loop"); diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp index cafa7a52..977d304b 100644 --- a/tests/lbpm_color_simulator.cpp +++ b/tests/lbpm_color_simulator.cpp @@ -66,7 +66,8 @@ int main( int argc, char **argv ) // structure and allocate variables ColorModel.Initialize(); // initializing the model will set initial conditions for variables - if (SimulationMode == "legacy"){ + auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "default" ); + if (PROTOCOL == "sw_steady") { ColorModel.Run(); } else { @@ -75,7 +76,6 @@ int main( int argc, char **argv ) bool ContinueSimulation = true; /* Variables for simulation protocols */ - auto PROTOCOL = ColorModel.color_db->getWithDefault( "protocol", "default" ); /* image sequence protocol */ int IMAGE_INDEX = 0; int IMAGE_COUNT = 0;