added membrane properties to input db

common/Membrane.cpp

@@ -1351,46 +1351,11 @@ void Membrane::IonTransport(double *dist, double *den){
 }
 
 //	std::shared_ptr<Database> db){
-void Membrane::AssignCoefficients(int *Map, double *Psi, string method){
+void Membrane::AssignCoefficients(int *Map, double *Psi, double Threshold,
+		double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn,
+		double ThresholdMassFractionOut){
 	/* Assign mass transfer coefficients to the membrane data structure */
 	
-	double Threshold;
-	double MassFractionIn,MassFractionOut,ThresholdMassFractionIn,ThresholdMassFractionOut;
-	
-	Threshold = -55.0;
-	MassFractionIn = 0.0;
-	MassFractionOut = 0.0;
-	ThresholdMassFractionOut = 0.0;				
-	ThresholdMassFractionIn = 0.0;
-	
-	if (method == "ones"){
-		/* Initializing  */
-		//printf(".... initialize permeable membrane \n");
-		MassFractionIn = 1.0;
-		MassFractionOut = 1.0;
-		ThresholdMassFractionOut = 1.0;				
-		ThresholdMassFractionIn = 1.0;
-	}
-	if (method == "Voltage Gated Potassium"){
-		MassFractionIn = 0.0;
-		MassFractionOut = 0.0;
-		ThresholdMassFractionOut = 0.0;				
-		ThresholdMassFractionIn = 1.0;
-	}
-	if (method == "impermeable"){
-		//printf(".... impermeable membrane \n");
-		MassFractionIn = 0.001;
-		MassFractionOut = 0.001;
-		ThresholdMassFractionOut = 0.001;				
-		ThresholdMassFractionIn = 0.0001;
-	}
-	if (method == "Na+"){
-		//printf(".... Na+ permeable membrane \n");
-		MassFractionIn = 0.05;
-		MassFractionOut = 0.05;
-		ThresholdMassFractionOut = 0.05;				
-		ThresholdMassFractionIn = 0.05;
-	}
 	
 	ScaLBL_D3Q7_Membrane_AssignLinkCoef(MembraneLinks, Map, MembraneDistance, Psi, MembraneCoef,
 			Threshold, MassFractionIn, MassFractionOut, ThresholdMassFractionIn, ThresholdMassFractionOut,

common/Membrane.h

@@ -97,7 +97,9 @@ public:
 	void RecvD3Q19AA(double *dist);
 	void SendD3Q7AA(double *dist);
 	void RecvD3Q7AA(double *dist);
-	void AssignCoefficients(int *Map, double *Psi, std::string method);
+	void AssignCoefficients(int *Map, double *Psi, double Threshold,
+			double MassFractionIn, double MassFractionOut, double ThresholdMassFractionIn,
+			double ThresholdMassFractionOut);
 	void IonTransport(double *dist, double *den);
 	//......................................................................................
 	// Buffers to store data sent and recieved by this MPI process

models/IonModel.cpp

@@ -17,7 +17,6 @@ ScaLBL_IonModel::~ScaLBL_IonModel() {}
 
 void ScaLBL_IonModel::ReadParams(string filename, vector<int> &num_iter) {
 
-	USE_MEMBRANE = true;
     // read the input database
     db = std::make_shared<Database>(filename);
     domain_db = db->getDatabase("Domain");
@@ -283,7 +282,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 = true;
     // read the input database
     db = std::make_shared<Database>(filename);
     domain_db = db->getDatabase("Domain");
@@ -321,7 +320,9 @@ void ScaLBL_IonModel::ReadParams(string filename) {
     if (domain_db->keyExists("voxel_length")) { //default unit: um/lu
         h = domain_db->getScalar<double>("voxel_length");
     }
-
+    if (ion_db->keyExists("use_membrane")) {
+    	USE_MEMBRANE = ion_db->getScalar<bool>("use_membrane");
+    }
     // LB-Ion Model parameters
     //if (ion_db->keyExists( "timestepMax" )){
     //	timestepMax = ion_db->getScalar<int>( "timestepMax" );
@@ -423,23 +424,104 @@ void ScaLBL_IonModel::ReadParams(string filename) {
         }
     }
     
+    if (USE_MEMBRANE){
+        membrane_db = db->getDatabase("Membrane");
 
-    if (ion_db->keyExists("MembraneIonConcentrationList")) {
+    	/* get membrane permeability parameters*/
-      if (rank == 0) printf(".... Read MembraneIonConcentrationList \n");
+    	if (membrane_db->keyExists("MassFractionIn")) {
-    	MembraneIonConcentration.clear();
+    		if (rank == 0) printf(".... Read membrane permeability (MassFractionIn) \n");
-    	MembraneIonConcentration = ion_db->getVector<double>("MembraneIonConcentrationList");
+    		MassFractionIn.clear();
-    	if (MembraneIonConcentration.size() != number_ion_species) {
+    		MassFractionIn = membrane_db->getVector<double>("MassFractionIn");
-    		ERROR("Error: number_ion_species and MembraneIonConcentrationList must be "
+    		if (MassFractionIn.size() != number_ion_species) {
-    				"the same length! \n");
+    			ERROR("Error: number_ion_species and membrane permeability (MassFractionIn) must be "
-    	} 
+    					"the same length! \n");
-    	else {
-    		for (size_t i = 0; i < MembraneIonConcentration.size(); i++) {
-    			MembraneIonConcentration[i] =
-    					MembraneIonConcentration[i] *
-    					(h * h * h *
-    							1.0e-18); //LB ion concentration has unit [mol/lu^3]
     		} 
     	}	
+    	else{
+    		MassFractionIn.resize(IonConcentration.size());
+            for (size_t i = 0; i < IonConcentration.size(); i++) {
+            	MassFractionIn[i] = 0.0;
+            }
+    	}
+    	if (membrane_db->keyExists("MassFractionOut")) {
+    		if (rank == 0) printf(".... Read membrane permeability (MassFractionOut) \n");
+    		MassFractionOut.clear();
+    		MassFractionOut = membrane_db->getVector<double>("MassFractionOut");
+    		if (MassFractionIn.size() != number_ion_species) {
+    			ERROR("Error: number_ion_species and membrane permeability (MassFractionOut) must be "
+    					"the same length! \n");
+    		} 
+    	}	
+    	else{
+    		MassFractionOut.resize(IonConcentration.size());
+            for (size_t i = 0; i < IonConcentration.size(); i++) {
+            	MassFractionOut[i] = 0.0;
+            }
+    	}
+    	if (membrane_db->keyExists("ThresholdMassFractionIn")) {
+    		if (rank == 0) printf(".... Read membrane permeability (ThresholdMassFractionIn) \n");
+    		ThresholdMassFractionIn.clear();
+    		ThresholdMassFractionIn = membrane_db->getVector<double>("ThresholdMassFractionIn");
+    		if (ThresholdMassFractionIn.size() != number_ion_species) {
+    			ERROR("Error: number_ion_species and membrane permeability (ThresholdMassFractionIn) must be "
+    					"the same length! \n");
+    		} 
+    	}	
+    	else{
+    		ThresholdMassFractionIn.resize(IonConcentration.size());
+            for (size_t i = 0; i < IonConcentration.size(); i++) {
+            	ThresholdMassFractionIn[i] = 0.0;
+            }
+    	}	
+    	if (membrane_db->keyExists("ThresholdMassFractionOut")) {
+    		if (rank == 0) printf(".... Read membrane permeability (ThresholdMassFractionOut) \n");
+    		ThresholdMassFractionOut.clear();
+    		ThresholdMassFractionOut = membrane_db->getVector<double>("ThresholdMassFractionOut");
+    		if (ThresholdMassFractionOut.size() != number_ion_species) {
+    			ERROR("Error: number_ion_species and membrane permeability (ThresholdMassFractionOut) must be "
+    					"the same length! \n");
+    		} 
+    	}	
+    	else{
+    		ThresholdMassFractionOut.resize(IonConcentration.size());
+            for (size_t i = 0; i < IonConcentration.size(); i++) {
+            	ThresholdMassFractionOut[i] = 0.0;
+            }
+    	}
+    	if (membrane_db->keyExists("ThresholdVoltage")) {
+    		if (rank == 0) printf(".... Read membrane threshold (ThresholdVoltage) \n");
+    		ThresholdVoltage.clear();
+    		ThresholdVoltage = membrane_db->getVector<double>("ThresholdVoltage");
+    		if (ThresholdVoltage.size() != number_ion_species) {
+    			ERROR("Error: number_ion_species and membrane voltage threshold (ThresholdVoltage) must be "
+    					"the same length! \n");
+    		} 
+    	}	
+    	else{
+    		ThresholdVoltage.resize(IonConcentration.size());
+            for (size_t i = 0; i < IonConcentration.size(); i++) {
+            	ThresholdVoltage[i] = 0.0;
+            }
+    	}
+    	    	
+    	if (ion_db->keyExists("MembraneIonConcentrationList")) {
+    		if (rank == 0) printf(".... Read MembraneIonConcentrationList \n");
+    		MembraneIonConcentration.clear();
+    		MembraneIonConcentration = ion_db->getVector<double>("MembraneIonConcentrationList");
+    		if (MembraneIonConcentration.size() != number_ion_species) {
+    			ERROR("Error: number_ion_species and MembraneIonConcentrationList must be "
+    					"the same length! \n");
+    		} 
+    		else {
+    			for (size_t i = 0; i < MembraneIonConcentration.size(); i++) {
+    				MembraneIonConcentration[i] =
+    						MembraneIonConcentration[i] *
+    						(h * h * h *
+    								1.0e-18); //LB ion concentration has unit [mol/lu^3]
+    			}
+    		}
+    	}
+
     }
     //Read solid boundary condition specific to Ion model
     BoundaryConditionSolid = 0;
@@ -1362,11 +1444,9 @@ void ScaLBL_IonModel::RunMembrane(double *Velocity, double *ElectricField, doubl
 
     for (size_t ic = 0; ic < number_ion_species; ic++) {
         /* set the mass transfer coefficients for the membrane */
-    	if (ic == 0)
+		IonMembrane->AssignCoefficients(dvcMap, Psi, ThresholdVoltage[ic],MassFractionIn[ic],
-    		IonMembrane->AssignCoefficients(dvcMap, Psi, "Na+");
+				MassFractionOut[ic],ThresholdMassFractionIn[ic],ThresholdMassFractionOut[ic]);
-    	else {
+
-    		IonMembrane->AssignCoefficients(dvcMap, Psi, "impermeable");
-    	}
         timestep = 0;
         while (timestep < timestepMax[ic]) {
             //************************************************************************/

models/IonModel.h

@@ -69,10 +69,13 @@ public:
     vector<int> IonValence;
     vector<double> IonConcentration; //unit [mol/m^3]
     vector<double> MembraneIonConcentration; //unit [mol/m^3]
-    vector<double>
+    vector<double> ThresholdVoltage;
-        Cin; //inlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec]
+    vector<double> MassFractionIn;
-    vector<double>
+    vector<double> MassFractionOut;
-        Cout; //outlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec]
+    vector<double> ThresholdMassFractionIn;
+    vector<double> ThresholdMassFractionOut;
+    vector<double> Cin; //inlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec]
+    vector<double> Cout; //outlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec]
     vector<double> tau;
     vector<double> time_conv;
 
@@ -83,7 +86,7 @@ public:
     std::shared_ptr<Domain> Dm;   // this domain is for analysis
     std::shared_ptr<Domain> Mask; // this domain is for lbm
     std::shared_ptr<ScaLBL_Communicator> ScaLBL_Comm;
-    // input database
+    // input databaseF
     std::shared_ptr<Database> db;
     std::shared_ptr<Database> domain_db;
     std::shared_ptr<Database> ion_db;

tests/TestMembrane.cpp

@@ -200,7 +200,7 @@ int main(int argc, char **argv)
         ScaLBL_D3Q19_AAodd_Compact(M.NeighborList, gq, Np);
 
         /* explicit mass transfer step with the membrane*/
-        M.AssignCoefficients(dvcMap, Psi, "ones");
+        M.AssignCoefficients(dvcMap, Psi, 0.0, 1.0, 1.0, 1.0, 1.0);
         M.IonTransport(gq, Cj);
         ScaLBL_CopyToHost(Ci_host, Cj, sizeof(double) * Np);
         

tests/lbpm_cell_simulator.cpp

@@ -115,7 +115,7 @@ int main(int argc, char **argv)
             PoissonSolver.Run(IonModel.ChargeDensity,SlipBC,timestep);//solve Poisson equtaion to get steady-state electrical potental
 	    comm.barrier();
 	    //if (rank == 0) printf("    Poisson step %i \n",timestep);
-            StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
+           StokesModel.Run_Lite(IonModel.ChargeDensity, PoissonSolver.ElectricField);// Solve the N-S equations to get velocity
 	    //fflush(stdout);
 
             IonModel.RunMembrane(StokesModel.Velocity,PoissonSolver.ElectricField,PoissonSolver.Psi); //solve for ion transport with membrane