Merge pull request #253 from ComputationalPhysiology/fix-Ta-bug

Fix Ta bug

demos/simple_demo.py

@@ -27,7 +27,7 @@
     geometry_schema_path=geometry_schema_path,
     loglevel=10,
     coupling_type="fully_coupled_ORdmm_Land",
-    T=20,
+    T=100,
 )
 
 # To see all different configuration options you can visit https://computationalphysiology.github.io/simcardems/api.html#module-simcardems.config

src/simcardems/mechanics_model.py

@@ -101,6 +101,9 @@ def __init__(self, *args, **kwargs):
         self.old_states = []
         self.old_controls = []
 
+    def _init_forms(self):
+        raise NotImplementedError
+
     def solve(self):
         self._init_forms()
         return super().solve()
@@ -151,6 +154,10 @@ def _init_spaces(self):
         )
         self._init_functions()
 
+    @property
+    def strong_coupling(self):
+        return hasattr(self.material.active, "Ta")
+
     @property
     def u_subspace_index(self) -> int:
         return 0
@@ -178,6 +185,13 @@ def _init_forms(self, init_solver: bool = True):
             self.state_test,
         )
 
+        if self.strong_coupling:
+            f0 = self.material.active.f0
+            f = self._F * f0
+            lmbda = dolfin.sqrt(f**2)
+            Pa = self.material.active.Ta(lmbda) * dolfin.outer(f, f0)
+            self._virtual_work += dolfin.inner(Pa, dolfin.grad(v)) * dx
+
         external_work = self._external_work(u, v)
         if external_work is not None:
             self._virtual_work += external_work
@@ -214,6 +228,25 @@ def solve(self):
         self._init_forms(init_solver=False)
         newton_iteration, newton_converged = self.solver.solve()
         getattr(self.solver, "check_overloads_called", None)
+
+        if self.strong_coupling:
+            u, p = self.state.split(deepcopy=True)
+
+            F = dolfin.grad(u) + dolfin.Identity(3)
+            f = F * self.material.active.f0
+            lmbda = dolfin.sqrt(f**2)
+
+            self.material.active._projector.project(self.material.active.lmbda, lmbda)
+            if self.material.active.dt > 0:
+                self.material.active._projector.project(
+                    self.material.active._dLambda,
+                    (lmbda - self.material.active.lmbda_prev) / self.material.active.dt,
+                )
+
+            self.material.active._projector.project(self.material.active.Ta_current, self.material.active.Ta(lmbda))
+            self.material.active.update_current(lmbda=lmbda)
+            self.material.active.update_prev()
+
         return newton_iteration, newton_converged
 
 
@@ -274,6 +307,12 @@ def _init_forms(self, init_solver: bool = True):
             self.state_test,
         )
 
+        f0 = self.material.active.f0
+        f = self._F * f0
+        lmbda = dolfin.sqrt(f**2)
+        Pa = self.material.active.Ta(lmbda) * dolfin.outer(f, f0)
+        self._virtual_work += dolfin.inner(Pa, dolfin.grad(v)) * dx
+
         self._jacobian = dolfin.derivative(
             self._virtual_work,
             self.state,

src/simcardems/models/fully_coupled_ORdmm_Land/active_model.py

@@ -3,7 +3,6 @@
 from typing import Optional
 
 import dolfin
-import numpy as np
 import pulse
 
 try:
@@ -90,13 +89,6 @@ def __init__(
         self._dLambda_tol = dLambda_tol
         self._t_prev = 0.0
 
-    @property
-    def dLambda(self):
-        logger.debug("Evaluate dLambda")
-        self._dLambda.vector()[:] = self.lmbda.vector() - self.lmbda_prev.vector()
-        self._dLambda.vector()[np.where(np.abs(self._dLambda.vector().get_local()) < self._dLambda_tol)[0]] = 0.0
-        return self._dLambda
-
     @property
     def Aw(self):
         Tot_A = self._parameters["Tot_A"]
@@ -131,93 +123,106 @@ def cs(self):
         scale_popu_rs = self._parameters["scale_popu_rs"]
         return kws * scale_popu_kws * phi * rw * scale_popu_rw * (1.0 - (rs * scale_popu_rs)) / (rs * scale_popu_rs)
 
-    def update_Zetas(self):
+    def register_time_stepper(self, time_stepper: TimeStepper) -> None:
+        self.time_stepper = time_stepper
+        self._t_prev = time_stepper.t
+
+    @property
+    def dt(self) -> float:
+        return TimeStepper.ns2ms(self.t - self._t_prev)
+
+    @property
+    def t(self) -> float:
+        if not hasattr(self, "time_stepper"):
+            return 0.0
+        return self.time_stepper.t
+
+    def dLambda(self, lmbda):
+        logger.debug("Evaluate dLambda")
+        if self.dt == 0:
+            return self._dLambda
+        else:
+            return (lmbda - self.lmbda_prev) / self.dt
+
+    def update_Zetas(self, lmbda):
         logger.debug("update Zetas")
-        self._Zetas.vector()[:] = _Zeta(
-            self.Zetas_prev.vector(),
+        self._projector(
+            self._Zetas,
+            _Zeta(
+                self.Zetas_prev,
+                self.As,
+                self.cs,
+                self.dLambda(lmbda),
+                self.dt,
+                self._scheme,
+            ),
+        )
+
+    def Zetas(self, lmbda):
+        # return self._Zetas
+        return _Zeta(
+            self.Zetas_prev,
             self.As,
             self.cs,
-            self.dLambda.vector(),
+            self.dLambda(lmbda),
             self.dt,
             self._scheme,
         )
 
-    @property
-    def Zetas(self):
-        return self._Zetas
-
-    def update_Zetaw(self):
+    def update_Zetaw(self, lmbda):
         logger.debug("update Zetaw")
-        self._Zetaw.vector()[:] = _Zeta(
-            self.Zetaw_prev.vector(),
+        self._projector(
+            self._Zetaw,
+            _Zeta(
+                self.Zetaw_prev,
+                self.Aw,
+                self.cw,
+                self.dLambda(lmbda),
+                self.dt,
+                self._scheme,
+            ),
+        )
+
+    def Zetaw(self, lmbda):
+        return _Zeta(
+            self.Zetaw_prev,
             self.Aw,
             self.cw,
-            self.dLambda.vector(),
+            self.dLambda(lmbda),
             self.dt,
             self._scheme,
         )
 
-    @property
-    def Zetaw(self):
-        return self._Zetaw
-
-    def register_time_stepper(self, time_stepper: TimeStepper) -> None:
-        self.time_stepper = time_stepper
-        self._t_prev = time_stepper.t
-
-    @property
-    def dt(self) -> float:
-        return TimeStepper.ns2ms(self.t - self._t_prev)
-
-    @property
-    def t(self) -> float:
-        if not hasattr(self, "time_stepper"):
-            return 0.0
-        return self.time_stepper.t
+    def update_current(self, lmbda):
+        self.update_Zetas(lmbda=lmbda)
+        self.update_Zetaw(lmbda=lmbda)
 
     def update_prev(self):
         logger.debug("update previous")
-        self.Zetas_prev.vector()[:] = self.Zetas.vector()
-        self.Zetaw_prev.vector()[:] = self.Zetaw.vector()
+        self.Zetas_prev.vector()[:] = self._Zetas.vector()
+        self.Zetaw_prev.vector()[:] = self._Zetaw.vector()
         self.lmbda_prev.vector()[:] = self.lmbda.vector()
-        self._projector.project(self.Ta_current, self.Ta)
+        # self.u_prev_prev.vector()[:] = self.u_prev.vector()
         self._t_prev = self.t
 
-    @property
-    def Ta(self):
+    def Ta(self, lmbda):
         logger.debug("Evaluate Ta")
         Tref = self._parameters["Tref"]
         rs = self._parameters["rs"]
-        scale_popu_Tref = self._parameters["scale_popu_Tref"]
-        scale_popu_rs = self._parameters["scale_popu_rs"]
+        scale_popu_Tref = 1.0  # self._parameters["scale_popu_Tref"]
+        scale_popu_rs = 1.0  # self._parameters["scale_popu_rs"]
         Beta0 = self._parameters["Beta0"]
 
         _min = ufl.min_value
         _max = ufl.max_value
-        if isinstance(self.lmbda, (int, float)):
+        if isinstance(lmbda, (int, float)):
             _min = min
             _max = max
-        lmbda = _min(1.2, self.lmbda)
+        lmbda = _min(1.2, lmbda)
         h_lambda_prima = 1.0 + Beta0 * (lmbda + _min(lmbda, 0.87) - 1.87)
         h_lambda = _max(0, h_lambda_prima)
 
-        return (
-            h_lambda
-            * (Tref * scale_popu_Tref / (rs * scale_popu_rs))
-            * (self.XS * (self.Zetas + 1.0) + self.XW * self.Zetaw)
-        )
+        Zetas = self.Zetas(lmbda)
+        Zetaw = self.Zetaw(lmbda)
 
-    def Wactive(self, F, **kwargs):
-        """Active stress energy"""
-        logger.debug("Compute active stress energy")
-        C = F.T * F
-        f = F * self.f0
-        self._projector.project(self.lmbda, dolfin.sqrt(f**2))
-        self.update_Zetas()
-        self.update_Zetaw()
-        return pulse.material.active_model.Wactive_transversally(
-            Ta=self.Ta,
-            C=C,
-            f0=self.f0,
-            eta=self.eta,
-        )
+        return h_lambda * (Tref * scale_popu_Tref / (rs * scale_popu_rs)) * (self.XS * (Zetas + 1.0) + self.XW * Zetaw)