New tutorial: Partitioned Burgers eq. 1D

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/clean-tutorial.sh

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

partitioned-burgers-1d/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/dirichlet-scipy/run.sh

@@ -0,0 +1,4 @@
+#!/usr/bin/env bash
+set -e -u
+
+python3 ../solver-scipy-fvolumes/solver.py Dirichlet

partitioned-burgers-1d/generate_ic.py

@@ -0,0 +1,67 @@
+import numpy as np
+import json
+import os
+import argparse
+import matplotlib.pyplot as plt
+
+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, domain_max, 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
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Generate initial conditions for the Burgers equation simulation.")
+    parser.add_argument("--epoch", type=int, default=0, help="Seed for the random number generator to ensure reproducibility.")
+    args = parser.parse_args()
+
+    script_dir = os.path.dirname(os.path.abspath(__file__))
+    with open(os.path.join(script_dir, "ic_params.json"), 'r') as f:
+        config = json.load(f)
+
+    ic_config = config["initial_conditions"]
+    domain_config = config["domain"]
+
+    # full domain
+    full_domain_min = domain_config["full_domain_min"]
+    full_domain_max = domain_config["full_domain_max"]
+    nelems_total = domain_config["nelems_total"]
+
+    # Generate IC
+    initial_condition = project_initial_condition(full_domain_min, full_domain_max, nelems_total, ic_config, args.epoch)
+
+    output_path = os.path.join(script_dir, "initial_condition.npz")
+    np.savez(output_path, initial_condition=initial_condition)
+
+    plt.figure(figsize=(8, 4))
+    x_coords = np.linspace(full_domain_min, full_domain_max, nelems_total, endpoint=False)
+    plt.figure(figsize=(10, 5))
+    plt.plot(x_coords, initial_condition, marker='.', linestyle='-')
+    plt.xlabel('Spatial Coordinate (x)')
+    plt.ylabel('Solution Value (u)')
+    plt.grid(True)
+    plt.savefig(os.path.join(script_dir, "initial_condition.png"))
+    print(f"Initial condition and plot saved to {output_path}")
+    plt.close()

partitioned-burgers-1d/ic_params.json

@@ -0,0 +1,12 @@
+{
+  "initial_conditions": {
+    "type": "sinusoidal",
+    "wavenumber_range": [1, 7],
+    "num_modes": 4
+  },
+  "domain": {
+    "nelems_total": 256,
+    "full_domain_min": 0.0,
+    "full_domain_max": 2.0
+  }
+}

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

@@ -0,0 +1,4 @@
+#!/usr/bin/env bash
+set -e -u
+
+python3 ../solver-scipy-fvolumes/solver.py Neumann

partitioned-burgers-1d/precice-config.xml

@@ -0,0 +1,56 @@
+<?xml version="1.0"?>
+<precice-configuration>
+  <log>
+    <sink filter="%Severity% > info" format="---[ %Severity% ] %Message%" enabled="true" />
+  </log>
+
+  <data:scalar name="Gradient" />
+  <data:scalar name="Velocity" />
+
+  <mesh name="Dirichlet-Mesh" dimensions="2">
+    <use-data name="Velocity" />
+    <use-data name="Gradient" />
+  </mesh>
+
+  <mesh name="Neumann-Mesh" dimensions="2">
+    <use-data name="Velocity" />
+    <use-data name="Gradient" />
+  </mesh>
+
+  <participant name="Dirichlet">
+    <provide-mesh name="Dirichlet-Mesh" />
+    <receive-mesh name="Neumann-Mesh" from="Neumann" />
+    <write-data name="Gradient" mesh="Dirichlet-Mesh" />
+    <read-data name="Velocity" mesh="Dirichlet-Mesh" />
+    <mapping:nearest-neighbor direction="read" from="Neumann-Mesh" to="Dirichlet-Mesh" constraint="consistent" />
+  </participant>
+
+  <participant name="Neumann">
+    <provide-mesh name="Neumann-Mesh" />
+    <receive-mesh name="Dirichlet-Mesh" from="Dirichlet" />
+    <write-data name="Velocity" mesh="Neumann-Mesh" />
+    <read-data name="Gradient" mesh="Neumann-Mesh" />
+    <mapping:nearest-neighbor direction="read" from="Dirichlet-Mesh" to="Neumann-Mesh" constraint="consistent" />
+  </participant>
+
+  <m2n:sockets acceptor="Dirichlet" connector="Neumann" exchange-directory=".."/>
+
+  <coupling-scheme:serial-implicit>
+    <participants first="Dirichlet" second="Neumann" />
+    <max-time value="0.1" />
+    <time-window-size value="0.001" />
+    <max-iterations value="10" />
+    <exchange data="Gradient" mesh="Dirichlet-Mesh" from="Dirichlet" to="Neumann" initialize="true" substeps="false" />
+    <exchange data="Velocity" mesh="Neumann-Mesh" from="Neumann" to="Dirichlet" initialize="true" substeps="false" />
+    <relative-convergence-measure data="Gradient" mesh="Dirichlet-Mesh" limit="1e-3" />
+    <relative-convergence-measure data="Velocity" mesh="Neumann-Mesh" limit="1e-3" />
+    <acceleration:aitken>
+      <data name="Velocity" mesh="Neumann-Mesh" />
+      <initial-relaxation value="0.5" />
+    </acceleration:aitken>
+    <!-- <acceleration:constant>
+    <relaxation value="0.01" />
+    </acceleration:constant> -->
+  </coupling-scheme:serial-implicit>
+
+</precice-configuration>

partitioned-burgers-1d/requirements.txt

@@ -0,0 +1,8 @@
+numpy
+scipy
+matplotlib
+pyprecice
+
+# Optional: PyTorch for neural surrogate solver
+# Uncomment if using surrogate-burgers/solver.py
+# torch

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

@@ -0,0 +1,15 @@
+#!/usr/bin/env bash
+set -e -u
+
+rm precice-run -rf
+
+python3 generate_ic.py --epoch ${1:-0}
+
+cd dirichlet-scipy; pwd; ./run.sh &
+cd ../neumann-scipy; pwd; ./run.sh && cd ..
+
+# full domain reference solution
+cd solver-scipy-fvolumes; python3 solver.py None; cd ..
+
+
+python3 visualize_partitioned_domain.py