Merge pull request #11 from LIBRA-project/sensitivity_analysis

Merge all changes implemented during sensitivity analysis: graph view of the input, fix source normalization and T leakage due to dispersion, several sensitivity analysis and UQ scripts along with data processing notebooks

Author: hmonroedd

.github/workflows/ci.yml

@@ -0,0 +1,36 @@
+name: CI
+on: [pull_request, push]
+
+jobs:
+  run-tests:
+    name: run-tests
+    runs-on: ubuntu-latest
+
+    steps:
+    - name: Checkout code
+      uses: actions/checkout@v5
+
+    - name: Set up Conda
+      uses: conda-incubator/setup-miniconda@v3
+      with:
+        activate-environment: libra_sparging
+        environment-file: environment.yml
+        miniforge-version: latest
+        use-mamba: false
+        channels: conda-forge
+
+    - name: Install package
+      shell: bash -l {0}
+      run: |
+        python -m pip install -e .[dev]
+
+    - name: Run tests
+      shell: bash -l {0}
+      run: |
+        python -m pytest test --cov src/sparging --cov-report xml
+    
+    - name: Upload to codecov
+      uses: codecov/codecov-action@v4
+      with:
+        token: ${{ secrets.CODECOV_TOKEN }}
+        files: ./coverage.xml

.gitignore

@@ -56,7 +56,6 @@ cover/
 *.pot
 
 # Django stuff:
-*.log
 local_settings.py
 db.sqlite3
 db.sqlite3-journal
@@ -207,6 +206,16 @@ marimo/_lsp/
 __marimo__/
 
 .vscode/
+misc_scripts/
+datasets/
+figures/
 *.yaml
 *.csv
-mwe.py
+*.json
+!standard_input.json
+*.ipynb
+!sensitivity.ipynb
+!LL_sensitivity.ipynb
+*.code-workspace
+*.pkl
+*.html

README.md

@@ -10,12 +10,24 @@ To run the model:
 > [!NOTE]
 > This uses `dolfinx` which doesn't run on Linux. For windows users, consider using Windows Subsystem for Linux (WSL)
 
+To simply use the library, use the libra_sparging environment:
 ```
 conda env create -f environment.yml
 conda activate libra_sparging
 ```
 
+To also train a surrogate model (for sensitivity analysis, for example), you need autoemulate, use the autoemulate environment:
+```
+conda env create -f autoemulate_env.yml
+conda activate autoemulate_env
+```
+
+```
+python -m pip install -e .[dev]
+```
+
+## How to run tests
 
 ```
-python model.py
+python -m pytest test
 ```

autoemulate_env.yml

@@ -0,0 +1,17 @@
+name: autoemulate_env
+channels:
+  - conda-forge
+dependencies:
+  - python=3.12
+  - fenics-dolfinx
+  - matplotlib
+  - pyvista
+  - pyyaml
+  - numpy
+  - scipy
+  - pandas
+  - pint
+  - networkx
+  - pip:
+    - autoemulate
+    - pyvis

data_processing/LL_sensitivity.ipynb

@@ -2,10 +2,13 @@ name: libra_sparging
 channels:
   - conda-forge
 dependencies:
+  - python>=3.12
   - fenics-dolfinx
   - matplotlib
   - pyvista
   - pyyaml
   - numpy
   - scipy
-  - pandas
+  - pandas
+  - pint
+  - networkx

data_processing/sensitivity.ipynb

@@ -0,0 +1,156 @@
+from datetime import datetime
+from pathlib import Path
+from sparging import (
+    LIBRA_PI_GEOM,
+    LIBRA_PI_MAT,
+    LIBRA_PI_OPERATING_PARAMS,
+    LIBRA_PI_SPARGING_PARAMS,
+    SimulationInput,
+    ureg,
+    all_correlations,
+    CorrelationType,
+)
+import logging
+from autoemulate.simulations.base import Simulator
+from autoemulate import AutoEmulate
+from autoemulate.core.sensitivity_analysis import SensitivityAnalysis
+import torch
+import pandas as pd
+import json
+import networkx as nx
+
+COMPUTE_SOBOL = False
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.WARNING)
+
+FOLDER = Path("datasets") / datetime.now().strftime("%Y%m%d_%H%M%S")
+FOLDER.mkdir(exist_ok=True, parents=True)
+FOLDER_SAMPLES = FOLDER / "samples"
+FOLDER_SAMPLES.mkdir(exist_ok=True, parents=True)
+FOLDER_PP = FOLDER / "postprocessing"
+FOLDER_PP.mkdir(exist_ok=True)
+
+outputs = []
+h_l_corrs = all_correlations.get_list(CorrelationType.MASS_TRANSFER_COEFF)
+D_l_corrs = all_correlations.get_list(CorrelationType.DIFFUSIVITY)
+
+
+class SpargingProblem(Simulator):
+    def __init__(self, parameters_range, output_names):
+        self.counter = 0
+        super().__init__(parameters_range, output_names)
+
+    def _forward(self, x: torch.Tensor) -> torch.Tensor:
+        # construct simulation input
+        LIBRA_PI_OPERATING_PARAMS.temperature = x[0, 0].item() * ureg.celsius
+        LIBRA_PI_OPERATING_PARAMS.P_top = x[0, 1].item() * ureg.bar
+        LIBRA_PI_OPERATING_PARAMS.flow_g_mol = x[0, 2].item() * ureg.sccm
+        LIBRA_PI_GEOM.nozzle_diameter = x[0, 3].item() * ureg.m
+        # LIBRA_PI_SPARGING_PARAMS.h_l = h_l_corrs[int(x[0, 4].item())]
+        LIBRA_PI_MAT.D_l = D_l_corrs[int(x[0, 4].item())]
+        # breakpoint()
+
+        graph = nx.Graph()
+        sim_input = SimulationInput.from_parameters(
+            LIBRA_PI_GEOM,
+            LIBRA_PI_MAT,
+            LIBRA_PI_OPERATING_PARAMS,
+            LIBRA_PI_SPARGING_PARAMS,
+            graph=graph,
+        )
+        tau = sim_input.get_tau()
+        h_l = sim_input.h_l
+        a = sim_input.a
+        eps_g = sim_input.eps_g
+        sim_input.to_json(
+            FOLDER_SAMPLES / f"sample_{self.counter}.json"
+        )  # for debugging and postprocessing
+        self.counter += 1
+
+        # for post processing
+        Pi = sim_input.get_Pi_number().to("dimensionless").magnitude
+        PP_numbers.append(
+            [
+                Pi,
+                tau.to("s").magnitude,
+                h_l.to("m/s").magnitude,
+                a.to("1/m").magnitude,
+                eps_g.to("dimensionless").magnitude,
+                graph.nodes["d_b"]["value"].to("m").magnitude,
+                graph.nodes["u_g0"]["value"].to("m/s").magnitude,
+                graph.nodes["h_l"]["origin"],
+                graph.nodes["D_l"]["origin"],
+            ]
+        )
+        y = torch.tensor(
+            [[tau.to("s").magnitude]],
+            dtype=torch.float64,
+        )
+        return y
+
+
+simulator = SpargingProblem(
+    parameters_range={  # realistic (wide) parameters range for LIBRA
+        "temperature": (500, 600),  # celsius
+        "P_top": (1, 2),  # bar
+        "flow_g_mol": (100, 500),  # sccm
+        "nozzle_diameter": (1e-3, 5e-3),  # m
+        "D_l_corr": (0, 0),  # index for selecting diffusivity correlation
+    },
+    output_names=["tau"],
+)
+
+n_samples = 5000
+
+X = simulator.sample_inputs(n_samples)
+
+PP_numbers = []
+Y, _ = simulator.forward_batch(X, allow_failures=False)
+
+# save training data
+pd.DataFrame(Y, columns=simulator.output_names).to_csv(
+    FOLDER / "simulator_outputs.csv", index=False
+)
+pd.DataFrame(X, columns=simulator.param_names).to_csv(
+    FOLDER / "simulator_inputs.csv", index=False
+)
+pd.DataFrame(
+    PP_numbers,
+    columns=["Pi", "tau", "h_l", "a", "eps_g", "d_b", "u_g0", "h_l_corr", "D_l_corr"],
+).to_csv(FOLDER / "PP_data.csv", index=False)
+
+# sensitivity analysis problem
+problem = {
+    "num_vars": simulator.in_dim,
+    "names": simulator.param_names,
+    "bounds": simulator.param_bounds,
+    "output_names": simulator.output_names,
+}
+
+with open(FOLDER / "problem.json", "w") as f:
+    json.dump(problem, f, indent=4)
+
+if COMPUTE_SOBOL:
+    # Run AutoEmulate with default settings
+    ae = AutoEmulate(X, Y, log_level="WARNING", models=["GaussianProcessRBF"])
+    ae.summarise()
+
+    # pick best model
+    emulator = ae.best_result()
+    print(f"Selected model: {emulator.model_name} with id: {emulator.id}")
+
+    # The use_timestamp paramater ensures a new result is saved each time the save method is called
+    best_result_filepath = ae.save(emulator, FOLDER, use_timestamp=False)
+    print("Model and metadata saved to: ", best_result_filepath)
+
+    ae.plot_preds(
+        emulator,
+        output_names=simulator.output_names,
+        fname=FOLDER_PP / "predictions.png",
+    )
+
+    # === Sensitivity analysis ===
+    sa = SensitivityAnalysis(emulator.model, problem=problem)
+    sobol_df = sa.run("sobol")
+    sa.plot_sobol(sobol_df, index="ST", fname=FOLDER_PP / "sobolTot.png")

environment.yml

@@ -0,0 +1,116 @@
+from datetime import datetime
+from pathlib import Path
+
+from sparging import (
+    LIBRA_PI_GEOM,
+    LIBRA_PI_MAT,
+    LIBRA_PI_OPERATING_PARAMS,
+    LIBRA_PI_SPARGING_PARAMS,
+    SimulationInput,
+    ureg,
+    all_correlations,
+    CorrelationType,
+)
+
+import logging
+import pandas as pd
+import json
+import networkx as nx
+import numpy as np
+
+logger = logging.getLogger(__name__)
+logging.basicConfig(level=logging.WARNING)
+
+FOLDER = Path("datasets") / datetime.now().strftime("%Y%m%d_%H%M%S")
+FOLDER.mkdir(exist_ok=True, parents=True)
+FOLDER_SAMPLES = FOLDER / "samples"
+FOLDER_SAMPLES.mkdir(exist_ok=True, parents=True)
+FOLDER_PP = FOLDER / "postprocessing"
+FOLDER_PP.mkdir(exist_ok=True)
+
+outputs = []
+h_l_corrs = all_correlations.get_list(CorrelationType.MASS_TRANSFER_COEFF)
+D_l_corrs = all_correlations.get_list(CorrelationType.DIFFUSIVITY)
+
+
+def forward(count, x: dict) -> list:
+    """return list of values for post processing"""
+    # construct simulation input
+    LIBRA_PI_OPERATING_PARAMS.temperature = x["temperature"] * ureg.celsius
+    LIBRA_PI_OPERATING_PARAMS.P_top = x["P_top"] * ureg.bar
+    LIBRA_PI_OPERATING_PARAMS.flow_g_mol = x["flow_g_mol"] * ureg.sccm
+    LIBRA_PI_GEOM.nozzle_diameter = x["nozzle_diameter"] * ureg.m
+    # LIBRA_PI_SPARGING_PARAMS.h_l = x["h_l_corr"]
+    LIBRA_PI_MAT.D_l = D_l_corrs[x["D_l_corr_number"]]
+    # breakpoint()
+
+    graph = nx.Graph()
+    sim_input = SimulationInput.from_parameters(
+        LIBRA_PI_GEOM,
+        LIBRA_PI_MAT,
+        LIBRA_PI_OPERATING_PARAMS,
+        LIBRA_PI_SPARGING_PARAMS,
+        graph=graph,
+    )
+    tau = sim_input.get_tau()
+    h_l = sim_input.h_l
+    a = sim_input.a
+    eps_g = sim_input.eps_g
+    # sim_input.to_json(
+    #     FOLDER_SAMPLES / f"sample_{count}.json"
+    # )  # for debugging and postprocessing
+
+    # for post processing
+    Pi = sim_input.get_Pi_number().to("dimensionless").magnitude
+    return [
+        Pi,
+        tau.to("s").magnitude,
+        h_l.to("m/s").magnitude,
+        a.to("1/m").magnitude,
+        eps_g.to("dimensionless").magnitude,
+        graph.nodes["d_b"]["value"].to("m").magnitude,
+        graph.nodes["u_g0"]["value"].to("m/s").magnitude,
+        graph.nodes["h_l"]["origin"],
+        graph.nodes["D_l"]["origin"],
+    ]
+
+
+n_samples = 1000
+inputs = []
+postprocess = []
+
+for i in range(n_samples):
+    input = {
+        "temperature": np.random.normal(550, 25),  # celsius
+        "P_top": 1.2,  # bar
+        "flow_g_mol": 400,  # sccm
+        "nozzle_diameter": 1.5e-3,  # m
+        "D_l_corr_number": np.random.randint(
+            0, 2
+        ),  # index for selecting diffusivity correlation
+    }
+    inputs.append(input)
+    print(f"Running simulation {i + 1}/{n_samples}")
+    postprocess.append(forward(i, input))
+
+# save training data
+PP_frame = pd.DataFrame(
+    postprocess,
+    columns=["Pi", "tau", "h_l", "a", "eps_g", "d_b", "u_g0", "h_l_corr", "D_l_corr"],
+)
+PP_frame.to_csv(FOLDER / "PP_data.csv", index=False)
+
+# for compatibility with postprocessing code
+pd.DataFrame(PP_frame, columns=["tau"]).to_csv(
+    FOLDER / "simulator_outputs.csv", index=False
+)
+input_frame = pd.DataFrame(inputs)
+input_frame.to_csv(FOLDER / "simulator_inputs.csv", index=False)
+
+problem = {
+    "names": [x for x in inputs[0].keys()],
+    "bounds": ["550 +- 25", "1.2", "400", "1.5e-3", "oishi or calderoni"],
+}
+
+with open(FOLDER / "problem.json", "w") as f:
+    json.dump(problem, f, indent=4)