Skip to content

Commit

Permalink
Merge pull request #29 from TS-CUBED/physiologicalWK-params
Browse files Browse the repository at this point in the history 
Physiological wk params
  • Loading branch information
@TS-CUBED
TS-CUBED authored Jul 23, 2022
2 parents 6c41c70 + 0f5893a commit 5de03a0
Show file tree
Hide file tree
Showing 3 changed files with 206 additions and 29 deletions.
25 changes: 17 additions & 8 deletions haemoPimpleFoam/WKFunctions.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ typedef struct
double R; /* Resistance */
double C; /* Compliance */
double Z; /* Impedance */
bool physUnits; /* Physiological or SI units for R, C, Z */
int order;
double time; // store time in windkessel struct to reset values
// should the time step be run multiple times
Expand Down Expand Up @@ -48,9 +49,10 @@ void initialise(const dictionary& windkesselProperties)

const dictionary& subDict = windkesselProperties.subDict(outletName);

scalar Z = readScalar(subDict.lookup("Z"));
scalar C = readScalar(subDict.lookup("C"));
scalar R = readScalar(subDict.lookup("R"));
scalar Z = readScalar(subDict.lookup("Z"));
bool physUnits = readBool(subDict.lookup("physiologicalUnits"));
scalar real_index = readScalar(subDict.lookup("outIndex"));
scalar Qpre1 = readScalar(subDict.lookup("Flowrate_oneStepBefore"));
scalar Qpre2 = readScalar(subDict.lookup("Flowrate_twoStepBefore"));
Expand All @@ -64,11 +66,21 @@ void initialise(const dictionary& windkesselProperties)
int order = FDMorder;


Info << "C, R, Z and index are " << C << ", " << R << ", " << Z << ", " << out_index << "." <<endl;
/* e.g., last saved time: 80.00, running from 80.01 ---> 0: 80.01, -1: 80.00, -2: 79.99, -3: 79.98*/
wk[out_index].R = R;
wk[out_index].C = C;
if (physUnits == true)
{
Z = Z * 133.322365e6;
C = C * 1e-6 / 133.322365;
R = R * 133.322365e6;
Info << "Converting Z, C, R from physioligical units to SI units!" <<endl;;
}

Info << "Z, C, R and index are " << Z << ", " << C << ", " << R << ", " << out_index << "." <<endl;

wk[out_index].Z = Z;
wk[out_index].C = C;
wk[out_index].R = R;
wk[out_index].physUnits = physUnits;
wk[out_index].id = out_index;
wk[out_index].Q_3 = Qpre3;
wk[out_index].Q_2 = Qpre2;
Expand Down Expand Up @@ -225,17 +237,14 @@ void execute_at_end(fvMesh& mesh, surfaceScalarField& phi, scalarIOList& store,
WK_dict.set("Pressure_start", wk[i].P1);

/* update */
// if (wk[i].time != time)
// {
// Info << "time step advanced, updating WK properties" << endl;
wk[i].P_2 = wk[i].P_1;
wk[i].P_1 = wk[i].P_0;
wk[i].P_0 = wk[i].P1;
wk[i].Q_3 = wk[i].Q_2;
wk[i].Q_2 = wk[i].Q_1;
wk[i].Q_1 = wk[i].Q_0;
wk[i].time = time;
// }

Info << "P(n) = " << wk[i].P1 << endl;


Expand Down
24 changes: 15 additions & 9 deletions ...ansient_coarse_tet/0/windkesselProperties → ..._coarse_tet/constant/windkesselProperties
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,30 +16,36 @@ FoamFile {
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ICA {
C 10e-10;
R 15e+08; // commonly R_2 or R_distal
Z 6e+07; // commonly R_1 or R_proximal
outIndex 0; // must equal 'index' value in 0/p
FDM_order 1;
Z 0.1; // commonly R_1 or R_proximal
C 0.5;
R 2; // commonly R_2 or R_distal
physiologicalUnits true;

// NOTE: ALL WK PARAMETERS ARE IN SI UNITS! Not in physiological units.
// NOTE: Set this to true to use physiological units [mmHg.ml^-1.s] and [ml.mmHg^-1]
// if this is false, then SI units are used. Physiological units are converted
// at runtime.

outIndex 0; // must equal 'index' value in 0/p
FDM_order 1;
// backward difference order: up to 3rd order

// Initialise WK parameters
// also useful to set initial condidions to speed up pressure development
Flowrate_threeStepBefore 0;
Flowrate_twoStepBefore 0;
Flowrate_oneStepBefore 0;
// NOTE: these pressures are in [Pa]!
Pressure_twoStepBefore 9000;
Pressure_oneStepBefore 9000;
Pressure_start 9000;
}

ECA {
C 10e-10;
R 15e+08;
Z 6e+07;
Z 0.2;
C 1.5;
R 0.6;
physiologicalUnits true;

outIndex 1;
FDM_order 1;
Flowrate_threeStepBefore 0;
Expand Down
186 changes: 174 additions & 12 deletions tutorial_cases/carotid_transient_coarse_tet/system/controlDict
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ deltaT 1e-3;

writeControl runTime;

writeInterval 1e-2;
writeInterval 1e-3;

purgeWrite 0;

Expand All @@ -56,8 +56,110 @@ runTimeModifiable true;

functions
{
//#include "surfaces"
CCA
wallShearStress
{
// Mandatory entries (unmodifiable)
type wallShearStress;
libs (fieldFunctionObjects);

// Optional entries (runtime modifiable)
// patches (WALL APEX SINUS); // (wall1 "(wall2|wall3)");

// Optional (inherited) entries
// writePrecision 8;
// writeToFile true;
// useUserTime true;
// region region0;
// enabled true;
// log true;
// timeStart 0;
// timeEnd 1000;
// executeControl timeStep;
// executeInterval 1;
writeControl outputTime;
}

wallShearStressMag
{
// Mandatory entries (unmodifiable)
type mag;
libs (fieldFunctionObjects);

// Mandatory (inherited) entries (runtime modifiable)
field wallShearStress;

// Optional (inherited) entries
result wallShearStressMag;
region region0;
enabled true;
log false;
timeStart 0;
timeEnd 1000;
executeControl timeStep;
executeInterval 1;
writeControl none;
writeInterval 1;
}

averagesWALL
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
log yes;
writeFields no;
regionType patch;
name WALL;
operation average;
// operation weightedAverage;
// weightField phi;

fields
(
wallShearStress
wallShearStressMag
);
}
averagesAPEX
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
log yes;
writeFields no;
regionType patch;
name APEX;
operation average;
// operation weightedAverage;
// weightField phi;

fields
(
wallShearStress
wallShearStressMag
);
}
averagesSINUS
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
log yes;
writeFields no;
regionType patch;
name SINUS;
operation average;
// operation weightedAverage;
// weightField phi;

fields
(
wallShearStress
wallShearStressMag
);
}

averagesCCA
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
Expand All @@ -77,7 +179,7 @@ functions
H
);
}
ICA
averagesICA
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
Expand All @@ -97,7 +199,7 @@ functions
H
);
}
ECA
averagesECA
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
Expand All @@ -117,12 +219,72 @@ functions
H
);
}
// residuals
// {
// type residuals;
// libs ("libutilityFunctionObjects.so");
// fields (U p);
// writeResidualFields yes;
// }

flowRateCCA
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
enabled yes;
writeControl timeStep;
writeInterval 1;
log yes;
writeFields no;
regionType patch;
name CCA;
operation sum;

fields
(
phi
);
}

flowRateICA
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
enabled yes;
writeControl timeStep;
writeInterval 1;
log yes;
writeFields no;
regionType patch;
name ICA;
operation sum;

fields
(
phi
);
}

flowRateECA
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
enabled yes;
writeControl timeStep;
writeInterval 1;
log yes;
writeFields no;
regionType patch;
name ECA;
operation sum;

fields
(
phi
);
}

solverInfo
{
type solverInfo;
libs ("libutilityFunctionObjects.so");
fields (U H p);
timeStart 3.01;
writeResidualFields no;
}

}
// ************************************************************************* //

0 comments on commit 5de03a0

Please sign in to comment.