up script for DSA

CITATION.cff

@@ -14,7 +14,7 @@ title: "Model order reduction and sensitivity analysis for complex heat transfer
 version: 1.0.0
 identifiers:
   - type: doi
-    value: 
+    value:
 date-released: 2024-09-06
 url: "https://github.com/feelpp/article.eye-heat-fom-rom-sa.ijnmbe24"
 preferred-citation:
@@ -30,7 +30,7 @@ preferred-citation:
     given-names: "Marcela"
     orcid:  https://orcid.org/0000-0002-7300-3267
   doi: 10.1002/cnm.3864
-  journal: "International Journal for Numerical Methods in Biomedical Engineering"  
+  journal: "International Journal for Numerical Methods in Biomedical Engineering"
   title: "Model order reduction and sensitivity analysis for complex heat transfer simulations inside the human eyeball"
   year: 2024
   month: 9

README.adoc

@@ -99,16 +99,20 @@ To use another mesh, update the content of the section `Meshes/heat/Import` in t
 NOTE: When running in parallel, make sure that the mesh is partitioned correctly. In the JSON file, the macro `$np` will be replaced by the number of MPI processes.
 
 
-=== Results #2
+=== Run the deterministic sensitivity analysis
 
-how to reproduce ...
+To reproduce the results of the deterministic sensitivity analysis, presented in section 4.1 of the pulication, run
+
+[source, bash]
+----
+[mpirun -np 12] python3 run-SA.py
+----
 
 
 == References
 
 * [[fpp24]] "Heat Transfer Toolbox Application Manual." Feel++ Documentation, 2024, https://docs.feelpp.org/toolboxes/latest/heat/index.html.
-* [[[pp]]] Andy Hunt & Dave Thomas. The Pragmatic Programmer:
-From Journeyman to Master. Addison-Wesley. 1999.
+
 
 == Acknowledgements
 

data/crb_param.json

@@ -17,15 +17,15 @@
         },
         "T_bl":    // [K]         \cite{https://hypertextbook.com/facts/1997/LenaWong.shtml}
         {
-            "min":"308.3",
-            "max":"312",
+            "min":"308.15",
+            "max":"312.15",
             "description":"Blood tempetature [K]",
             "source": "pm 5% of 310.15K = 37°C"
         },
         "T_amb":    // [K]         \cite{MAPSTONE1968237}
         {
-            "min":"283.15",
-            "max":"303.15",
+            "min":"273.15",
+            "max":"308.15",
             "description":"Ambiant temperature [K]",
             "source": "de 10°C à 30°C"
         },

data/deterministic-sensitivity-analysis/points_coord.csv

@@ -0,0 +1,10 @@
+Point,x,y,z
+O,-0.0136,0,0
+A,-0.0126,0,0
+B,-0.00979266,0,0
+B1,-0.009524,0.00194062,0
+C,-0.00615819,0,0
+D,0.0106,0,0
+D1,0.011,0,0
+F,0.0115,0,0
+G,0.0125,0,0

data/deterministic-sensitivity-analysis/results/.gitignore

@@ -0,0 +1 @@
+*.csv

data/deterministic-sensitivity-analysis/run-SA.py

@@ -0,0 +1,103 @@
+import sys, os
+import feelpp.core as fppc
+import feelpp.toolboxes.heat as heat
+import feelpp.mor as mor
+import pandas as pd
+
+PWD = os.getcwd()
+
+o = heat.toolboxes_options("heat")
+config = fppc.globalRepository("SA_det")
+e = fppc.Environment(sys.argv, opts=o, config=config)
+
+DIM = 3
+assert DIM in [2,3]
+
+
+fppc.Environment.setConfigFile(os.path.join(PWD, "..", "eye.cfg"))
+
+MEASURES_PATH = f"{e.rootRepository()}/SA_det/np_{e.expand('$np')}/heat.measures/values.csv"
+
+
+def updateParameters(tb, mu):
+    for i in range(0,mu.size()):
+        tb.addParameterInModelProperties(mu.parameterName(i),mu(i))
+    tb.updateParameterValues()
+
+def value(tb, param):
+    return tb.modelProperties().parameters().at(param).value()
+
+def printParam(tb, names):
+    param = tb.modelProperties().parameters()
+    for name in names:
+        print(name, param.at(name).value())
+
+
+def to_kelvin(T_C): return T_C + 273.15
+def to_celcius(T):  return T - 273.15
+
+def res_of_meas(meas, points):
+    res = {}
+    for point in points:
+        res[point] = meas[f'Points_{point}_expr_T{point}_C']
+    return res
+
+heatBox = heat.heat(dim=DIM, order=2)
+heatBox.init()
+
+
+crb_model_properties = mor.CRBModelProperties(worldComm=fppc.Environment.worldCommPtr())
+crb_model_properties.setup(PWD + '/../crb_param.json')
+crb_model_parameters = crb_model_properties.parameters()
+
+# modelParameters = heatBox.modelProperties().parameters()
+D = mor._mor.ParameterSpace.New(crb_model_parameters, fppc.Environment.worldCommPtr())
+
+control_values = {"k_lens":0.40, "E":40, "h_bl":65, "h_amb":10, "T_amb":to_kelvin(20), "T_bl":to_kelvin(37)}
+mu = D.element()
+mu.setParameters(control_values)
+
+if e.isMasterRank():
+    print(mu.parameterNames())
+
+
+params = D.parameterNames()
+values = {
+    "E": [20, 40, 70, 100, 320],
+    "T_amb": [273.15, 278, to_kelvin(20), to_kelvin(25), to_kelvin(30), 308],
+    "T_bl": [to_kelvin(35), to_kelvin(37), to_kelvin(37.7), to_kelvin(38), to_kelvin(38.5), to_kelvin(39)],
+    "h_amb": [8, 10, 12, 15, 100],
+    "h_bl": [65, 90, 110],
+    "k_lens": [0.21, 0.30, 0.40, 0.544]
+}
+
+df_points = pd.read_csv(os.path.join(PWD, "points_coord.csv"))
+points_name = list(df_points['Point'])
+
+
+for param in params:
+    if e.isMasterRank():
+        print("=================================")
+        print("Running for the parameter", param)
+        print("param will vary in ", values[param])
+    mu.setParameters(control_values)
+
+    df_res = pd.DataFrame(columns=points_name, dtype=float)
+    for value in values[param]:
+        try:
+            mu.setParameterNamed(param, value)
+        except RuntimeError as e:
+            print(e)
+            sys.exit(1)
+        if e.isMasterRank():
+            print(f"  Running with {param} = {value}")
+        updateParameters(heatBox, mu)
+        heatBox.solve()
+        heatBox.exportResults()
+        meas = heatBox.postProcessMeasures().values()
+
+        if param[0] == 'T':
+            value = to_celcius(value)
+        df_res.loc[value] = res_of_meas(meas, points_name)
+
+    df_res.to_csv(os.path.join(PWD, "results", f"{param}_feel.csv"), index_label=param)

data/deterministic-sensitivity-analysis/run-SA.slurm

@@ -0,0 +1,9 @@
+#!/bin/bash
+#SBATCH -J DSA                  # name of the job
+#SBATCH -N 1                    # number of nodes
+#SBATCH -n 12                   # number of cores
+#SBATCH -t 06:45:00             # job duration
+#SBATCH -o log/DSA%j_out.log    # standard output
+#SBATCH -e log/DSA%j_err.log    # standard error
+
+mpiexec -bind-to core python3 run-SA.py

data/eye.cfg

@@ -13,9 +13,9 @@ json.filename=$cfgdir/boundary-conditions.json
 json.filename=$cfgdir/post-process.json
 solver=Linear #Newton #Oseen,Picard,Newton
 ksp-rtol=1e-15
-ksp-monitor=1
+ksp-monitor=0
 pc-type=gamg
-snes-monitor=1
+snes-monitor=0
 snes-ksp-rtol=1e-10