diff --git a/two-scale-heat-conduction/micro-dumux-mada/mada_switcher.py b/two-scale-heat-conduction/micro-dumux-mada/mada_switcher.py new file mode 100644 index 000000000..84fd05713 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-mada/mada_switcher.py @@ -0,0 +1,20 @@ +import numpy as np + + +def switching_function(resolution, location, t, input, prev_output): + result = 0 + # in the beginning we only want FOM + if t == 0.0: + if resolution > 0: result = -1 + else: result = 0 + # after small init phase, we want dynamic model selection + # for test purposes we say ROM is accurate in range + else: + concentration = input['concentration'] + is_valid_range = 0.45 < concentration < 0.55 + is_fom = resolution == 0 + + if is_fom and is_valid_range: result = 1 + if not is_fom and not is_valid_range: result = -1 + + return result diff --git a/two-scale-heat-conduction/micro-dumux-mada/micro-manager-mada-config.json b/two-scale-heat-conduction/micro-dumux-mada/micro-manager-mada-config.json new file mode 100644 index 000000000..7ccf5e012 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-mada/micro-manager-mada-config.json @@ -0,0 +1,35 @@ +{ + "micro_file_name": "micro_sim", + "coupling_params": { + "participant_name": "Micro-Manager", + "precice_config_file_name": "../precice-config.xml", + "macro_mesh_name": "macro-mesh", + "write_data_names": ["k_00", "k_01", "k_10", "k_11", "porosity"], + "read_data_names": ["concentration"] + }, + "simulation_params": { + "micro_dt": 0.01, + "macro_domain_bounds": [0.0, 1.0, 0.0, 0.5], + "decomposition": [1, 1], + "adaptivity": "True", + "adaptivity_settings": { + "type": "global", + "data": ["k_00", "k_11", "porosity", "concentration"], + "history_param": 0.1, + "coarsening_constant": 0.2, + "refining_constant": 0.05, + "every_implicit_iteration": "False", + "similarity_measure": "L2rel" + }, + "model_adaptivity": true, + "model_adaptivity_settings": { + "micro_file_names": ["micro_sim", "micro_sim_sur"], + "data": [], + "thresholds": [1.0, 1.0], + "switching_function": "mada_switcher" + } + }, + "diagnostics": { + "data_from_micro_sims": ["grain_size"] + } +} diff --git a/two-scale-heat-conduction/micro-dumux-mada/micro_sim.pc.in b/two-scale-heat-conduction/micro-dumux-mada/micro_sim.pc.in new file mode 100644 index 000000000..566aba5e4 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-mada/micro_sim.pc.in @@ -0,0 +1,15 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ +CXX=@CXX@ +CC=@CC@ +DEPENDENCIES=@REQUIRES@ + +Name: @PACKAGE_NAME@ +Version: @VERSION@ +Description: micro_sim module +URL: http://dune-project.org/ +Requires: dumux-phasefield dumux-precice +Libs: -L${libdir} +Cflags: -I${includedir} diff --git a/two-scale-heat-conduction/micro-dumux-mada/micro_sim_sur.py b/two-scale-heat-conduction/micro-dumux-mada/micro_sim_sur.py new file mode 100644 index 000000000..d2502c60e --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-mada/micro_sim_sur.py @@ -0,0 +1,59 @@ +""" +Micro simulation Surrogate, requrie previous computation of surrogate model +""" +import os +import subprocess +from bayesvalidrox import PyLinkForwardModel, Input, PCE, ExpDesigns, Engine +import h5py +import joblib +import numpy as np +import math + + +class MicroSimulation: + + def __init__(self, sim_id): + """ + Constructor of MicroSimulation class. + """ + self._sim_id = sim_id + self._state = None + + self._model = None + with open('micro-dumux-surrogate.pkl', 'rb') as input: + self._model = joblib.load(input) + if self._model is None: + raise RuntimeError("Failed to load model.") + + def initialize(self): + output_data = dict() + output_data["k_00"] = 0.4912490635619572 + output_data["k_11"] = 0.4912490635989945 + output_data["porosity"] = 0.4933482661391027 + + if self._sim_id == 0: + output_data["k_00"] = 0.4912490640081466 + output_data["k_11"] = 0.4912490640081367 + + return output_data + + def get_state(self): + return self._state; + + def set_state(self, state): + self._state = state + + def solve(self, macro_data, dt): + model_eval, _ = self._model.eval_metamodel(np.array([macro_data["concentration"]])[:, np.newaxis]) + output_data = dict() + output_data["k_00"] = model_eval["k_00"][0][0] + output_data["k_01"] = model_eval["k_01"][0][0] + output_data["k_10"] = model_eval["k_10"][0][0] + output_data["k_11"] = model_eval["k_11"][0][0] + output_data["porosity"] = model_eval["porosity"][0][0] + output_data["grain_size"] = math.sqrt((1 - model_eval["porosity"][0][0]) / math.pi) + + return output_data + + + diff --git a/two-scale-heat-conduction/micro-dumux-mada/params.input b/two-scale-heat-conduction/micro-dumux-mada/params.input new file mode 100644 index 000000000..964a6896d --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-mada/params.input @@ -0,0 +1,25 @@ +[Assembly] +Multithreading = false + +[TimeLoop] +TEnd = 0.25 # end time of the simulation +DtInitial = 0.01 # initial time step size +MaxTimeStepSize = 0.01 # maximal time step size + +[Grid] +LowerLeft = 0.0 0.0 # lower left (front) corner of the domain (keep this fixed at 0 0!) +UpperRight = 1.0 1.0 # upper right (back) corner of the domain +Cells = 80 80 # grid resolution in each coordinate direction +Periodic = 1 1 # Periodic Boundary conditions in both dimensions + +[Problem] +xi = 0.08 # phasefield parameter (lambda, set to around 4/Ncells) +omega = 0.01 # phasefield diffusivity/surface tension parameter (gamma) +kt = 1.0 # constant deciding speed of expansion/contraction +eqconc = 0.5 # equilibrium concentration +ks = 1.0 # conductivity of sand material +kg = 0.0 # conductivity of void material +Name = cell_phase # base name for VTK output files +Radius = 0.4 # initial radius of the grain +PhasefieldICScaling = 4.0 # factor in initial phasefield function +MaxPorosity = 0.9686 # porosity cap diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/clean.sh b/two-scale-heat-conduction/micro-dumux-surrogate/clean.sh new file mode 100755 index 000000000..3e9fd4e51 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_dumux . diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/create-snapshots.sh b/two-scale-heat-conduction/micro-dumux-surrogate/create-snapshots.sh new file mode 100755 index 000000000..8d616e8ae --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/create-snapshots.sh @@ -0,0 +1,29 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +usage() { echo "Usage: cmd [-s] [-p n]" 1>&2; exit 1; } + +# Check if no input argument was provided +if [ -z "$*" ] ; then + echo "No input argument provided. Micro Manager for snapshot computation is launched in serial" + micro-manager-precice --snapshot micro-manager-snapshot-config.json +fi + +while getopts ":sp" opt; do + case ${opt} in + s) + micro-manager-precice --snapshot micro-manager-snapshot-config.json + ;; + p) + mpiexec -n "$2" micro-manager-precice --snapshot micro-manager-snapshot-config.json + ;; + *) + usage + ;; + esac +done + +close_log diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-snapshot-config.json b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-snapshot-config.json new file mode 100644 index 000000000..3eabcff11 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-snapshot-config.json @@ -0,0 +1,15 @@ +{ + "micro_file_name": "micro_sim", + "coupling_params": { + "parameter_file_name": "input_samples.hdf5", + "read_data_names": ["concentration"], + "write_data_names": ["k_00", "k_01", "k_10", "k_11", "porosity"] + }, + "simulation_params": { + "micro_dt": 0.01 + }, + "snapshot_params": { + "initialize_once": false + }, + "output_directory": "output" +} diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-surrogate-config.json b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-surrogate-config.json new file mode 100644 index 000000000..198abbffd --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro-manager-surrogate-config.json @@ -0,0 +1,28 @@ +{ + "micro_file_name": "micro_sim_sur", + "coupling_params": { + "participant_name": "Micro-Manager", + "precice_config_file_name": "../precice-config.xml", + "macro_mesh_name": "macro-mesh", + "write_data_names": ["k_00", "k_01", "k_10", "k_11", "porosity"], + "read_data_names": ["concentration"] + }, + "simulation_params": { + "micro_dt": 0.01, + "macro_domain_bounds": [0.0, 1.0, 0.0, 0.5], + "decomposition": [1, 1], + "adaptivity": "True", + "adaptivity_settings": { + "type": "global", + "data": ["k_00", "k_11", "porosity", "concentration"], + "history_param": 0.1, + "coarsening_constant": 0.2, + "refining_constant": 0.05, + "every_implicit_iteration": "False", + "similarity_measure": "L2rel" + } + }, + "diagnostics": { + "data_from_micro_sims": ["grain_size"] + } +} diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim.pc.in b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim.pc.in new file mode 100644 index 000000000..566aba5e4 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim.pc.in @@ -0,0 +1,15 @@ +prefix=@prefix@ +exec_prefix=@exec_prefix@ +libdir=@libdir@ +includedir=@includedir@ +CXX=@CXX@ +CC=@CC@ +DEPENDENCIES=@REQUIRES@ + +Name: @PACKAGE_NAME@ +Version: @VERSION@ +Description: micro_sim module +URL: http://dune-project.org/ +Requires: dumux-phasefield dumux-precice +Libs: -L${libdir} +Cflags: -I${includedir} diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim_sur.py b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim_sur.py new file mode 100644 index 000000000..d2502c60e --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/micro_sim_sur.py @@ -0,0 +1,59 @@ +""" +Micro simulation Surrogate, requrie previous computation of surrogate model +""" +import os +import subprocess +from bayesvalidrox import PyLinkForwardModel, Input, PCE, ExpDesigns, Engine +import h5py +import joblib +import numpy as np +import math + + +class MicroSimulation: + + def __init__(self, sim_id): + """ + Constructor of MicroSimulation class. + """ + self._sim_id = sim_id + self._state = None + + self._model = None + with open('micro-dumux-surrogate.pkl', 'rb') as input: + self._model = joblib.load(input) + if self._model is None: + raise RuntimeError("Failed to load model.") + + def initialize(self): + output_data = dict() + output_data["k_00"] = 0.4912490635619572 + output_data["k_11"] = 0.4912490635989945 + output_data["porosity"] = 0.4933482661391027 + + if self._sim_id == 0: + output_data["k_00"] = 0.4912490640081466 + output_data["k_11"] = 0.4912490640081367 + + return output_data + + def get_state(self): + return self._state; + + def set_state(self, state): + self._state = state + + def solve(self, macro_data, dt): + model_eval, _ = self._model.eval_metamodel(np.array([macro_data["concentration"]])[:, np.newaxis]) + output_data = dict() + output_data["k_00"] = model_eval["k_00"][0][0] + output_data["k_01"] = model_eval["k_01"][0][0] + output_data["k_10"] = model_eval["k_10"][0][0] + output_data["k_11"] = model_eval["k_11"][0][0] + output_data["porosity"] = model_eval["porosity"][0][0] + output_data["grain_size"] = math.sqrt((1 - model_eval["porosity"][0][0]) / math.pi) + + return output_data + + + diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/model.py b/two-scale-heat-conduction/micro-dumux-surrogate/model.py new file mode 100644 index 000000000..5284f7485 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/model.py @@ -0,0 +1,5 @@ +# Dummy model function for the micro-dumux surrogate. +# We do not wrap the original DuMuX model because we will directly provide +# snapshots (computed by the Micro Manager) to BayesValidRox. +def model(samples): + return None diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/params.input b/two-scale-heat-conduction/micro-dumux-surrogate/params.input new file mode 100644 index 000000000..964a6896d --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/params.input @@ -0,0 +1,25 @@ +[Assembly] +Multithreading = false + +[TimeLoop] +TEnd = 0.25 # end time of the simulation +DtInitial = 0.01 # initial time step size +MaxTimeStepSize = 0.01 # maximal time step size + +[Grid] +LowerLeft = 0.0 0.0 # lower left (front) corner of the domain (keep this fixed at 0 0!) +UpperRight = 1.0 1.0 # upper right (back) corner of the domain +Cells = 80 80 # grid resolution in each coordinate direction +Periodic = 1 1 # Periodic Boundary conditions in both dimensions + +[Problem] +xi = 0.08 # phasefield parameter (lambda, set to around 4/Ncells) +omega = 0.01 # phasefield diffusivity/surface tension parameter (gamma) +kt = 1.0 # constant deciding speed of expansion/contraction +eqconc = 0.5 # equilibrium concentration +ks = 1.0 # conductivity of sand material +kg = 0.0 # conductivity of void material +Name = cell_phase # base name for VTK output files +Radius = 0.4 # initial radius of the grain +PhasefieldICScaling = 4.0 # factor in initial phasefield function +MaxPorosity = 0.9686 # porosity cap diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/requirements.txt b/two-scale-heat-conduction/micro-dumux-surrogate/requirements.txt new file mode 100644 index 000000000..f167b1a92 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/requirements.txt @@ -0,0 +1,3 @@ +numpy +bayesvalidrox +micro-manager-precice diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/run-surrogate-workflow.sh b/two-scale-heat-conduction/micro-dumux-surrogate/run-surrogate-workflow.sh new file mode 100755 index 000000000..ddf636773 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/run-surrogate-workflow.sh @@ -0,0 +1,9 @@ +#!/usr/bin/env bash +set -e -u + +python3 -m venv .venv +. .venv/bin/activate + +pip install -r requirements.txt + +python surrogate_workflow.py diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/run.sh b/two-scale-heat-conduction/micro-dumux-surrogate/run.sh new file mode 100755 index 000000000..f28acb260 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/run.sh @@ -0,0 +1,29 @@ +#!/usr/bin/env bash +set -e -u + +. ../../tools/log.sh +exec > >(tee --append "$LOGFILE") 2>&1 + +usage() { echo "Usage: cmd [-s] [-p n]" 1>&2; exit 1; } + +# Check if no input argument was provided +if [ -z "$*" ] ; then + echo "No input argument provided. Micro Manager is launched in serial" + micro-manager-precice micro-manager-surrogate-config.json +fi + +while getopts ":sp" opt; do + case ${opt} in + s) + micro-manager-precice micro-manager-surrogate-config.json + ;; + p) + mpiexec -n "$2" micro-manager-precice micro-manager-surrogate-config.json + ;; + *) + usage + ;; + esac +done + +close_log diff --git a/two-scale-heat-conduction/micro-dumux-surrogate/surrogate_workflow.py b/two-scale-heat-conduction/micro-dumux-surrogate/surrogate_workflow.py new file mode 100644 index 000000000..cb7cac941 --- /dev/null +++ b/two-scale-heat-conduction/micro-dumux-surrogate/surrogate_workflow.py @@ -0,0 +1,197 @@ +import os +import subprocess +from bayesvalidrox import PyLinkForwardModel, Input, PCE, ExpDesigns, Engine +import h5py +import joblib +import numpy as np +import matplotlib.pyplot as plt + + +def create_snapshots() -> None: + """ + Create snapshots of the DuMuX model in micro-dumux/ using the Micro Manager. + The snapshots are saved in the specified directory. + """ + # What effect do uniformly distributed concentration samples have on the quality of the surrogate model? + concentration_samples = np.linspace(0.0, 0.9, 500) + + with h5py.File("input_samples.hdf5", "w") as f: + f.create_dataset("concentration", data=concentration_samples) + + # Run the Micro Manager to create snapshot database + subprocess.call('micro-manager-precice --snapshot micro-manager-snapshot-config.json', shell=True) + + +def read_snapshots(snapshots_dir: str) -> tuple: + """ + Read the snapshots created by the Micro Manager and return the inputs and outputs. + The inputs are the concentration samples, and the outputs are the porosity and conductivity matrix values. + + Parameters + ---------- + snapshots_dir : str + The directory where the snapshots are stored. + + Returns + ------- + tuple + A tuple containing the inputs (concentration samples) and outputs (porosity and conductivity values). + """ + with h5py.File(os.path.join(snapshots_dir, "snapshot_data.hdf5"), "r") as f: + concentration_data = f["concentration"][:] + porosity_data = f["porosity"][:] + k_00_data = f["k_00"][:] + k_01_data = f["k_01"][:] + k_10_data = f["k_10"][:] + k_11_data = f["k_11"][:] + + concentration = np.swapaxes(np.array([concentration_data]), 0, 1) + porosity = np.swapaxes(np.array([porosity_data]), 0, 1) + k_00 = np.swapaxes(np.array([k_00_data]), 0, 1) + k_01 = np.swapaxes(np.array([k_01_data]), 0, 1) + k_10 = np.swapaxes(np.array([k_10_data]), 0, 1) + k_11 = np.swapaxes(np.array([k_11_data]), 0, 1) + + outputs = {"porosity": porosity, "k_00": k_00, "k_01": k_01, "k_10": k_10, "k_11": k_11, "x_values": np.array([0])} + return concentration, outputs + + +def split_samples(X, y, n_valid): + """ + Split the samples and evaluations into training and validation/test data. + The split is performed randomly. + + Parameters + ---------- + X : np.ndarray + Samples, shape (#samples, #parameters) + y : dict + Corresponding model evaluations. Expected to match BVR output format. + n_valid : int + Number of samples to keep for validation. + + Returns + ------- + X_train, y_train, + X_valid, y_valid + """ + n_samples = X.shape[0] + if n_valid >= n_samples: + raise AttributeError('The set number of validation points is invalid.') + + # Random split + n_train = n_samples - n_valid + choice = np.random.choice( + range(n_samples), size=(n_train,), replace=False + ) + ind = np.zeros(n_samples, dtype=bool) + ind[choice] = True + + # Split samples + X_train = X[ind] + X_valid = X[~ind] + + # Split outputs + y_train = {} + y_valid = {} + for key in y: + if key != "x_values": + y_train[key] = y[key][ind] + y_valid[key] = y[key][~ind] + + return X_train, y_train, X_valid, y_valid + + +def create_surrogate(snapshots_dir: str) -> tuple: + """ + Create a surrogate model from the Micro Manager snapshots. + + """ + # We create a fake model from model.py because we directly provide the + # input and outputs from the Micro Manager snapshots. + model = PyLinkForwardModel() + model.py_file = "model" + model.name = "micro-dumux-surrogate" + model.link_type = "function" + model.output_names = ["porosity", "k_00", "k_01", "k_10", "k_11"] + + x, y = read_snapshots(snapshots_dir) + + # Split the samples into training and validation sets + n_valid = 200 + x_train, y_train, x_valid, y_valid = split_samples(x, y, n_valid) + + inputs = Input() + inputs.add_marginals(name="concentration", dist_type="unif", parameters=[0, 0.5]) + + exp_design = ExpDesigns(inputs) + exp_design.x = x_train + exp_design.y = y_train + + # Create the surrogate model + meta_model = PCE(inputs) + meta_model.meta_model_type = "aPCE" + meta_model.pce_reg_method = "FastARD" + meta_model.pce_deg = 5 + + # Train the surrogate model + engine = Engine(meta_model, model, exp_design) + engine.train_normal() + + with open(f'{model.name}.pkl', 'wb') as output: + joblib.dump(engine, output, 2) + + return x_valid, y_valid + + +def validate_surrogate(x_valid, y_valid, model_name="micro-dumux-surrogate.pkl"): + """ + Validate the surrogate model using the validation samples and outputs. + + Parameters + ---------- + x_valid : np.ndarray + Validation samples. + y_valid : dict + Corresponding model evaluations for validation samples. + model_name : str + Name of the surrogate model file. + """ + with open(model_name, 'rb') as input: + engine = joblib.load(input) + + y_metamod, _ = engine.eval_metamodel(x_valid) + + # engine.plot_adapt(y_valid, y_metamod, y_metamod_std, x_valid) + + # Compare predictions with true values + plt.figure() + plt.scatter(y_valid["porosity"], y_metamod["porosity"]) + plt.xlabel("True Values") + plt.ylabel("Predictions") + plt.title(f"Validation: porosity") + plt.plot([0.4, 1.1], [0.4, 1.1], "k--") + plt.xlim(0.4, 1.1) + plt.ylim(0.4, 1.1) + plt.show() + + plt.figure() + plt.scatter(x_valid, y_valid["porosity"]) + plt.xlabel("valid concentration") + plt.ylabel("valid porosity") + plt.title("IO: con to poro") + plt.show() + +def main(): + snapshots_dir = "output" + if not os.path.exists(snapshots_dir): + os.makedirs(snapshots_dir) + + # os.chdir(snapshots_dir) + create_snapshots() + x_valid, y_valid = create_surrogate(snapshots_dir) + validate_surrogate(x_valid, y_valid) + + +if __name__ == "__main__": + main()