diff --git a/applications/svZeroDGUI/static/index.html b/applications/svZeroDGUI/static/index.html
index e1422bc53..9bd1eeda2 100644
--- a/applications/svZeroDGUI/static/index.html
+++ b/applications/svZeroDGUI/static/index.html
@@ -4,9 +4,8 @@
Vasculature Simulation
-
+
-
diff --git a/scripts/ChamberSphere_expmat.yaml b/scripts/ChamberSphere_expmat.yaml
new file mode 100644
index 000000000..57cdecc34
--- /dev/null
+++ b/scripts/ChamberSphere_expmat.yaml
@@ -0,0 +1,61 @@
+variables:
+ - Pin
+ - Qin
+ - Pout
+ - Qout
+ - radius
+ - velo
+ - stress
+ - tau
+ - volume
+
+derivatives:
+ - dPin_dt
+ - dQin_dt
+ - dPout_dt
+ - dQout_dt
+ - dradius_dt
+ - dvelo_dt
+ - dstress_dt
+ - dtau_dt
+ - dvolume_dt
+
+constants:
+ - rho
+ - thick0
+ - radius0
+ - C0
+ - C1
+ - C2
+ - C3
+ - eta
+ - act
+ - act_plus
+ - sigma_max
+
+time_dependent:
+ - act
+ - act_plus
+
+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
+ - dtau_dt + act * tau - sigma_max * act_plus
+ - dradius_dt - velo
+ - Qin - Qout - dvolume_dt
+ - Pin - Pout
\ No newline at end of file
diff --git a/scripts/jacobian.py b/scripts/jacobian.py
index 17cfd1694..aa8d8ece2 100755
--- a/scripts/jacobian.py
+++ b/scripts/jacobian.py
@@ -1,4 +1,4 @@
-from sympy import symbols, Matrix, simplify, pi, Abs
+from sympy import symbols, Matrix, simplify, pi, Abs, exp
from sympy.printing import ccode
import re
import pdb
@@ -16,6 +16,7 @@ def load_model(filepath):
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)
diff --git a/src/model/Block.h b/src/model/Block.h
index 752a5c0c0..6861efe92 100644
--- a/src/model/Block.h
+++ b/src/model/Block.h
@@ -15,6 +15,7 @@
#include "ActivationFunction.h"
#include "BlockType.h"
+#include "SphereMaterial.h"
#include "DOFHandler.h"
#include "Parameter.h"
#include "SparseSystem.h"
@@ -297,6 +298,18 @@ class Block {
virtual void set_activation_function(std::unique_ptr af) {
(void)af; // Included to avoid unused parameter warning
}
+
+ /**
+ * @brief Set wall material (for ChamberSphere blocks).
+ *
+ * Default no-op. Overridden by ChamberSphere to take ownership of the
+ * material.
+ *
+ * @param m Unique pointer to the material (caller transfers ownership)
+ */
+ virtual void set_material(std::unique_ptr m) {
+ (void)m; // Included to avoid unused parameter warning
+ }
};
#endif
diff --git a/src/model/CMakeLists.txt b/src/model/CMakeLists.txt
index 88ba98d5f..579d6f22c 100644
--- a/src/model/CMakeLists.txt
+++ b/src/model/CMakeLists.txt
@@ -5,8 +5,9 @@
set(lib svzero_model_library)
-set(CXXSRCS
+set(CXXSRCS
ActivationFunction.cpp
+ SphereMaterial.cpp
BloodVesselRC.cpp
Block.cpp
BloodVessel.cpp
@@ -37,8 +38,9 @@ set(CXXSRCS
PiecewiseValve.cpp
)
-set(HDRS
+set(HDRS
ActivationFunction.h
+ SphereMaterial.h
BloodVesselRC.h
Block.h
BlockType.h
diff --git a/src/model/ChamberSphere.cpp b/src/model/ChamberSphere.cpp
index 645a7089d..30edf99a7 100644
--- a/src/model/ChamberSphere.cpp
+++ b/src/model/ChamberSphere.cpp
@@ -53,13 +53,11 @@ void ChamberSphere::update_solution(
SparseSystem& system, std::vector& parameters,
const Eigen::Matrix& y,
const Eigen::Matrix& dy) {
- const double W1 = parameters[global_param_ids[ParamId::W1]];
- const double W2 = parameters[global_param_ids[ParamId::W2]];
- const double eta = parameters[global_param_ids[ParamId::eta]];
const double thick0 = parameters[global_param_ids[ParamId::thick0]];
const double sigma_max = parameters[global_param_ids[ParamId::sigma_max]];
-
const double radius0 = parameters[global_param_ids[ParamId::radius0]];
+ const double eta = parameters[global_param_ids[ParamId::eta]];
+
const double velo = y[global_var_ids[5]];
const double Pout = y[global_var_ids[2]];
const double stress = y[global_var_ids[6]];
@@ -77,24 +75,25 @@ void ChamberSphere::update_solution(
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[6]) =
thick0 * (radius + radius0) / pow(radius0, 2);
- // spherical stress
- system.C.coeffRef(global_eqn_ids[1]) =
- 2 *
- (dradius_dt * eta * (2 * pow(radius0, 12) + pow(radius + radius0, 12)) +
- 2 * pow(radius + radius0, 5) *
- (-pow(radius0, 6) + pow(radius + radius0, 6)) *
- (W1 * pow(radius0, 2) + W2 * pow(radius + radius0, 2))) /
+ // spherical stress: material-dependent elastic term plus viscous damping
+ const auto mat = material_->compute(radius, radius0);
+
+ const double C_val_visc =
+ 2 * dradius_dt * eta * (2 * pow(radius0, 12) + pow(radius + radius0, 12)) /
(pow(radius0, 2) * pow(radius + radius0, 11));
- system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[4]) =
- 24 * W1 * pow(radius0, 6) / pow(radius + radius0, 7) +
- 8 * W2 * radius / pow(radius0, 2) +
- 16 * W2 * pow(radius0, 4) / pow(radius + radius0, 5) + 8 * W2 / radius0 -
- 44 * dradius_dt * eta * pow(radius0, 10) / pow(radius + radius0, 12) +
- 2 * dradius_dt * eta / pow(radius0, 2);
- system.dC_dydot.coeffRef(global_eqn_ids[1], global_var_ids[4]) =
+ const double dC_dy_visc =
+ 2 * dradius_dt * eta / pow(radius0, 2) -
+ 44 * dradius_dt * eta * pow(radius0, 10) / pow(radius + radius0, 12);
+ const double dC_dydot_visc =
2 * eta * (2 * pow(radius0, 12) + pow(radius + radius0, 12)) /
(pow(radius0, 2) * pow(radius + radius0, 11));
+ system.C.coeffRef(global_eqn_ids[1]) = mat.C_val + C_val_visc;
+ system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[4]) =
+ mat.dC_dy_radius + dC_dy_visc;
+ system.dC_dydot.coeffRef(global_eqn_ids[1], global_var_ids[4]) =
+ dC_dydot_visc;
+
// volume change
system.C.coeffRef(global_eqn_ids[2]) =
4 * M_PI * velo * pow(radius + radius0, 2);
@@ -107,6 +106,10 @@ void ChamberSphere::update_solution(
system.C.coeffRef(global_eqn_ids[3]) = -act_plus * sigma_max;
}
+void ChamberSphere::set_material(std::unique_ptr m) {
+ material_ = std::move(m);
+}
+
void ChamberSphere::get_elastance_values(std::vector& parameters) {
const double alpha_max = parameters[global_param_ids[ParamId::alpha_max]];
const double alpha_min = parameters[global_param_ids[ParamId::alpha_min]];
diff --git a/src/model/ChamberSphere.h b/src/model/ChamberSphere.h
index 4039288a9..ea53614e3 100644
--- a/src/model/ChamberSphere.h
+++ b/src/model/ChamberSphere.h
@@ -9,7 +9,10 @@
#include
+#include
+
#include "Block.h"
+#include "SphereMaterial.h"
#include "SparseSystem.h"
/**
@@ -41,9 +44,11 @@
*
* 2. Spherical stress:
* \f[
- * -S + \tau + 4 (1 - C^{-3}) (W_1 + C W_2) + \eta \dot{C}
+ * -S + \tau + S_\text{el}(r, r_0) + \eta \dot{C}
* (1 + 2 C^{-6}) = 0
* \f]
+ * where \f$S_\text{el}\f$ is the elastic stress term supplied by the
+ * chamber's \ref SphereMaterial.
*
* 3. Volume change:
* \f[
@@ -83,8 +88,6 @@
* * `rho` - Density \f$\rho\f$
* * `thick0` - Wall thickness \f$d_0\f$
* * `radius0` - Reference radius \f$r_0\f$
- * * `W1` - Material constant \f$W_1\f$
- * * `W2` - Material constant \f$W_2\f$
* * `eta` - Viscosity parameter \f$\eta\f$
* * `sigma_max` - Maximum active stress \f$\sigma_\text{max}\f$
* * `alpha_max` - Maximum activation parameter \f$\alpha_\text{max}\f$
@@ -106,8 +109,6 @@
* "rho" : 1e3,
* "thick0" : 0.01,
* "radius0" : 0.05,
- * "W1" : 10e3,
- * "W2" : 40,
* "eta" : 10.0,
* "sigma_max" : 185e3,
* "alpha_max": 30.0,
@@ -116,8 +117,11 @@
* "tdias": 0.484,
* "steepness": 0.005
* }
+ * "material": {}
* }
* ]
+ *
+ * Material definition can be found in \ref SphereMaterial.
*
* ### Internal variables
*
@@ -140,15 +144,13 @@ class ChamberSphere : public Block {
rho = 0,
thick0 = 1,
radius0 = 2,
- W1 = 3,
- W2 = 4,
- eta = 5,
- sigma_max = 6,
- alpha_max = 7,
- alpha_min = 8,
- tsys = 9,
- tdias = 10,
- steepness = 11
+ eta = 3,
+ sigma_max = 4,
+ alpha_max = 5,
+ alpha_min = 6,
+ tsys = 7,
+ tdias = 8,
+ steepness = 9
};
/**
@@ -162,8 +164,6 @@ class ChamberSphere : public Block {
{{"rho", InputParameter()},
{"thick0", InputParameter()},
{"radius0", InputParameter()},
- {"W1", InputParameter()},
- {"W2", InputParameter()},
{"eta", InputParameter()},
{"sigma_max", InputParameter()},
{"alpha_max", InputParameter()},
@@ -225,10 +225,14 @@ class ChamberSphere : public Block {
*/
void get_elastance_values(std::vector& parameters);
+ void set_material(std::unique_ptr m) override;
+
private:
double act = 0.0; // activation function
double act_plus = 0.0; // act_plus = max(act, 0)
+ std::unique_ptr material_;
+
/**
* @brief Number of triplets of element
*
diff --git a/src/model/SphereMaterial.cpp b/src/model/SphereMaterial.cpp
new file mode 100644
index 000000000..09176c6ae
--- /dev/null
+++ b/src/model/SphereMaterial.cpp
@@ -0,0 +1,103 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#include "SphereMaterial.h"
+
+#include
+#include
+
+SphereMaterial::SphereMaterial(
+ const std::vector>& props)
+ : input_param_properties(props) {
+ for (const auto& p : props) {
+ if (p.second.is_number) {
+ params_[p.first] = p.second.is_optional ? p.second.default_val : 0.0;
+ }
+ }
+}
+
+void SphereMaterial::set_param(const std::string& name, double value) {
+ params_[name] = value;
+}
+
+std::unique_ptr SphereMaterial::create(
+ const std::string& type_str) {
+ if (type_str == "mooney_rivlin") {
+ return std::make_unique();
+ }
+ if (type_str == "exponential") {
+ return std::make_unique();
+ }
+ throw std::runtime_error(
+ "Unknown material type '" + type_str +
+ "'. Must be one of: mooney_rivlin, exponential");
+}
+
+SphericalStressResult MooneyRivlinMaterial::compute(double radius,
+ double radius0) const {
+ const double W1 = params_.at("W1");
+ const double W2 = params_.at("W2");
+
+ SphericalStressResult res;
+ res.C_val =
+ 2 *
+ (2 * pow(radius + radius0, 5) *
+ (-pow(radius0, 6) + pow(radius + radius0, 6)) *
+ (W1 * pow(radius0, 2) + W2 * pow(radius + radius0, 2))) /
+ (pow(radius0, 2) * pow(radius + radius0, 11));
+
+ res.dC_dy_radius =
+ 24 * W1 * pow(radius0, 6) / pow(radius + radius0, 7) +
+ 8 * W2 * radius / pow(radius0, 2) +
+ 16 * W2 * pow(radius0, 4) / pow(radius + radius0, 5) +
+ 8 * W2 / radius0;
+
+ return res;
+}
+
+SphericalStressResult ExponentialMaterial::compute(double radius,
+ double radius0) const {
+ const double C0 = params_.at("C0");
+ const double C1 = params_.at("C1");
+ const double C2 = params_.at("C2");
+ const double C3 = params_.at("C3");
+
+ // S2, S6: differences of the deformed/reference radius raised to the
+ // matching power. A: the I1,iso-like bracket driving the C0/C1 term.
+ const double S2 = pow(radius + radius0, 2) - pow(radius0, 2);
+ const double S6 = pow(radius + radius0, 6) - pow(radius0, 6);
+ const double A = pow(radius0, 6) -
+ 3 * pow(radius0, 2) * pow(radius + radius0, 4) +
+ 2 * pow(radius + radius0, 6);
+
+ // Each exponential is evaluated identically several times in the
+ // unexpanded expression, so it is computed once here.
+ const double E1 =
+ exp(C1 * A * A / (pow(radius0, 4) * pow(radius + radius0, 8)));
+ const double E2 = exp(C3 * S2 * S2 / pow(radius0, 4));
+
+ SphericalStressResult res;
+ res.C_val =
+ 2 *
+ (4 * C0 * C1 * (radius + radius0) * S6 * A * E1 +
+ 2 * C2 * C3 * pow(radius + radius0, 11) * S2 * E2) /
+ (pow(radius0, 2) * pow(radius + radius0, 11));
+
+ res.dC_dy_radius =
+ 2 *
+ (-32 * C0 * pow(C1, 2) * S6 * (A * A) *
+ (A - 3 * pow(radius + radius0, 4) * S2) * E1 +
+ 8 * C2 * pow(C3, 2) * pow(radius + radius0, 20) * (S2 * S2) * E2 +
+ 2 * pow(radius0, 4) * pow(radius + radius0, 8) *
+ (12 * C0 * C1 * pow(radius + radius0, 6) * A * E1 +
+ 24 * C0 * C1 * pow(radius + radius0, 4) * S2 * S6 * E1 +
+ 2 * C0 * C1 * S6 * A * E1 +
+ 2 * C2 * C3 * pow(radius + radius0, 12) * E2 +
+ 11 * C2 * C3 * pow(radius + radius0, 10) * S2 * E2) +
+ 11 * pow(radius0, 4) * pow(radius + radius0, 7) *
+ (-4 * C0 * C1 * (radius + radius0) * S6 * A * E1 -
+ 2 * C2 * C3 * pow(radius + radius0, 11) * S2 * E2)) /
+ (pow(radius0, 6) * pow(radius + radius0, 19));
+
+ return res;
+}
diff --git a/src/model/SphereMaterial.h b/src/model/SphereMaterial.h
new file mode 100644
index 000000000..4e77af438
--- /dev/null
+++ b/src/model/SphereMaterial.h
@@ -0,0 +1,133 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+/**
+ * @file SphereMaterial.h
+ * @brief Material models for the ChamberSphere block
+ */
+
+#ifndef SVZERODSOLVER_MODEL_SPHEREMATERIAL_HPP_
+#define SVZERODSOLVER_MODEL_SPHEREMATERIAL_HPP_
+
+#include