update docs for cell model

docs/source/examples/membrane/cellModel.rst

@@ -9,6 +9,14 @@ file ``Bacterium.swc``, which specifies an oblong cell shape, relying on the ``.
 is commonly used to approximate neuron structures. The case considered is the four ion membrane transport
 problem considered in Figure 4 from McClure & Li 
 
+The cell simulation is performed by the executable ``lbpm_nernst_planck_cell_simulator``, which is launched
+in the same way as other parallel tools
+
+.. code:: bash
+
+	  mpirun -np 2 $LBPM_BIN/lbpm_nernst_planck_cell_simulator Bacterium.db
+
+
 The input file ``Bacterium.db`` specifies the following
 
 .. code:: c
@@ -85,16 +93,13 @@ The input file ``Bacterium.db`` specifies the following
 	     ThresholdMassFractionIn = 1e-1, 1.0, 5e-3, 0.0
 	     ThresholdMassFractionOut = 1e-1, 1.0, 5e-3, 0.0
 	  }
-
+	  
-Once this has been set, we launch ``lbpm_nernst_planck_simulator`` in the same way as other parallel tools
+*******************
-
+Example Output
-.. code:: bash
+*******************
-
-	  mpirun -np 2 $LBPM_BIN/lbpm_nernst_planck_simulator bacterium.db
 
 Successful output looks like the following
 
-
 .. code:: bash
 
 	  ********************************************************

docs/source/userGuide/models/cell/cell.rst

@@ -3,6 +3,7 @@ Cell model
 =============================================
 
 LBPM includes a whole-cell simulator based on a coupled solution of the Nernst-Planck equations with Gauss's law. 
+The lattice Boltzmann formulation is described below.
 
 *********************
 Nernst-Planck model
@@ -157,6 +158,27 @@ The local value of the potential is then updated based on a relaxation scheme, w
    
 The algorithm can then proceed to the next timestep.
 
+***************************
+Membrane Model
+***************************
+
+The LBPM membrane model provides the basis to model cellular dynamics. There are currently two supported ways
+to specify the membrane location:
+
+1. provide a segemented image that is labeled to differentiate the cell
+interior and exterior. See the script ``NaCl-cell.py`` and input file ``NaCl.db`` as a reference for how to use labeled images.
+
+- ``IonConcentrationFile`` -- list of files that specify the initial concentration for each ion
+- ``Filename`` -- 8-bit binary file provided in the ``Domain`` section of the input database
+- ``ReadType`` -- this should be ``"8bit"`` (this is the default)
+
+2. provide a ``.swc`` file that specifies the geometry (see example input file below).
+ 
+- ``Filename`` -- swc file name should be provided in the ``Domain`` section of the input database
+- ``ReadType`` -- this should be ``"swc"`` (required since ``"8bit"`` is the internal default)
+
+Both examples are stored within the LBPM repository, located at ``example/SingleCell/``
+
 ****************************
 Example Input File
 ****************************
@@ -170,13 +192,9 @@ Example Input File
 	 tolerance = 1.0e-9
 	 visualization_interval = 1000        // Frequency to write visualization data
      }
-     Stokes {
-	 tau = 1.0
-	 F = 0, 0, 0
-	 ElectricField = 0, 0, 0 //body electric field; user-input unit: [V/m]
-	 nu_phys = 0.889e-6      //fluid kinematic viscosity; user-input unit: [m^2/sec]
-     }
      Ions {
+         use_membrane = true
+         Restart = false
 	 MembraneIonConcentrationList = 150.0e-3, 10.0e-3, 15.0e-3, 155.0e-3 //user-input unit: [mol/m^3]
 	 temperature = 293.15 //unit [K]
 	 number_ion_species = 4  //number of ions

docs/source/userGuide/models/index.rst

@@ -12,9 +12,7 @@ Currently supported lattice Boltzmann models
 
    mrt/*
 
-   nernstPlanck/*
+   cell/*
-
-   PoissonBoltzmann/*
    
    greyscale/*
 

example/SingleCell/NaCl-cell.py

@@ -1,41 +1,38 @@
 import numpy as np
 import matplotlib.pylab as plt
 
-Nx = 64
+D=np.ones((40,40,40),dtype="uint8")
-Ny = 64
-Nz = 64
-cx = Nx/2
-cy = Ny/2
-cz = Nz/2
-radius = 12
 
-D=np.ones((Nx,Ny,Nz),dtype="uint8")
+cx = 20
+cy = 20
+cz = 20
+radius = 8
 
-for i in range(0, Nx):
+for i in range(0, 40):
-    for j in range (0, Ny):
+    for j in range (0, 40):
-        for k in range (0,Nz):
+        for k in range (0,40):
             dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
             if (dist < radius ) :
                 D[i,j,k] = 2
 
-D.tofile("cell_64x64x64.raw")
+D.tofile("cell_40x40x40.raw")
 
                 
-C1=np.zeros((Nx,Ny,Nz),dtype="double")
+C1=np.zeros((40,40,40),dtype="double")
-C2=np.zeros((Nx,Ny,Nz),dtype="double")
+C2=np.zeros((40,40,40),dtype="double")
 
-for i in range(0, Nx):
+for i in range(0, 40):
-    for j in range (0, Ny):
+    for j in range (0, 40):
-        for k in range (0,Nz):
+        for k in range (0,40):
             #outside the cell
             C1[i,j,k] = 125.0e-6  # Na
             C2[i,j,k] = 125.0e-6  # Cl
             dist = np.sqrt((i-cx)*(i-cx) + (j-cx)*(j-cx) + (k-cz)*(k-cz))
             # inside the cell
             if (dist < radius ) :
-                C1[i,j,k] = 110.0e-6
+                C1[i,j,k] = 5.0e-6
-                C2[i,j,k] = 110.0e-6
+                C2[i,j,k] = 5.0e-6
 
-C1.tofile("cell_concentration_Na_64x64x64.raw")
+C1.tofile("cell_concentration_Na_40x40x40.raw")
-C2.tofile("cell_concentration_Cl_64x64x64.raw")
+C2.tofile("cell_concentration_Cl_40x40x40.raw")
 

example/SingleCell/NaCl.db

@@ -1,9 +1,10 @@
 MultiphysController {
-    timestepMax = 2000
+    timestepMax = 20
     num_iter_Ion_List = 2
-    analysis_interval  = 40
+    analysis_interval  = 50
     tolerance = 1.0e-9
-    visualization_interval = 40        // Frequency to write visualization data
+    visualization_interval = 100        // Frequency to write visualization data
+    analysis_interval = 50    // Frequency to perform analysis
 }
 Stokes {
     tau = 1.0
@@ -12,7 +13,9 @@ Stokes {
     nu_phys = 0.889e-6      //fluid kinematic viscosity; user-input unit: [m^2/sec]
 }
 Ions {
-    IonConcentrationFile = "cell_concentration_Na_64x64x64.raw", "double", "cell_concentration_Cl_64x64x64.raw", "double"
+    use_membrane = true
+    Restart = false
+    IonConcentrationFile = "cell_concentration_Na_40x40x40.raw", "double", "cell_concentration_Cl_40x40x40.raw", "double"
     temperature = 293.15 //unit [K]
     number_ion_species = 2  //number of ions
     tauList = 1.0, 1.0
@@ -25,6 +28,8 @@ Ions {
     FluidVelDummy = 0.0, 0.0, 0.0 // dummy fluid velocity for debugging
 }
 Poisson {
+    lattice_scheme = "D3Q19"
+    Restart = false
     epsilonR = 78.5 //fluid dielectric constant [dimensionless]
     BC_Inlet  = 0  // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
     BC_Outlet = 0  // ->1: fixed electric potential; ->2: sine/cosine periodic electric potential
@@ -41,16 +46,19 @@ Poisson {
     TestPeriodicTimeConv = 0.01 //unit:[sec]
     TestPeriodicSaveInterval = 0.2 //unit:[sec]
     //------------------------------ advanced setting ------------------------------------
-    timestepMax = 4000 //max timestep for obtaining steady-state electrical potential
+    timestepMax = 100000 //max timestep for obtaining steady-state electrical potential
-    analysis_interval  = 25 //timestep checking steady-state convergence
+    analysis_interval  = 200 //timestep checking steady-state convergence
-    tolerance = 1.0e-10  //stopping criterion for steady-state solution
+    tolerance = 1.0e-6  //stopping criterion for steady-state solution
+    InitialValueLabels = 1, 2
+    InitialValues = 0.0, 0.0
+
 }
 Domain {
-    Filename = "cell_64x64x64.raw"
+    Filename = "cell_40x40x40.raw"
     nproc = 1, 1, 1     // Number of processors (Npx,Npy,Npz)
-    n = 64, 64, 64      // Size of local domain (Nx,Ny,Nz)
+    n = 40, 40, 40      // Size of local domain (Nx,Ny,Nz)
-    N = 64, 64, 64         // size of the input image
+    N = 40, 40, 40         // size of the input image
-    voxel_length = 0.1   //resolution; user-input unit: [um]
+    voxel_length = 1.0   //resolution; user-input unit: [um]
     BC = 0              // Boundary condition type
     ReadType = "8bit"
     ReadValues  = 0, 1, 2
@@ -72,9 +80,9 @@ Visualization {
 Membrane {
     MembraneLabels = 2
     VoltageThreshold = 0.0, 0.0
-    MassFractionIn = 1e-2, 1e-8
+    MassFractionIn = 0e-2, 5e-2
-    MassFractionOut = 1e-2, 1e-8
+    MassFractionOut = 0e-2, 5e-2
-    ThresholdMassFractionIn = 1e-2, 1e-8
+    ThresholdMassFractionIn = 0e-2, 5e-2
-    ThresholdMassFractionOut = 1e-2, 1e-8
+    ThresholdMassFractionOut = 0e-2, 5e-2
 
 }

models/IonModel.cpp

@@ -300,7 +300,7 @@ void ScaLBL_IonModel::ReadParams(string filename, vector<int> &num_iter) {
 void ScaLBL_IonModel::ReadParams(string filename) {
     //NOTE: the maximum iteration timesteps for ions are left unspecified
     //      it relies on the multiphys controller to compute the max timestep
-	USE_MEMBRANE = false;
+	USE_MEMBRANE = true;
 	Restart = false;
     // read the input database
     db = std::make_shared<Database>(filename);
@@ -1081,20 +1081,42 @@ void ScaLBL_IonModel::Initialize() {
         printf("LB Ion Solver: initializing D3Q7 distributions\n");
     //USE_MEMBRANE = true; 
     if (USE_MEMBRANE){
-            double *Ci_host;
+    	double *Ci_host;
-            if (rank == 0)
+    	Ci_host = new double[number_ion_species * Np];
-                printf("   ...initializing based on membrane list \n");
+
-            Ci_host = new double[number_ion_species * Np];
+    	if (ion_db->keyExists("IonConcentrationFile")) {
-            for (size_t ic = 0; ic < number_ion_species; ic++) {
+    		//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
-            	AssignIonConcentrationMembrane( &Ci_host[ic * Np],  ic);
+    		auto File_ion = ion_db->getVector<std::string>("IonConcentrationFile");
-            }
+    		if (File_ion.size() == 2*number_ion_species) {
-            ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species * sizeof(double) * Np);
+    			for (size_t ic = 0; ic < number_ion_species; ic++) {
-            comm.barrier();
+    				AssignIonConcentration_FromFile(&Ci_host[ic * Np], File_ion,
-            for (size_t ic = 0; ic < number_ion_species; ic++) {
+    						ic);
-                ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np], Np);
+    			}
-            }
+    			ScaLBL_CopyToDevice(Ci, Ci_host,
-            delete[] Ci_host;
+    					number_ion_species * sizeof(double) * Np);
+    			comm.barrier();
+    			for (size_t ic = 0; ic < number_ion_species; ic++) {
+    				ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np],
+    						Np);
+    			}
+    		}
+    	}
+    	else{
+    		if (rank == 0)
+    			printf("   ...initializing based on membrane list \n");
+    		for (size_t ic = 0; ic < number_ion_species; ic++) {
+    			AssignIonConcentrationMembrane( &Ci_host[ic * Np],  ic);
+    		}
+    		ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species * sizeof(double) * Np);
+    		comm.barrier();
+    		for (size_t ic = 0; ic < number_ion_species; ic++) {
+    			ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic * Np * 7], &Ci[ic * Np], Np);
+    		}
+    	}
+    	delete[] Ci_host;
+
     }
+
     else if (ion_db->keyExists("IonConcentrationFile")) {
         //NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
         auto File_ion = ion_db->getVector<std::string>("IonConcentrationFile");