Skip to content

Commit d40189e

Browse files
rjrios915ncdorn
rjrios915
and
ncdorn
authored
CRL Vessel (SimVascular#202)
Co-authored-by: ncdorn <ndorn@stanford.edu>
1 parent a1f5eca commit d40189e

13 files changed

+2121
-1826
lines changed

‎src/model/BlockType.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ enum class BlockType {
2929
valve_tanh = 13,
3030
chamber_elastance_inductor = 14,
3131
chamber_sphere = 15,
32+
blood_vessel_CRL = 16,
3233
piecewise_cosine_chamber = 17,
3334
piecewise_valve = 18
3435
};

‎src/model/BloodVesselCRL.cpp‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
2+
// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
3+
4+
#include "BloodVesselCRL.h"
5+
6+
void BloodVesselCRL::setup_dofs(DOFHandler& dofhandler) {
7+
Block::setup_dofs_(dofhandler, 2, {});
8+
}
9+
10+
void BloodVesselCRL::update_constant(SparseSystem& system,
11+
std::vector<double>& parameters) {
12+
// Get parameters
13+
double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]];
14+
double inductance = parameters[global_param_ids[ParamId::INDUCTANCE]];
15+
double resistance = parameters[global_param_ids[ParamId::RESISTANCE]];
16+
17+
// Set element contributions
18+
// coeffRef args are the indices (i,j) of the matrix
19+
// global_eqn_ids: number of rows in the matrix, set in setup_dofs
20+
// global_var_ids: number of columns, organized as pressure and flow of all
21+
// inlets and then all outlets of the block
22+
system.E.coeffRef(global_eqn_ids[0], global_var_ids[3]) = -inductance;
23+
system.E.coeffRef(global_eqn_ids[1], global_var_ids[0]) = -capacitance;
24+
system.F.coeffRef(global_eqn_ids[0], global_var_ids[0]) = 1.0;
25+
system.F.coeffRef(global_eqn_ids[0], global_var_ids[3]) = -resistance;
26+
system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -1.0;
27+
system.F.coeffRef(global_eqn_ids[1], global_var_ids[1]) = 1.0;
28+
system.F.coeffRef(global_eqn_ids[1], global_var_ids[3]) = -1.0;
29+
}
30+
31+
void BloodVesselCRL::update_solution(
32+
SparseSystem& system, std::vector<double>& parameters,
33+
const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
34+
const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
35+
// Get parameters
36+
double capacitance = parameters[global_param_ids[ParamId::CAPACITANCE]];
37+
double stenosis_coeff =
38+
parameters[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
39+
double q_out = y[global_var_ids[3]];
40+
double dq_out = dy[global_var_ids[3]];
41+
double stenosis_resistance = stenosis_coeff * fabs(q_out);
42+
43+
// Set element contributions
44+
system.C(global_eqn_ids[0]) = stenosis_resistance * -q_out;
45+
46+
double sgn_q_out = (0.0 < q_out) - (q_out < 0.0);
47+
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[1]) =
48+
stenosis_coeff * sgn_q_out * -2.0 * q_out;
49+
}
50+
51+
void BloodVesselCRL::update_gradient(
52+
Eigen::SparseMatrix<double>& jacobian,
53+
Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
54+
Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, std::vector<double>& y,
55+
std::vector<double>& dy) {
56+
auto y0 = y[global_var_ids[0]];
57+
auto y1 = y[global_var_ids[1]];
58+
auto y2 = y[global_var_ids[2]];
59+
auto y3 = y[global_var_ids[3]];
60+
61+
auto dy0 = dy[global_var_ids[0]];
62+
auto dy1 = dy[global_var_ids[1]];
63+
auto dy3 = dy[global_var_ids[3]];
64+
65+
auto resistance = alpha[global_param_ids[ParamId::RESISTANCE]];
66+
auto capacitance = alpha[global_param_ids[ParamId::CAPACITANCE]];
67+
auto inductance = alpha[global_param_ids[ParamId::INDUCTANCE]];
68+
double stenosis_coeff = 0.0;
69+
70+
if (global_param_ids.size() > 3) {
71+
stenosis_coeff = alpha[global_param_ids[ParamId::STENOSIS_COEFFICIENT]];
72+
}
73+
auto stenosis_resistance = stenosis_coeff * fabs(y3);
74+
75+
jacobian.coeffRef(global_eqn_ids[0], global_param_ids[0]) = -y3;
76+
jacobian.coeffRef(global_eqn_ids[0], global_param_ids[2]) = -dy3;
77+
78+
if (global_param_ids.size() > 3) {
79+
jacobian.coeffRef(global_eqn_ids[0], global_param_ids[3]) = -fabs(y3) * y3;
80+
}
81+
82+
jacobian.coeffRef(global_eqn_ids[1], global_param_ids[1]) = -dy0;
83+
84+
residual(global_eqn_ids[0]) =
85+
y0 - (resistance + stenosis_resistance) * y3 - y2 - inductance * dy3;
86+
residual(global_eqn_ids[1]) = y1 - y3 - capacitance * dy0;
87+
}

‎src/model/BloodVesselCRL.h‎

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
2+
// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
3+
4+
/**
5+
* @file BloodVesselCRL.h
6+
* @brief model::BloodVesselCRL source file
7+
*/
8+
#ifndef SVZERODSOLVER_MODEL_BLOODVESSELCRL_HPP_
9+
#define SVZERODSOLVER_MODEL_BLOODVESSELCRL_HPP_
10+
11+
#include <math.h>
12+
13+
#include "Block.h"
14+
#include "SparseSystem.h"
15+
16+
/**
17+
* @brief Capacitor-resistor-inductor blood vessel with optional stenosis
18+
*
19+
* Models the mechanical behavior of a bloodvesselCRL with optional stenosis.
20+
*
21+
* \f[
22+
* \begin{circuitikz} \draw
23+
* node[left] {$Q_{in}$} [-latex] (0,0) -- (0.8,0);
24+
* \draw (1,0) node[anchor=south]{$P_{in}$}
25+
* to [R, l=$R$, *-] (3,0)
26+
* to [R, l=$S$, -] (5,0)
27+
* (5,0) to [L, l=$L$, -*] (7,0)
28+
* node[anchor=south]{$P_{out}$}
29+
* (1,0) to [C, l=$C$, -] (1,-1.5)
30+
* node[ground]{};
31+
* \draw [-latex] (7.2,0) -- (8,0) node[right] {$Q_{out}$};
32+
* \end{circuitikz}
33+
* \f]
34+
*
35+
* ### Governing equations
36+
*
37+
* \f[
38+
* P_\text{in}-P_\text{out} - (R + S|Q_\text{out}|) Q_\text{out}-L
39+
* \dot{Q}_\text{out}=0 \f]
40+
*
41+
* \f[
42+
* Q_\text{in}-Q_\text{out} - C \dot{P}_\text{in}=0 \f]
43+
*
44+
* ### Local contributions
45+
*
46+
* \f[
47+
* \mathbf{y}^{e}=\left[\begin{array}{llll}P_{i n} & Q_{in} &
48+
* P_{out} & Q_{out}\end{array}\right]^\text{T} \f]
49+
*
50+
* \f[
51+
* \mathbf{F}^{e}=\left[\begin{array}{cccc}
52+
* 1 & 0 & -1 & -R \\
53+
* 0 & 1 & 0 & -1
54+
* \end{array}\right]
55+
* \f]
56+
*
57+
* \f[
58+
* \mathbf{E}^{e}=\left[\begin{array}{cccc}
59+
* 0 & 0 & 0 & -L \\
60+
* -C & 0 & 0 & 0
61+
* \end{array}\right]
62+
* \f]
63+
*
64+
* \f[
65+
* \mathbf{c}^{e} = S|Q_\text{out}|
66+
* \left[\begin{array}{c}
67+
* -Q_\text{out} \\
68+
* 0
69+
* \end{array}\right]
70+
* \f]
71+
*
72+
* \f[
73+
* \left(\frac{\partial\mathbf{c}}{\partial\mathbf{y}}\right)^{e} =
74+
* S \text{sgn} (Q_\text{in})
75+
* \left[\begin{array}{cccc}
76+
* 0 & -2Q_\text{out} & 0 & 0 \\
77+
* 0 & 0 & 0 & 0
78+
* \end{array}\right]
79+
* \f]
80+
*
81+
* \f[
82+
* \left(\frac{\partial\mathbf{c}}{\partial\dot{\mathbf{y}}}\right)^{e} =
83+
* S|Q_\text{out}|
84+
* \left[\begin{array}{cccc}
85+
* 0 & 0 & 0 & 0 \\
86+
* 0 & 0 & 0 & 0
87+
* \end{array}\right]
88+
* \f]
89+
*
90+
* with the stenosis resistance \f$ S=K_{t} \frac{\rho}{2
91+
* A_{o}^{2}}\left(\frac{A_{o}}{A_{s}}-1\right)^{2} \f$.
92+
* \f$R\f$, \f$C\f$, and \f$L\f$ refer to
93+
* Poisieuille resistance, capacitance and inductance, respectively.
94+
*
95+
* ### Gradient
96+
*
97+
* Gradient of the equations with respect to the parameters:
98+
*
99+
* \f[
100+
* \mathbf{J}^{e} = \left[\begin{array}{cccc}
101+
* -y_3 & 0 & -\dot{y}_3 & -|y_3|y_3 \\
102+
* 0 & 0 & -\dot{y}_0 & 0 \\
103+
* \end{array}\right]
104+
* \f]
105+
*
106+
* ### Parameters
107+
*
108+
* Parameter sequence for constructing this block
109+
*
110+
* * `0` Poiseuille resistance
111+
* * `1` Capacitance
112+
* * `2` Inductance
113+
* * `3` Stenosis coefficient
114+
*
115+
*/
116+
class BloodVesselCRL : public Block {
117+
public:
118+
/**
119+
* @brief Local IDs of the parameters
120+
*
121+
*/
122+
enum ParamId {
123+
RESISTANCE = 0,
124+
CAPACITANCE = 1,
125+
INDUCTANCE = 2,
126+
STENOSIS_COEFFICIENT = 3,
127+
};
128+
129+
/**
130+
* @brief Construct a new BloodVesselCRL object
131+
*
132+
* @param id Global ID of the block
133+
* @param model The model to which the block belongs
134+
*/
135+
BloodVesselCRL(int id, Model* model)
136+
: Block(id, model, BlockType::blood_vessel_CRL, BlockClass::vessel,
137+
{{"R_poiseuille", InputParameter()},
138+
{"C", InputParameter(true)},
139+
{"L", InputParameter(true)},
140+
{"stenosis_coefficient", InputParameter(true)}}) {}
141+
142+
/**
143+
* @brief Set up the degrees of freedom (DOF) of the block
144+
*
145+
* Set \ref global_var_ids and \ref global_eqn_ids of the element based on the
146+
* number of equations and the number of internal variables of the
147+
* element.
148+
*
149+
* @param dofhandler Degree-of-freedom handler to register variables and
150+
* equations at
151+
*/
152+
void setup_dofs(DOFHandler& dofhandler);
153+
154+
/**
155+
* @brief Update the constant contributions of the element in a sparse
156+
system
157+
*
158+
* @param system System to update contributions at
159+
* @param parameters Parameters of the model
160+
*/
161+
void update_constant(SparseSystem& system, std::vector<double>& parameters);
162+
163+
/**
164+
* @brief Update the solution-dependent contributions of the element in a
165+
* sparse system
166+
*
167+
* @param system System to update contributions at
168+
* @param parameters Parameters of the model
169+
* @param y Current solution
170+
* @param dy Current derivate of the solution
171+
*/
172+
void update_solution(SparseSystem& system, std::vector<double>& parameters,
173+
const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
174+
const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy);
175+
176+
/**
177+
* @brief Set the gradient of the block contributions with respect to the
178+
* parameters
179+
*
180+
* @param jacobian Jacobian with respect to the parameters
181+
* @param alpha Current parameter vector
182+
* @param residual Residual with respect to the parameters
183+
* @param y Current solution
184+
* @param dy Time-derivative of the current solution
185+
*/
186+
void update_gradient(Eigen::SparseMatrix<double>& jacobian,
187+
Eigen::Matrix<double, Eigen::Dynamic, 1>& residual,
188+
Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha,
189+
std::vector<double>& y, std::vector<double>& dy);
190+
191+
/**
192+
* @brief Number of triplets of element
193+
*
194+
* Number of triplets that the element contributes to the global system
195+
* (relevant for sparse memory reservation)
196+
*/
197+
TripletsContributions num_triplets{5, 3, 2};
198+
};
199+
200+
#endif // SVZERODSOLVER_MODEL_BLOODVESSELCRL_HPP_

‎src/model/BloodVesselJunction.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99

1010
#include "Block.h"
1111
#include "BloodVessel.h"
12+
#include "BloodVesselCRL.h"
1213
#include "SparseSystem.h"
1314

1415
/**

‎src/model/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ set(lib svzero_model_library)
88
set(CXXSRCS
99
Block.cpp
1010
BloodVessel.cpp
11+
BloodVesselCRL.cpp
1112
BloodVesselJunction.cpp
1213
ChamberSphere.cpp
1314
ChamberElastanceInductor.cpp
@@ -36,6 +37,7 @@ set(HDRS
3637
Block.h
3738
BlockType.h
3839
BloodVessel.h
40+
BloodVesselCRL.h
3941
BloodVesselJunction.h
4042
ChamberSphere.h
4143
ChamberElastanceInductor.h

‎src/model/Model.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ Model::Model() {
2929
{"resistive_junction", block_factory<ResistiveJunction>()},
3030
{"ValveTanh", block_factory<ValveTanh>()},
3131
{"ChamberElastanceInductor", block_factory<ChamberElastanceInductor>()},
32+
{"BloodVesselCRL", block_factory<BloodVesselCRL>()},
3233
{"PiecewiseValve", block_factory<PiecewiseValve>()},
3334
{"PiecewiseCosineChamber", block_factory<PiecewiseCosineChamber>()}};
3435
}

‎src/model/Model.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@
1818
#include "Block.h"
1919
#include "BlockFactory.h"
2020
#include "BloodVessel.h"
21+
#include "BloodVesselCRL.h"
2122
#include "BloodVesselJunction.h"
2223
#include "ChamberElastanceInductor.h"
2324
#include "ChamberSphere.h"

‎src/solve/csv_writer.cpp‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -46,7 +46,8 @@ std::string to_vessel_csv(const std::vector<double>& times,
4646
if (dynamic_cast<const BloodVessel*>(block) == nullptr &&
4747
dynamic_cast<const ChamberSphere*>(block) == nullptr &&
4848
dynamic_cast<const PiecewiseCosineChamber*>(block) == nullptr &&
49-
dynamic_cast<const PiecewiseValve*>(block) == nullptr) {
49+
dynamic_cast<const PiecewiseValve*>(block) == nullptr &&
50+
dynamic_cast<const BloodVesselCRL*>(block) == nullptr) {
5051
continue;
5152
}
5253

tests/cases/double_pulsatileFlow_CRL.json renamed to tests/cases/pulsatileFlow_CRL.json

Lines changed: 9 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
"bc_type": "PRESSURE",
1919
"bc_values": {
2020
"P": [
21-
0.01047121, 0.020937827, 0.03139526, 0.041838922, 0.052264232,
21+
0.0, 0.01047121, 0.020937827, 0.03139526, 0.041838922, 0.052264232,
2222
0.062666617, 0.073041514, 0.083384373, 0.093690657, 0.103955845,
2323
0.114175435, 0.124344944, 0.13445991, 0.144515898, 0.154508497,
2424
0.164433323, 0.174286024, 0.184062276, 0.193757793, 0.203368322,
@@ -80,7 +80,7 @@
8080
-0.010471211, -5.89793e-10
8181
],
8282
"t": [
83-
0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
83+
0.0, 0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
8484
0.125663706, 0.146607657, 0.167551608, 0.188495559, 0.20943951,
8585
0.230383461, 0.251327412, 0.272271363, 0.293215314, 0.314159265,
8686
0.335103216, 0.356047167, 0.376991118, 0.397935069, 0.41887902,
@@ -148,7 +148,7 @@
148148
"bc_type": "FLOW",
149149
"bc_values": {
150150
"Q": [
151-
0.99912283, 0.996492859, 0.992114701, 0.985996037, 0.978147601,
151+
1.0, 0.99912283, 0.996492859, 0.992114701, 0.985996037, 0.978147601,
152152
0.968583161, 0.957319498, 0.94437637, 0.929776486, 0.913545458,
153153
0.89571176, 0.87630668, 0.85536426, 0.832921241, 0.809016994,
154154
0.783693457, 0.756995056, 0.728968628, 0.699663341, 0.669130606,
@@ -207,10 +207,10 @@
207207
0.756995054, 0.783693456, 0.809016993, 0.832921239, 0.855364259,
208208
0.876306679, 0.895711759, 0.913545457, 0.929776485, 0.944376369,
209209
0.957319497, 0.968583161, 0.9781476, 0.985996037, 0.992114701,
210-
0.996492859, 0.99912283, 1
210+
0.996492859, 0.99912283, 1.0
211211
],
212212
"t": [
213-
0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
213+
0.0, 0.020943951, 0.041887902, 0.062831853, 0.083775804, 0.104719755,
214214
0.125663706, 0.146607657, 0.167551608, 0.188495559, 0.20943951,
215215
0.230383461, 0.251327412, 0.272271363, 0.293215314, 0.314159265,
216216
0.335103216, 0.356047167, 0.376991118, 0.397935069, 0.41887902,
@@ -277,8 +277,10 @@
277277

278278
"simulation_parameters": {
279279
"number_of_cardiac_cycles": 10,
280-
"number_of_time_pts_per_cardiac_cycle": 100,
281-
"output_variable_based": false
280+
"number_of_time_pts_per_cardiac_cycle": 300,
281+
"output_variable_based": false,
282+
"output_all_cycles": false,
283+
"cardiac_period": 6.283185306
282284
},
283285
"vessels": [
284286
{

0 commit comments

Comments
 (0)