'sw_steady' protocol for color model

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
+   }

models/ColorModel.cpp

@@ -172,7 +172,7 @@ void ScaLBL_ColorModel::ReadParams(string filename) {
                        "periodic boundary condition \n");
         }
         domain_db->putScalar<int>("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");

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<std::string>( "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<std::string>( "protocol", "default" );
 			/* image sequence protocol */
 			int IMAGE_INDEX = 0;
 			int IMAGE_COUNT = 0;