Skip to content

Commit 3cfd49c

Browse files
mrp089claude
mrp089
and
claude
committed
Set ChamberSphere shape factor n=1 and remove viscosity
Bake n=1 and eta=0 into the ChamberSphere block: substitute n=1 in the symbolic stretch helper, drop the gamma_eta viscosity term from the stress residual, and regenerate the C++ from scripts/ChamberSphere.yaml. Remove n and gamma_eta from the ParamId enum and input parameter list, fix the triplet reservation counts, and update the documented equations/parameters. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
1 parent 7a1b728 commit 3cfd49c

3 files changed

Lines changed: 36 additions & 67 deletions

File tree

‎scripts/ChamberSphere.yaml‎

Lines changed: 2 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -22,10 +22,8 @@ derivatives:
2222
# stress here represents gamma * stress
2323
# tau here represents gamma * tau
2424
constants:
25-
- n
2625
- volume0
2726
- gamma_W1
28-
- gamma_eta
2927
- gamma_sigma_max # gamma * sigma_max
3028
- act
3129
- act_plus
@@ -41,24 +39,15 @@ helper_functions: |
4139
4240
def stretch(volume):
4341
V_val = V(volume)
44-
return (V_val / volume0) ** (1 / (3 * n))
42+
return (V_val / volume0) ** (1 / 3)
4543
4644
def CG(volume):
4745
stretch_val = stretch(volume)
4846
return (stretch_val) ** 2
4947
50-
def dstretch_dt(volume, dvolume_dt):
51-
V_val = V(volume)
52-
return (1 / (3 * n)) * (V_val / volume0) ** ((1 / (3 * n)) - 1) * (1 / volume0) * dvolume_dt
53-
54-
def dCG(volume, dvolume_dt):
55-
stretch_val = stretch(volume)
56-
dstretch = dstretch_dt(volume, dvolume_dt)
57-
return 2 * stretch_val * dstretch
58-
5948
residuals:
6049
- stretch(volume) * stress - Pout * CG(volume)
61-
- -stress + tau + 4 * gamma_W1 * (1 - CG(volume) ** (-3)) + prestress + gamma_eta * dCG(volume, dvolume_dt) * (1 + 2 * CG(volume) ** (-6))
50+
- -stress + tau + 4 * gamma_W1 * (1 - CG(volume) ** (-3)) + prestress
6251
- dtau_dt + act * tau - gamma_sigma_max * act_plus
6352
- Qin - Qout - dvolume_dt
6453
- Pin - Pout

‎src/model/ChamberSphere.cpp‎

Lines changed: 12 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -12,10 +12,6 @@ void ChamberSphere::setup_dofs(DOFHandler& dofhandler) {
1212

1313
void ChamberSphere::update_constant(SparseSystem& system,
1414
std::vector<double>& parameters) {
15-
const double gamma_eta = parameters[global_param_ids[ParamId::gamma_eta]];
16-
const double volume0 = parameters[global_param_ids[ParamId::volume0]];
17-
const double n = parameters[global_param_ids[ParamId::n]];
18-
system.E.coeffRef(global_eqn_ids[1], global_var_ids[6]) = 2*gamma_eta/(n*volume0);
1915
system.E.coeffRef(global_eqn_ids[2], global_var_ids[5]) = 1;
2016
system.E.coeffRef(global_eqn_ids[3], global_var_ids[6]) = -1;
2117
system.F.coeffRef(global_eqn_ids[0], global_var_ids[2]) = -1;
@@ -30,35 +26,30 @@ void ChamberSphere::update_constant(SparseSystem& system,
3026

3127
void ChamberSphere::update_time(SparseSystem& system,
3228
std::vector<double>& parameters) {
29+
const double gamma_sigma_max = parameters[global_param_ids[ParamId::gamma_sigma_max]];
30+
3331
// active stress
3432
get_elastance_values(parameters);
3533
system.F.coeffRef(global_eqn_ids[2], global_var_ids[5]) = act;
34+
system.C.coeffRef(global_eqn_ids[2]) = -act_plus*gamma_sigma_max;
3635
}
3736

3837
void ChamberSphere::update_solution(
3938
SparseSystem& system, std::vector<double>& parameters,
4039
const Eigen::Matrix<double, Eigen::Dynamic, 1>& y,
4140
const Eigen::Matrix<double, Eigen::Dynamic, 1>& dy) {
4241
const double prestress = parameters[global_param_ids[ParamId::prestress]];
43-
const double gamma_W1 = parameters[global_param_ids[ParamId::gamma_W1]];
44-
const double gamma_eta = parameters[global_param_ids[ParamId::gamma_eta]];
4542
const double volume0 = parameters[global_param_ids[ParamId::volume0]];
46-
const double gamma_sigma_max = parameters[global_param_ids[ParamId::gamma_sigma_max]];
47-
const double n = parameters[global_param_ids[ParamId::n]];
48-
const double volume = y[global_var_ids[6]];
49-
const double dvolume_dt = dy[global_var_ids[6]];
50-
const double stress = y[global_var_ids[4]];
43+
const double gamma_W1 = parameters[global_param_ids[ParamId::gamma_W1]];
5144
const double Pout = y[global_var_ids[2]];
52-
system.C.coeffRef(global_eqn_ids[0]) = -Pout*pow((volume + volume0)/volume0, (2.0/3.0)/n) + Pout + stress*pow((volume + volume0)/volume0, (1.0/3.0)/n) - stress;
53-
system.C.coeffRef(global_eqn_ids[1]) = -2*dvolume_dt*gamma_eta/(n*volume0) + (2.0/3.0)*dvolume_dt*gamma_eta*pow(volume/volume0 + 1, -(5.0/3.0)/n)*pow(volume/volume0 + 1, -1 + (7.0/3.0)/n)/(n*volume0) + (4.0/3.0)*dvolume_dt*gamma_eta*pow(volume/volume0 + 1, -(17.0/3.0)/n)*pow(volume/volume0 + 1, -1 + (7.0/3.0)/n)/(n*volume0) + 4*gamma_W1 - 4*gamma_W1*pow(volume/volume0 + 1, -2/n) + prestress;
54-
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[2]) = 1 - pow((volume + volume0)/volume0, (2.0/3.0)/n);
55-
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[4]) = pow((volume + volume0)/volume0, (1.0/3.0)/n) - 1;
56-
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[6]) = (1.0/3.0)*pow((volume + volume0)/volume0, (1.0/3.0)/n)*(-2*Pout*pow((volume + volume0)/volume0, (1.0/3.0)/n) + stress)/(n*(volume + volume0));
57-
system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[6]) = (2.0/9.0)*pow((volume + volume0)/volume0, -(17.0/3.0)/n)*(dvolume_dt*gamma_eta*pow((volume + volume0)/volume0, (19.0/3.0 - n)/n)*(2 - 3*n) - 2*dvolume_dt*gamma_eta*pow((volume + volume0)/volume0, (1.0/3.0)*(7 - 3*n)/n)*(3*n + 10) + 36*gamma_W1*n*volume0*pow((volume + volume0)/volume0, (11.0/3.0)/n))/(pow(n, 2)*volume0*(volume + volume0));
58-
system.dC_dydot.coeffRef(global_eqn_ids[1], global_var_ids[6]) = (2.0/3.0)*gamma_eta*pow((volume + volume0)/volume0, -(22.0/3.0)/n)*(-3*pow((volume + volume0)/volume0, (22.0/3.0)/n) + 2*pow((volume + volume0)/volume0, (4 - n)/n) + pow((volume + volume0)/volume0, (8 - n)/n))/(n*volume0);
59-
60-
// active stress
61-
system.C.coeffRef(global_eqn_ids[2]) = -act_plus*gamma_sigma_max;
45+
const double stress = y[global_var_ids[4]];
46+
const double volume = y[global_var_ids[6]];
47+
system.C.coeffRef(global_eqn_ids[0]) = -Pout*pow((volume + volume0)/volume0, 0.66666666666666663) + Pout + stress*pow((volume + volume0)/volume0, 0.33333333333333331) - stress;
48+
system.C.coeffRef(global_eqn_ids[1]) = -4*gamma_W1*pow(volume/volume0 + 1, -2.0) + 4*gamma_W1 + prestress;
49+
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[2]) = 1 - pow((volume + volume0)/volume0, 0.66666666666666663);
50+
system.dC_dy.coeffRef(global_eqn_ids[0], global_var_ids[4]) = pow((volume + volume0)/volume0, 0.33333333333333331) - 1;
51+
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);
52+
system.dC_dy.coeffRef(global_eqn_ids[1], global_var_ids[6]) = 8.0*gamma_W1*pow((volume + volume0)/volume0, -3.0)/volume0;
6253
}
6354

6455
void ChamberSphere::get_elastance_values(std::vector<double>& parameters) {

‎src/model/ChamberSphere.h‎

Lines changed: 22 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -39,10 +39,9 @@
3939
P_\text{out} C = 0
4040
* \f]
4141
*
42-
* 2. Spherical stress:
42+
* 2. Spherical stress (shape correction factor \f$n = 1\f$, no viscosity):
4343
* \f[
44-
* -S + \tau + 4 (1 - C^{-3}) (W_1 + C W_2) + \eta \dot{C}
45-
* (1 + 2 C^{-6}) = 0
44+
* -S + \tau + 4 (1 - C^{-3}) W_1 = 0
4645
* \f]
4746
*
4847
* 3. Volume change:
@@ -80,18 +79,15 @@
8079
*
8180
* Parameter sequence for constructing this block:
8281
*
83-
* * `rho` - Density \f$\rho\f$
84-
* * `thick0` - Wall thickness \f$d_0\f$
85-
* * `radius0` - Reference radius \f$r_0\f$
86-
* * `W1` - Material constant \f$W_1\f$
87-
* * `W2` - Material constant \f$W_2\f$
88-
* * `eta` - Viscosity parameter \f$\eta\f$
89-
* * `sigma_max` - Maximum active stress \f$\sigma_\text{max}\f$
82+
* * `volume0` - Reference (unloaded) chamber volume \f$V_0\f$
83+
* * `gamma_W1` - Scaled material constant \f$\gamma W_1\f$
84+
* * `gamma_sigma_max` - Scaled maximum active stress \f$\gamma \sigma_\text{max}\f$
85+
* * `prestress` - Prestress
9086
* * `alpha_max` - Maximum activation parameter \f$\alpha_\text{max}\f$
9187
* * `alpha_min` - Minimum activation parameter \f$\alpha_\text{min}\f$
9288
* * `tsys` - Systole timing parameter \f$t_\text{sys}\f$
9389
* * `tdias` - Diastole timing parameter \f$t_\text{dias}\f$
94-
* * `steepness` - Activation steepness parameter \f$\gamma\f$
90+
* * `steepness` - Activation steepness parameter
9591
*
9692
* ### Usage in json configuration file
9793
*
@@ -103,13 +99,10 @@
10399
* "vessel_name": "ventricle",
104100
* "zero_d_element_type": "ChamberSphere",
105101
* "zero_d_element_values": {
106-
* "rho" : 1e3,
107-
* "thick0" : 0.01,
108-
* "radius0" : 0.05,
109-
* "W1" : 10e3,
110-
* "W2" : 40,
111-
* "eta" : 10.0,
112-
* "sigma_max" : 185e3,
102+
* "volume0" : 1e-4,
103+
* "gamma_W1" : 10e3,
104+
* "gamma_sigma_max" : 185e3,
105+
* "prestress" : 0.0,
113106
* "alpha_max": 30.0,
114107
* "alpha_min": -30.0,
115108
* "tsys": 0.170,
@@ -137,17 +130,15 @@ class ChamberSphere : public Block {
137130
*
138131
*/
139132
enum ParamId {
140-
n = 0,
141-
volume0 = 1,
142-
gamma_W1 = 2,
143-
gamma_eta = 3,
144-
gamma_sigma_max = 4,
145-
prestress = 5,
146-
alpha_max = 6,
147-
alpha_min = 7,
148-
tsys = 8,
149-
tdias = 9,
150-
steepness = 10
133+
volume0 = 0,
134+
gamma_W1 = 1,
135+
gamma_sigma_max = 2,
136+
prestress = 3,
137+
alpha_max = 4,
138+
alpha_min = 5,
139+
tsys = 6,
140+
tdias = 7,
141+
steepness = 8
151142
};
152143

153144
/**
@@ -158,10 +149,8 @@ class ChamberSphere : public Block {
158149
*/
159150
ChamberSphere(int id, Model* model)
160151
: Block(id, model, BlockType::chamber_sphere, BlockClass::vessel,
161-
{{"n", InputParameter()},
162-
{"volume0", InputParameter()},
152+
{{"volume0", InputParameter()},
163153
{"gamma_W1", InputParameter()},
164-
{"gamma_eta", InputParameter()},
165154
{"gamma_sigma_max", InputParameter()},
166155
{"prestress", InputParameter()},
167156
{"alpha_max", InputParameter()},
@@ -230,7 +219,7 @@ class ChamberSphere : public Block {
230219
* Number of triplets that the element contributes to the global system
231220
* (relevant for sparse memory reservation)
232221
*/
233-
TripletsContributions num_triplets{0, 0, 18};
222+
TripletsContributions num_triplets{9, 2, 4};
234223
};
235224

236225
#endif // SVZERODSOLVER_MODEL_ChamberSphere_HPP_

0 commit comments

Comments
 (0)