Add scripts for running sensitivity analysis

demos/stiffness_sensitivity.py

@@ -0,0 +1,76 @@
+# # Simple demo
+#
+# In this demo we show the most simple usage of the `simcardems` library using the python API.
+#
+# Import the necessary libraries
+#
+
+from pathlib import Path
+import simcardems
+import matplotlib.pyplot as plt
+import numpy as np
+
+here = Path(__file__).absolute().parent
+geometry_path = here / "geometries/slab.h5"
+geometry_schema_path = here / "geometries/slab.json"
+
+
+def run(outdir: Path, a: float = 2.28):
+    outdir.mkdir(exist_ok=True, parents=True)
+    results_file = outdir.joinpath("results.h5")
+    if results_file.exists():
+        values = np.load(outdir.joinpath("values.npy"), allow_pickle=True).item()
+
+        t = values["time"]
+        lmbda = values["mechanics"]["lambda"]
+        return t, lmbda
+    config = simcardems.Config(
+        outdir=outdir,
+        geometry_path=geometry_path,
+        geometry_schema_path=geometry_schema_path,
+        coupling_type="fully_coupled_ORdmm_Land",
+        T=1000,
+        material_parameter_a=a,
+    )
+    runner = simcardems.Runner(config)
+    runner.solve(T=config.T, save_freq=config.save_freq, show_progress_bar=True)
+
+    loader = simcardems.datacollector.DataLoader(results_file)
+    values = simcardems.postprocess.extract_traces(loader, reduction="center")
+    np.save(outdir.joinpath("values.npy"), values)
+
+    t = values["time"]
+    lmbda = values["mechanics"]["lambda"]
+    fig, ax = plt.subplots()
+    ax.plot(t, lmbda)
+    ax.set_xlabel("Time [ms]")
+    ax.set_ylabel("Lambda")
+    ax.set_title(f"Material parameter a = {a}")
+    fig.savefig(outdir.joinpath(f"lambda_{a}.png"))
+
+    return t, lmbda
+    # simcardems.postprocess.plot_state_traces(outdir.joinpath("results.h5"), "center")
+
+
+def main():
+    outdir = here / "results_stiffness_sensitivity"
+    fig, ax = plt.subplots()
+    a_default = 2.28
+    for a in [
+        a_default / 4.0,
+        a_default / 2.0,
+        a_default,
+        a_default * 2.0,
+        a_default * 4.0,
+    ]:
+        t, lmbda = run(outdir=outdir / f"a_{a}", a=a)
+        ax.plot(t, lmbda, label=f"a = {a}")
+    ax.legend()
+    ax.set_xlabel("Time [ms]")
+    ax.set_ylabel("Lambda")
+    ax.set_title("Stiffness sensitivity")
+    fig.savefig(outdir / "stiffness_sensitivity.png")
+
+
+if __name__ == "__main__":
+    main()

demos/timestep_sensitivity.py

@@ -0,0 +1,76 @@
+# # Simple demo
+#
+# In this demo we show the most simple usage of the `simcardems` library using the python API.
+#
+# Import the necessary libraries
+#
+
+from pathlib import Path
+import simcardems
+import matplotlib.pyplot as plt
+import numpy as np
+
+here = Path(__file__).absolute().parent
+geometry_path = here / "geometries/slab.h5"
+geometry_schema_path = here / "geometries/slab.json"
+
+
+def run(outdir: Path, dt_mech: float = 2.28):
+    outdir.mkdir(exist_ok=True, parents=True)
+    results_file = outdir.joinpath("results.h5")
+    if results_file.exists():
+        values = np.load(outdir.joinpath("values.npy"), allow_pickle=True).item()
+
+        t = values["time"]
+        lmbda = values["mechanics"]["lambda"]
+        return t, lmbda
+    config = simcardems.Config(
+        outdir=outdir,
+        geometry_path=geometry_path,
+        geometry_schema_path=geometry_schema_path,
+        coupling_type="fully_coupled_ORdmm_Land",
+        T=1000,
+        dt_mech=dt_mech,
+    )
+    runner = simcardems.Runner(config)
+    runner.solve(T=config.T, save_freq=config.save_freq, show_progress_bar=True)
+
+    loader = simcardems.datacollector.DataLoader(results_file)
+    values = simcardems.postprocess.extract_traces(loader, reduction="center")
+    np.save(outdir.joinpath("values.npy"), values)
+
+    t = values["time"]
+    lmbda = values["mechanics"]["lambda"]
+    fig, ax = plt.subplots()
+    ax.plot(t, lmbda)
+    ax.set_xlabel("Time [ms]")
+    ax.set_ylabel("Lambda")
+    ax.set_title(f"Timestep = {dt_mech}")
+    fig.savefig(outdir.joinpath(f"lambda_{dt_mech}.png"))
+
+    return t, lmbda
+    # simcardems.postprocess.plot_state_traces(outdir.joinpath("results.h5"), "center")
+
+
+def main():
+    outdir = here / "results_timestep_sensitivity"
+    fig, ax = plt.subplots()
+    dt_mech_default = 1.0
+    for dt in [
+        dt_mech_default / 4.0,
+        dt_mech_default / 2.0,
+        dt_mech_default,
+        dt_mech_default * 2.0,
+        dt_mech_default * 4.0,
+    ]:
+        t, lmbda = run(outdir=outdir / f"dt_{dt}", dt_mech=dt)
+        ax.plot(t, lmbda, label=f"dt = {dt}")
+    ax.legend()
+    ax.set_xlabel("Time [ms]")
+    ax.set_ylabel("Lambda")
+    ax.set_title("Timestep sensitivity")
+    fig.savefig(outdir / "timestep_sensitivity.png")
+
+
+if __name__ == "__main__":
+    main()

demos/timestep_threshold_sensitivity.py

@@ -0,0 +1,76 @@
+# # Simple demo
+#
+# In this demo we show the most simple usage of the `simcardems` library using the python API.
+#
+# Import the necessary libraries
+#
+
+from pathlib import Path
+import simcardems
+import matplotlib.pyplot as plt
+import numpy as np
+
+here = Path(__file__).absolute().parent
+geometry_path = here / "geometries/slab.h5"
+geometry_schema_path = here / "geometries/slab.json"
+
+
+def run(outdir: Path, mech_threshold: float = 0.05):
+    outdir.mkdir(exist_ok=True, parents=True)
+    results_file = outdir.joinpath("results.h5")
+    if results_file.exists():
+        values = np.load(outdir.joinpath("values.npy"), allow_pickle=True).item()
+
+        t = values["time"]
+        lmbda = values["mechanics"]["lambda"]
+        return t, lmbda
+    config = simcardems.Config(
+        outdir=outdir,
+        geometry_path=geometry_path,
+        geometry_schema_path=geometry_schema_path,
+        coupling_type="fully_coupled_ORdmm_Land",
+        T=1000,
+        mech_threshold=mech_threshold,
+    )
+    runner = simcardems.Runner(config)
+    runner.solve(T=config.T, save_freq=config.save_freq, show_progress_bar=True)
+
+    loader = simcardems.datacollector.DataLoader(results_file)
+    values = simcardems.postprocess.extract_traces(loader, reduction="center")
+    np.save(outdir.joinpath("values.npy"), values)
+
+    t = values["time"]
+    lmbda = values["mechanics"]["lambda"]
+    fig, ax = plt.subplots()
+    ax.plot(t, lmbda)
+    ax.set_xlabel("Time [ms]")
+    ax.set_ylabel("Lambda")
+    ax.set_title(f"Threshold = {mech_threshold}")
+    fig.savefig(outdir.joinpath(f"lambda_{mech_threshold}.png"))
+
+    return t, lmbda
+    # simcardems.postprocess.plot_state_traces(outdir.joinpath("results.h5"), "center")
+
+
+def main():
+    outdir = here / "results_timestep_threshold_sensitivity"
+    fig, ax = plt.subplots()
+    threshold_default = 0.05
+    for thr in [
+        threshold_default / 4.0,
+        threshold_default / 2.0,
+        threshold_default,
+        threshold_default * 2.0,
+        threshold_default * 4.0,
+    ]:
+        t, lmbda = run(outdir=outdir / f"thr_{thr}", mech_threshold=thr)
+        ax.plot(t, lmbda, label=f"thr = {thr}")
+    ax.legend()
+    ax.set_xlabel("Time [ms]")
+    ax.set_ylabel("Lambda")
+    ax.set_title("Timestep sensitivity")
+    fig.savefig(outdir / "timestep_threshold_sensitivity.png")
+
+
+if __name__ == "__main__":
+    main()

src/simcardems/config.py

@@ -54,6 +54,7 @@ class Config:
         "explicit_ORdmm_Land",
         "pureEP_ORdmm_Land",
     ] = "fully_coupled_ORdmm_Land"
+    material_parameter_a: float = 2.28
 
     def as_dict(self):
         return {k: v for k, v in self.__dict__.items()}

src/simcardems/mechanics_model.py

@@ -34,6 +34,7 @@ def setup_solver(
     linear_solver="mumps",
     use_custom_newton_solver: bool = config.Config.mechanics_use_custom_newton_solver,
     state_prev=None,
+    material_parameter_a: float = 2.28,
 ):
     """Setup mechanics model with dirichlet boundary conditions or rigid motion."""
 
@@ -42,8 +43,9 @@ def setup_solver(
     logger.info("Set up mechanics model")
 
     # Use parameters from Biaxial test in Holzapfel 2019 (Table 1).
+
     material_parameters = dict(
-        a=2.28,
+        a=material_parameter_a,
         a_f=1.686,
         b=9.726,
         b_f=15.779,

src/simcardems/models/em_model.py

@@ -70,6 +70,7 @@ def setup_EM_model(
         use_custom_newton_solver=config.mechanics_use_custom_newton_solver,
         debug_mode=config.debug_mode,
         ActiveModel=cls_ActiveModel,
+        material_parameter_a=config.material_parameter_a,
     )
     if mech_state_init is not None:
         mech_heart.state.assign(mech_state_init)