Strain-dependent active stress model

scripts/ChamberSphere_StrainDepActStress.yaml

@@ -0,0 +1,83 @@
+variables:
+  - Pin
+  - Qin
+  - Pout
+  - Qout
+  - radius
+  - velo
+  - stress
+  - tau
+  - volume
+  - e_c
+  - tau_c
+  - k_c
+  - omega
+
+derivatives:
+  - dPin_dt
+  - dQin_dt
+  - dPout_dt
+  - dQout_dt
+  - dradius_dt
+  - dvelo_dt
+  - dstress_dt
+  - dtau_dt
+  - dvolume_dt
+  - de_c_dt
+  - dtau_c_dt
+  - dk_c_dt
+  - domega_dt
+
+constants:
+  - rho
+  - thick0
+  - radius0
+  - C0
+  - C1
+  - C2
+  - C3
+  - eta
+  - E_s
+  - mu
+  - alpha_r
+  - alpha
+  - k_0
+  - sigma_0
+  - m_0
+  - n_0
+  - u_plus
+  - u_minus
+
+time_dependent:
+  - u_plus 
+  - u_minus
+
+strain_dependent:
+  - m_0
+  - n_0
+
+helper_functions: |
+  def CG(radius):
+      return (1 + (radius / radius0)) ** 2
+
+  def dCG(radius, dradius_dt):
+      return 2 * (1 + (radius / radius0)) * (1 / radius0) * dradius_dt
+
+  def dW1(radius): 
+      return C0 * exp(C1 * (((1 + (radius / radius0))**(-4) + 2*(1 + (radius / radius0))**2) - 3)**2) * 2 * C1 * (((1 + (radius / radius0))**(-4) + 2*(1 + (radius / radius0))**2) - 3)
+
+  def dW4(radius):
+      return C2 * exp(C3 * ((((1 + (radius / radius0)) ** 2)-1)**2)) * 2 * C3 * (((1 + (radius / radius0)) ** 2) - 1)
+
+residuals:
+  - rho * thick0 * dvelo_dt + (thick0 / radius0) * (1 + (radius / radius0)) * stress - Pout * CG(radius)
+  - -stress + 4 * dW1(radius) *(1 - CG(radius)**(-3)) + 2 * dW4(radius)  + eta * dCG(radius,dradius_dt) * (1 + 2  * (CG(radius)**(-6))) + tau 
+  - 4 * pi * radius0 ** 2 * CG(radius) * velo - dvolume_dt
+  - dradius_dt - velo
+  - Qin - Qout - dvolume_dt
+  - Pin - Pout
+  - -tau + E_s * (radius/radius0 +0.5 * (radius / radius0) ** 2 - e_c) / ((1 + 2 * e_c) ** 2) 
+  - -alpha_r * domega_dt + m_0 - omega
+  - -tau_c - mu * de_c_dt + E_s * ((radius/radius0 +0.5 * (radius / radius0) ** 2 - e_c) * (1 + 2 * radius/radius0 + (radius/ radius0) ** 2)) / ((1 + 2 * e_c) ** 3)
+  - -dk_c_dt - (u_plus + omega * abs(u_minus) + alpha * abs(de_c_dt)) * k_c + n_0 * k_0 * u_plus
+  - -dtau_c_dt - (u_plus + omega * abs(u_minus) + alpha * abs(de_c_dt)) * tau_c + n_0 * sigma_0 * u_plus + k_c * de_c_dt

scripts/jacobian_strainDependence.py

@@ -0,0 +1,148 @@
+from sympy import symbols, Matrix, simplify, pi, Abs, exp
+from sympy.printing import ccode
+import re
+import pdb
+import yaml
+import argparse
+
+def load_model(filepath):
+    with open(filepath, 'r') as file:
+        data = yaml.safe_load(file)
+
+    variables = Matrix(symbols(' '.join(data['variables']), real=True))
+    derivatives = Matrix(symbols(' '.join(data['derivatives']), real=True))
+    constants = symbols(' '.join(data['constants']), real=True)
+
+    context = {str(s): s for s in list(variables) + list(derivatives) + list(constants)}
+    context['pi'] = pi
+    context['abs'] = Abs
+    context['exp'] = exp
+
+    if 'helper_functions' in data:
+        exec(data['helper_functions'], context, context)
+    residual_exprs = [eval(res, context, context) for res in data['residuals']]
+
+    residuals = Matrix(residual_exprs)
+    time_dependent = {context[s] for s in data.get('time_dependent', [])}
+    strain_dependent = {context[s] for s in data.get('strain_dependent', [])}
+
+    assert len(variables) == len(derivatives), f"Number of variables must be equal to number of derivatives"
+    assert len(variables) - 2 == len(residuals), f"Number of residuals must be number of unknowns minus 2"
+
+    return variables, derivatives, constants, residuals, time_dependent, strain_dependent
+
+def extract_linear(residuals, y, dy):
+    def is_constant_coeff(coeff):
+        return (coeff and coeff.free_symbols.isdisjoint(y.free_symbols | dy.free_symbols))
+    nr = residuals.shape[0]
+    ny = y.shape[0]
+    E = Matrix.zeros(nr, ny)
+    F = Matrix.zeros(nr, ny)
+    for i in range(nr):
+        for j in range(ny):
+            for mat, dd in zip([E, F], [dy, y]):
+                coeff = residuals[i].coeff(dd[j])
+                if is_constant_coeff(coeff):
+                    mat[i, j] = coeff
+    return E, F
+
+def extract_nonlinear(residuals, E, F, y, dy):
+    C = simplify(residuals - E * dy - F * y)
+    dC_dy = simplify(C.jacobian(y))
+    dC_dydot = simplify(C.jacobian(dy))
+    return C, dC_dy, dC_dydot
+
+def depends_on(expr, symbols_set):
+    return any(sym in expr.free_symbols for sym in symbols_set)
+
+def partition_terms(E, F, C, dC_dy, dC_dydot, time_dependent_symbols, strain_dependent_symbols):
+    parts = {"constant": [], "time": [], "solution": [], "strain": []}
+    for i in range(E.shape[0]):
+        parts["solution"].append(("C", i, C[i]))
+        for j in range(E.shape[1]):
+            for label, mat in [("E", E), ("F", F), ("dC_dy", dC_dy), ("dC_dydot", dC_dydot)]:
+                if label.startswith('d'):
+                    target = "solution"
+                elif depends_on(mat[i, j], time_dependent_symbols):  
+                    target = "time" 
+                elif depends_on(mat[i, j], strain_dependent_symbols):  
+                    target = "strain"   
+                else:
+                    target = "constant"
+                parts[target].append((label, i, j, mat[i, j]))
+    return parts
+
+def get_expressions(parts, type):
+    expressions = [part[-1] for part in parts]
+    deps = set().union(*(expr.free_symbols for expr in expressions))
+    return deps.intersection(type)
+
+def replace_symbolic_indices(expr, base):
+    pattern = re.compile(rf'\b{base}(\d+)\b')
+    return pattern.sub(rf'{base}[global_var_ids[\1]]', expr)
+
+def format_cpp_expr(expr):
+    try:
+        cpp_expr = ccode(expr).replace('3.141592653589793', 'M_PI')
+    except:
+        pdb.set_trace()
+    cpp_expr = replace_symbolic_indices(cpp_expr, 'y')
+    cpp_expr = replace_symbolic_indices(cpp_expr, 'dy')
+    return cpp_expr
+
+def print_index(i, j=None):
+    print(f".coeffRef(global_eqn_ids[{i}]", end='')
+    if j is not None:
+        print(f", global_var_ids[{j}]", end='')
+    print(f")", end='')
+
+def print_system(parts):
+    for part in parts:
+        if part[-1] != 0:
+            print(f"  system.{part[0]}", end='')
+            if len(part) == 3:
+                print_index(part[1])
+            elif len(part) == 4:
+                print_index(part[1], part[2])
+            print(f" = {format_cpp_expr(part[-1])};")
+
+def print_constants(parts, constants, time_dependent, strain_dependent):
+    for out in get_expressions(parts, constants):
+        if out in time_dependent:
+            print(f"  // compute time dependent constant {out}")
+        elif out in strain_dependent:
+            print(f"  // compute strain dependent constant {out}")
+        else:
+            print(f"  const double {out} = parameters[global_param_ids[ParamId::{out}]];")
+
+def print_variables(parts, y, dy):
+    for out in get_expressions(parts, y.free_symbols | dy.free_symbols):
+        for name, vec in zip(['y', 'dy'], [y, dy]):
+            if out in vec:
+                index = next(i for i, sym in enumerate(vec) if sym == out)
+                print(f"  const double {out} = {name}[global_var_ids[{index}]];")
+
+def main(yaml_path):
+    # read model from yaml file
+    y, dy, constants, residuals, time_dependent, strain_dependent = load_model(yaml_path)
+
+    # extract linear and nonlinear terms
+    E, F = extract_linear(residuals, y, dy)
+    C, dC_dy, dC_dydot = extract_nonlinear(residuals, E, F, y, dy)
+
+    # split into constant, time dependent and solution parts
+    parts = partition_terms(E, F, C, dC_dy, dC_dydot, time_dependent, strain_dependent) 
+
+    # print C++ code
+    for section in ['constant', 'time', 'strain', 'solution']:
+        print(f'update_{section}')
+        print_constants(parts[section], constants, time_dependent, strain_dependent)
+        print_variables(parts[section], y, dy)
+        print_system(parts[section])
+        print()
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description='Process YAML model file to generate C++ code for svZeroDSolver')
+    parser.add_argument('yaml_file', help='Path to YAML model file')
+    args = parser.parse_args()
+    main(args.yaml_file)

src/model/ActivationFunction.cpp

@@ -50,10 +50,13 @@ std::unique_ptr<ActivationFunction> ActivationFunction::create_default(
   if (type_str == "wrapping_cosine") {
     return std::make_unique<WrappingCosineActivation>(cardiac_period);
   }
+  if (type_str == "piecewise_rate") {
+    return std::make_unique<PiecewiseRateActivation>(cardiac_period);
+  }
   throw std::runtime_error(
       "Unknown activation_function type '" + type_str +
       "'. Must be one of: half_cosine, piecewise_cosine, two_hill, "
-      "double_tanh, fourier, wrapping_cosine");
+      "double_tanh, fourier, wrapping_cosine, piecewise_rate");
 }
 
 double HalfCosineActivation::compute(double time) {
@@ -251,3 +254,31 @@ double FourierActivation::compute(double time) {
   double t_in_cycle = std::fmod(time, cardiac_period_);
   return (compute_raw(t_in_cycle) - norm_min_) / norm_range_;
 }
+
+// ============================================================
+// PiecewiseRateActivation — activation rate signal u(t) for the
+// strain-dependent active-stress model (Caruel 2013)
+// ============================================================
+
+double PiecewiseRateActivation::compute(double time) {
+  const double t_in_cycle = std::fmod(time, cardiac_period_);
+  double u;
+  if (t_in_cycle <= 0.1) {
+    u = -3.2;
+  } else if (t_in_cycle > 0.1 && t_in_cycle <= 0.27) {
+    u = -3.2 + 38.2 / 0.17 * (t_in_cycle - 0.1);
+  } else if (t_in_cycle > 0.27 && t_in_cycle <= 0.32) {
+    u = 35;
+  } else if (t_in_cycle > 0.32 && t_in_cycle <= 0.34) {
+    u = 35 - (35 / 0.02) * (t_in_cycle - 0.32);
+  } else if (t_in_cycle > 0.34 && t_in_cycle <= 0.46) {
+    u = 24 - 12 / 0.12 * (t_in_cycle - 0.1);
+  } else if (t_in_cycle > 0.46 && t_in_cycle <= 0.59) {
+    u = -12;
+  } else if (t_in_cycle > 0.59 && t_in_cycle <= 0.6) {
+    u = -443.2 + 8.8 / 0.01 * (t_in_cycle - 0.1);
+  } else {
+    u = -3.2;
+  }
+  return u;
+}

src/model/ActivationFunction.h

@@ -60,7 +60,7 @@ class ActivationFunction {
    * @brief Create a default activation function from activation function type
    *
    * @param type_str One of: "half_cosine", "piecewise_cosine", "two_hill",
-   * "double_tanh", "wrapping_cosine", "fourier"
+   * "double_tanh", "wrapping_cosine", "fourier", "piecewise_rate"
    * @param cardiac_period Cardiac cycle period
    * @return Unique pointer to the created activation function
    */
@@ -304,4 +304,25 @@ class FourierActivation : public ActivationFunction {
   bool normalization_initialized_;
 };
 
+/**
+ * @brief Piecewise-linear activation rate activation signal
+ *
+ * Reproduces the activation rate signal \f$u(t)\f$ driving the
+ * \ref StrainDependentActiveStress model (Caruel et al. 2013). It is a piecewise-linear function of time within the
+ * cardiac cycle whose sign is meaningful (split into positive/negative parts
+ * by the strain-dependent active stress model).
+ */
+class PiecewiseRateActivation : public ActivationFunction {
+ public:
+  /**
+   * @brief Construct a new PiecewiseRateActivation object
+   *
+   * @param cardiac_period Cardiac cycle period
+   */
+  explicit PiecewiseRateActivation(double cardiac_period)
+      : ActivationFunction(cardiac_period, {}) {}
+
+  double compute(double time) override;
+};
+
 #endif  // SVZERODSOLVER_MODEL_ACTIVATIONFUNCTION_HPP_