Valve status output

src/model/ValveTanh.cpp

@@ -4,8 +4,9 @@
 
 void ValveTanh::setup_dofs(DOFHandler& dofhandler) {
   // set_up_dofs args: dofhandler (passed in), num equations, list of internal
-  // variable names (strings) 2 eqns, one for Pressure, one for Flow
-  Block::setup_dofs_(dofhandler, 2, {});
+  // variable names (strings) 3 eqns, one for Pressure, one for Flow, one for
+  // the valve status output
+  Block::setup_dofs_(dofhandler, 3, {"valve_status"});
 }
 
 // update_constant updates matrices E and F from E(y,t)*y_dot + F(y,t)*y +
@@ -25,6 +26,7 @@ void ValveTanh::update_constant(SparseSystem& system,
       -0.5 * (Rmax + Rmin);
   system.F.coeffRef(global_eqn_ids[1], global_var_ids[1]) = 1.0;
   system.F.coeffRef(global_eqn_ids[1], global_var_ids[3]) = -1.0;
+  system.F.coeffRef(global_eqn_ids[2], global_var_ids[4]) = 1.0;
 }
 
 // update_solution updates matrices E and F from E(y,t)*y_dot + F(y,t)*y +
@@ -43,17 +45,21 @@ void ValveTanh::update_solution(
   double Rmax = parameters[global_param_ids[ParamId::RMAX]];
   double steep = parameters[global_param_ids[ParamId::STEEPNESS]];
 
-  // Nonlinear term
-  system.C(global_eqn_ids[0]) =
-      -0.5 * q_in * (Rmax - Rmin) * tanh(steep * (p_out - p_in));
+  // Helper functions
+  double fun_tanh = tanh(steep * (p_out - p_in));
+  double fun_cosh = 0.5 * steep / pow(cosh(steep * (p_in - p_out)), 2);
 
-  // Derivatives of non-linear term
+  // Nonlinear terms
+  system.C(global_eqn_ids[0]) = -0.5 * q_in * (Rmax - Rmin) * fun_tanh;
+  system.C(global_eqn_ids[2]) = -0.5 * (1 + fun_tanh);
+
+  // Derivatives of non-linear terms
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[0]) =
-      0.5 * q_in * (Rmax - Rmin) * steep *
-      (1.0 - tanh(steep * (p_out - p_in)) * tanh(steep * (p_out - p_in)));
+      0.5 * q_in * (Rmax - Rmin) * steep * (1.0 - pow(fun_tanh, 2));
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[1]) =
-      -0.5 * (Rmax - Rmin) * tanh(steep * (p_out - p_in));
+      -0.5 * (Rmax - Rmin) * fun_tanh;
   system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[2]) =
-      -0.5 * q_in * (Rmax - Rmin) * steep *
-      (1.0 - tanh(steep * (p_out - p_in)) * tanh(steep * (p_out - p_in)));
+      -0.5 * q_in * (Rmax - Rmin) * steep * (1.0 - pow(fun_tanh, 2));
+  system.dC_dy.coeffRef(global_eqn_ids[2], global_var_ids[0]) = fun_cosh;
+  system.dC_dy.coeffRef(global_eqn_ids[2], global_var_ids[2]) = -fun_cosh;
 }

src/model/ValveTanh.h

@@ -40,38 +40,43 @@
  * Q_{in}-Q_{out}=0
  * \f]
  *
+ * \f[
+ * \text{valve\_status} = \frac{1}{2}\left(1+tanh\{k(P_{out}-P{in})\}\right)
+ * \f]
+ *
  * ### Local contributions
  *
  * \f[
- * \mathbf{y}^{e}=\left[\begin{array}{llll}P_{in} & Q_{in} &
- * P_{out} & Q_{out}\end{array}\right]^{T} \f]
+ * \mathbf{y}^{e}=\left[\begin{array}{lllll}P_{in} & Q_{in} &
+ * P_{out} & Q_{out} & \text{valve\_status} \end{array}\right]^{T} \f]
  *
  * \f[
- * \mathbf{E}^{e}=\left[\begin{array}{cccc}
- * 0 & 0 & 0 & 0 \\
- * 0 & 0 & 0 & 0
- * \end{array}\right]
+ * \mathbf{E}^{e}=\mathbf{0}
  * \f]
  *
  * \f[
- * \mathbf{F}^{e}=\left[\begin{array}{cccc}
- * 1 & -(R_{max}+R_{min})/2.0 & -1 & 0 \\
- * 0 &      1                 &  0 & -1
+ * \mathbf{F}^{e}=\left[\begin{array}{ccccc}
+ * 1 & -(R_{max}+R_{min})/2.0 & -1 & 0  & 0\\
+ * 0 &      1                 &  0 & -1 & 0\\
+ * 0 &      0                 &  0 & 0  & 1
  * \end{array}\right]
  * \f]
  *
  * \f[
  * \mathbf{c}^{e}=\left[\begin{array}{c}
  * -\frac{1}{2}Q_{in}(R_{max}-R_{min})tanh\{k(P_{out}-P_{in})\} \\
- * 0
+ * 0 \\
+ * -\frac{1}{2}\left[1+tanh\{k(P_{out}-P_{in})\}\right]
  * \end{array}\right]
  * \f]
  *
  * \f[
  * \left(\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right)^{e} =
- * \left[\begin{array}{cccc}
- * A & B & C & 0 \\
- * 0 & 0 & 0 & 0 \end{array}\right] \f]
+ * \left[\begin{array}{ccccc}
+ * A & B &  C & 0 & 0\\
+ * 0 & 0 &  0 & 0 & 0\\
+ * D & 0 & -D & 0 & 0
+ * \end{array}\right] \f]
  * where,
  * \f[
  * A = \frac{1}{2} k Q_{in}
@@ -82,14 +87,15 @@
  * \f]
  * \f[
  * C = -\frac{1}{2} k Q_{in}
- * (R_{max}-R_{min})\left[1-tanh^2\{k(P_{out}-P_{in})\}\right] \f]
+ * (R_{max}-R_{min})\left[1-tanh^2\{k(P_{out}-P_{in})\}\right]
+ * \f]
+ * \f[
+ * D = \frac{1}{2} \frac{k}{cosh^2\{k(P_{in}-P_{out})\} }
+ * \f]
  *
  * \f[
  * \left(\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)^{e} =
- * \left[\begin{array}{cccc}
- * 0 & 0 & 0 & 0 \\
- * 0 & 0 & 0 & 0
- * \end{array}\right]
+ *  \mathbf{E}^{e}=\mathbf{0}
  * \f]
  *
  * ### Parameters