Use Guccione passive law in the ChamberSphere block

Replace the neo-Hookean passive stress (4*gamma_W1*(1-CG^-3)) with the Guccione
spherical wall stress derived for a thin-walled incompressible sphere, so the
forward 0D model and the Yale calibration use the same material:
  S_pas = C exp(Q) [ (b_f+b_t)/2 lam^2 E_p - b_t lam^-4 E_r ]
Rename gamma_W1 -> guccione_C and add the dimensionless exponents b_f, b_t as
input parameters (b_fs drops out: no shear in the equibiaxial deformation).
Regenerate update_solution/update_gradient from scripts/ChamberSphere.yaml.

The synthetic full-state identifiability test recovers all 11 parameters
(including guccione_C, b_f, b_t) to machine precision; all calibrator tests pass.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: mrp089

scripts/ChamberSphere.yaml

@@ -23,16 +23,19 @@ derivatives:
 # below, which are exactly the quantities provided in the solver input file.
 #   stress  here represents  gamma * stress
 #   tau     here represents  gamma * tau
+#   guccione_C represents    gamma * C  (the Guccione passive scaling)
 constants:
   - volume0
-  - gamma_W1
+  - guccione_C # gamma * C, scaling in front of the Guccione exponential
   - gamma_sigma_max # gamma * sigma_max
   - prestress
   - alpha_max
   - alpha_min
   - tsys
   - tdias
   - steepness
+  - b_f # Guccione fiber exponent (dimensionless)
+  - b_t # Guccione transverse exponent (dimensionless)
 
 # In-cycle time of the activation. Not a state variable and not a calibration
 # parameter (no Jacobian column), but the activation depends on it. The solver
@@ -43,7 +46,7 @@ constants:
 time_symbol: t
 
 helper_functions: |
-  from sympy import Abs, Max, tanh
+  from sympy import Abs, Max, tanh, exp
 
   def V(volume):
       return volume + volume0
@@ -56,6 +59,27 @@ helper_functions: |
       stretch_val = stretch(volume)
       return (stretch_val) ** 2
 
+  # Guccione passive spherical wall stress for a thin-walled incompressible
+  # sphere. The in-plane stretch is lam = stretch, the radial stretch lam^-2, so
+  # the Green-Lagrange strains are E_p = (lam^2 - 1)/2 (in-plane fiber/sheet) and
+  # E_r = (lam^-4 - 1)/2 (radial). With no shear, Q = (b_f + b_t) E_p^2 + b_t E_r^2
+  # and the mean in-plane Cauchy stress (sigma_radial = 0 for the thin wall) is
+  #   S = C exp(Q) [ (b_f + b_t)/2 lam^2 E_p - b_t lam^-4 E_r ].
+  # b_fs drops out: it only multiplies shear terms, which vanish here.
+  def E_p(volume):
+      return 0.5 * (CG(volume) - 1)
+
+  def E_r(volume):
+      return 0.5 * (CG(volume) ** (-2) - 1)
+
+  def Q_guc(volume):
+      return (b_f + b_t) * E_p(volume) ** 2 + b_t * E_r(volume) ** 2
+
+  def passive(volume):
+      cg = CG(volume)
+      return guccione_C * exp(Q_guc(volume)) * (
+          0.5 * (b_f + b_t) * cg * E_p(volume) - b_t * cg ** (-2) * E_r(volume))
+
   # Active contraction activation, identical to
   # ChamberSphere::get_elastance_values (t is the in-cycle time).
   S_plus = 0.5 * (1.0 + tanh((t - tsys) / steepness))
@@ -67,7 +91,7 @@ helper_functions: |
 
 residuals:
   - stretch(volume) * stress - Pout * CG(volume)
-  - -stress + tau + 4 * gamma_W1 * (1 - CG(volume) ** (-3)) + prestress
+  - -stress + tau + passive(volume) + prestress
   - dtau_dt + act * tau - gamma_sigma_max * act_plus
   - Qin - Qout - dvolume_dt
   - Pin - Pout

scripts/make_chamber_sphere_calibration.py

@@ -35,15 +35,26 @@
 # meaningful). Order matches the ParamId enum in ChamberSphere.h.
 TRUE = {
     "volume0": 1.0e-4,
-    "gamma_W1": 1.0e4,
+    "guccione_C": 1.0e3,
     "gamma_sigma_max": 1.85e5,
     "prestress": 1.0e3,
     "alpha_max": 30.0,
     "alpha_min": -30.0,
     "tsys": 0.17,
     "tdias": 0.484,
     "steepness": 0.05,
+    "b_f": 8.0,
+    "b_t": 3.0,
 }
+
+
+def passive_guccione(lam, C, b_f, b_t):
+    """Guccione passive spherical wall stress (matches ChamberSphere)."""
+    Ep = 0.5 * (lam ** 2 - 1.0)
+    Er = 0.5 * (lam ** (-4) - 1.0)
+    Q = (b_f + b_t) * Ep ** 2 + b_t * Er ** 2
+    return C * np.exp(Q) * (
+        0.5 * (b_f + b_t) * lam ** 2 * Ep - b_t * lam ** (-4) * Er)
 PERIOD = 1.0
 NUM_OBS = 200
 PERTURB = 1.20  # 20% perturbed start
@@ -98,7 +109,8 @@ def main():
     # Reconstruct the consistent full state (every residual vanishes exactly).
     act, act_plus = activation(t, p)
     dtau_dt = -act * tau + p["gamma_sigma_max"] * act_plus          # residual 2
-    stress = tau + 4.0 * p["gamma_W1"] * (1.0 - CG ** (-3)) + p["prestress"]  # residual 1
+    stress = tau + passive_guccione(stretch, p["guccione_C"], p["b_f"], p["b_t"]) \
+        + p["prestress"]                                            # residual 1
     Pout = stress / stretch                                          # residual 0
     Pin = Pout                                                       # residual 4
 
@@ -179,7 +191,7 @@ def main():
         rel = abs(calibrated[k] - TRUE[k]) / abs(TRUE[k])
         ok = ok and rel < 1e-6
         print(f"{k:<18}{TRUE[k]:>14.6g}{start[k]:>14.6g}{calibrated[k]:>16.8g}{rel:>12.2e}")
-    print(f"\nAll 9 parameters recovered to rtol<1e-6: {ok}")
+    print(f"\nAll {len(names)} parameters recovered to rtol<1e-6: {ok}")
     return 0 if ok else 1
 
 

src/model/ChamberSphere.cpp

@@ -38,18 +38,20 @@ void ChamberSphere::update_solution(
     SparseSystem& system, std::vector<double>& parameters,
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
-  const double prestress = parameters[global_param_ids[ParamId::prestress]];
+  const double b_t = parameters[global_param_ids[ParamId::b_t]];
+  const double b_f = parameters[global_param_ids[ParamId::b_f]];
   const double volume0 = parameters[global_param_ids[ParamId::volume0]];
-  const double gamma_W1 = parameters[global_param_ids[ParamId::gamma_W1]];
+  const double guccione_C = parameters[global_param_ids[ParamId::guccione_C]];
+  const double prestress = parameters[global_param_ids[ParamId::prestress]];
   const double Pout = y[global_var_ids[2]];
-  const double stress = y[global_var_ids[4]];
   const double volume = y[global_var_ids[6]];
+  const double stress = y[global_var_ids[4]];
   system.C.coeffRef(global_eqn_ids[0]) = -Pout*pow((volume + volume0)/volume0, 0.66666666666666663) + Pout + stress*pow((volume + volume0)/volume0, 0.33333333333333331) - stress;
-  system.C.coeffRef(global_eqn_ids[1]) = -4*gamma_W1*pow(volume/volume0 + 1, -2.0) + 4*gamma_W1 + prestress;
+  system.C.coeffRef(global_eqn_ids[1]) = pow((volume + volume0)/volume0, -2.6666666666666665)*(-guccione_C*(0.5*b_t*(1 - pow((volume + volume0)/volume0, 1.3333333333333333)) - 0.25*pow((volume + volume0)/volume0, 3.333333333333333)*(b_f + b_t)*(pow((volume + volume0)/volume0, 0.66666666666666663) - 1))*exp(0.25*pow((volume + volume0)/volume0, -2.6666666666666665)*(b_t*pow(1 - pow((volume + volume0)/volume0, 1.3333333333333333), 2) + pow((volume + volume0)/volume0, 2.6666666666666665)*(b_f + b_t)*pow(pow((volume + volume0)/volume0, 0.66666666666666663) - 1, 2))) + prestress*pow((volume + volume0)/volume0, 2.6666666666666665));
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[2]) = 1 - pow((volume + volume0)/volume0, 0.66666666666666663);
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[4]) = pow((volume + volume0)/volume0, 0.33333333333333331) - 1;
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[6]) = (-0.66666666666666663*Pout*pow((volume + volume0)/volume0, 0.66666666666666663) + 0.33333333333333331*stress*pow((volume + volume0)/volume0, 0.33333333333333331))/(volume + volume0);
-  system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[6]) = 8.0*gamma_W1*pow((volume + volume0)/volume0, -3.0)/volume0;
+  system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[6]) = pow((volume + volume0)/volume0, -8.0)*(pow((volume + volume0)/volume0, 2.6666666666666665)*(guccione_C*(0.5*b_t*(pow((volume + volume0)/volume0, 1.3333333333333333) - 1) + 0.25*pow((volume + volume0)/volume0, 3.333333333333333)*(b_f + b_t)*(pow((volume + volume0)/volume0, 0.66666666666666663) - 1))*(0.66666666666666663*b_t*pow((volume + volume0)/volume0, 1.3333333333333333)*(pow((volume + volume0)/volume0, 1.3333333333333333) - 1) - 0.66666666666666663*b_t*pow(pow((volume + volume0)/volume0, 1.3333333333333333) - 1, 2) + 0.33333333333333331*pow((volume + volume0)/volume0, 3.333333333333333)*(b_f + b_t)*(pow((volume + volume0)/volume0, 0.66666666666666663) - 1))*exp(0.25*pow((volume + volume0)/volume0, -2.6666666666666665)*(b_t*pow(pow((volume + volume0)/volume0, 1.3333333333333333) - 1, 2) + pow((volume + volume0)/volume0, 2.6666666666666665)*(b_f + b_t)*pow(pow((volume + volume0)/volume0, 0.66666666666666663) - 1, 2))) + pow((volume + volume0)/volume0, 2.6666666666666665)*(guccione_C*(0.66666666666666663*b_t*pow((volume + volume0)/volume0, 1.3333333333333333) + 0.83333333333333326*pow((volume + volume0)/volume0, 3.333333333333333)*(b_f + b_t)*(pow((volume + volume0)/volume0, 0.66666666666666663) - 1) + 0.16666666666666666*pow((volume + volume0)/volume0, 3.9999999999999996)*(b_f + b_t))*exp(0.25*pow((volume + volume0)/volume0, -2.6666666666666665)*(b_t*pow(pow((volume + volume0)/volume0, 1.3333333333333333) - 1, 2) + pow((volume + volume0)/volume0, 2.6666666666666665)*(b_f + b_t)*pow(pow((volume + volume0)/volume0, 0.66666666666666663) - 1, 2))) + 2.6666666666666665*prestress*pow((volume + volume0)/volume0, 2.6666666666666665))) - 2.6666666666666665*pow((volume + volume0)/volume0, 5.333333333333333)*(guccione_C*(0.5*b_t*(pow((volume + volume0)/volume0, 1.3333333333333333) - 1) + 0.25*pow((volume + volume0)/volume0, 3.333333333333333)*(b_f + b_t)*(pow((volume + volume0)/volume0, 0.66666666666666663) - 1))*exp(0.25*pow((volume + volume0)/volume0, -2.6666666666666665)*(b_t*pow(pow((volume + volume0)/volume0, 1.3333333333333333) - 1, 2) + pow((volume + volume0)/volume0, 2.6666666666666665)*(b_f + b_t)*pow(pow((volume + volume0)/volume0, 0.66666666666666663) - 1, 2))) + prestress*pow((volume + volume0)/volume0, 2.6666666666666665)))/(volume + volume0);
 }
 
 void ChamberSphere::get_elastance_values(std::vector<double>& parameters) {
@@ -82,22 +84,18 @@ void ChamberSphere::update_gradient(
     Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
     Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
     std::vector<double>& dy) {
-  // Calibrates the chamber parameters from a full-state observation set. The
-  // time-independent parameters (volume0, gamma_W1, prestress) appear in the
-  // momentum and spherical-stress equations; the active-stress equation (eq 2)
-  // and its parameters (gamma_sigma_max and the activation timing) depend on the
-  // observation time, which the optimizer supplies via model->time when a time
-  // vector is given. Shared subexpressions are hoisted into temporaries (xN).
-  // Generated by scripts/jacobian.py from scripts/ChamberSphere.yaml.
+  // Guccione passive law; generated by scripts/jacobian.py from scripts/ChamberSphere.yaml.
   const double volume0 = alpha[global_param_ids[ParamId::volume0]];
-  const double gamma_W1 = alpha[global_param_ids[ParamId::gamma_W1]];
+  const double guccione_C = alpha[global_param_ids[ParamId::guccione_C]];
   const double gamma_sigma_max = alpha[global_param_ids[ParamId::gamma_sigma_max]];
   const double prestress = alpha[global_param_ids[ParamId::prestress]];
   const double alpha_max = alpha[global_param_ids[ParamId::alpha_max]];
   const double alpha_min = alpha[global_param_ids[ParamId::alpha_min]];
   const double tsys = alpha[global_param_ids[ParamId::tsys]];
   const double tdias = alpha[global_param_ids[ParamId::tdias]];
   const double steepness = alpha[global_param_ids[ParamId::steepness]];
+  const double b_f = alpha[global_param_ids[ParamId::b_f]];
+  const double b_t = alpha[global_param_ids[ParamId::b_t]];
   const double Pin = y[global_var_ids[0]];
   const double Qin = y[global_var_ids[1]];
   const double Pout = y[global_var_ids[2]];
@@ -112,54 +110,76 @@ void ChamberSphere::update_gradient(
   const double x0 = 1.0/volume0;
   const double x1 = volume + volume0;
   const double x2 = x0*x1;
-  const double x3 = Pout*pow(x2, 0.66666666666666663);
-  const double x4 = pow(x2, 0.33333333333333331);
-  const double x5 = pow(x2, -2.0);
-  const double x6 = 1.0/steepness;
-  const double x7 = t - tdias;
-  const double x8 = tanh(x6*x7);
-  const double x9 = 0.5*x8 - 0.5;
-  const double x10 = t - tsys;
-  const double x11 = tanh(x10*x6);
-  const double x12 = 0.5*x11 + 0.5;
-  const double x13 = x12*x9;
-  const double x14 = x13 + 1;
-  const double x15 = alpha_max*x13 - alpha_min*x14;
-  const double x16 = -x9;
-  const double x17 = x12*x16;
-  const double x18 = 1 - x17;
-  const double x19 = alpha_max*x17 + alpha_min*x18;
-  const double x20 = fmax(0, x19);
-  const double x21 = volume*x0;
-  const double x22 = x21/x1;
-  const double x23 = gamma_sigma_max*(x19 > 0 ? 1.0 : 0.0);
-  const double x24 = (((x15) > 0) - ((x15) < 0));
-  const double x25 = 1 - pow(x11, 2);
-  const double x26 = 0.5*x6;
-  const double x27 = x25*x26;
-  const double x28 = 1 - pow(x8, 2);
-  const double x29 = x12*x28;
-  const double x30 = alpha_max*x26*x29 - 0.5*alpha_min*x12*x28*x6;
-  const double x31 = 0.5/pow(steepness, 2);
-  const double x32 = x10*x25*x31;
-  const double x33 = x16*x32;
-  const double x34 = x29*x31*x7;
-  const double x35 = alpha_max*x34;
-  const double x36 = x32*x9;
-
-  residual(global_eqn_ids[0]) = stress*x4 - x3;
-  residual(global_eqn_ids[1]) = 4*gamma_W1*(1 - x5) + prestress - stress + tau;
-  residual(global_eqn_ids[2]) = dtau_dt - gamma_sigma_max*x20 + tau*fabs(x15);
+  const double x3 = pow(x2, 0.66666666666666663);
+  const double x4 = Pout*x3;
+  const double x5 = pow(x2, 0.33333333333333331);
+  const double x6 = pow(x2, 1.3333333333333333);
+  const double x7 = 1.0/x6;
+  const double x8 = 0.5*x7 - 0.5;
+  const double x9 = b_f + b_t;
+  const double x10 = 0.5*x3 - 0.5;
+  const double x11 = 0.5*b_t;
+  const double x12 = 1.0/steepness;
+  const double x13 = t - tdias;
+  const double x14 = tanh(x12*x13);
+  const double x15 = 0.5*x14 - 0.5;
+  const double x16 = t - tsys;
+  const double x17 = tanh(x12*x16);
+  const double x18 = 0.5*x17 + 0.5;
+  const double x19 = x15*x18;
+  const double x20 = x19 + 1;
+  const double x21 = alpha_max*x19 - alpha_min*x20;
+  const double x22 = -x15;
+  const double x23 = x18*x22;
+  const double x24 = 1 - x23;
+  const double x25 = alpha_max*x23 + alpha_min*x24;
+  const double x26 = fmax(0, x25);
+  const double x27 = volume*x0/x1;
+  const double x28 = pow(x2, 5.333333333333333);
+  const double x29 = 0.66666666666666663*b_t;
+  const double x30 = pow(x2, 2.6666666666666665);
+  const double x31 = x3 - 1;
+  const double x32 = 1 - x6;
+  const double x33 = pow(x2, 3.333333333333333);
+  const double x34 = x11*x32 - 0.25*x31*x33*x9;
+  const double x35 = x30*x34;
+  const double x36 = pow(x32, 2);
+  const double x37 = x30*pow(x31, 2);
+  const double x38 = 1.0/x30;
+  const double x39 = 0.25*x38;
+  const double x40 = exp(x39*(b_t*x36 + x37*x9));
+  const double x41 = guccione_C*x40;
+  const double x42 = pow(x2, -8.0)*x41;
+  const double x43 = gamma_sigma_max*(x25 > 0 ? 1.0 : 0.0);
+  const double x44 = (((x21) > 0) - ((x21) < 0));
+  const double x45 = 1 - pow(x17, 2);
+  const double x46 = 0.5*x12;
+  const double x47 = x45*x46;
+  const double x48 = 1 - pow(x14, 2);
+  const double x49 = x18*x48;
+  const double x50 = alpha_max*x46*x49 - 0.5*alpha_min*x12*x18*x48;
+  const double x51 = 0.5/pow(steepness, 2);
+  const double x52 = x16*x45*x51;
+  const double x53 = x22*x52;
+  const double x54 = x13*x49*x51;
+  const double x55 = alpha_max*x54;
+  const double x56 = x15*x52;
+
+  residual(global_eqn_ids[0]) = stress*x5 - x4;
+  residual(global_eqn_ids[1]) = guccione_C*(-b_t*x7*x8 + x10*x3*(0.5*b_f + x11))*exp(b_t*pow(x8, 2) + pow(x10, 2)*x9) + prestress - stress + tau;
+  residual(global_eqn_ids[2]) = dtau_dt - gamma_sigma_max*x26 + tau*fabs(x21);
   residual(global_eqn_ids[3]) = Qin - Qout - dvolume_dt;
   residual(global_eqn_ids[4]) = Pin - Pout;
-  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[ParamId::volume0]) = x22*(-0.33333333333333331*stress*x4 + 0.66666666666666663*x3);
-  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::volume0]) = -8.0*gamma_W1*x22*x5;
-  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::gamma_W1]) = 4 - 4*pow(x21 + 1, -2.0);
-  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::gamma_sigma_max]) = -x20;
+  jacobian.coeffRef(global_eqn_ids[0], global_param_ids[ParamId::volume0]) = x27*(-0.33333333333333331*stress*x5 + 0.66666666666666663*x4);
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::volume0]) = x27*x42*(x28*(x29*(x6 - 2) + 0.16666666666666666*x30*x9*(-x3*x31 - x6)) + x35*(-x29*x32 + 0.33333333333333331*x31*x33*x9));
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::guccione_C]) = -x34*x38*x40;
+  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::gamma_sigma_max]) = -x26;
   jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::prestress]) = 1;
-  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::alpha_max]) = tau*x12*x24*x9 - x17*x23;
-  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::alpha_min]) = -tau*x14*x24 - x18*x23;
-  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::tsys]) = tau*x24*(-alpha_max*x27*x9 + 0.5*alpha_min*x25*x6*x9) - x23*(-alpha_max*x16*x27 + 0.5*alpha_min*x16*x25*x6);
-  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::tdias]) = -tau*x24*x30 - x23*x30;
-  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::steepness]) = tau*x24*(-alpha_max*x36 - alpha_min*(-x34 - x36) - x35) - x23*(-alpha_max*x33 + alpha_min*(x33 - x34) + x35);
+  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::alpha_max]) = tau*x15*x18*x44 - x23*x43;
+  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::alpha_min]) = -tau*x20*x44 - x24*x43;
+  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::tsys]) = tau*x44*(-alpha_max*x15*x47 + 0.5*alpha_min*x12*x15*x45) - x43*(-alpha_max*x22*x47 + 0.5*alpha_min*x12*x22*x45);
+  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::tdias]) = -tau*x44*x50 - x43*x50;
+  jacobian.coeffRef(global_eqn_ids[2], global_param_ids[ParamId::steepness]) = tau*x44*(-alpha_max*x56 - alpha_min*(-x54 - x56) - x55) - x43*(-alpha_max*x53 + alpha_min*(x53 - x54) + x55);
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::b_f]) = x31*x39*x41*(-x31*x34 + x33);
+  jacobian.coeffRef(global_eqn_ids[1], global_param_ids[ParamId::b_t]) = x42*(x28*(x33*(0.25*x3 - 0.25) + 0.5*x6 - 0.5) - 0.25*x35*(x36 + x37));
 }