From 74ecbaa41b9d9418cd23b2ff58ea5e26f8b9d579 Mon Sep 17 00:00:00 2001
From: TS-CUBED <ts-cubed@t-cubed.org.uk>
Date: Sat, 23 Jul 2022 12:45:34 +0100
Subject: [PATCH 1/5] Use physiological parameters, incl. toggle for physUnits
 and SI

---
 haemoPimpleFoam/WKFunctions.C                 | 19 ++++++++++++-----
 .../{0 => constant}/windkesselProperties      | 21 +++++++++++++------
 2 files changed, 29 insertions(+), 11 deletions(-)
 rename tutorial_cases/carotid_transient_coarse_tet/{0 => constant}/windkesselProperties (80%)

diff --git a/haemoPimpleFoam/WKFunctions.C b/haemoPimpleFoam/WKFunctions.C
index 5764320..d19c274 100644
--- a/haemoPimpleFoam/WKFunctions.C
+++ b/haemoPimpleFoam/WKFunctions.C
@@ -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
@@ -51,6 +52,7 @@ void initialise(const dictionary& windkesselProperties)
         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"));
@@ -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*/
+        if (physUnits == true)
+            {
+            R = R * 133.322365e6;
+            Z = Z * 133.322365e6;
+            C = C * 1e-6 / 133.322365;
+            Info << "Converting C, R, Z from physioligical units to SI units!";
+            }
+
+        Info << "C, R, Z and index are " << C << ", " << R << ", " << Z << ", " << out_index << "." <<endl;
+
         wk[out_index].R = R;
         wk[out_index].C = C;
         wk[out_index].Z = Z;
+        wk[out_index].physUnits = physUnits;
         wk[out_index].id = out_index;
         wk[out_index].Q_3 = Qpre3;
         wk[out_index].Q_2 = Qpre2;
@@ -225,9 +237,6 @@ 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;
@@ -235,7 +244,7 @@ void execute_at_end(fvMesh& mesh, surfaceScalarField& phi, scalarIOList& store,
         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;
 
     
diff --git a/tutorial_cases/carotid_transient_coarse_tet/0/windkesselProperties b/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
similarity index 80%
rename from tutorial_cases/carotid_transient_coarse_tet/0/windkesselProperties
rename to tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
index 6e979b8..5a0cce1 100644
--- a/tutorial_cases/carotid_transient_coarse_tet/0/windkesselProperties
+++ b/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
@@ -16,9 +16,16 @@ FoamFile {
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 ICA {
-  C 10e-10;
-  R 15e+08;   // commonly R_2 or R_distal
-  Z 6e+07;    // commonly R_1 or R_proximal
+  C 0.5;
+  R 0.1;   // commonly R_2 or R_distal
+  Z 2.0;    // commonly R_1 or R_proximal
+
+  // 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.
+  physiologicalUnits true;
+
+
   outIndex 0; // must equal 'index' value in 0/p
   FDM_order 1;
 
@@ -37,9 +44,11 @@ ICA {
 }
 
 ECA {
-  C 10e-10;
-  R 15e+08;
-  Z 6e+07;
+  C 1.5;
+  R 0.2;
+  Z 0.6;
+  physiologicalUnits true;
+
   outIndex 1;
   FDM_order 1;
   Flowrate_threeStepBefore 0;

From 2b36e690de4700a8ab97e52f9f2c62d69e589ce3 Mon Sep 17 00:00:00 2001
From: TS-CUBED <ts-cubed@t-cubed.org.uk>
Date: Sat, 23 Jul 2022 12:46:20 +0100
Subject: [PATCH 2/5] Update windkesselProperties and move to constant

---
 .../constant/windkesselProperties                             | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties b/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
index 5a0cce1..4a1f806 100644
--- a/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
+++ b/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
@@ -28,9 +28,6 @@ ICA {
 
   outIndex 0; // must equal 'index' value in 0/p
   FDM_order 1;
-
-  // NOTE: ALL WK PARAMETERS ARE IN SI UNITS! Not in physiological units.
-
   // backward difference order: up to 3rd order
 
   // Initialise WK parameters
@@ -38,6 +35,7 @@ ICA {
   Flowrate_threeStepBefore 0;
   Flowrate_twoStepBefore 0;
   Flowrate_oneStepBefore 0;
+  // NOTE: these pressures are in [Pa]!
   Pressure_twoStepBefore 9000;
   Pressure_oneStepBefore 9000;
   Pressure_start 9000;

From f99814bfd68716295d60381bcec75d896022543d Mon Sep 17 00:00:00 2001
From: TS-CUBED <ts-cubed@t-cubed.org.uk>
Date: Sat, 23 Jul 2022 13:57:01 +0100
Subject: [PATCH 3/5] Update tutorial case to include more functionObjects

---
 .../system/controlDict                        | 133 ++++++++++++++++--
 1 file changed, 122 insertions(+), 11 deletions(-)

diff --git a/tutorial_cases/carotid_transient_coarse_tet/system/controlDict b/tutorial_cases/carotid_transient_coarse_tet/system/controlDict
index dcd14fb..90b849c 100644
--- a/tutorial_cases/carotid_transient_coarse_tet/system/controlDict
+++ b/tutorial_cases/carotid_transient_coarse_tet/system/controlDict
@@ -38,7 +38,7 @@ deltaT          1e-3;
 
 writeControl    runTime;
 
-writeInterval   1e-2;
+writeInterval   1e-3;
 
 purgeWrite      0;
 
@@ -56,8 +56,53 @@ runTimeModifiable true;
 
 functions
 {
+    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;
+    }
+
+    Mag
+    {
+        // Mandatory entries (unmodifiable)
+        type            mag;
+        libs            (fieldFunctionObjects);
+
+        // Mandatory (inherited) entries (runtime modifiable)
+        field           U; // wallShearStress;
+
+        // Optional (inherited) entries
+        result          wallShearStressMag;
+        region          region0;
+        enabled         true;
+        log             false;
+        timeStart       0;
+        timeEnd         1000;
+        executeControl  timeStep;
+        executeInterval 1;
+        writeControl    none;
+        writeInterval   1;
+    }
+
     //#include "surfaces"
-    CCA
+    averagesCCA
     {
         type            surfaceFieldValue;
         libs            ("libfieldFunctionObjects.so");
@@ -75,9 +120,11 @@ functions
                 p
                 U
                 H
+                wallShearStress
+                wallShearStressMag
             );
     }
-    ICA
+    averagesICA
     {
         type            surfaceFieldValue;
         libs            ("libfieldFunctionObjects.so");
@@ -95,9 +142,11 @@ functions
                 p
                 U
                 H
+                wallShearStress
+                wallShearStressMag
             );
     }
-    ECA
+    averagesECA
     {
         type            surfaceFieldValue;
         libs            ("libfieldFunctionObjects.so");
@@ -115,14 +164,76 @@ functions
                 p
                 U
                 H
+                wallShearStress
+                wallShearStressMag
             );
     }
-    //    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;
+   }
+
 }
 // ************************************************************************* //

From f668ab561bdbe20f83286c08cc724c80ddecb95a Mon Sep 17 00:00:00 2001
From: TS-CUBED <ts-cubed@t-cubed.org.uk>
Date: Sat, 23 Jul 2022 14:03:31 +0100
Subject: [PATCH 4/5] Fix WSS functionObjects

---
 .../system/controlDict                        | 69 ++++++++++++++++---
 1 file changed, 60 insertions(+), 9 deletions(-)

diff --git a/tutorial_cases/carotid_transient_coarse_tet/system/controlDict b/tutorial_cases/carotid_transient_coarse_tet/system/controlDict
index 90b849c..f838308 100644
--- a/tutorial_cases/carotid_transient_coarse_tet/system/controlDict
+++ b/tutorial_cases/carotid_transient_coarse_tet/system/controlDict
@@ -79,14 +79,14 @@ functions
     writeControl    outputTime;
     }
 
-    Mag
+    wallShearStressMag
     {
         // Mandatory entries (unmodifiable)
         type            mag;
         libs            (fieldFunctionObjects);
 
         // Mandatory (inherited) entries (runtime modifiable)
-        field           U; // wallShearStress;
+        field           wallShearStress;
 
         // Optional (inherited) entries
         result          wallShearStressMag;
@@ -101,7 +101,64 @@ functions
         writeInterval   1;
     }
 
-    //#include "surfaces"
+    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;
@@ -120,8 +177,6 @@ functions
                 p
                 U
                 H
-                wallShearStress
-                wallShearStressMag
             );
     }
     averagesICA
@@ -142,8 +197,6 @@ functions
                 p
                 U
                 H
-                wallShearStress
-                wallShearStressMag
             );
     }
     averagesECA
@@ -164,8 +217,6 @@ functions
                 p
                 U
                 H
-                wallShearStress
-                wallShearStressMag
             );
     }
 

From 31e64ccc2ac4c782227d88dea4026f76de030020 Mon Sep 17 00:00:00 2001
From: TS-CUBED <ts-cubed@t-cubed.org.uk>
Date: Sat, 23 Jul 2022 14:48:48 +0100
Subject: [PATCH 5/5] Rearrange Z, C, R to match 3EWK order (avoid confusion)

---
 haemoPimpleFoam/WKFunctions.C                        | 12 ++++++------
 .../constant/windkesselProperties                    | 11 +++++------
 2 files changed, 11 insertions(+), 12 deletions(-)

diff --git a/haemoPimpleFoam/WKFunctions.C b/haemoPimpleFoam/WKFunctions.C
index d19c274..0d48d00 100644
--- a/haemoPimpleFoam/WKFunctions.C
+++ b/haemoPimpleFoam/WKFunctions.C
@@ -49,9 +49,9 @@ 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"));
@@ -69,17 +69,17 @@ void initialise(const dictionary& windkesselProperties)
         /* e.g., last saved time: 80.00, running from 80.01 ---> 0: 80.01, -1: 80.00, -2: 79.99, -3: 79.98*/
         if (physUnits == true)
             {
-            R = R * 133.322365e6;
             Z = Z * 133.322365e6;
             C = C * 1e-6 / 133.322365;
-            Info << "Converting C, R, Z from physioligical units to SI units!";
+            R = R * 133.322365e6;
+            Info << "Converting Z, C, R from physioligical units to SI units!" <<endl;;
             }
 
-        Info << "C, R, Z and index are " << C << ", " << R << ", " << Z << ", " << out_index << "." <<endl;
+        Info << "Z, C, R and index are " << Z << ", " << C << ", " << R << ", " << out_index << "." <<endl;
 
-        wk[out_index].R = R;
-        wk[out_index].C = C;
         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;
diff --git a/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties b/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
index 4a1f806..c75ab65 100644
--- a/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
+++ b/tutorial_cases/carotid_transient_coarse_tet/constant/windkesselProperties
@@ -16,15 +16,14 @@ FoamFile {
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 ICA {
+  Z 0.1;    // commonly R_1 or R_proximal
   C 0.5;
-  R 0.1;   // commonly R_2 or R_distal
-  Z 2.0;    // commonly R_1 or R_proximal
+  R 2;   // commonly R_2 or R_distal
+  physiologicalUnits true;
 
   // 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.
-  physiologicalUnits true;
-
 
   outIndex 0; // must equal 'index' value in 0/p
   FDM_order 1;
@@ -42,9 +41,9 @@ ICA {
 }
 
 ECA {
+  Z 0.2;
   C 1.5;
-  R 0.2;
-  Z 0.6;
+  R 0.6;
   physiologicalUnits true;
 
   outIndex 1;