add solver files

src/knpemidg/__init__.py

@@ -10,7 +10,7 @@
 from knpemidg.utils import pcws_constant_project
 from knpemidg.utils import CellCenterDistance
 from knpemidg.solver import Solver
-#from knpemidg.solver_emi import SolverEMI
+from knpemidg.solver_emi import SolverEMI
 
 __all__ = ["Solver", "MembraneModel", "subdomain_marking_foo",
         "interface_normal", "plus", "minus", "pcws_constant_project",

src/knpemidg/solver.py

@@ -43,21 +43,20 @@ class bcolors:
 #      (step III: solve ODEs at interface, and update membrane potential)
 #
 # Membrane potential is defined as phi_i - phi_e, since we have
-# marked cell in ECS with 2 and cells in ICS with 1 we have an
-# interface normal pointing inwards
+# marked cell in ECS with 0 and cells in ICS with 1 we have an
+# interface normal pointing inwards (from lower to higher)
 #    ____________________
 #   |                    |
 #   |      ________      |
 #   |     |        |     |
 #   | ECS |   ICS  |     |
-#   |  2  |->  1   |     |
+#   |  0  |->  1   |     |
 #   | (+) |   (-)  |     |
 #   |     |________|     |
 #   |                    |
 #   |____________________|
 #
-# Normal will always point from higher to lower (e.g. from 2 -> 1)
-#
+# Normal will always point from lower to higher (e.g. from 0 -> 1)
 # NB! The code assumes that all interior facets are tagged with 0.
 
 class Solver:
@@ -220,7 +219,7 @@ def setup_membrane_model(self, stim_params, odes):
             ode_model.set_parameter_values({'Cm': lambda x: self.params.C_M})
 
             # dictionary for ion specific currents (i.e src terms PDEs)
-            I_ch_k = {} 
+            I_ch_k = {}
             # Initialize src terms for PDEs
             for ion in self.ion_list:
                 # function for src term pde
@@ -915,7 +914,7 @@ def solve_system_passive(self, Tstop, t, solver_params, membrane_params, filenam
         self.atol_knp = solver_params.atol_knp           # absolute tolerance knp
         self.threshold_knp = solver_params.threshold_knp # threshold knp
 
-        self.splitting_scheme = False               # no splitting scheme
+        self.splitting_scheme = False                    # no splitting scheme
 
         # Set filename for saving results
         self.filename = filename
@@ -940,22 +939,6 @@ def solve_system_passive(self, Tstop, t, solver_params, membrane_params, filenam
             # Start timer (ODE solve)
             ts = time.perf_counter()
 
-            if self.mms is None:
-                # if not mms, get src terms for PDEs from passive model
-                for mem_model in self.mem_models:
-
-                    ode_model = mem_model['ode']
-
-                    # Update membrane potential and Nernst potential in membrane model
-                    ode_model.set_membrane_potential(self.phi_M_prev_PDE)
-                    ode_model.set_parameter('E_K', self.ion_list[0]['E'])
-                    ode_model.set_parameter('E_Na', self.ion_list[2]['E'])
-
-                    # Get src terms for next PDE step
-                    for ion, I_ch_k in mem_model['I_ch_k'].items():
-                        # update src term for each ion species
-                        ode_model.get_parameter("I_ch_" + ion, I_ch_k)
-
             # End timer (ODE solve)
             te = time.perf_counter()
             res = te - ts
@@ -1036,6 +1019,14 @@ def solve_system_active(self, Tstop, t, solver_params, filename=None):
                 ode_model.set_parameter('E_K', self.ion_list[0]['E'])
                 ode_model.set_parameter('E_Na', self.ion_list[2]['E'])
 
+                # set extracellular trace of K concentration at membrane
+                K_e = plus(self.c_prev_k.split()[0], self.n_g)
+                ode_model.set_parameter('K_e', pcws_constant_project(K_e, self.Q))
+
+                # set intracellular trace of Na concentration at membrane
+                Na_i = minus(self.ion_list[-1]['c'], self.n_g)
+                ode_model.set_parameter('Na_i', pcws_constant_project(Na_i, self.Q))
+
                 # Solve ODEs
                 ode_model.step_lsoda(dt=dt_ode, \
                     stimulus=stimulus, stimulus_locator=stimulus_locator)
@@ -1127,10 +1118,6 @@ def initialize_solver_savefile(self, path_timings):
         self.f_Na = File('results/active_knp/Na.pvd')
         self.f_K = File('results/active_knp/K.pvd')
         self.f_Cl = File('results/active_knp/Cl.pvd')
-
-        self.f_grad_Na = File('results/active_knp/Na_grad.pvd')
-        self.f_grad_K = File('results/active_knp/K_grad.pvd')
-
         self.f_kappa = File('results/active_knp/Kappa.pvd')
 
         return
@@ -1147,17 +1134,11 @@ def save_solver(self, k):
             Cl_ = project(self.c.split()[1], VDG0)
             Na_ = project(self.ion_list[-1]['c'], VDG0)
 
-            grad_Na_ = project(grad(self.c.split()[0]), VDG1)
-            grad_K_ = project(grad(self.c.split()[1]), VDG1)
-            self.f_pot_grad << (phi_, k)
             self.f_pot << (phi, k)
             self.f_Na << (Na_, k)
             self.f_K << (K_, k)
             self.f_Cl << (Cl_, k)
 
-            self.f_grad_Na << (grad_Na_, k)
-            self.f_grad_K << (grad_K_, k)
-
             self.f_kappa << (project(self.kappa, VDG0), k)
 
             return