New tutorial: Partitioned Burgers eq. 1D

partitioned-burgers-1d/.solver-nutils-dgalerkin/compare_burgers.py

@@ -0,0 +1,54 @@
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+DATA_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__)))
+
+data_dir_dgalerkin = os.path.join(DATA_DIR, "solver-nutils-dgalerkin", "data-training")
+data_dir_diffuse = os.path.join(DATA_DIR, "solver-scipy-fvolumes", "data-training")
+
+# file names are the same, folders are different
+files_ = os.listdir(data_dir_dgalerkin)
+files_ = sorted(f for f in files_ if f.endswith(".npz"))
+
+
+file_num = 0
+timestep = 0
+
+file_dgalerkin = os.path.join(data_dir_dgalerkin, files_[file_num])
+file_diffuse = os.path.join(data_dir_diffuse, files_[file_num])
+
+data_dgalerkin = np.load(file_dgalerkin)["Solver-Mesh-1D-Internal"]
+data_diffuse = np.load(file_diffuse)["Solver-Mesh-1D-Internal"]
+
+print(f"DG data shape: {data_dgalerkin.shape}")
+print(f"FV data shape: {data_diffuse.shape}")
+
+plt.figure(figsize=(12, 5))
+plt.plot(data_dgalerkin[timestep, :], label='DG Solver', color='blue')
+plt.plot(data_diffuse[timestep, :], label='FV Solver', color='orange', linestyle='--')
+plt.title(f'Comparison at Timestep {timestep}')
+plt.xlabel('Spatial Position')
+plt.ylabel('Solution Value')
+plt.legend()
+plt.savefig(os.path.join(DATA_DIR, f'comparison_timestep_{timestep}.png'))
+
+# plot the imshow with unified colormap
+vmin = min(data_dgalerkin.min(), data_diffuse.min())
+vmax = max(data_dgalerkin.max(), data_diffuse.max())
+
+plt.figure(figsize=(12, 5))
+plt.subplot(1, 2, 1)
+plt.imshow(data_dgalerkin.T, aspect='auto', cmap='viridis', vmin=vmin, vmax=vmax)
+plt.title('DG Solver Evolution')
+plt.ylabel('Spatial Position')
+plt.xlabel('Timestep')
+plt.colorbar(label='Solution Value')
+plt.subplot(1, 2, 2)
+plt.imshow(data_diffuse.T, aspect='auto', cmap='viridis', vmin=vmin, vmax=vmax)
+plt.title('FV Solver Evolution')
+plt.ylabel('Spatial Position')
+plt.xlabel('Timestep')
+plt.colorbar(label='Solution Value')
+plt.tight_layout()
+plt.show()

partitioned-burgers-1d/.solver-nutils-dgalerkin/solver.py

@@ -0,0 +1,148 @@
+# Original source: https://github.com/evalf/nutils/blob/d73749ff7d64c9ccafdbb88cd442f80b9448c118/examples/burgers.py
+
+from nutils import mesh, function, export, testing
+from nutils.solver import System
+from nutils.expression_v2 import Namespace
+import treelog as log
+import numpy as np
+import itertools
+import precice
+import json
+import os
+import argparse
+
+def _generate_initial_condition(x_coords, ic_config, epoch):
+    np.random.seed(epoch)
+    ic_values = np.zeros(len(x_coords))
+    if ic_config["type"] == "sinusoidal":
+        num_modes = ic_config.get("num_modes", 1)
+        superpositions = np.random.randint(2, num_modes + 1)
+        for _ in range(superpositions):
+            amp = np.random.uniform(0.1, 2)
+            k = np.random.randint(ic_config["wavenumber_range"][0], ic_config["wavenumber_range"][1] + 1)
+            phase_shift = np.random.uniform(0, 2 * np.pi)
+            ic_values += amp * np.sin(2 * np.pi * k * x_coords + phase_shift)
+    return ic_values
+
+def project_initial_condition(domain_min, domain_max, nelems, ic_config, epoch):
+    # 1. Generate a high-resolution "truth" on a fine grid
+    fine_res = nelems * 10
+    fine_x = np.linspace(domain_min[0], domain_max[0], fine_res, endpoint=False)
+    fine_u = _generate_initial_condition(fine_x, ic_config, epoch)
+
+    # 2. Average the high-resolution truth over each coarse cell
+    u_projected = np.zeros(nelems)
+    for i in range(nelems):
+        cell_start = i * 10
+        cell_end = (i + 1) * 10
+        u_projected[i] = np.mean(fine_u[cell_start:cell_end])
+
+    return u_projected
+
+def main(dim: int,
+         epoch: int = 0,
+         btype: str = 'discont',
+         degree: int = 1,
+         newtontol: float = 1e-5,
+         config_file: str = "precice-config.xml"):
+
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    with open(os.path.join(script_dir, "..", "python_participant", "config.json"), 'r') as f:
+        config = json.load(f)["solver"]
+    with open(os.path.join(script_dir, "ic_params.json"), 'r') as f:
+        ic_config = json.load(f)["initial_conditions"]
+
+    config_path = os.path.join(script_dir, "..", f"{dim}d", config_file)
+    participant = precice.Participant("Solver", config_path, 0, 1)
+
+    mesh_internal_name = f"Solver-Mesh-{dim}D-Internal"
+    mesh_boundaries_name = f"Solver-Mesh-{dim}D-Boundaries"
+    data_name = f"Data_{dim}D"
+    res = config[f"{dim}d_resolution"]
+    domain_min = config[f"{dim}d_domain_min"]
+    domain_max = config[f"{dim}d_domain_max"]
+    nelems = res[0]
+
+    domain, geom = mesh.line(np.linspace(domain_min[0], domain_max[0], nelems + 1), periodic=True)
+
+    # Define all nelems +1 nodes for evaluation
+    eval_coords_x = np.linspace(domain_min[0], domain_max[0], nelems + 1)
+
+    # Define the nelems vertices for saving (all but the last)
+    trunc_coords_x = eval_coords_x[:-1]
+    internal_coords = np.array([trunc_coords_x, np.full(len(trunc_coords_x), domain_min[1])]).T
+    boundary_coords = np.array([[domain_min[0], domain_min[1]], [domain_max[0], domain_max[1]]])
+
+    internal_vertex_ids = participant.set_mesh_vertices(mesh_internal_name, internal_coords)
+    boundary_vertex_ids = participant.set_mesh_vertices(mesh_boundaries_name, boundary_coords)
+
+    sample = domain.locate(geom, eval_coords_x, tol=1e-5)
+
+    ns = Namespace()
+    ns.x = geom
+    ns.define_for('x', gradient='∇', normal='n', jacobians=('dV', 'dS'))
+    ns.u = domain.field('u', btype=btype, degree=degree)
+    ns.du = ns.u - function.replace_arguments(ns.u, 'u:u0')
+    ns.v = domain.field('v', btype=btype, degree=degree)
+    ns.t = function.field('t')
+    ns.dt = ns.t - function.field('t0')
+    ns.f = '.5 u^2'
+    ns.C = 1
+
+    res_pde = domain.integral('(v du / dt - ∇(v) f) dV' @ ns, degree=degree*2)
+    res_pde -= domain.interfaces.integral('[v] n ({f} - .5 C [u] n) dS' @ ns, degree=degree*2)
+    system = System(res_pde, trial='u', test='v')
+
+    # Project the initial condition
+    u_averaged = project_initial_condition(domain_min, domain_max, nelems, ic_config, epoch)
+    ns.uic = domain.basis('discont', degree=0).dot(u_averaged)
+
+    sqr = domain.integral('(u - uic)^2 dV' @ ns, degree=max(degree*2, 5))
+    args = System(sqr, trial='u').solve()
+
+    if participant.requires_initial_data():
+        # Evaluate at all nodes
+        all_data = sample.eval(ns.u, arguments=args)
+        # Truncate last element
+        trunc_data = all_data[:-1]
+        boundary_data_values = np.array([trunc_data[0], trunc_data[-1]])
+
+        participant.write_data(mesh_internal_name, data_name, internal_vertex_ids, trunc_data)
+        participant.write_data(mesh_boundaries_name, data_name, boundary_vertex_ids, boundary_data_values)
+
+    participant.initialize()
+
+    args['t'] = 0.
+
+    with log.iter.plain('timestep', itertools.count()) as steps:
+        for _ in steps:
+            if not participant.is_coupling_ongoing():
+                break
+
+            timestep = participant.get_max_time_step_size()
+
+            args = system.step(timestep=timestep, arguments=args, timearg='t', suffix='0', tol=newtontol)
+
+            all_data = sample.eval(ns.u, arguments=args)
+            trunc_data = all_data[:-1]
+            boundary_data_values = np.array([trunc_data[0], trunc_data[-1]])
+
+            participant.write_data(mesh_internal_name, data_name, internal_vertex_ids, trunc_data)
+            participant.write_data(mesh_boundaries_name, data_name, boundary_vertex_ids, boundary_data_values)
+
+            participant.advance(timestep)
+
+    participant.finalize()
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument("dim", type=int, choices=[1, 2], help="Dimension of the simulation")
+    parser.add_argument('--config_file', type=str, default="precice-config.xml")
+    parser.add_argument('--btype', type=str, default='discont')
+    parser.add_argument('--degree', type=int, default=1)
+    parser.add_argument('--newtontol', type=float, default=1e-5)
+    parser.add_argument("--epoch", type=int, default=0, help="Current epoch number")
+
+    args_cli = parser.parse_args()
+
+    main(dim=args_cli.dim, epoch=args_cli.epoch, btype=args_cli.btype, degree=args_cli.degree, newtontol=args_cli.newtontol, config_file=args_cli.config_file)

partitioned-burgers-1d/README.md

@@ -0,0 +1,142 @@
+---
+title: Partitioned 1D Burgers' Equation
+permalink: tutorials-partitioned-burgers-1d.html
+keywords: Python, Neural Network, Surrogate, Burgers Equation, Finite Volume, CFD
+summary: This tutorial demonstrates the partitioned solution of the 1D Burgers' equation using preCICE and a neural network surrogate solver.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/partitioned-burgers-1d). Read how in the [tutorials introduction](https://precice.org/tutorials.html).
+{% endnote %}
+
+## Setup
+
+We solve the 1D viscous Burgers' equation on the domain $[0,2]$:
+
+
+$$
+\frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2} - u \frac{\partial u}{\partial x},
+$$
+
+where $u(x,t)$ is the scalar velocity field and $\nu$ is the viscosity. In this tutorial by default $\nu$ is very small ($10^{-12}$), but can be changed in the solver.
+
+
+The domain is partitioned into participants at $x=1$:
+
+- **Dirichlet**: Solves the left half $[0,1]$ and receives Dirichlet boundary conditions at the interface.
+- **Neumann**: Solves the right half $[1,2]$ and receives Neumann boundary conditions at the interface.
+
+Both outer boundaries use zero-gradient conditions $\frac{\partial u}{\partial x} = 0$. The problem is solved for different initial conditions of superimposed sine waves, which can be generated using the provided script.
+
+<p align="center">
+  <img src="images/tutorials-partitioned-burgers-1d-initial-condition.png" alt="Initial Condition" width="400"/>
+  <br><em>Initial condition used for the simulation (seed 0, see <code>generate_ic.py</code>)</em>
+</p>
+
+The conservative formulation of the Burgers' equation `solver-scipy` is implemented in a first-order finite volume code using Lax-Friedrichs fluxes and implicit Euler time stepping.
+
+This tutorial includes two versions for the Neumann participant:
+- A standard finite volume solver (`neumann-scipy`).
+- A pre-trained neural network surrogate that approximates the solver (`neumann-surrogate`).
+
+
+## Configuration
+
+preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)):
+
+<p align="center">
+  <img src="images/tutorials-partitioned-burgers-1d-precice-config.png" alt="preCICE configuration visualization" width="600"/>
+</p>
+
+
+## Running the tutorial
+
+### Initial condition 
+
+Before running the participants, you must generate the initial condition file. From the root of the tutorial folder, run the command and replace `<seed>` with any integer.:
+
+```bash
+python3 generate_ic.py --epoch <seed>
+```
+
+This script requires the Python libraries `numpy` and `matplotlib`. It generates the initial condition file `initial_condition.npz` used by both participants. If you do not have the dependencies, then you can run either of the participant run scripts to create a Python virtual environment with the required packages and try again. 
+
+Or, simply, run the helper scripts (see below), which will also generate the initial condition file if it does not exist.
+
+---
+
+#### Helper scripts
+
+There are helper scripts that automate runs and visualization of both participants. They also accept an integer seed argument to specify the initial condition.
+
+```bash
+run-partitioned-scipy.sh
+```
+
+and
+
+```bash
+run-partitioned-surrogate.sh
+```
+
+### Running the participants
+
+To run the partitioned simulation, open two separate terminals and start each participant individually:
+
+You can find the corresponding `run.sh` script for running the case in the folders corresponding to the participant you want to use:
+
+```bash
+cd dirichlet-scipy
+./run.sh
+```
+
+and
+
+```bash
+cd neumann-scipy
+./run.sh
+```
+
+or, to use the pretrained neural network surrogate participant:
+
+```bash
+cd neumann-surrogate
+./run.sh
+```
+
+**Note:** The surrogate participant requires PyTorch and related dependencies, which requires several gigabytes of disk space.
+
+### Monolithic solution (reference)
+
+You can run the whole domain using the monolithic solver for comparison:
+
+```bash
+cd solver-scipy
+./run.sh
+```
+
+
+
+## Visualization
+
+After both participants (and/or monolithic simulation) have finished, you can run the visualization script.
+`visualize_partitioned_domain.py` generates plots comparing the partitioned and monolithic solutions. You can specify which timestep to plot:
+
+```bash
+python3 visualize_partitioned_domain.py --neumann neumann-surrogate/surrogate.npz [timestep]
+```
+
+The script will produce the following output files in the `images/` directory:
+- `full-domain-timestep-slice.png`: Solution $u$ at a selected timestep
+
+<p align="left">
+  <img src="images/tutorials-partitioned-burgers-1d-full-domain-timestep-slice.png" alt="Full Domain Timestep Slice" width="400"/>
+</p>
+
+- `gradient-timestep-slice.png`: Gradient $du/dx$ at a selected timestep
+
+- `full-domain-evolution.png`: Time evolution of the solution
+
+<p align="left">
+  <img src="images/tutorials-partitioned-burgers-1d-full-domain-evolution.png" alt="Full Domain Evolution" width="400"/>
+</p>

partitioned-burgers-1d/clean-tutorial.sh

@@ -0,0 +1 @@
+../tools/clean-tutorial-base.sh

partitioned-burgers-1d/clean.sh

@@ -0,0 +1,6 @@
+#!/usr/bin/env sh
+set -e -u
+
+rm -f solver-scipy/full_domain.npz
+rm -f images/full_domain_evolution.png images/full_domain_timestep_slice.png images/gradient_timestep_slice.png images/initial_condition.png
+rm -f initial_condition.npz # comment out?

partitioned-burgers-1d/dirichlet-scipy/clean.sh

@@ -0,0 +1,6 @@
+#!/usr/bin/env sh
+set -e -u
+
+rm -rf precice-profiling
+rm -f dirichlet-scipy.log precice-Dirichlet-convergence.log precice-Dirichlet-iterations.log
+rm -f dirichlet.npz

partitioned-burgers-1d/dirichlet-scipy/run.sh

@@ -0,0 +1,30 @@
+#!/usr/bin/env bash
+set -e -u
+
+solver_path="../solver-scipy"
+
+if [ -d "$solver_path/.venv" ]; then
+	echo "Using existing virtual environment"
+	. "$solver_path/.venv/bin/activate"
+else
+	echo "Creating new virtual environment"
+	python3 -m venv "$solver_path/.venv"
+	. "$solver_path/.venv/bin/activate"
+	pip install -r $solver_path/requirements.txt && pip freeze > $solver_path/pip-installed-packages.log
+fi
+
+if [ -f "../initial_condition.npz" ]; then
+	:
+else
+	echo "Error, missing initial condition file."
+	echo "Run 'python3 ../generate_ic.py --epoch <seed>' to create one."
+	exit 1
+fi
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+echo "[preCICE] Waiting for Neumann participant..."
+python3 "$solver_path/solver.py" Dirichlet
+
+close_log