From 7e286f16f16ef4c4f5198f25adf3aaa25e318d4e Mon Sep 17 00:00:00 2001 From: James McClure Date: Mon, 16 Jan 2023 13:29:58 -0500 Subject: [PATCH] update docs for cell model --- docs/source/examples/membrane/cellModel.rst | 19 +++++--- docs/source/userGuide/models/cell/cell.rst | 30 ++++++++++--- docs/source/userGuide/models/index.rst | 4 +- example/SingleCell/NaCl-cell.py | 39 ++++++++-------- example/SingleCell/NaCl.db | 38 +++++++++------- models/IonModel.cpp | 50 +++++++++++++++------ 6 files changed, 114 insertions(+), 66 deletions(-) diff --git a/docs/source/examples/membrane/cellModel.rst b/docs/source/examples/membrane/cellModel.rst index 54876065..e6e7048c 100644 --- a/docs/source/examples/membrane/cellModel.rst +++ b/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 - -.. code:: bash - - mpirun -np 2 $LBPM_BIN/lbpm_nernst_planck_simulator bacterium.db + +******************* +Example Output +******************* Successful output looks like the following - .. code:: bash ******************************************************** diff --git a/docs/source/userGuide/models/cell/cell.rst b/docs/source/userGuide/models/cell/cell.rst index 1e0d987a..262cf669 100644 --- a/docs/source/userGuide/models/cell/cell.rst +++ b/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 diff --git a/docs/source/userGuide/models/index.rst b/docs/source/userGuide/models/index.rst index dbae0914..39a526c0 100644 --- a/docs/source/userGuide/models/index.rst +++ b/docs/source/userGuide/models/index.rst @@ -12,9 +12,7 @@ Currently supported lattice Boltzmann models mrt/* - nernstPlanck/* - - PoissonBoltzmann/* + cell/* greyscale/* diff --git a/example/SingleCell/NaCl-cell.py b/example/SingleCell/NaCl-cell.py index dcb0adbf..d7828924 100644 --- a/example/SingleCell/NaCl-cell.py +++ b/example/SingleCell/NaCl-cell.py @@ -1,41 +1,38 @@ import numpy as np import matplotlib.pylab as plt -Nx = 64 -Ny = 64 -Nz = 64 -cx = Nx/2 -cy = Ny/2 -cz = Nz/2 -radius = 12 +D=np.ones((40,40,40),dtype="uint8") -D=np.ones((Nx,Ny,Nz),dtype="uint8") +cx = 20 +cy = 20 +cz = 20 +radius = 8 -for i in range(0, Nx): - for j in range (0, Ny): - for k in range (0,Nz): +for i in range(0, 40): + for j in range (0, 40): + 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") -C2=np.zeros((Nx,Ny,Nz),dtype="double") +C1=np.zeros((40,40,40),dtype="double") +C2=np.zeros((40,40,40),dtype="double") -for i in range(0, Nx): - for j in range (0, Ny): - for k in range (0,Nz): +for i in range(0, 40): + for j in range (0, 40): + 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 - C2[i,j,k] = 110.0e-6 + C1[i,j,k] = 5.0e-6 + C2[i,j,k] = 5.0e-6 -C1.tofile("cell_concentration_Na_64x64x64.raw") -C2.tofile("cell_concentration_Cl_64x64x64.raw") +C1.tofile("cell_concentration_Na_40x40x40.raw") +C2.tofile("cell_concentration_Cl_40x40x40.raw") diff --git a/example/SingleCell/NaCl.db b/example/SingleCell/NaCl.db index 95d222e0..c297a898 100644 --- a/example/SingleCell/NaCl.db +++ b/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 - analysis_interval = 25 //timestep checking steady-state convergence - tolerance = 1.0e-10 //stopping criterion for steady-state solution + timestepMax = 100000 //max timestep for obtaining steady-state electrical potential + analysis_interval = 200 //timestep checking steady-state convergence + 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 = 64, 64, 64 // size of the input image - voxel_length = 0.1 //resolution; user-input unit: [um] + n = 40, 40, 40 // Size of local domain (Nx,Ny,Nz) + N = 40, 40, 40 // size of the input image + 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 - MassFractionOut = 1e-2, 1e-8 - ThresholdMassFractionIn = 1e-2, 1e-8 - ThresholdMassFractionOut = 1e-2, 1e-8 + MassFractionIn = 0e-2, 5e-2 + MassFractionOut = 0e-2, 5e-2 + ThresholdMassFractionIn = 0e-2, 5e-2 + ThresholdMassFractionOut = 0e-2, 5e-2 } diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 16679884..099b034e 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -300,7 +300,7 @@ void ScaLBL_IonModel::ReadParams(string filename, vector &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(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; - if (rank == 0) - printf(" ...initializing based on membrane list \n"); - Ci_host = new double[number_ion_species * Np]; - 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; + double *Ci_host; + Ci_host = new double[number_ion_species * Np]; + + if (ion_db->keyExists("IonConcentrationFile")) { + //NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype" + auto File_ion = ion_db->getVector("IonConcentrationFile"); + if (File_ion.size() == 2*number_ion_species) { + for (size_t ic = 0; ic < number_ion_species; ic++) { + AssignIonConcentration_FromFile(&Ci_host[ic * Np], File_ion, + 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); + } + } + } + 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("IonConcentrationFile");