Merge branch 'master' of github.com:OPM/LBPM

Add perm converter

.github/workflows/c-cpp.yml

@@ -29,7 +29,7 @@ jobs:
         sudo apt-get update -y
              
         wget https://bitbucket.org/AdvancedMultiPhysics/tpl-builder/downloads/silo-4.10.2.tar.gz
-        wget https://www.zlib.net/zlib-1.2.11.tar.gz
+        wget https://www.zlib.net/fossils/zlib-1.2.11.tar.gz
         wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.12/src/hdf5-1.8.12.tar.gz
         #wget https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.10/src/hdf5-1.8.10.tar.gz 
        

IO/PackData.cpp

@@ -1,6 +1,6 @@
 #include "IO/PackData.h"
 
-#include <string.h>
+#include <string>
 
 
 /********************************************************

analysis/FlowAdaptor.cpp

@@ -17,7 +17,9 @@ FlowAdaptor::FlowAdaptor(ScaLBL_ColorModel &M) {
     phi_t.fill(0); // time derivative for the phase indicator field
 }
 
-FlowAdaptor::~FlowAdaptor() {}
+FlowAdaptor::~FlowAdaptor() {
+	
+}
 
 double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename) {
     int rank = M.rank;
@@ -71,6 +73,7 @@ double FlowAdaptor::ImageInit(ScaLBL_ColorModel &M, std::string Filename) {
     ScaLBL_CopyToHost(M.Averages->Phi.data(), M.Phi,
                       Nx * Ny * Nz * sizeof(double));
 
+    delete PhaseLabel;
     double saturation = Count / PoreCount;
     return saturation;
 }
@@ -234,6 +237,13 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M) {
     //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
     ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7 * Np * sizeof(double));
     ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7 * Np * sizeof(double));
+    
+    delete Aq_tmp; 
+    delete Bq_tmp;
+    delete Vel_x; 
+    delete Vel_y; 
+    delete Vel_z; 
+    delete Phase;
 
     return (TOTAL_MASS_CHANGE);
 }
@@ -403,7 +413,6 @@ double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M,
             }
         }
     }
-
     if (rank == 0)
         printf("Pathway volume / next largest ganglion %f \n",
                volume_connected / second_biggest);
@@ -585,6 +594,8 @@ double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M,
     //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
     ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7 * Np * sizeof(double));
     ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7 * Np * sizeof(double));
-
+    
+    delete Aq_tmp; 
+    delete Bq_tmp;
     return (mass_loss);
 }

analysis/TwoPhase.cpp

@@ -14,22 +14,6 @@
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
-/*
-  Copyright 2013--2018 James E. McClure, Virginia Polytechnic & State University
-  Copyright Equnior ASA
-
-  This file is part of the Open Porous Media project (OPM).
-  OPM is free software: you can redistribute it and/or modify
-  it under the terms of the GNU General Public License as published by
-  the Free Software Foundation, either version 3 of the License, or
-  (at your option) any later version.
-  OPM is distributed in the hope that it will be useful,
-  but WITHOUT ANY WARRANTY; without even the implied warranty of
-  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-  GNU General Public License for more details.
-  You should have received a copy of the GNU General Public License
-  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-*/
 #include "analysis/TwoPhase.h"
 
 #include "analysis/pmmc.h"
@@ -41,6 +25,8 @@
 #include "IO/MeshDatabase.h"
 #include "IO/Reader.h"
 #include "IO/Writer.h"
+#include "analysis/filters.h"
+
 
 #include <memory>
 
@@ -375,6 +361,7 @@ void TwoPhase::Initialize() {
     wwndnw = 0.0;
     wwnsdnwn = 0.0;
     Jwnwwndnw = 0.0;
+    Xwn = Xns = Xws = 0.0;
 }
 
 void TwoPhase::SetParams(double rhoA, double rhoB, double tauA, double tauB,
@@ -430,28 +417,38 @@ void TwoPhase::UpdateSolid() {
 
 void TwoPhase::UpdateMeshValues() {
     int i, j, k, n;
+    fillHalo<double> fillData(Dm->Comm, Dm->rank_info, {Nx-2,Ny-2,Nz-2}, {1, 1, 1}, 0, 1);
+
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDn);
+    //Dm->CommunicateMeshHalo(SDn);
+    fillData.fill(SDn);
     //...........................................................................
     // Compute the gradients of the phase indicator and signed distance fields
     pmmc_MeshGradient(SDn, SDn_x, SDn_y, SDn_z, Nx, Ny, Nz);
     //...........................................................................
     // Gradient of the phase indicator field
+    fillData.fill(SDn_x);
+    fillData.fill(SDn_y);
+    fillData.fill(SDn_z);
+    fillData.fill(SDs);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDn_x);
+    //Dm->CommunicateMeshHalo(SDn_x);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDn_y);
+    //Dm->CommunicateMeshHalo(SDn_y);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDn_z);
+    //Dm->CommunicateMeshHalo(SDn_z);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDs);
+    //Dm->CommunicateMeshHalo(SDs);
     pmmc_MeshGradient(SDs, SDs_x, SDs_y, SDs_z, Nx, Ny, Nz);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDs_x);
+    fillData.fill(SDs_x);
+    fillData.fill(SDs_y);
+    fillData.fill(SDs_z);
+    //Dm->CommunicateMeshHalo(SDs_x);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDs_y);
+    //Dm->CommunicateMeshHalo(SDs_y);
     //...........................................................................
-    Dm->CommunicateMeshHalo(SDs_z);
+    //Dm->CommunicateMeshHalo(SDs_z);
     //...........................................................................
     // Compute the mesh curvature of the phase indicator field
     pmmc_MeshCurvature(SDn, MeanCurvature, GaussCurvature, Nx, Ny, Nz);
@@ -579,6 +576,7 @@ void TwoPhase::ComputeLocal() {
                     Kwn += pmmc_CubeSurfaceInterpValue(
                         CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i,
                         j, k, n_nw_pts, n_nw_tris);
+                    
                     Jwn += pmmc_CubeSurfaceInterpValue(
                         CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
                         j, k, n_nw_pts, n_nw_tris);
@@ -609,7 +607,7 @@ void TwoPhase::ComputeLocal() {
                     efawns += pmmc_CubeContactAngle(
                         CubeValues, Values, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
                         SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
-
+                    
                     wwnsdnwn += pmmc_CommonCurveSpeed(
                         CubeValues, dPdt, vawns, SDn_x, SDn_y, SDn_z, SDs_x,
                         SDs_y, SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
@@ -715,6 +713,217 @@ void TwoPhase::ComputeLocal() {
     nonwet_morph->ComputeScalar(phase_distance, 0.f);
     //printf("rank=%i completed \n",Dm->rank());
 }
+void TwoPhase::ComputeStatic() {
+    int i, j, k, n, imin, jmin, kmin, kmax;
+    int cube[8][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {1, 1, 0},
+                      {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, {1, 1, 1}};
+
+    kmin = 1;
+    kmax = Nz - 1;
+    imin = jmin = 1;
+    
+    /* set fluid isovalue to "grow" NWP for contact angle measurement */
+    fluid_isovalue = -1.0;
+    
+    string FILENAME = "ContactAngle";
+    
+    char LocalRankString[8];
+    char LocalRankFilename[40];
+    sprintf(LocalRankString, "%05d", Dm->rank());
+    sprintf(LocalRankFilename, "%s%s%s", "ContactAngle.", LocalRankString,".csv");
+    FILE *ANGLES = fopen(LocalRankFilename, "a+");
+    fprintf(ANGLES,"x y z angle\n");
+
+    for (k = kmin; k < kmax; k++) {
+        for (j = jmin; j < Ny - 1; j++) {
+            for (i = imin; i < Nx - 1; i++) {
+                //...........................................................................
+                n_nw_pts = n_ns_pts = n_ws_pts = n_nws_pts = n_local_sol_pts =
+                    n_local_nws_pts = 0;
+                n_nw_tris = n_ns_tris = n_ws_tris = n_nws_seg =
+                    n_local_sol_tris = 0;
+                //...........................................................................
+                // Compute volume averages
+                for (int p = 0; p < 8; p++) {
+                    n = i + cube[p][0] + (j + cube[p][1]) * Nx +
+                        (k + cube[p][2]) * Nx * Ny;
+                    if (Dm->id[n] > 0) {
+                        // 1-D index for this cube corner
+                        // compute the norm of the gradient of the phase indicator field
+                        // Compute the non-wetting phase volume contribution
+                        if (Phase(i + cube[p][0], j + cube[p][1],
+                                  k + cube[p][2]) > 0) {
+                            nwp_volume += 0.125;
+                        } else {
+                            wp_volume += 0.125;
+                        }
+                    }
+                }
+
+                //...........................................................................
+                // Construct the interfaces and common curve
+                pmmc_ConstructLocalCube(
+                    SDs, SDn, solid_isovalue, fluid_isovalue, nw_pts, nw_tris,
+                    Values, ns_pts, ns_tris, ws_pts, ws_tris, local_nws_pts,
+                    nws_pts, nws_seg, local_sol_pts, local_sol_tris,
+                    n_local_sol_tris, n_local_sol_pts, n_nw_pts, n_nw_tris,
+                    n_ws_pts, n_ws_tris, n_ns_tris, n_ns_pts, n_local_nws_pts,
+                    n_nws_pts, n_nws_seg, i, j, k, Nx, Ny, Nz);
+
+                // wn interface averages
+                if (n_nw_pts > 0) {
+                    awn += pmmc_CubeSurfaceOrientation(Gwn, nw_pts, nw_tris,
+                                                       n_nw_tris);
+                    Kwn += pmmc_CubeSurfaceInterpValue(
+                        CubeValues, GaussCurvature, nw_pts, nw_tris, Values, i,
+                        j, k, n_nw_pts, n_nw_tris);
+                    
+                    Jwn += pmmc_CubeSurfaceInterpValue(
+                        CubeValues, MeanCurvature, nw_pts, nw_tris, Values, i,
+                        j, k, n_nw_pts, n_nw_tris);
+                    
+                    Xwn += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
+                                                         n_nw_tris, i, j, k);
+
+                    // Integrate the trimmed mean curvature (hard-coded to use a distance of 4 pixels)
+                    pmmc_CubeTrimSurfaceInterpValues(
+                        CubeValues, MeanCurvature, SDs, nw_pts, nw_tris, Values,
+                        DistanceValues, i, j, k, n_nw_pts, n_nw_tris, trimdist,
+                        dummy, trJwn);
+
+                    pmmc_CubeTrimSurfaceInterpInverseValues(
+                        CubeValues, MeanCurvature, SDs, nw_pts, nw_tris, Values,
+                        DistanceValues, i, j, k, n_nw_pts, n_nw_tris, trimdist,
+                        dummy, trRwn);
+                }
+                // wns common curve averages
+                if (n_local_nws_pts > 0) {
+                    efawns += pmmc_CubeContactAngle(
+                        CubeValues, Values, SDn_x, SDn_y, SDn_z, SDs_x, SDs_y,
+                        SDs_z, local_nws_pts, i, j, k, n_local_nws_pts);
+                    
+                    for (int p = 0; p < n_local_nws_pts; p++) {
+                        // Extract the line segment
+                        Point A = local_nws_pts(p);
+                        double value = Values(p);
+                        fprintf(ANGLES, "%.8g %.8g %.8g %.8g\n", A.x, A.y, A.z, value);
+                    }
+
+                    pmmc_CurveCurvature(SDn, SDs, SDn_x, SDn_y, SDn_z, SDs_x,
+                                        SDs_y, SDs_z, KNwns_values,
+                                        KGwns_values, KNwns, KGwns, nws_pts,
+                                        n_nws_pts, i, j, k);
+
+                    lwns +=
+                        pmmc_CubeCurveLength(local_nws_pts, n_local_nws_pts);
+                    
+                    /* half contribution for vertices / edges at the common line 
+                     *   each cube with contact line has a net of undercounting vertices 
+                     *   each cube is undercounting edges due to internal counts
+                    */
+                    Xwn += 0.25*n_local_nws_pts - 0.5;
+                    Xws += 0.25*n_local_nws_pts - 0.5;
+                    Xns += 0.25*n_local_nws_pts - 0.5;
+                }
+
+                // Solid interface averagees
+                if (n_local_sol_tris > 0) {
+                    As += pmmc_CubeSurfaceArea(local_sol_pts, local_sol_tris,
+                                               n_local_sol_tris);
+
+                    // Compute the surface orientation and the interfacial area
+                    ans += pmmc_CubeSurfaceOrientation(Gns, ns_pts, ns_tris,
+                                                       n_ns_tris);
+                    aws += pmmc_CubeSurfaceOrientation(Gws, ws_pts, ws_tris,
+                                                       n_ws_tris);
+                    
+                    Xws += geomavg_EulerCharacteristic(ws_pts, ws_tris, n_ws_pts,
+                                                         n_ws_tris, i, j, k);
+                    
+                    Xns += geomavg_EulerCharacteristic(ns_pts, ns_tris, n_ns_pts,
+                                                         n_ns_tris, i, j, k);
+                }
+                //...........................................................................
+                // Compute the integral curvature of the non-wetting phase
+
+                n_nw_pts = n_nw_tris = 0;
+                // Compute the non-wetting phase surface and associated area
+                An +=
+                    geomavg_MarchingCubes(SDn, fluid_isovalue, i, j, k, nw_pts,
+                                          n_nw_pts, nw_tris, n_nw_tris);
+                // Compute the integral of mean curvature
+                if (n_nw_pts > 0) {
+                    pmmc_CubeTrimSurfaceInterpValues(
+                        CubeValues, MeanCurvature, SDs, nw_pts, nw_tris, Values,
+                        DistanceValues, i, j, k, n_nw_pts, n_nw_tris, trimdist,
+                        trawn, dummy);
+                }
+
+                Jn += pmmc_CubeSurfaceInterpValue(CubeValues, MeanCurvature,
+                                                  nw_pts, nw_tris, Values, i, j,
+                                                  k, n_nw_pts, n_nw_tris);
+                // Compute Euler characteristic from integral of gaussian curvature
+                Kn += pmmc_CubeSurfaceInterpValue(CubeValues, GaussCurvature,
+                                                  nw_pts, nw_tris, Values, i, j,
+                                                  k, n_nw_pts, n_nw_tris);
+                
+
+                euler += geomavg_EulerCharacteristic(nw_pts, nw_tris, n_nw_pts,
+                                                     n_nw_tris, i, j, k);
+
+            }
+        }
+    }
+    fclose(ANGLES);
+
+    Array<char> phase_label(Nx, Ny, Nz);
+    Array<double> phase_distance(Nx, Ny, Nz);
+    // Analyze the wetting fluid
+    for (k = 0; k < Nz; k++) {
+        for (j = 0; j < Ny; j++) {
+            for (i = 0; i < Nx; i++) {
+                n = k * Nx * Ny + j * Nx + i;
+                if (!(Dm->id[n] > 0)) {
+                    // Solid phase
+                    phase_label(i, j, k) = 1;
+                } else if (SDn(i, j, k) < 0.0) {
+                    // wetting phase
+                    phase_label(i, j, k) = 0;
+                } else {
+                    // non-wetting phase
+                    phase_label(i, j, k) = 1;
+                }
+                phase_distance(i, j, k) =
+                    2.0 * double(phase_label(i, j, k)) - 1.0;
+            }
+        }
+    }
+    CalcDist(phase_distance, phase_label, *Dm);
+    wet_morph->ComputeScalar(phase_distance, 0.f);
+    //printf("generating distance at rank=%i \n",Dm->rank());
+    // Analyze the wetting fluid
+    for (k = 0; k < Nz; k++) {
+        for (j = 0; j < Ny; j++) {
+            for (i = 0; i < Nx; i++) {
+                n = k * Nx * Ny + j * Nx + i;
+                if (!(Dm->id[n] > 0)) {
+                    // Solid phase
+                    phase_label(i, j, k) = 1;
+                } else if (SDn(i, j, k) < 0.0) {
+                    // wetting phase
+                    phase_label(i, j, k) = 1;
+                } else {
+                    // non-wetting phase
+                    phase_label(i, j, k) = 0;
+                }
+                phase_distance(i, j, k) =
+                    2.0 * double(phase_label(i, j, k)) - 1.0;
+            }
+        }
+    }
+    CalcDist(phase_distance, phase_label, *Dm);
+    nonwet_morph->ComputeScalar(phase_distance, 0.f);
+}
 
 void TwoPhase::AssignComponentLabels() {
     //int LabelNWP=1;
@@ -1219,7 +1428,7 @@ void TwoPhase::ComponentAverages() {
 
 void TwoPhase::Reduce() {
     int i;
-    double iVol_global = 1.0 / Volume;
+    //double iVol_global = 1.0 / Volume;
     //...........................................................................
     Dm->Comm.barrier();
     nwp_volume_global = Dm->Comm.sumReduce(nwp_volume);
@@ -1258,10 +1467,14 @@ void TwoPhase::Reduce() {
     trawn_global = Dm->Comm.sumReduce(trawn);
     trJwn_global = Dm->Comm.sumReduce(trJwn);
     trRwn_global = Dm->Comm.sumReduce(trRwn);
-    euler_global = Dm->Comm.sumReduce(euler);
+    Xwn_global = Dm->Comm.sumReduce(Xwn);
+    Xws_global = Dm->Comm.sumReduce(Xws);
+    Xns_global = Dm->Comm.sumReduce(Xns);
     An_global = Dm->Comm.sumReduce(An);
     Jn_global = Dm->Comm.sumReduce(Jn);
     Kn_global = Dm->Comm.sumReduce(Kn);
+    euler_global = Dm->Comm.sumReduce(euler);
+
     Dm->Comm.barrier();
 
     // Normalize the phase averages
@@ -1285,17 +1498,18 @@ void TwoPhase::Reduce() {
     // Normalize surface averages by the interfacial area
     if (awn_global > 0.0) {
         Jwn_global /= awn_global;
-        Kwn_global /= awn_global;
+        //Kwn_global /= awn_global;
         wwndnw_global /= awn_global;
         for (i = 0; i < 3; i++)
             vawn_global(i) /= awn_global;
         for (i = 0; i < 6; i++)
             Gwn_global(i) /= awn_global;
     }
+
     if (lwns_global > 0.0) {
         efawns_global /= lwns_global;
-        KNwns_global /= lwns_global;
-        KGwns_global /= lwns_global;
+        //KNwns_global /= lwns_global;
+        //KGwns_global /= lwns_global;
         for (i = 0; i < 3; i++)
             vawns_global(i) /= lwns_global;
     }
@@ -1314,15 +1528,17 @@ void TwoPhase::Reduce() {
             Gws_global(i) /= aws_global;
 
     euler_global /= (2 * PI);
-
-    //sat_w = 1.0 - nwp_volume_global*iVol_global/porosity;
     sat_w = 1.0 - nwp_volume_global / (nwp_volume_global + wp_volume_global);
+
     // Compute the specific interfacial areas and common line length (dimensionless per unit volume)
+    /* 
     awn_global = awn_global * iVol_global;
     ans_global = ans_global * iVol_global;
     aws_global = aws_global * iVol_global;
     dEs = dEs * iVol_global;
     lwns_global = lwns_global * iVol_global;
+    */
+
 }
 
 void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) {
@@ -1334,6 +1550,55 @@ void TwoPhase::NonDimensionalize(double D, double viscosity, double IFT) {
     lwns_global *= D * D;
 }
 
+void TwoPhase::PrintStatic() {
+    if (Dm->rank() == 0) {
+    	FILE *STATIC;
+        STATIC = fopen("geometry.csv", "a+");
+        if (fseek(STATIC, 0, SEEK_SET) == fseek(STATIC, 0, SEEK_CUR)) {
+            // If timelog is empty, write a short header to list the averages
+            fprintf(STATIC, "sw awn ans aws Jwn Kwn lwns cwns KGws "
+                             "KGwn Xwn Xws Xns "); // Scalar averages
+            fprintf(STATIC,
+                "Gwnxx Gwnyy Gwnzz Gwnxy Gwnxz Gwnyz "); // Orientation tensors
+            fprintf(STATIC, "Gwsxx Gwsyy Gwszz Gwsxy Gwsxz Gwsyz ");
+            fprintf(STATIC, "Gnsxx Gnsyy Gnszz Gnsxy Gnsxz Gnsyz ");
+            fprintf(STATIC, "trawn trJwn trRwn "); //trimmed curvature,
+            fprintf(STATIC, "Vw Aw Jw Xw ");      //miknowski measures,
+            fprintf(STATIC, "Vn An Jn Xn\n");     //miknowski measures,
+            //fprintf(STATIC,"Euler Kn2 Jn2 An2\n"); 			//miknowski measures,
+        }
+
+        fprintf(STATIC, "%.5g ", sat_w); // saturation 
+        fprintf(STATIC, "%.5g %.5g %.5g ", awn_global, ans_global,
+                aws_global); // interfacial areas
+        fprintf(STATIC, "%.5g %.5g ", Jwn_global,
+                Kwn_global);                      // curvature of wn interface
+        fprintf(STATIC, "%.5g ", lwns_global);   // common curve length
+        fprintf(STATIC, "%.5g ", efawns_global); // average contact angle
+        fprintf(STATIC, "%.5g %.5g ", KNwns_global,
+                KGwns_global); // total curvature contributions of common line
+        fprintf(STATIC, "%.5g %.5g %.5g ", Xwn_global, Xns_global,
+                Xws_global); // Euler characteristic
+        fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gwn_global(0),
+                Gwn_global(1), Gwn_global(2), Gwn_global(3), Gwn_global(4),
+                Gwn_global(5)); // orientation of wn interface
+        fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gns_global(0),
+                Gns_global(1), Gns_global(2), Gns_global(3), Gns_global(4),
+                Gns_global(5)); // orientation of ns interface
+        fprintf(STATIC, "%.5g %.5g %.5g %.5g %.5g %.5g ", Gws_global(0),
+                Gws_global(1), Gws_global(2), Gws_global(3), Gws_global(4),
+                Gws_global(5)); // orientation of ws interface
+        fprintf(STATIC, "%.5g %.5g %.5g ", trawn_global, trJwn_global,
+                trRwn_global); // Trimmed curvature
+        fprintf(STATIC, "%.5g %.5g %.5g %.5g ", wet_morph->V(), wet_morph->A(),
+                wet_morph->H(), wet_morph->X());
+       fprintf(STATIC, "%.5g %.5g %.5g %.5g\n", nonwet_morph->V(),
+                nonwet_morph->A(), nonwet_morph->H(), nonwet_morph->X());
+        //fprintf(STATIC,"%.5g %.5g %.5g %.5g\n",euler_global, Kn_global, Jn_global, An_global);			// minkowski measures
+        fclose(STATIC);
+    }
+}
+
 void TwoPhase::PrintAll(int timestep) {
     if (Dm->rank() == 0) {
         fprintf(TIMELOG, "%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g ",

analysis/TwoPhase.h

@@ -103,7 +103,9 @@ public:
     double lwns_global;
     double efawns, efawns_global; // averaged contact angle
     double euler, Kn, Jn, An;
+    double Xwn, Xns, Xws;
     double euler_global, Kn_global, Jn_global, An_global;
+    double Xwn_global, Xns_global, Xws_global;
 
     double rho_n, rho_w;
     double nu_n, nu_w;
@@ -184,6 +186,8 @@ public:
     void ColorToSignedDistance(double Beta, DoubleArray &ColorData,
                                DoubleArray &DistData);
     void ComputeLocal();
+    void ComputeStatic();
+    void PrintStatic();
     void AssignComponentLabels();
     void ComponentAverages();
     void Reduce();

analysis/morphology.cpp

@@ -137,8 +137,7 @@ void Morphology::Initialize(std::shared_ptr<Domain> Dm, DoubleArray &Distance) {
     morphRadius.resize(recvLoc);
     //..............................
     /* send the morphological radius */
-    Dm->Comm.Irecv(&morphRadius[recvOffset_x], recvCount, Dm->rank_x(),
-                   recvtag + 0);
+    Dm->Comm.Irecv(&morphRadius[recvOffset_x], recvCount, Dm->rank_x(), recvtag + 0);
     Dm->Comm.send(&tmpDistance[0], sendCount, Dm->rank_X(), sendtag + 0);
     /* send the shift values */
     Dm->Comm.Irecv(&xShift[recvOffset_x], recvCount, Dm->rank_x(), recvtag + 1);
@@ -502,7 +501,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
     if (rank == 0)
         printf("Maximum pore size: %f \n", maxdistGlobal);
     final_void_fraction = volume_fraction; //initialize
-
+    
     int ii, jj, kk;
     int imin, jmin, kmin, imax, jmax, kmax;
     int Nx = nx;
@@ -524,27 +523,28 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
 
     int numTry = 0;
     int maxTry = 100;
-    while (void_fraction_new > VoidFraction && numTry < maxTry) {
+    while ( !(void_fraction_new < VoidFraction) && numTry < maxTry) {
         numTry++;
         void_fraction_diff_old = void_fraction_diff_new;
         void_fraction_old = void_fraction_new;
         Rcrit_old = Rcrit_new;
         Rcrit_new -= deltaR * Rcrit_old;
+        if (rank==0) printf("Try %i with radius %f \n", numTry, Rcrit_new);
         if (Rcrit_new < 0.5) {
             numTry = maxTry;
         }
         int Window = round(Rcrit_new);
         if (Window == 0)
-            Window =
-                1; // If Window = 0 at the begining, after the following process will have sw=1.0
+            Window =  1; // If Window = 0 at the begining, after the following process will have sw=1.0
         // and sw<Sw will be immediately broken
         double LocalNumber = 0.f;
-        for (int k = 0; k < Nz; k++) {
-            for (int j = 0; j < Ny; j++) {
-                for (int i = 0; i < Nx; i++) {
+        for (int k = 1; k < Nz-1; k++) {
+            for (int j = 1; j < Ny-1; j++) {
+                for (int i = 1; i < Nx-1; i++) {
                     n = k * nx * ny + j * nx + i;
                     if (SignDist(i, j, k) > Rcrit_new) {
                         // loop over the window and update
+                    	//printf("Distance(%i %i %i) = %f \n",i,j,k, SignDist(i,j,k));
                         imin = max(1, i - Window);
                         jmin = max(1, j - Window);
                         kmin = max(1, k - Window);
@@ -571,7 +571,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
                 }
             }
         }
-        LocalNumber += Structure.GetOverlaps(Dm, id, ErodeLabel, NewLabel);
+        //LocalNumber += Structure.GetOverlaps(Dm, id, ErodeLabel, NewLabel);
 
         count = 0.f;
         for (int k = 1; k < Nz - 1; k++) {
@@ -611,7 +611,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
 
 //***************************************************************************************
 double MorphDrain(DoubleArray &SignDist, signed char *id,
-                  std::shared_ptr<Domain> Dm, double VoidFraction) {
+                  std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius) {
     // SignDist is the distance to the object that you want to constaing the morphological opening
     // VoidFraction is the the empty space where the object inst
     // id is a labeled map
@@ -688,6 +688,11 @@ double MorphDrain(DoubleArray &SignDist, signed char *id,
     double deltaR = 0.05; // amount to change the radius in voxel units
     double Rcrit_old = maxdistGlobal;
     double Rcrit_new = maxdistGlobal;
+    
+    if (InitialRadius < maxdistGlobal){
+    	Rcrit_old = InitialRadius;
+    	Rcrit_new = InitialRadius;
+    }
     //if (argc>2){
     //	Rcrit_new = strtod(argv[2],NULL);
     //	if (rank==0) printf("Max. distance =%f, Initial critical radius = %f \n",maxdistGlobal,Rcrit_new);

analysis/morphology.h

@@ -7,7 +7,7 @@ double MorphOpen(DoubleArray &SignDist, signed char *id,
                  std::shared_ptr<Domain> Dm, double VoidFraction,
                  signed char ErodeLabel, signed char ReplaceLabel);
 double MorphDrain(DoubleArray &SignDist, signed char *id,
-                  std::shared_ptr<Domain> Dm, double VoidFraction);
+                  std::shared_ptr<Domain> Dm, double VoidFraction, double InitialRadius);
 double MorphGrow(DoubleArray &BoundaryDist, DoubleArray &Dist, Array<char> &id,
                  std::shared_ptr<Domain> Dm, double TargetVol,
                  double WallFactor);

analysis/pmmc.h

@@ -4057,7 +4057,7 @@ inline double pmmc_CubeContactAngle(DoubleArray &CubeValues,
                       (A.z - B.z) * (A.z - B.z));
         integral += 0.5 * length * (vA + vB);
     }
-
+    
     return integral;
 }
 //--------------------------------------------------------------------------------------------------------
@@ -4420,7 +4420,9 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
     double twnsx, twnsy, twnsz, nwnsx, nwnsy, nwnsz,
         K; // tangent,normal and curvature
     double nsx, nsy, nsz, norm;
+    double nwx, nwy, nwz;
     double nwsx, nwsy, nwsz;
+    double nwnx, nwny, nwnz;
 
     Point P, A, B;
     // Local trilinear approximation for tangent and normal vector
@@ -4430,47 +4432,18 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
     for (k = kc; k < kc + 2; k++) {
         for (j = jc; j < jc + 2; j++) {
             for (i = ic; i < ic + 2; i++) {
-
-                // Compute all of the derivatives using finite differences
-                // fluid phase indicator field
-                //	fx = 0.5*(f(i+1,j,k) - f(i-1,j,k));
-                //    fy = 0.5*(f(i,j+1,k) - f(i,j-1,k));
-                //	fz = 0.5*(f(i,j,k+1) - f(i,j,k-1));
-                /*fxx = f(i+1,j,k) - 2.0*f(i,j,k) + f(i-1,j,k);
-				fyy = f(i,j+1,k) - 2.0*f(i,j,k) + f(i,j-1,k);
-				fzz = f(i,j,k+1) - 2.0*f(i,j,k) + f(i,j,k-1);
-				fxy = 0.25*(f(i+1,j+1,k) - f(i+1,j-1,k) - f(i-1,j+1,k) + f(i-1,j-1,k));
-				fxz = 0.25*(f(i+1,j,k+1) - f(i+1,j,k-1) - f(i-1,j,k+1) + f(i-1,j,k-1));
-				fyz = 0.25*(f(i,j+1,k+1) - f(i,j+1,k-1) - f(i,j-1,k+1) + f(i,j-1,k-1));
-				*/
-                // solid distance function
-                //	sx = 0.5*(s(i+1,j,k) - s(i-1,j,k));
-                //	sy = 0.5*(s(i,j+1,k) - s(i,j-1,k));
-                //	sz = 0.5*(s(i,j,k+1) - s(i,j,k-1));
-                /*	sxx = s(i+1,j,k) - 2.0*s(i,j,k) + s(i-1,j,k);
-				syy = s(i,j+1,k) - 2.0*s(i,j,k) + s(i,j-1,k);
-				szz = s(i,j,k+1) - 2.0*s(i,j,k) + s(i,j,k-1);
-				sxy = 0.25*(s(i+1,j+1,k) - s(i+1,j-1,k) - s(i-1,j+1,k) + s(i-1,j-1,k));
-				sxz = 0.25*(s(i+1,j,k+1) - s(i+1,j,k-1) - s(i-1,j,k+1) + s(i-1,j,k-1));
-				syz = 0.25*(s(i,j+1,k+1) - s(i,j+1,k-1) - s(i,j-1,k+1) + s(i,j-1,k-1));
-				*/
-                /*				// Compute the Jacobean matrix for tangent vector
-				Axx = sxy*fz + sy*fxz - sxz*fy - sz*fxy;
-				Axy = sxz*fx + sz*fxx - sxx*fz - sx*fxz;
-				Axz = sxx*fy + sx*fxy - sxy*fx - sy*fxx;
-				Ayx = syy*fz + sy*fyz - syz*fy - sz*fyy;
-				Ayy = syz*fx + sz*fxy - sxy*fz - sx*fyz;
-				Ayz = sxy*fy + sx*fyy - syy*fx - sy*fxy;
-				Azx = syz*fz + sy*fzz - szz*fy - sz*fyz;
-				Azy = szz*fx + sz*fxz - sxz*fz - sx*fzz;
-				Azz = sxz*fy + sx*fyz - syz*fx - sy*fxz;
-*/
                 sx = s_x(i, j, k);
                 sy = s_y(i, j, k);
                 sz = s_z(i, j, k);
                 fx = f_x(i, j, k);
                 fy = f_y(i, j, k);
                 fz = f_z(i, j, k);
+                
+                // Normal to fluid surface
+                Nx.Corners(i - ic, j - jc, k - kc) = fx;
+                Ny.Corners(i - ic, j - jc, k - kc) = fy;
+                Nz.Corners(i - ic, j - jc, k - kc) = fz;
+                
                 // Normal to solid surface
                 Sx.Corners(i - ic, j - jc, k - kc) = sx;
                 Sy.Corners(i - ic, j - jc, k - kc) = sy;
@@ -4577,8 +4550,19 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
         nsx /= norm;
         nsy /= norm;
         nsz /= norm;
+        
+        // Normal vector to the fluid surface
+        nwx = Nx.eval(P);
+        nwy = Ny.eval(P);
+        nwz = Nz.eval(P);
+        norm = sqrt(nwx * nwx + nwy * nwy + nwz * nwz);
+        if (norm == 0.0)
+            norm = 1.0;
+        nwx /= norm;
+        nwy /= norm;
+        nwz /= norm;
 
-        // normal in the surface tangent plane (rel. geodesic curvature)
+        // common curve normal in the solid surface tangent plane (rel. geodesic curvature)
         nwsx = twnsy * nsz - twnsz * nsy;
         nwsy = twnsz * nsx - twnsx * nsz;
         nwsz = twnsx * nsy - twnsy * nsx;
@@ -4588,16 +4572,35 @@ inline void pmmc_CurveCurvature(DoubleArray &f, DoubleArray &s,
         nwsx /= norm;
         nwsy /= norm;
         nwsz /= norm;
-
-        if (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz < 0.0) {
-            nwnsx = -nwnsx;
-            nwnsy = -nwnsy;
-            nwnsz = -nwnsz;
+        /* normal to ws interface boundary should point into fluid (same direction as gradient) */
+        if (nwx * nwsx + nwy * nwsy + nwz * nwsz < 0.0) {
+            nwsx = -nwsx;
+            nwsy = -nwsy;
+            nwsz = -nwsz;
         }
+        
+        // common curve normal in the fluid surface tangent plane (rel. geodesic curvature)
+        nwnx = twnsy * nwz - twnsz * nwy;
+        nwny = twnsz * nwx - twnsx * nwz;
+        nwnz = twnsx * nwy - twnsy * nwx;
+        norm = sqrt(nwnx * nwnx + nwny * nwny + nwnz * nwnz);
+        if (norm == 0.0)
+            norm = 1.0;
+        nwnx /= norm;
+        nwny /= norm;
+        nwnz /= norm;
+        /* normal to wn interface boundary should point into the solid */
+        if (nsx * nwnx + nsy * nwny + nsz * nwnz > 0.0) {
+            nwnx = -nwnx;
+            nwny = -nwny;
+            nwnz = -nwnz;
+        }
+
 
         if (length > 0.0) {
             // normal curvature component in the direction of the solid surface
-            KNavg += K * (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz) * length;
+            //KNavg += K * (nsx * nwnsx + nsy * nwnsy + nsz * nwnsz) * length;
+            KNavg += K * (nwnx * nwnsx + nwny * nwnsy + nwnz * nwnsz) * length;
             //geodesic curvature
             KGavg += K * (nwsx * nwnsx + nwsy * nwnsy + nwsz * nwnsz) * length;
         }

cpu/D3Q7BC.cpp

@@ -35,13 +35,31 @@ extern "C" void ScaLBL_Solid_Neumann_D3Q7(double *dist, double *BoundaryValue,
     }
 }
 
-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, double *lattice_weight,
-    float *lattice_cx, float *lattice_cy, float *lattice_cz, int count,
-    int Np) {
+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; idx<N; idx++){
+		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;
+        }
+	}
+}
+
+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,
+                                                      double *lattice_weight, float *lattice_cx, float *lattice_cy, float *lattice_cz,
+                                                      int count, int Np){
 
     int idx;
     int iq, ib, ifluidBC;

cpu/Poisson.cpp

@@ -95,7 +95,7 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(
 extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
                                           double *dist, double *Den_charge,
                                           double *Psi, double *ElectricField,
-                                          double tau, double epsilon_LB,
+                                          double tau, double epsilon_LB, bool UseSlippingVelBC,
                                           int start, int finish, int Np) {
     int n;
     double psi;        //electric potential
@@ -109,8 +109,9 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
     for (n = start; n < finish; n++) {
 
         //Load data
-        rho_e = Den_charge[n];
-        rho_e = rho_e / epsilon_LB;
+        //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+        //and thus the net space charge density is zero. 
+        rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
         idx = Map[n];
         psi = Psi[idx];
 
@@ -175,8 +176,8 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map,
 extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
                                            double *Den_charge, double *Psi,
                                            double *ElectricField, double tau,
-                                           double epsilon_LB, int start,
-                                           int finish, int Np) {
+                                           double epsilon_LB, bool UseSlippingVelBC,
+                                           int start, int finish, int Np) {
     int n;
     double psi;        //electric potential
     double Ex, Ey, Ez; //electric field
@@ -188,8 +189,9 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist,
     for (n = start; n < finish; n++) {
 
         //Load data
-        rho_e = Den_charge[n];
-        rho_e = rho_e / epsilon_LB;
+        //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+        //and thus the net space charge density is zero. 
+        rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
         idx = Map[n];
         psi = Psi[idx];
 

cpu/Stokes.cpp

@@ -4,7 +4,7 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
     double *dist, double *Velocity, double *ChargeDensity,
     double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
     double Gy, double Gz, double rho0, double den_scale, double h,
-    double time_conv, int start, int finish, int Np) {
+    double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
     double fq;
     // conserved momemnts
     double rho, jx, jy, jz;
@@ -38,13 +38,11 @@ extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(
         Ey = ElectricField[n + 1 * Np];
         Ez = ElectricField[n + 2 * Np];
         //compute total body force, including input body force (Gx,Gy,Gz)
-        Fx =
-            Gx +
-            rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                den_scale; //the extra factors at the end necessarily convert unit from phys to LB
-        Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                      den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+        Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
                       den_scale;
-        Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
                       den_scale;
 
         // q=0
@@ -479,7 +477,7 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
     int *neighborList, double *dist, double *Velocity, double *ChargeDensity,
     double *ElectricField, double rlx_setA, double rlx_setB, double Gx,
     double Gy, double Gz, double rho0, double den_scale, double h,
-    double time_conv, int start, int finish, int Np) {
+    double time_conv, bool UseSlippingVelBC, int start, int finish, int Np) {
     double fq;
     // conserved momemnts
     double rho, jx, jy, jz;
@@ -513,12 +511,21 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(
         Ex = ElectricField[n + 0 * Np];
         Ey = ElectricField[n + 1 * Np];
         Ez = ElectricField[n + 2 * Np];
+
         //compute total body force, including input body force (Gx,Gy,Gz)
-        Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        //Fx = Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        //              den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+        //Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        //              den_scale;
+        //Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        //              den_scale;
+        //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+        //and body force induced by external efectric field is reduced to slipping velocity BC. 
+        Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                      den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+        Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
                       den_scale;
-        Fy = Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
-                      den_scale;
-        Fz = Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+        Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
                       den_scale;
 
         // q=0

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<<<GRID,512>>>(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,

cuda/Poisson.cu

@@ -104,7 +104,7 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	int n;
 	double psi;//electric potential
@@ -122,8 +122,9 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
 		if (n<finish) {
 
             //Load data
-            rho_e = Den_charge[n];
-            rho_e = rho_e/epsilon_LB;
+            //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+            //and thus the net space charge density is zero. 
+            rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
             idx=Map[n];
             psi = Psi[idx];
             
@@ -184,7 +185,7 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	int n;
 	double psi;//electric potential
@@ -201,8 +202,9 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
 		if (n<finish) {
 
             //Load data
-            rho_e = Den_charge[n];
-            rho_e = rho_e/epsilon_LB;
+            //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+            //and thus the net space charge density is zero. 
+            rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
             idx=Map[n];
             psi = Psi[idx];
 
@@ -293,10 +295,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
+	dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
 
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
@@ -305,10 +307,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
+	dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
 
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){

cuda/Stokes.cu

@@ -5,7 +5,7 @@
 #define NBLOCKS 1024
 #define NTHREADS 256
 
-__global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
 
     int n;
     double fq;
@@ -46,9 +46,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
             Ey = ElectricField[n+1*Np];
             Ez = ElectricField[n+2*Np];
             //compute total body force, including input body force (Gx,Gy,Gz)
-            Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+            Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+            Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
+            Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
 
             // q=0
             fq = dist[n];
@@ -510,7 +513,7 @@ __global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
 
     int n;
     double fq;
@@ -550,9 +553,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
             Ey = ElectricField[n+1*Np];
             Ez = ElectricField[n+2*Np];
             //compute total body force, including input body force (Gx,Gy,Gz)
-            Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
-            Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+            Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+            Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
+            Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
 
             // q=0
             fq = dist[n];
@@ -969,10 +975,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
                 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
+	dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np);
 
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){
@@ -981,10 +987,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
                 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
+	dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC,start,finish,Np);
 
 	cudaError_t err = cudaGetLastError();
 	if (cudaSuccess != err){

docs/source/examples/analysis/TwoPhase.rst

@@ -0,0 +1,137 @@
+*************************
+Measuring Contact Angles
+*************************
+
+LBPM includes specialized data analysis capabilities for two-fluid systems. While these
+components are generally designed for in situ analysis of simulation data, they can also
+be applied independently to analyze 3D image data. In this example we consider applying
+the analysis tools implemented in ``lbpm_TwoPhase_analysis``, which are designed to
+analyze two-fluid configurations in porous media. The numerical implementation used to
+construct the common line are described in ( https://doi.org/10.1016/j.advwatres.2006.06.010 ).
+Methods used to measure the contact angle are described in ( https://doi.org/10.1017/jfm.2016.212 ).
+
+Source files for the example are included in the LBPM repository
+in the directory ``examples/Droplet``. A simple python code is included
+to set up a fluid droplet on a flat surface
+
+.. code:: python
+	  
+	  import numpy as np
+	  import matplotlib.pylab as plt
+
+	  D=np.ones((80,80,40),dtype="uint8")
+
+	  cx = 40
+	  cy = 40
+	  cz = 0
+
+	  for i in range(0,80):
+	      for j in range (0, 80):
+		  for k in range (0,40):
+		      dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
+		      if (dist < 25) :
+			  D[i,j,k] = 2
+		      if k < 4 : 
+			  D[i,j,k] = 0
+
+	  D.tofile("droplet_40x80x80.raw")
+
+
+The input file provided below will specify how the analysis should be performed. The name and the dimensions of the input
+file are provided in the ``Domain`` section, as with other LBPM simulators. For large images, additional processors can be
+used to speed up the analysis or take advantage of distributed memory. The ``ReadValues`` list should be used to specify the
+labels to use for analysis. The first label will be taken to be the solid.  The second label will be taken to be the fluid
+to analyze, which in this case will be the droplet labeled with ``2`` above. 
+	  
+.. code:: bash
+	  
+       Domain {
+	  Filename = "droplet_40x80x80.raw"
+	  nproc = 1, 1, 1     // Number of processors (Npx,Npy,Npz)
+	  n = 40, 80, 80      // Size of local domain (Nx,Ny,Nz)
+	  N = 40, 80, 80      // size of the input image
+	  voxel_length = 1.0 
+	  BC = 0              // Boundary condition type
+	  ReadType = "8bit"
+	  ReadValues = 0, 2, 1
+	  WriteValues = 0, 2, 1
+      }
+      Visualization {
+      }
+
+
+The analysis can be launched as ``mpirun -np 1 $LBPM_DIR/lbpm_TwoPhase_analysis input.db``. Output should appear as
+follows:
+
+.. code:: bash
+ 
+	   Input data file: input.db
+	   voxel length = 1.000000 micron 
+	   Input media: droplet_40x80x80.raw
+	   Relabeling 3 values
+	   oldvalue=0, newvalue =0 
+	   oldvalue=2, newvalue =2 
+	   oldvalue=1, newvalue =1 
+	   Dimensions of segmented image: 40 x 80 x 80 
+	   Reading 8-bit input data 
+	   Read segmented data from droplet_40x80x80.raw 
+	   Label=0, Count=25600 
+	   Label=2, Count=25773 
+	   Label=1, Count=204627 
+	   Distributing subdomains across 1 processors 
+	   Process grid: 1 x 1 x 1 
+	   Subdomain size: 40 x 80 x 80 
+	   Size of transition region: 0 
+	   Media porosity = 0.900000 
+	   Initialized solid phase -- Converting to Signed Distance function 
+	   Initialized fluid phase -- Converting to Signed Distance function 
+	   Computing Minkowski functionals 
+
+The ``TwoPhase`` analysis class will generate signed distance functions for the solid and fluid surfaces.
+Using the distance functions, the interfaces and common line will be constructed. Contact angles are logged
+for each processor, e.g. ``ContactAngle.00000.csv``, which specifies the x, y, z coordinates for each measurement
+along with the cosine of the contact angle. Averaged measures (determined over the entire input image)
+are logged to ``geometry.csv``
+
+* ``sw`` -- water saturation
+* ``awn`` -- surface area of meniscus between wn fluids
+* ``ans`` -- surface area between fluid n and solid
+* ``aws`` -- surface area between fluid w and solid
+* ``Jwn`` -- integral of mean curvature of meniscus
+* ``Kwn`` -- integral of Gaussian curvature of meniscus
+* ``lwns`` -- length of common line
+* ``cwns`` -- average contact angle
+* ``KGws`` -- geodesic curvature of common line relative to ws surface
+* ``KGwn`` -- geodesic curvature of common line relative to wn surface
+* ``Gwnxx`` -- orientation tensor component for wn surface
+* ``Gwnyy`` -- orientation tensor component for wn surface
+* ``Gwnzz`` -- orientation tensor component for wn surface
+* ``Gwnxy`` -- orientation tensor component for wn surface
+* ``Gwnxz`` -- orientation tensor component for wn surface
+* ``Gwnyz`` -- orientation tensor component for wn surface
+* ``Gwsxx`` -- orientation tensor component for ws surface
+* ``Gwsyy`` -- orientation tensor component for ws surface
+* ``Gwszz`` -- orientation tensor component for ws surface
+* ``Gwsxy`` -- orientation tensor component for ws surface
+* ``Gwsxz`` -- orientation tensor component for ws surface
+* ``Gwsyz`` -- orientation tensor component for ws surface
+* ``Gnsxx`` -- orientation tensor component for ns surface
+* ``Gnsyy`` -- orientation tensor component for ns surface
+* ``Gnszz`` -- orientation tensor component for ns surface
+* ``Gnsxy`` -- orientation tensor component for ns surface
+* ``Gnsxz`` -- orientation tensor component for ns surface
+* ``Gnsyz`` -- orientation tensor component for ns surface
+* ``trawn`` -- trimmed surface area for meniscus (one voxel from solid)
+* ``trJwn`` -- mean curvature for trimmed meniscus
+* ``trRwn`` -- radius of curvature for trimmed meniscus
+* ``Vw`` -- volume of fluid w
+* ``Aw`` -- boundary surface area for fluid w
+* ``Jw`` -- integral of mean curvature for fluid w
+* ``Xw`` -- Euler characteristic for fluid w
+* ``Vn`` -- volume of fluid n
+* ``An``  -- boundary surface area for fluid n
+* ``Jn`` -- integral of mean curvature for fluid n
+* ``Xn`` -- Euler characteristic for fluid n
+
+
+  

docs/source/examples/color/fractionalFlow.rst

@@ -0,0 +1,125 @@
+********************************
+Steady-state fractional flow 
+********************************
+
+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 {
+	     protocol = "fractional flow" 
+	     capillary_number = -1e-5           // capillary number for the displacement, positive="oil injection"
+	     timestepMax = 200000               // 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 {
+	     max_steady_timesteps = 25000       // maximum number of timesteps per steady point
+	     min_steady_timesteps = 25000       // minimum number of timesteps per steady point
+	     fractional_flow_increment = 0.0003 // parameter that controls rate of mass seeding
+	     skip_timesteps = 10000             // number of timesteps to spend in flow adaptor
+	  }
+
+	  
+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)  
+

docs/source/examples/color/steadyState.rst

@@ -1,5 +1,126 @@
-*******************
-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 {
+	   protocol = "fractional flow"
+	   capillary_number = 1e-5            // capillary number for the displacement, positive="oil injection"
+	   timestepMax = 500000               // 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, 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 = 1000000    // 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 {
+	   min_steady_timesteps = 250000       // minimum number of timesteps per steady point
+	   max_steady_timesteps = 300000       // maximum number of timesteps per steady point
+	   fractional_flow_increment = 0.1   // parameter that controls rate of mass seeding
+	   skip_timesteps = 10000             // number of timesteps to spend in flow adaptor
+	   endpoint_threshold = 0.1           // endpoint exit criterion
+	}
+
+	  
+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.

docs/source/examples/running.rst

@@ -18,4 +18,5 @@ a basic introduction to working with LBPM.
    morphology/*	      
    
    color/*
-   
+
+   analysis/*

docs/source/publications/publications.rst

@@ -1,13 +1,15 @@
-============
+\============
 Publications
 ============
 
-* James E McClure, Zhe Li, Mark Berrill, Thomas Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9
+* J.E. McClure, Z. Li, M. Berrill, T. Ramstad, "The LBPM software package for simulating multiphase flow on digital images of porous rocks" Computational Geosciences (25) 871–895 (2021) https://doi.org/10.1007/s10596-020-10028-9
 
 
-* James E. McClure, Zhe Li, Adrian P. Sheppard, Cass T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670
+* J. E. McClure, Z. Li, A.P. Sheppard, C.T. Miller, "An adaptive volumetric flux boundary condition for lattice Boltzmann methods" Computers & Fluids (210) (2020) https://doi.org/10.1016/j.compfluid.2020.104670
 
 
 * Y.D. Wang, T. Chung, R.T. Armstrong, J. McClure, T. Ramstad, P. Mostaghimi, "Accelerated Computation of Relative Permeability by Coupled Morphological and Direct Multiphase Flow Simulation" Journal of Computational Physics (401) (2020) https://doi.org/10.1016/j.jcp.2019.108966
 
- 
+* J.E. McClure, M. Berrill, W. Gray, C.T. Miller, C.T. "Tracking interface and common curve dynamics for two-fluid flow in porous media. Journal of Fluid Mechanics, 796, 211-232 (2016) https://doi.org/10.1017/jfm.2016.212 
+
+* J.E.McClure, D.Adalsteinsson, C.Pan, W.G.Gray, C.T.Miller "Approximation of interfacial properties in multiphase porous medium systems" Advances in Water Resources, 30 (3): 354-365 (2007) https://doi.org/10.1016/j.advwatres.2006.06.010

docs/source/userGuide/models/greyscaleColor/greyscaleColor.rst

@@ -120,7 +120,7 @@ two fluids are permitted to freely mix between the endpoints. Beyond the endpoin
 term is used to drive spontaneous imbibition into the grey voxels
 
 
-..math::
+.. math::
    :nowrap:
 
       $$

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

example/CMakeLists.txt

@@ -3,6 +3,8 @@
 INSTALL_EXAMPLE( Bubble )
 INSTALL_EXAMPLE( ConstrainedBubble )
 INSTALL_EXAMPLE( Piston )
+INSTALL_EXAMPLE( Droplet )
+INSTALL_EXAMPLE( DropletCoalescence )
 INSTALL_EXAMPLE( Plates )
 INSTALL_EXAMPLE( SquareTube )
 INSTALL_EXAMPLE( InkBottle )

example/Droplet/CreateDroplet.py

@@ -0,0 +1,20 @@
+import numpy as np
+import matplotlib.pylab as plt
+
+D=np.ones((80,80,40),dtype="uint8")
+
+cx = 40
+cy = 40
+cz = 0
+
+for i in range(0,80):
+    for j in range (0, 80):
+        for k in range (0,40):
+            dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
+            if (dist < 25) :
+                D[i,j,k] = 2
+            if k < 4 : 
+                D[i,j,k] = 0
+                
+D.tofile("droplet_40x80x80.raw")
+

example/Droplet/input.db

@@ -0,0 +1,14 @@
+Domain {
+    Filename = "droplet_40x80x80.raw"
+    nproc = 1, 1, 1     // Number of processors (Npx,Npy,Npz)
+    n = 40, 80, 80      // Size of local domain (Nx,Ny,Nz)
+    N = 40, 80, 80      // size of the input image
+    voxel_length = 1.0 
+    BC = 0              // Boundary condition type
+    ReadType = "8bit"
+    ReadValues = 0, 2, 1
+    WriteValues = 0, 2, 1
+}
+
+Visualization {
+}

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<<<GRID,512>>>(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<<<GRID,512>>>(list, dist, Vin, count, Np);

hip/Poisson.cu

@@ -104,7 +104,7 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, doub
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	int n;
 	double psi;//electric potential
@@ -122,8 +122,9 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
 		if (n<finish) {
 
             //Load data
-            rho_e = Den_charge[n];
-            rho_e = rho_e/epsilon_LB;
+            //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+            //and thus the net space charge density is zero. 
+            rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
             idx=Map[n];
             psi = Psi[idx];
             
@@ -184,7 +185,7 @@ __global__  void dvc_ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, doub
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	int n;
 	double psi;//electric potential
@@ -201,8 +202,9 @@ __global__  void dvc_ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *
 		if (n<finish) {
 
             //Load data
-            rho_e = Den_charge[n];
-            rho_e = rho_e/epsilon_LB;
+            //When Helmholtz-Smoluchowski slipping velocity BC is used, the bulk fluid is considered as electroneutral
+            //and thus the net space charge density is zero. 
+            rho_e = (UseSlippingVelBC==1) ? 0.0 : Den_charge[n] / epsilon_LB;
             idx=Map[n];
             psi = Psi[idx];
 
@@ -293,10 +295,10 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
+	dvc_ScaLBL_D3Q7_AAodd_Poisson<<<NBLOCKS,NTHREADS >>>(neighborList,Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
@@ -305,10 +307,10 @@ extern "C" void ScaLBL_D3Q7_AAodd_Poisson(int *neighborList, int *Map, double *d
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q7_AAeven_Poisson(int *Map, double *dist, double *Den_charge, double *Psi, double *ElectricField, double tau, double epsilon_LB,bool UseSlippingVelBC,int start, int finish, int Np){
 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,start,finish,Np);
+	dvc_ScaLBL_D3Q7_AAeven_Poisson<<<NBLOCKS,NTHREADS >>>(Map,dist,Den_charge,Psi,ElectricField,tau,epsilon_LB,UseSlippingVelBC,start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){

hip/Stokes.cu

@@ -6,7 +6,7 @@
 #define NTHREADS 256
 
 
-__global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz, double rho0, double den_scale, double h, double time_conv,bool UseSlippingVelBC,int start, int finish, int Np){
 
     int n;
     double fq;
@@ -47,9 +47,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
             Ey = ElectricField[n+1*Np];
             Ez = ElectricField[n+2*Np];
             //compute total body force, including input body force (Gx,Gy,Gz)
-            Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+            Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+            Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
+            Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
 
             // q=0
             fq = dist[n];
@@ -511,7 +514,7 @@ __global__  void dvc_ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dis
 	}
 }
 
-__global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+__global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC,int start, int finish, int Np){
 
     int n;
     double fq;
@@ -551,9 +554,12 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
             Ey = ElectricField[n+1*Np];
             Ez = ElectricField[n+2*Np];
             //compute total body force, including input body force (Gx,Gy,Gz)
-            Fx = Gx + rhoE*Ex*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;//the extra factors at the end necessarily convert unit from phys to LB
-            Fy = Gy + rhoE*Ey*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
-            Fz = Gz + rhoE*Ez*(time_conv*time_conv)/(h*h*1.0e-12)/den_scale;
+            Fx = (UseSlippingVelBC==1) ? Gx : Gx + rhoE * Ex * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale; //the extra factors at the end necessarily convert unit from phys to LB
+            Fy = (UseSlippingVelBC==1) ? Gy : Gy + rhoE * Ey * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
+            Fz = (UseSlippingVelBC==1) ? Gz : Gz + rhoE * Ez * (time_conv * time_conv) / (h * h * 1.0e-12) /
+                          den_scale;
 
             // q=0
             fq = dist[n];
@@ -970,10 +976,10 @@ __global__  void dvc_ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocit
 	}
 }
 
-extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
                 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
+	dvc_ScaLBL_D3Q19_AAodd_StokesMRT<<<NBLOCKS,NTHREADS >>>(neighborList,dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){
@@ -982,10 +988,10 @@ extern "C" void ScaLBL_D3Q19_AAodd_StokesMRT(int *neighborList, double *dist, do
 	//cudaProfilerStop();
 }
 
-extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, int start, int finish, int Np){
+extern "C" void ScaLBL_D3Q19_AAeven_StokesMRT(double *dist, double *Velocity, double *ChargeDensity, double *ElectricField, double rlx_setA, double rlx_setB, double Gx, double Gy, double Gz,double rho0, double den_scale, double h, double time_conv, bool UseSlippingVelBC, int start, int finish, int Np){
                 
 	//cudaProfilerStart();
-	dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,start,finish,Np);
+	dvc_ScaLBL_D3Q19_AAeven_StokesMRT<<<NBLOCKS,NTHREADS >>>(dist,Velocity,ChargeDensity,ElectricField,rlx_setA,rlx_setB,Gx,Gy,Gz,rho0,den_scale,h,time_conv,UseSlippingVelBC, start,finish,Np);
 
 	hipError_t err = hipGetLastError();
 	if (hipSuccess != err){

models/ColorModel.cpp

@@ -618,6 +618,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
     bool SET_CAPILLARY_NUMBER = false;
     bool TRIGGER_FORCE_RESCALE = false;
     double tolerance = 0.01;
+    auto WettingConvention = color_db->getWithDefault<std::string>( "WettingConvention", "none" );
     auto current_db = db->cloneDatabase();
     auto flow_db = db->getDatabase("FlowAdaptor");
     int MIN_STEADY_TIMESTEPS =
@@ -645,7 +646,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
     if (analysis_db->keyExists("tolerance")) {
         tolerance = analysis_db->getScalar<double>("tolerance");
     }
-
+    
     runAnalysis analysis(current_db, rank_info, ScaLBL_Comm, Dm, Np, Regular,
                          Map);
     auto t1 = std::chrono::system_clock::now();
@@ -757,7 +758,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
             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;
@@ -875,6 +876,14 @@ double ScaLBL_ColorModel::Run(int returntime) {
                     double vBd_x = Averages->gwd.Px / Mass_w;
                     double vBd_y = Averages->gwd.Py / Mass_w;
                     double vBd_z = Averages->gwd.Pz / Mass_w;
+                    
+                    /* contribution from water films */
+                    double water_film_flow_rate =
+                    		Averages->gwb.V * (Averages->giwn.Pwx * dir_x + Averages->giwn.Pwy * dir_y + Averages->giwn.Pwz * dir_z) /
+                    		Averages->gwb.M / Dm->Volume;
+                    double not_water_film_flow_rate =
+                    		Averages->gnb.V * (Averages->giwn.Pnx * dir_x + Averages->giwn.Pny * dir_y + Averages->giwn.Pnz * dir_z) /
+                    		Averages->gnb.M / Dm->Volume;
 
                     double flow_rate_A_connected =
                         Vol_nc *
@@ -911,6 +920,15 @@ double ScaLBL_ColorModel::Run(int returntime) {
                     double kAeff = h * h * muA * (flow_rate_A) / (force_mag);
                     double kBeff = h * h * muB * (flow_rate_B) / (force_mag);
 
+                    
+                    /* not counting films */
+                    double krnf = kAeff - h * h * muA * not_water_film_flow_rate / force_mag;
+                    double krwf = kBeff - h * h * muB * water_film_flow_rate / force_mag;
+
+		    /* not counting films */
+                    double krnf_low = (1.0 - current_saturation) * krnf;
+		    double krwf_low = current_saturation * krwf;
+
                     // Saturation normalized effective permeability to account for decoupled phases and
                     // effective porosity.
                     double kAeff_low = (1.0 - current_saturation) * h * h *
@@ -935,6 +953,10 @@ double ScaLBL_ColorModel::Run(int returntime) {
                         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.film "
+                                             "eff.perm.water.film ");
+			fprintf(kr_log_file, "eff.perm.oil.film.lower.bound "
+                                             "eff.perm.water.film.lower.bound ");
                         fprintf(kr_log_file, "eff.perm.oil.lower.bound "
                                              "eff.perm.water.lower.bound ");
                         fprintf(kr_log_file,
@@ -952,6 +974,8 @@ double ScaLBL_ColorModel::Run(int returntime) {
                     fprintf(kr_log_file, "%i %.5g ", CURRENT_TIMESTEP,
                             current_saturation);
                     fprintf(kr_log_file, "%.5g %.5g ", kAeff, kBeff);
+                    fprintf(kr_log_file, "%.5g %.5g ", krnf, krwf);
+                    fprintf(kr_log_file, "%.5g %.5g ", krnf_low, krwf_low);
                     fprintf(kr_log_file, "%.5g %.5g ", kAeff_low, kBeff_low);
                     fprintf(kr_log_file, "%.5g %.5g ", kAeff_connected,
                             kBeff_connected);
@@ -964,6 +988,53 @@ double ScaLBL_ColorModel::Run(int returntime) {
                     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.disconnected "
+					"eff.perm.water.disconnected ");
+				fprintf(scal_log_file,
+					"eff.perm.oil.film.lower.bound "
+					"eff.perm.water.film.lower.bound ");
+                    		fprintf(scal_log_file,
+					"cap.pressure "
+					"cap.pressure.connected "
+					"Ca " 
+					"eff.pressure\n");
+                    	}
+
+			// Adjust the effperms with porosity term to make the nominal
+			// values closer to measured data
+                    	fprintf(scal_log_file, "%i %.5g ", CURRENT_TIMESTEP,
+                    			current_saturation);
+			fprintf(scal_log_file, "%.5g %.5g ", kAeff_low * Mask->Porosity() * mDarcy_converter,
+				kBeff_low * Mask->Porosity() * mDarcy_converter);
+			fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected_low * Mask->Porosity() * mDarcy_converter,
+				kBeff_connected_low * Mask->Porosity() * mDarcy_converter);
+			fprintf(scal_log_file, "%.5g %.5g ", kAeff_disconnected * Mask->Porosity() * mDarcy_converter,
+				kBeff_disconnected * Mask->Porosity() * mDarcy_converter);
+			fprintf(scal_log_file, "%.5g %.5g ", krnf_low * Mask->Porosity() * mDarcy_converter,
+				krwf_low * Mask->Porosity() * mDarcy_converter);
+			fprintf(scal_log_file, "%.5g %.5g %.5g ", pAB,
+				pAB_connected, Ca);
+                    	fprintf(scal_log_file, "%.5g\n", eff_pres);
+                    	fclose(scal_log_file);
+
+                    }
+
                     printf("  Measured capillary number %f \n ", Ca);
                 }
                 if (SET_CAPILLARY_NUMBER) {

models/ColorModel.h

@@ -117,6 +117,7 @@ public:
     double tauA, tauB, rhoA, rhoB, alpha, beta;
     double Fx, Fy, Fz, flux;
     double din, dout, inletA, inletB, outletA, outletB;
+    const double mDarcy_converter = 1013.0;
 
     int Nx, Ny, Nz, N, Np;
     int rank, nprocx, nprocy, nprocz, nprocs;

models/MRTModel.cpp

@@ -20,13 +20,18 @@
 #include "models/MRTModel.h"
 #include "analysis/distance.h"
 #include "common/ReadMicroCT.h"
+
+
 ScaLBL_MRTModel::ScaLBL_MRTModel(int RANK, int NP, const Utilities::MPI &COMM)
     : rank(RANK), nprocs(NP), Restart(0), timestep(0), timestepMax(0), tau(0),
       Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), mu(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), comm(COMM) {}
+
+
 ScaLBL_MRTModel::~ScaLBL_MRTModel() {}
 
+
 void ScaLBL_MRTModel::ReadParams(string filename) {
     // read the input database
     db = std::make_shared<Database>(filename);
@@ -254,7 +259,7 @@ void ScaLBL_MRTModel::Run() {
 
         if (WriteHeader) {
             log_file = fopen("Permeability.csv", "a+");
-            fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n");
+            fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz absperm absperm_adjusted\n");
             fclose(log_file);
         }
     }
@@ -384,22 +389,25 @@ void ScaLBL_MRTModel::Run() {
             double h = Dm->voxel_length;
             double absperm =
                 h * h * mu * Mask->Porosity() * flow_rate / force_mag;
+            double absperm_adjusted =
+                h * h * mu * Mask->Porosity() * Mask->Porosity() * flow_rate / force_mag;
+	    absperm_adjusted *= 1013.0; // Convert to mDarcy
+
             if (rank == 0) {
                 printf("     %f\n", absperm);
                 FILE *log_file = fopen("Permeability.csv", "a");
                 fprintf(log_file,
                         "%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g "
-                        "%.8g %.8g\n",
+                        "%.8g %.8g %.8g\n",
                         timestep, Fx, Fy, Fz, mu, h * h * h * Vs, h * h * As,
-                        h * Hs, Xs, vax, vay, vaz, absperm);
+                        h * Hs, Xs, vax, vay, vaz, absperm, absperm_adjusted);
                 fclose(log_file);
             }
         }
     }
     //************************************************************************/
     if (rank == 0)
-        printf("---------------------------------------------------------------"
-               "----\n");
+        printf("--------------------------------------------------------\n");
     // Compute the walltime per timestep
     auto t2 = std::chrono::system_clock::now();
     double cputime = std::chrono::duration<double>(t2 - t1).count() / timestep;

models/PoissonSolver.h

@@ -34,7 +34,7 @@ public:
     void ReadInput();
     void Create();
     void Initialize(double time_conv_from_Study);
-    void Run(double *ChargeDensity, int timestep_from_Study);
+    void Run(double *ChargeDensity, bool UseSlippingVelBC, int timestep_from_Study);
     void getElectricPotential(DoubleArray &ReturnValues);
     void getElectricPotential_debug(int timestep);
     void getElectricField(DoubleArray &Values_x, DoubleArray &Values_y,
@@ -45,22 +45,22 @@ public:
     //bool Restart,pBC;
     int timestep, timestepMax;
     int analysis_interval;
-    int BoundaryConditionInlet;
-    int BoundaryConditionOutlet;
-    int BoundaryConditionSolid;
-    double tau;
-    double tolerance;
+	int BoundaryConditionInlet;
+	int BoundaryConditionOutlet;
+    vector<int> BoundaryConditionSolidList;
+	double tau;
+	double tolerance;
     std::string tolerance_method;
     double k2_inv;
     double epsilon0, epsilon0_LB, epsilonR, epsilon_LB;
     double Vin, Vout;
     double chargeDen_dummy; //for debugging
     bool WriteLog;
-    double Vin0, freqIn, t0_In, Vin_Type;
-    double Vout0, freqOut, t0_Out, Vout_Type;
-    bool TestPeriodic;
-    double TestPeriodicTime;         //unit: [sec]
-    double TestPeriodicTimeConv;     //unit [sec/lt]
+    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]
     double TestPeriodicSaveInterval; //unit [sec]
 
     int Nx, Ny, Nz, N, Np;
@@ -86,7 +86,8 @@ public:
     int *dvcMap;
     //signed char *dvcID;
     double *fq;
-    double *Psi;
+    double *Psi; 
+    int *Psi_BCLabel;
     double *ElectricField;
     double *ChargeDensityDummy; // for debugging
     double *ResidualError;
@@ -103,17 +104,17 @@ private:
     char OutputFilename[200];
 
     //int rank,nprocs;
-    void LoadParams(std::shared_ptr<Database> db0);
-    void AssignSolidBoundary(double *poisson_solid);
+    void LoadParams(std::shared_ptr<Database> db0);    	
+    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);
     void SolveElectricPotentialAAeven(int timestep_from_Study);
     //void SolveElectricField();
-    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);
+    void SolvePoissonAAodd(double *ChargeDensity, bool UseSlippingVelBC);
+    void SolvePoissonAAeven(double *ChargeDensity, bool UseSlippingVelBC);
+    void getConvergenceLog(int timestep,double error);
+    double getBoundaryVoltagefromPeriodicBC(double V0,double freq,double t0,int time_step);
+    
 };
 #endif

models/StokesModel.cpp

@@ -593,7 +593,7 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
         ScaLBL_Comm->SendD3Q19AA(fq); //READ FROM NORMAL
         ScaLBL_D3Q19_AAodd_StokesMRT(
             NeighborList, fq, Velocity, ChargeDensity, ElectricField, rlx_setA,
-            rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv,
+            rlx_setB, Fx, Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC,
             ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
         // Set boundary conditions
@@ -610,8 +610,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
         }
         ScaLBL_D3Q19_AAodd_StokesMRT(NeighborList, fq, Velocity, ChargeDensity,
                                      ElectricField, rlx_setA, rlx_setB, Fx, Fy,
-                                     Fz, rho0, den_scale, h, time_conv, 0,
-                                     ScaLBL_Comm->LastExterior(), Np);
+                                     Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, 
+                                     0, ScaLBL_Comm->LastExterior(), Np);
 
         if (UseSlippingVelBC == true) {
             ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
@@ -626,8 +626,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
         ScaLBL_Comm->SendD3Q19AA(fq); //READ FORM NORMAL
         ScaLBL_D3Q19_AAeven_StokesMRT(
             fq, Velocity, ChargeDensity, ElectricField, rlx_setA, rlx_setB, Fx,
-            Fy, Fz, rho0, den_scale, h, time_conv, ScaLBL_Comm->FirstInterior(),
-            ScaLBL_Comm->LastInterior(), Np);
+            Fy, Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, 
+            ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior(), Np);
         ScaLBL_Comm->RecvD3Q19AA(fq); //WRITE INTO OPPOSITE
         // Set boundary conditions
         if (BoundaryCondition == 3) {
@@ -643,8 +643,8 @@ void ScaLBL_StokesModel::Run_Lite(double *ChargeDensity,
         }
         ScaLBL_D3Q19_AAeven_StokesMRT(fq, Velocity, ChargeDensity,
                                       ElectricField, rlx_setA, rlx_setB, Fx, Fy,
-                                      Fz, rho0, den_scale, h, time_conv, 0,
-                                      ScaLBL_Comm->LastExterior(), Np);
+                                      Fz, rho0, den_scale, h, time_conv, UseSlippingVelBC, 
+                                      0, ScaLBL_Comm->LastExterior(), Np);
         if (UseSlippingVelBC == true) {
             ScaLBL_Comm->SolidSlippingVelocityBCD3Q19(
                 fq, ZetaPotentialSolid, ElectricField, SolidGrad, epsilon_LB,

sample_scripts/Hipify.sh

@@ -0,0 +1,9 @@
+#!/bin/bash
+
+echo "Hipify cuda file $1"
+
+sed -i 's/cudaError_t/hipError_t/g' $1
+sed -i 's/cudaGetLastError/hipGetLastError/g' $1
+sed -i 's/cudaSuccess/hipSuccess/g' $1
+sed -i 's/cudaGetErrorString/hipGetErrorString/g' $1
+sed -i 's/CUDA error/hip error/g' $1

sample_scripts/configure_arden

@@ -5,19 +5,16 @@ cmake -D CMAKE_C_COMPILER:PATH=/opt/arden/openmpi/3.1.2/bin/mpicc          \
     -D CMAKE_C_FLAGS="-O3 -fPIC"         \
     -D CMAKE_CXX_FLAGS="-O3 -fPIC "      \
     -D CMAKE_CXX_STANDARD=14    \
-    -D MPIEXEC=mpirun                     \
+    -D MPIEXEC=mpirun  \
+    -D CMAKE_BUILD_PARALLEL_LEVEL=4      \
     -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
     -D CMAKE_BUILD_TYPE:STRING=Release     \
     -D CUDA_FLAGS="-arch sm_35"          \
     -D CUDA_HOST_COMPILER="/usr/bin/gcc" \
+    -D USE_HDF5=1 \
     -D HDF5_DIRECTORY="/opt/arden/hdf5/1.8.12" \
-    -D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\
-    -D USE_SILO=1 \
-    -D SILO_LIB="/opt/arden/silo/4.10.2/lib/libsiloh5.a" \
-    -D SILO_DIRECTORY="/opt/arden/silo/4.10.2" \
-    -D USE_NETCDF=0  \
-    -D NETCDF_DIRECTORY="/opt/arden/netcdf/4.6.1" \
     -D USE_CUDA=0                        \
     -D USE_TIMER=0 \
     ~/Programs/LBPM-WIA
 
+#    -D HDF5_LIB="/opt/arden/hdf5/1.8.12/lib/libhdf5.a"\

sample_scripts/configure_summit

@@ -2,19 +2,23 @@
 #module load cmake
 
 #source /gpfs/gpfs_stage1/b6p315aa/setup/setup-mpi.sh
-module load cmake gcc
+module load cmake gcc/7.5.0
 module load cuda
-#/ccs/proj/csc380/mcclurej
-export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
-export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
-export NETCDF_DIR=/ccs/proj/geo136/install/netcdf/4.6.1
+module load hdf5
 
+#/ccs/proj/csc380/mcclurej
+#export HDF5_DIR=/ccs/proj/csc380/mcclurej/install/hdf5/1.8.12/
+#export SILO_DIR=/ccs/proj/csc380/mcclurej/install/silo/4.10.2/
+#export NETCDF_DIR=/ccs/proj/geo136/install/netcdf/4.6.1
+export HDF5_DIR="$OLCF_HDF5_ROOT"
 # configure
+
 rm -rf CMake*
 cmake                                       \
     -D CMAKE_BUILD_TYPE:STRING=Release      \
     -D CMAKE_C_COMPILER:PATH=mpicc          \
     -D CMAKE_CXX_COMPILER:PATH=mpicxx       \
+    -D CMAKE_CXX_FLAGS="-lstdc++"	    \
     -D CMAKE_CXX_STANDARD=14                \
     -D USE_CUDA=1                           \
         -D CMAKE_CUDA_FLAGS="-arch sm_70 -Xptxas=-v -Xptxas -dlcm=cg -lineinfo" \
@@ -25,7 +29,7 @@ cmake                                       \
     -D USE_HDF5=1                        \
         -D HDF5_DIRECTORY="$HDF5_DIR"		   \
     	-D HDF5_LIB="$HDF5_DIR/lib/libhdf5.a"	   \
-     -D USE_SILO=1                        \
+     -D USE_SILO=0                        \
         -D SILO_LIB="$SILO_DIR/lib/libsiloh5.a"    \
    	-D SILO_DIRECTORY="$SILO_DIR"		     \
     -D USE_NETCDF=0				     \

sample_scripts/configure_tinker_A100

@@ -0,0 +1,33 @@
+# load the module for cmake
+module reset
+module load CMake
+#module load CUDA
+#module load OpenMPI
+#module load cuda11.2/toolkit
+ 
+module load OpenMPI/4.0.5-gcccuda-2020b
+module load HDF5
+
+#export GCC_BIN=/apps/easybuild/software/tinkercliffs-rome_a100/GCCcore/10.2.0/bin
+
+cmake                                    \
+    -D CMAKE_C_COMPILER:PATH=mpicc          \
+    -D CMAKE_CXX_COMPILER:PATH=mpicxx        \
+    -D CMAKE_C_FLAGS="-fPIC"            \
+    -D CMAKE_CXX_FLAGS="-fPIC"          \
+    -D MPIEXEC=mpirun                     \
+    -D USE_EXT_MPI_FOR_SERIAL_TESTS:BOOL=TRUE \
+    -D CMAKE_BUILD_TYPE:STRING=Release     \
+    -D USE_CUDA=1                        \
+    -D CMAKE_CUDA_FLAGS="-arch sm_80"   \
+    -D CMAKE_CUDA_HOST_COMPILER="$GCC_BIN/gcc" \
+    -D USE_HDF5=1					 \
+       -D HDF5_DIRECTORY=$HDF5_DIR \
+    -D USE_SILO=0				 \
+       -D SILO_DIRECTORY=$SILO_DIR \
+    -D USE_TIMER=0				 \
+    ~/LBPM-WIA
+
+make VERBOSE=1 -j4 && make install
+
+ 

tests/CMakeLists.txt

@@ -39,6 +39,7 @@ ADD_LBPM_EXECUTABLE( convertIO )
 ADD_LBPM_EXECUTABLE( DataAggregator )
 #ADD_LBPM_EXECUTABLE( BlobAnalyzeParallel )(
 ADD_LBPM_EXECUTABLE( lbpm_minkowski_scalar )
+ADD_LBPM_EXECUTABLE( lbpm_TwoPhase_analysis )
 ADD_LBPM_EXECUTABLE( TestPoissonSolver )
 ADD_LBPM_EXECUTABLE( TestIonModel )
 ADD_LBPM_EXECUTABLE( TestNernstPlanck )

tests/TestInterfaceSpeed.cpp

@@ -107,8 +107,8 @@ int main (int argc, char *argv[])
 	printf("Area ws = %f, Analytical = %f \n", Averages->aws, 4*PI*RADIUS*HEIGHT);
 	printf("Area s = %f, Analytical = %f \n", Averages->As, 2*PI*RADIUS*(Nz-2));
 	printf("Length wns = %f, Analytical = %f \n", Averages->lwns, 4*PI*RADIUS);
-	printf("Geodesic curvature (wns) = %f, Analytical = %f \n", Averages->KGwns_global, 0.0);
-	printf("Normal curvature (wns) = %f, Analytical = %f \n", Averages->KNwns_global, 1.0/RADIUS);
+	printf("Geodesic curvature (wn) = %f, Analytical = %f \n", Averages->KGwns_global, 0.0);
+	printf("Geodesic curvature (ws) = %f, Analytical = %f \n", Averages->KNwns_global, 4*PI);
 //	printf("Cos(theta_wns) = %f, Analytical = %f \n",efawns/lwns,1.0*RADIUS/CAPRAD);
 	printf("Interface Velocity = %f,%f,%f \n",Averages->vawn_global(0),Averages->vawn_global(1),Averages->vawn_global(2));
 	printf("Common Curve Velocity = %f,%f,%f \n",Averages->vawns_global(0),Averages->vawns_global(1),Averages->vawns_global(2));

tests/TestNernstPlanck.cpp

@@ -74,7 +74,7 @@ int main(int argc, char **argv)
         while (timestep < Study.timestepMax && error > Study.tolerance){
             
             timestep++;
-            PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental
+            PoissonSolver.Run(IonModel.ChargeDensity,false,0);//solve Poisson equtaion to get steady-state electrical potental
             IonModel.Run(IonModel.FluidVelocityDummy,PoissonSolver.ElectricField); //solve for ion transport and electric potential
             
             timestep++;//AA operations

tests/TestPNP_Stokes.cpp

@@ -94,7 +94,7 @@ int main(int argc, char **argv)
         while (timestep < Study.timestepMax && error > Study.tolerance){
             
             timestep++;
-            PoissonSolver.Run(IonModel.ChargeDensity,0);//solve Poisson equtaion to get steady-state electrical potental
+            PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,0);//solve Poisson equtaion to get steady-state electrical potental
             StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
             IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
             

tests/TestPoissonSolver.cpp

@@ -69,7 +69,7 @@ int main(int argc, char **argv)
             int timeSave = int(PoissonSolver.TestPeriodicSaveInterval/PoissonSolver.TestPeriodicTimeConv);
             while (timestep<timeMax){
                 timestep++;
-                PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,timestep);
+                PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,timestep);
                 if (timestep%timeSave==0){
                     if (rank==0) printf("   Time = %.3g[s]; saving electric potential and field\n",timestep*PoissonSolver.TestPeriodicTimeConv);
                     PoissonSolver.getElectricPotential_debug(timestep);
@@ -78,7 +78,7 @@ int main(int argc, char **argv)
             }
         }
         else {
-            PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,1);
+            PoissonSolver.Run(PoissonSolver.ChargeDensityDummy,false,1);
             PoissonSolver.getElectricPotential_debug(1);
             PoissonSolver.getElectricField_debug(1);
         }

tests/lbpm_TwoPhase_analysis.cpp

@@ -0,0 +1,239 @@
+/*
+ * Compute porous media marching cubes analysis (PMMC) based on input image data 
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <functional>
+
+#include "common/Array.h"
+#include "common/Domain.h"
+#include "common/Communication.h"
+#include "common/MPI.h"
+#include "IO/MeshDatabase.h"
+#include "IO/Mesh.h"
+#include "IO/Writer.h"
+#include "IO/netcdf.h"
+#include "analysis/analysis.h"
+#include "analysis/filters.h"
+#include "analysis/distance.h"
+#include "analysis/Minkowski.h"
+#include "analysis/TwoPhase.h"
+
+#include "ProfilerApp.h"
+
+int main(int argc, char **argv)
+{
+
+	// Initialize MPI
+        Utilities::startup( argc, argv );
+	Utilities::MPI comm( MPI_COMM_WORLD );
+        int rank = comm.getRank();
+        //int nprocs = comm.getSize();
+	{
+		Utilities::setErrorHandlers();
+		PROFILE_START("Main");
+
+		//std::vector<std::string> filenames;
+		if ( argc<2 ) {
+			if ( rank == 0 ){
+				printf("At least one filename must be specified\n");
+			}
+			return 1;
+		}
+		std::string filename = std::string(argv[1]);
+		if ( rank == 0 ){
+			printf("Input data file: %s\n",filename.c_str());
+		}
+
+		auto db = std::make_shared<Database>( filename );
+		auto domain_db = db->getDatabase( "Domain" );
+		auto vis_db = db->getDatabase( "Visualization" );
+
+		// Read domain parameters
+		auto Filename = domain_db->getScalar<std::string>( "Filename" );
+		auto size = domain_db->getVector<int>( "n" );
+		auto SIZE = domain_db->getVector<int>( "N" );
+		auto nproc = domain_db->getVector<int>( "nproc" );
+		auto ReadValues = domain_db->getVector<char>( "ReadValues" );
+		auto WriteValues = domain_db->getVector<char>( "WriteValues" );
+
+		auto nx = size[0];
+		auto ny = size[1];
+		auto nz = size[2];
+		int i,j,k,n;
+
+		std::shared_ptr<Domain> Dm  = std::shared_ptr<Domain>(new Domain(domain_db,comm));      // full domain for analysis
+		comm.barrier();
+		for (k=0;k<nz;k++){
+			for (j=0;j<ny;j++){
+				for (i=0;i<nx;i++){
+					n = k*nx*ny+j*nx+i;
+					// Initialize the object
+					Dm->id[n] = 1;
+				}
+			}
+		}
+		Dm->CommInit();
+		
+		/* read the data */
+		if (domain_db->keyExists( "Filename" )){
+			auto Filename = domain_db->getScalar<std::string>( "Filename" );
+			Dm->Decomp(Filename);
+		}
+		else{
+			Dm->ReadIDs();
+		}
+		
+	    fillHalo<double> fillDouble(Dm->Comm, Dm->rank_info, {nx,ny,nz}, {1, 1, 1}, 0, 1);
+	    		
+		// Compute the Minkowski functionals
+		comm.barrier();				
+		std::shared_ptr<TwoPhase> Averages( new TwoPhase(Dm) );
+		std::shared_ptr<Minkowski> Geometry(new Minkowski(Dm));
+
+		// Calculate the distance		
+		// Initialize the domain and communication
+		nx+=2; ny+=2; nz+=2;
+		Array<char> id(nx,ny,nz);
+		DoubleArray Distance(nx,ny,nz);
+
+		/* * * * * * * * * * * * * * * * * * * * * */
+		// Solve for the position of the solid phase
+		for (k=0;k<nz;k++){
+			for (j=0;j<ny;j++){
+				for (i=0;i<nx;i++){
+					n = k*nx*ny+j*nx+i;
+					// Initialize the object
+					if (Dm->id[n] == ReadValues[0])		id(i,j,k) = 0;
+					else		      					id(i,j,k) = 1;
+				}
+			}
+		}
+		for (k=0;k<nz;k++){
+			for (j=0;j<ny;j++){
+				for (i=0;i<nx;i++){
+					n=k*nx*ny+j*nx+i;
+					// Initialize distance to +/- 1
+					Averages->SDs(i,j,k) = 2.0*double(id(i,j,k))-1.0;
+				}
+			}
+		}
+		fillDouble.fill(Averages->SDs);
+		if (rank==0) printf("Initialized solid phase -- Converting to Signed Distance function \n");
+		CalcDist(Averages->SDs,id,*Dm);
+		
+		/* move the solid surface by a voxel to improve contact angle measurement
+		for (k=0;k<nz;k++){
+			for (j=0;j<ny;j++){
+				for (i=0;i<nx;i++){
+					n=k*nx*ny+j*nx+i;
+					Averages->SDs(i,j,k) -= 1.0;
+				}
+			}
+		}*/
+		/* * * * * * * * * * * * * * * * * * * * * */
+		// Solve for the position of the fluid phase
+		for (k=0;k<nz;k++){
+			for (j=0;j<ny;j++){
+				for (i=0;i<nx;i++){
+					n = k*nx*ny+j*nx+i;
+					// Initialize the object
+					if (Dm->id[n] == ReadValues[1]){
+						id(i,j,k) = 1;
+						Averages->Phase(i,j,k) = 1.0;
+					}
+					else	{
+						id(i,j,k) = 0;
+						Averages->Phase(i,j,k) = -1.0;
+					}
+				}
+			}
+		}
+		for (k=0;k<nz;k++){
+			for (j=0;j<ny;j++){
+				for (i=0;i<nx;i++){
+					n=k*nx*ny+j*nx+i;
+					// Initialize distance to +/- 1
+					Distance(i,j,k) = 2.0*double(id(i,j,k))-1.0;
+				}
+			}
+		}
+		fillDouble.fill(Averages->Phase);
+		fillDouble.fill(Distance);
+
+		if (rank==0) printf("Initialized fluid phase -- Converting to Signed Distance function \n");
+		CalcDist(Distance,id,*Dm);
+		Mean3D(Distance,Averages->SDn);
+		/* * * * * * * * * * * * * * * * * * * * * */
+
+		if (rank==0) printf("Computing Minkowski functionals \n");
+		Geometry->ComputeScalar(Distance,0.f);
+		Geometry->PrintAll();
+
+		Averages->Initialize();
+		Averages->UpdateMeshValues();
+		Averages->ComputeStatic();
+		Averages->Reduce();
+		Averages->PrintStatic();
+		/*        
+        Averages.Initialize();
+        Averages.ComponentAverages();
+        Averages.SortBlobs();
+        Averages.PrintComponents(timestep);
+		 */      
+
+		auto format = vis_db->getWithDefault<string>("format", "hdf5");
+
+		vis_db = db->getDatabase("Visualization");
+
+		std::vector<IO::MeshDataStruct> visData;
+		fillHalo<double> fillData(Dm->Comm, Dm->rank_info,
+				{Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2},
+				{1, 1, 1}, 0, 1);
+
+		auto PhaseVar = std::make_shared<IO::Variable>();
+		auto SignDistVar = std::make_shared<IO::Variable>();
+
+		IO::initialize("", format, "false");
+		// Create the MeshDataStruct
+		visData.resize(1);
+		visData[0].meshName = "domain";
+		visData[0].mesh = std::make_shared<IO::DomainMesh>(
+				Dm->rank_info, Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2, Dm->Lx, Dm->Ly,
+				Dm->Lz);
+		SignDistVar->name = "SignDist";
+		SignDistVar->type = IO::VariableType::VolumeVariable;
+		SignDistVar->dim = 1;
+		SignDistVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2);
+		visData[0].vars.push_back(SignDistVar);
+
+		PhaseVar->name = "Phase";
+		PhaseVar->type = IO::VariableType::VolumeVariable;
+		PhaseVar->dim = 1;
+		PhaseVar->data.resize(Dm->Nx - 2, Dm->Ny - 2, Dm->Nz - 2);
+		visData[0].vars.push_back(PhaseVar);
+
+		Array<double> &SignData = visData[0].vars[0]->data;
+		Array<double> &PhaseData = visData[0].vars[1]->data;
+
+		ASSERT(visData[0].vars[0]->name == "SignDist");
+		ASSERT(visData[0].vars[1]->name == "Phase");
+
+		fillData.copy(Averages->SDs, SignData);
+		fillData.copy(Averages->SDn, PhaseData);
+
+		int timestep = 0;
+		IO::writeData(timestep, visData, Dm->Comm);
+	}
+	PROFILE_STOP("Main");
+	PROFILE_SAVE("TwoPhase",true);
+	Utilities::shutdown();
+	return 0;
+}
+
+

tests/lbpm_color_simulator.cpp

@@ -32,10 +32,10 @@ int main( int argc, char **argv )
 		int nprocs = comm.getSize();
 		std::string SimulationMode = "production";
 		// Load the input database
-		auto db = std::make_shared<Database>( argv[1] );
-		if (argc > 2) {
-			SimulationMode = "legacy";
-		}
+		//auto db = std::make_shared<Database>( argv[1] );
+		//if (argc > 2) {
+		//	SimulationMode = "legacy";
+		//}
 
 		if ( rank == 0 ) {
 			printf( "********************************************************\n" );
@@ -187,12 +187,13 @@ int main( int argc, char **argv )
 				}
 			}
 		}
-		
+		/*
 		PROFILE_STOP( "Main" );
 		auto file  = db->getWithDefault<std::string>( "TimerFile", "lbpm_color_simulator" );
 		auto level = db->getWithDefault<int>( "TimerLevel", 1 );
 		NULL_USE(level);
 		PROFILE_SAVE( file, level );
+		*/
 		// ****************************************************
 
 

tests/lbpm_electrokinetic_SingleFluid_simulator.cpp

@@ -94,7 +94,7 @@ int main(int argc, char **argv)
         while (timestep < Study.timestepMax){
             
             timestep++;
-            PoissonSolver.Run(IonModel.ChargeDensity,timestep);//solve Poisson equtaion to get steady-state electrical potental
+            PoissonSolver.Run(IonModel.ChargeDensity,StokesModel.UseSlippingVelBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
             StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
             IonModel.Run(StokesModel.Velocity,PoissonSolver.ElectricField); //solve for ion transport and electric potential
             

tests/lbpm_morphdrain_pp.cpp

@@ -53,6 +53,7 @@ int main(int argc, char **argv)
 		auto WriteValues = domain_db->getVector<int>( "WriteValues" );
 		SW = domain_db->getScalar<double>("Sw");
 		auto READFILE = domain_db->getScalar<std::string>( "Filename" );
+		auto MORPH_RADIUS = domain_db->getWithDefault<double>( "MorphRadius", 100000.0);
 
 		// Generate the NWP configuration
 		//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
@@ -122,7 +123,7 @@ int main(int argc, char **argv)
 		comm.barrier();
 
 		// Run the morphological opening
-		MorphDrain(SignDist, id, Dm, SW);
+		MorphDrain(SignDist, id, Dm, SW, MORPH_RADIUS);
 	
 		// calculate distance to non-wetting fluid
 		if (domain_db->keyExists( "HistoryLabels" )){

tests/lbpm_morphopen_pp.cpp

@@ -81,7 +81,7 @@ int main(int argc, char **argv)
 		id = new signed char [N];
 		Mask->Decomp(READFILE);
 		Mask->CommInit();
-
+		
 		// Generate the NWP configuration
 		//if (rank==0) printf("Initializing morphological distribution with critical radius %f \n", Rcrit);
 		if (rank==0) printf("Performing morphological opening with target saturation %f \n", SW);

tests/test_morph.cpp

@@ -94,10 +94,11 @@ int main(int argc, char **argv)
 		printf("   Measure the opening \n");
 		sphere.MeasureObject();
 		//sphere.ComputeScalar(Distance,0.f);
-		printf("   Volume  = %f (analytical = %f) \n", sphere.Vi,0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
-		double error = fabs(sphere.Vi - 0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)))/ (0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
+		/* Note 0.856 = (0.95)^3 */
+		printf("   Volume  = %f (analytical = %f) \n", sphere.Vi,0.256*0.33333333333333*0.856*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
+		double error = fabs(sphere.Vi - 0.256*0.856*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)))/ (0.256*0.33333333333333*3.14159*double((Nx-2)*(Nx-2)*(Nx-2)));
 		printf("   Relative error %f \n", error);
-		if (error  > 0.03){
+		if (error  > 0.05){
 			toReturn = 10;
 			printf("ERROR (test_morph): difference from analytical volume is too large\n");