Reuse forward-solver assembly for calibrator residual

src/model/Block.cpp

@@ -62,11 +62,10 @@ void Block::update_solution(
 void Block::post_solve(Eigen::Matrix<double, Eigen::Dynamic, 1>& y) {}
 
 void Block::update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                            Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                             Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                             std::vector<double>& y, std::vector<double>& dy) {
-  throw std::runtime_error("Gradient calculation not implemented for block " +
-                           get_name());
+  throw std::runtime_error(
+      "Calibrator parameter Jacobian not implemented for block " + get_name());
 }
 
 TripletsContributions Block::get_num_triplets() { return num_triplets; }

src/model/Block.h

@@ -230,6 +230,19 @@ class Block {
   virtual void update_time(SparseSystem& system,
                            std::vector<double>& parameters);
 
+  /**
+   * @brief Whether this block contributes time-dependent terms to the
+   * assembly.
+   *
+   * The calibrator does not support models containing time-dependent blocks
+   * because the observations passed to it carry no time stamp. Blocks that
+   * implement non-trivial ``update_time`` logic must override this to return
+   * true; the calibrator throws on setup if any such block is present.
+   *
+   * @return true if ``update_time`` writes time-dependent contributions
+   */
+  virtual bool has_time_dependent_assembly() const { return false; }
+
   /**
    * @brief Update the solution-dependent contributions of the element in a
    * sparse system
@@ -252,20 +265,22 @@ class Block {
   virtual void post_solve(Eigen::Matrix<double, Eigen::Dynamic, 1>& y);
 
   /**
-   * @brief Set the gradient of the block contributions with respect to the
-   * parameters
+   * @brief Write this block's contributions to the parameter Jacobian
+   * \f$\partial r / \partial \alpha\f$.
+   *
+   * The residual is recomputed at the top level of the calibrator from the
+   * solver assembly (E, F, c via update_constant / update_time /
+   * update_solution); each block only needs to populate the columns of the
+   * sparse Jacobian that correspond to its own parameters.
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
-   * @param y Current solution
-   * @param dy Time-derivative of the current solution
+   * @param y Observed solution at the current data point
+   * @param dy Observed time derivative at the current data point
    */
-  virtual void update_gradient(
-      Eigen::SparseMatrix<double>& jacobian,
-      Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-      Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
-      std::vector<double>& dy);
+  virtual void update_gradient(Eigen::SparseMatrix<double>& jacobian,
+                               Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
+                               std::vector<double>& y, std::vector<double>& dy);
 
   /**
    * @brief Number of triplets of element

src/model/BloodVessel.cpp

@@ -58,47 +58,32 @@ void BloodVessel::update_solution(
 
 void BloodVessel::update_gradient(
     Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
     Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
     std::vector<double>& dy) {
-  auto y0 = y[global_var_ids[0]];
   auto y1 = y[global_var_ids[1]];
-  auto y2 = y[global_var_ids[2]];
-  auto y3 = y[global_var_ids[3]];
-
   auto dy0 = dy[global_var_ids[0]];
   auto dy1 = dy[global_var_ids[1]];
   auto dy3 = dy[global_var_ids[3]];
 
   auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
   auto capacitance = alpha[global_param_ids[ParamId::CAPACITANCE]];
-  auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
   double stenosis_coeff = 0.0;
-
   if (global_param_ids.size() > 3) {
     stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
   }
   auto stenosis_resistance = stenosis_coeff * fabs(y1);
 
   jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y1;
   jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -dy3;
-
   if (global_param_ids.size() > 3) {
     jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = -fabs(y1) * y1;
   }
 
   jacobian.coeffRef(global_eqn_ids[1], global_param_ids[0]) = capacitance * dy1;
   jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) =
       -dy0 + (resistance + 2 * stenosis_resistance) * dy1;
-
   if (global_param_ids.size() > 3) {
     jacobian.coeffRef(global_eqn_ids[1], global_param_ids[3]) =
         2.0 * capacitance * fabs(y1) * dy1;
   }
-
-  residual(global_eqn_ids[0]) =
-      y0 - (resistance + stenosis_resistance) * y1 - y2 - inductance * dy3;
-  residual(global_eqn_ids[1]) =
-      y1 - y3 - capacitance * dy0 +
-      capacitance * (resistance + 2.0 * stenosis_resistance) * dy1;
 }

src/model/BloodVessel.h

@@ -195,17 +195,14 @@ class BloodVessel : public Block {
                        const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy);
 
   /**
-   * @brief Set the gradient of the block contributions with respect to the
-   * parameters
+   * @brief Write this block's columns of the parameter Jacobian.
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
-   * @param y Current solution
-   * @param dy Time-derivative of the current solution
+   * @param y Observed solution
+   * @param dy Observed time derivative of the solution
    */
   void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                        Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                        std::vector<double>& y, std::vector<double>& dy);
 

src/model/BloodVesselCRL.cpp

@@ -50,38 +50,16 @@ void BloodVesselCRL::update_solution(
 
 void BloodVesselCRL::update_gradient(
     Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
     Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
     std::vector<double>& dy) {
-  auto y0 = y[global_var_ids[0]];
-  auto y1 = y[global_var_ids[1]];
-  auto y2 = y[global_var_ids[2]];
   auto y3 = y[global_var_ids[3]];
-
   auto dy0 = dy[global_var_ids[0]];
-  auto dy1 = dy[global_var_ids[1]];
   auto dy3 = dy[global_var_ids[3]];
 
-  auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
-  auto capacitance = alpha[global_param_ids[ParamId::CAPACITANCE]];
-  auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
-  double stenosis_coeff = 0.0;
-
-  if (global_param_ids.size() > 3) {
-    stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
-  }
-  auto stenosis_resistance = stenosis_coeff * fabs(y3);
-
   jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y3;
   jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -dy3;
-
   if (global_param_ids.size() > 3) {
     jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = -fabs(y3) * y3;
   }
-
   jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) = -dy0;
-
-  residual(global_eqn_ids[0]) =
-      y0 - (resistance + stenosis_resistance) * y3 - y2 - inductance * dy3;
-  residual(global_eqn_ids[1]) = y1 - y3 - capacitance * dy0;
 }

src/model/BloodVesselCRL.h

@@ -174,17 +174,14 @@ class BloodVesselCRL : public Block {
                        const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy);
 
   /**
-   * @brief Set the gradient of the block contributions with respect to the
-   * parameters
+   * @brief Write this block's columns of the parameter Jacobian.
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
-   * @param y Current solution
-   * @param dy Time-derivative of the current solution
+   * @param y Observed solution
+   * @param dy Observed time derivative of the solution
    */
   void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                        Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                        std::vector<double>& y, std::vector<double>& dy);
 

src/model/BloodVesselJunction.cpp

@@ -54,43 +54,22 @@ void BloodVesselJunction::update_solution(
 
 void BloodVesselJunction::update_gradient(
     Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
     Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
     std::vector<double>& dy) {
-  auto p_in = y[global_var_ids[0]];
-  auto q_in = y[global_var_ids[1]];
-
-  residual(global_eqn_ids[0]) = q_in;
   for (size_t i = 0; i < num_outlets; i++) {
-    // Get parameters
-    auto resistance = alpha[global_param_ids[i]];
-    auto inductance = alpha[global_param_ids[num_outlets + i]];
-    double stenosis_coeff = 0.0;
-    if (global_param_ids.size() / num_outlets > 2) {
-      stenosis_coeff = alpha[global_param_ids[2 * num_outlets + i]];
-    }
     auto q_out = y[global_var_ids[3 + 2 * i]];
-    auto p_out = y[global_var_ids[2 + 2 * i]];
     auto dq_out = dy[global_var_ids[3 + 2 * i]];
-    auto stenosis_resistance = stenosis_coeff * fabs(q_out);
 
     // Resistance
     jacobian.coeffRef(global_eqn_ids[i + 1], global_param_ids[i]) = -q_out;
-
     // Inductance
     jacobian.coeffRef(global_eqn_ids[i + 1],
                       global_param_ids[num_outlets + i]) = -dq_out;
-
-    // Stenosis Coefficent
+    // Stenosis coefficient
     if (global_param_ids.size() / num_outlets > 2) {
       jacobian.coeffRef(global_eqn_ids[i + 1],
                         global_param_ids[2 * num_outlets + i]) =
           -fabs(q_out) * q_out;
     }
-
-    residual(global_eqn_ids[0]) -= q_out;
-    residual(global_eqn_ids[i + 1]) =
-        p_in - p_out - (resistance + stenosis_resistance) * q_out -
-        inductance * dq_out;
   }
 }

src/model/BloodVesselJunction.h

@@ -204,17 +204,14 @@ class BloodVesselJunction : public Block {
       const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy);
 
   /**
-   * @brief Set the gradient of the block contributions with respect to the
-   * parameters
+   * @brief Write this block's columns of the parameter Jacobian.
    *
    * @param jacobian Jacobian with respect to the parameters
    * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
-   * @param y Current solution
-   * @param dy Time-derivative of the current solution
+   * @param y Observed solution
+   * @param dy Observed time derivative of the solution
    */
   void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
                        Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
                        std::vector<double>& y, std::vector<double>& dy);
 

src/model/ChamberElastanceInductor.h

@@ -170,6 +170,8 @@ class ChamberElastanceInductor : public Block {
   void update_time(SparseSystem& system,
                    std::vector<double>& parameters) override;
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /// @brief Number of triplets of element
   TripletsContributions num_triplets{6, 2, 0};
 

src/model/ChamberElastanceInductorExponential.h

@@ -66,6 +66,7 @@ class ChamberElastanceInductorExponential : public ChamberElastanceInductor {
   };
 
   void update_time(SparseSystem& system, std::vector<double>& parameters);
+  bool has_time_dependent_assembly() const override { return true; }
   void update_solution(SparseSystem& system, std::vector<double>& parameters,
                        const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
                        const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy);

src/model/ChamberSphere.h

@@ -202,6 +202,8 @@ class ChamberSphere : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Update the solution-dependent contributions of the element in a
    * sparse system

src/model/ClosedLoopHeartPulmonary.h

@@ -215,6 +215,8 @@ class ClosedLoopHeartPulmonary : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Update the solution-dependent contributions of the element in a
    * sparse system

src/model/FlowReferenceBC.h

@@ -116,6 +116,8 @@ class FlowReferenceBC : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Number of triplets of element
    *

src/model/Junction.cpp

@@ -31,14 +31,3 @@ void Junction::update_constant(SparseSystem& system,
                       global_var_ids[i]) = -1.0;
   }
 }
-
-void Junction::update_gradient(
-    Eigen::SparseMatrix<double>& jacobian,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-    Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
-    std::vector<double>& dy) {
-  // Pressure conservation
-  residual(global_eqn_ids[0]) = y[global_var_ids[0]] - y[global_var_ids[2]];
-
-  residual(global_eqn_ids[1]) = y[global_var_ids[1]] - y[global_var_ids[3]];
-}

src/model/Junction.h

@@ -119,21 +119,6 @@ class Junction : public Block {
    */
   void update_constant(SparseSystem& system, std::vector<double>& parameters);
 
-  /**
-   * @brief Set the gradient of the block contributions with respect to the
-   * parameters
-   *
-   * @param jacobian Jacobian with respect to the parameters
-   * @param alpha Current parameter vector
-   * @param residual Residual with respect to the parameters
-   * @param y Current solution
-   * @param dy Time-derivative of the current solution
-   */
-  void update_gradient(Eigen::SparseMatrix<double>& jacobian,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
-                       Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
-                       std::vector<double>& y, std::vector<double>& dy);
-
   /**
    * @brief Number of triplets of element
    *

src/model/LinearElastanceChamber.h

@@ -171,6 +171,8 @@ class LinearElastanceChamber : public Block {
   void update_time(SparseSystem& system,
                    std::vector<double>& parameters) override;
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Number of triplets of element
    *

src/model/OpenLoopCoronaryBC.h

@@ -171,6 +171,8 @@ class OpenLoopCoronaryBC : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Number of triplets of element
    *

src/model/PressureReferenceBC.h

@@ -117,6 +117,8 @@ class PressureReferenceBC : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Number of triplets of element
    *

src/model/ResistanceBC.h

@@ -110,6 +110,8 @@ class ResistanceBC : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Number of triplets of element
    *

src/model/WindkesselBC.h

@@ -144,6 +144,8 @@ class WindkesselBC : public Block {
    */
   void update_time(SparseSystem& system, std::vector<double>& parameters);
 
+  bool has_time_dependent_assembly() const override { return true; }
+
   /**
    * @brief Number of triplets of element
    *