diff --git a/common/Membrane.cpp b/common/Membrane.cpp index 70219f46..98fe25a8 100644 --- a/common/Membrane.cpp +++ b/common/Membrane.cpp @@ -1351,46 +1351,11 @@ void Membrane::IonTransport(double *dist, double *den){ } // std::shared_ptr 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, diff --git a/common/Membrane.h b/common/Membrane.h index bf6aa0b7..234f4917 100644 --- a/common/Membrane.h +++ b/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 diff --git a/models/IonModel.cpp b/models/IonModel.cpp index 0781eb4e..4eaa480b 100644 --- a/models/IonModel.cpp +++ b/models/IonModel.cpp @@ -17,7 +17,6 @@ ScaLBL_IonModel::~ScaLBL_IonModel() {} void ScaLBL_IonModel::ReadParams(string filename, vector &num_iter) { - USE_MEMBRANE = true; // read the input database db = std::make_shared(filename); domain_db = db->getDatabase("Domain"); @@ -283,7 +282,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 = true; // read the input database db = std::make_shared(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("voxel_length"); } - + if (ion_db->keyExists("use_membrane")) { + USE_MEMBRANE = ion_db->getScalar("use_membrane"); + } // LB-Ion Model parameters //if (ion_db->keyExists( "timestepMax" )){ // timestepMax = ion_db->getScalar( "timestepMax" ); @@ -423,23 +424,104 @@ void ScaLBL_IonModel::ReadParams(string filename) { } } - - if (ion_db->keyExists("MembraneIonConcentrationList")) { - if (rank == 0) printf(".... Read MembraneIonConcentrationList \n"); - MembraneIonConcentration.clear(); - MembraneIonConcentration = ion_db->getVector("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] + if (USE_MEMBRANE){ + membrane_db = db->getDatabase("Membrane"); + + /* get membrane permeability parameters*/ + if (membrane_db->keyExists("MassFractionIn")) { + if (rank == 0) printf(".... Read membrane permeability (MassFractionIn) \n"); + MassFractionIn.clear(); + MassFractionIn = membrane_db->getVector("MassFractionIn"); + if (MassFractionIn.size() != number_ion_species) { + ERROR("Error: number_ion_species and membrane permeability (MassFractionIn) must be " + "the same length! \n"); + } + } + 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("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("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("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("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("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, "Na+"); - else { - IonMembrane->AssignCoefficients(dvcMap, Psi, "impermeable"); - } + IonMembrane->AssignCoefficients(dvcMap, Psi, ThresholdVoltage[ic],MassFractionIn[ic], + MassFractionOut[ic],ThresholdMassFractionIn[ic],ThresholdMassFractionOut[ic]); + timestep = 0; while (timestep < timestepMax[ic]) { //************************************************************************/ diff --git a/models/IonModel.h b/models/IonModel.h index 0bf115b9..f3f844ec 100644 --- a/models/IonModel.h +++ b/models/IonModel.h @@ -69,10 +69,13 @@ public: vector IonValence; vector IonConcentration; //unit [mol/m^3] vector MembraneIonConcentration; //unit [mol/m^3] - vector - Cin; //inlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec] - vector - Cout; //outlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec] + vector ThresholdVoltage; + vector MassFractionIn; + vector MassFractionOut; + vector ThresholdMassFractionIn; + vector ThresholdMassFractionOut; + vector Cin; //inlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec] + vector Cout; //outlet boundary value, can be either concentration [mol/m^3] or flux [mol/m^2/sec] vector tau; vector time_conv; @@ -83,7 +86,7 @@ public: std::shared_ptr Dm; // this domain is for analysis std::shared_ptr Mask; // this domain is for lbm std::shared_ptr ScaLBL_Comm; - // input database + // input databaseF std::shared_ptr db; std::shared_ptr domain_db; std::shared_ptr ion_db; diff --git a/tests/TestMembrane.cpp b/tests/TestMembrane.cpp index 31e63ca9..39df505a 100644 --- a/tests/TestMembrane.cpp +++ b/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); diff --git a/tests/lbpm_cell_simulator.cpp b/tests/lbpm_cell_simulator.cpp index 18c719c8..ce2e7100 100644 --- a/tests/lbpm_cell_simulator.cpp +++ b/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