13_euler-with-amr.py

Go to the documentation of this file.

# This file is part of the Peano project. For conditions of distribution and
# use, please see the copyright notice at www.peano-framework.org
"""!
 
Create the PNG illustrations of the tutorial 13_euler-with-amr.ipynb
 
The script builds the tutorial's point explosion on a regular grid, runs
it and writes all pictures into the directory of this script. It then
builds the setup with the point-wise refinement criterion, runs it on the
resulting adaptive mesh, and draws the mesh, the solution, a comparison
with the regular mesh, plus the setup's task graph. Finally, it runs the
dynamically adaptive setup, which also erases cells while it runs time
steps, and draws its mesh adaptation task graph, the number of patches
over time and the mesh before and after it changes. Finally, it runs the
same dynamically adaptive setup with the gradient-based criterion and
compares the density profiles, the meshes and the numbers of patches of
the two criteria. The PDE terms are injected from Python, i.e. we do not
need any hand-written C++ file. Call it from anywhere:
 
     python3 tutorials/exahype2/illustrations/13_euler-with-amr.py
 
Pass --gradient-only to create only the pictures of the gradient section.
 
All builds and outputs go into a temporary directory (or --work-directory).
 
"""
import argparse
import os
import re
import types
import subprocess
import sys
import tempfile
 
Here = os.path.dirname(os.path.abspath(__file__))
PeanoRoot = os.path.abspath(os.path.join(Here, "..", "..", ".."))
sys.path.insert(0, os.path.join(PeanoRoot, "src"))
 
import matplotlib
 
matplotlib.use("Agg")
import matplotlib.patches
import matplotlib.pyplot as plt
import numpy as np
 
import peano4.api
import exahype2
import tarch.plotter.api.notebook
 
 
Flux = """
  const double irho  = 1.0 / Q[0];
  constexpr double gamma = 1.4;
  const double p = (gamma - 1.0) * (Q[3] - 0.5 * irho * (Q[1] * Q[1] + Q[2] * Q[2]));
  switch (normal) {
    case 0:
      F[0] = Q[1];
      F[1] = irho * Q[1] * Q[1] + p;
      F[2] = irho * Q[1] * Q[2];
      F[3] = irho * Q[1] * (Q[3] + p);
      break;
    case 1:
      F[0] = Q[2];
      F[1] = irho * Q[2] * Q[1];
      F[2] = irho * Q[2] * Q[2] + p;
      F[3] = irho * Q[2] * (Q[3] + p);
      break;
  }
"""
 
MaxEigenvalue = """
  const double irho  = 1.0 / Q[0];
  constexpr double gamma = 1.4;
  const double p = (gamma - 1.0) * (Q[3] - 0.5 * irho * (Q[1] * Q[1] + Q[2] * Q[2]));
  const double c = std::sqrt(gamma * p * irho);
  const double u = Q[1 + normal] * irho;
  lambdaMax[0] = std::max(std::abs(u + c), std::abs(u - c));
"""
 
InitPointWisely = """
  constexpr double gamma  = 1.4;
  constexpr double radius = 0.1;
  const bool   inside = x[0] * x[0] + x[1] * x[1] < radius * radius;
  const double rho    = inside ? 1.0 : 0.125;
  const double p      = inside ? 1.0 : 0.1;
  Q[0] = rho;
  Q[1] = 0.0;
  Q[2] = 0.0;
  Q[3] = p / (gamma - 1.0);
"""
 
VolumeSize = 0.005
AdaptiveVolumeSize = 0.02
AdaptiveMinVolumeSize = 0.002
RefineErasePointwisely = """
  if (Q[0] > 0.8) *refinementCommand = ::exahype2::RefinementCommand::Refine;
"""
DynamicRefineErasePointwisely = """
  if (Q[0] > 0.8) *refinementCommand = ::exahype2::RefinementCommand::Refine;
  else if (Q[0] < 0.1) *refinementCommand = ::exahype2::RefinementCommand::Erase;
"""
GradientRefinement = """
  const double densityGradient = std::abs(gradient[0]);
  if (densityGradient > 5.0)      *refinementCommand = ::exahype2::RefinementCommand::Refine;
  else if (densityGradient < 1.0) *refinementCommand = ::exahype2::RefinementCommand::Erase;
"""
EndTime = 0.2
DynamicEndTime = 0.3
TimeInBetweenPlots = 0.01
PatchFile = "solution-EulerFV.peano-patch-file"
 
FD4InitPointWisely = """
  constexpr double gamma = 1.4;
  constexpr double width = 0.05;
  const double r2  = x[0] * x[0] + x[1] * x[1];
  const double rho = 1.0 + 0.5 * std::exp(-r2 / (width * width));
  const double p   = std::pow(rho, gamma);
  Q[0] = rho;
  Q[1] = 0.0;
  Q[2] = 0.0;
  Q[3] = p / (gamma - 1.0);
"""
FD4PatchFile = "solution-EulerFD4.peano-patch-file"
 
 
def create_solver(adaptive, dynamic=False, gradient=False):
    """!
 
    The tutorial's solver over a regular, an adaptive or a dynamically
    adaptive mesh
 
    With gradient, the solver refines and erases with the gradient-based
    criterion instead of the point-wise one. The flag implies a dynamically
    adaptive mesh.
 
    """
    solver = exahype2.api.Solver("EulerFV", 4, 0)
    solver.add_user_solver_includes("#include <cmath>\n#include <algorithm>\n")
    solver.add_pde_term(exahype2.api.Flux(implementation=Flux, stateless=False))
    solver.add_pde_term(
        exahype2.api.MaxEigenvalue(implementation=MaxEigenvalue, stateless=False)
    )
    solver.add_pde_term(
        exahype2.api.InitPointWisely(implementation=InitPointWisely, stateless=False)
    )
    if adaptive and gradient:
        solver.add_pde_term(
            exahype2.api.GradientRefinement(implementation=GradientRefinement)
        )
    elif adaptive:
        solver.add_pde_term(
            exahype2.api.RefineErasePointwisely(
                implementation=DynamicRefineErasePointwisely
                if dynamic
                else RefineErasePointwisely
            )
        )
    if adaptive:
        solver.numerical_scheme = exahype2.api.fv.RusanovWithExplicitEuler(
            solver,
            9,
            min_volume_size=AdaptiveMinVolumeSize,
            max_volume_size=AdaptiveVolumeSize,
        )
    else:
        solver.numerical_scheme = exahype2.api.fv.RusanovWithExplicitEuler(solver, 9)
    if adaptive:
        solver.grid_construction = (
            exahype2.api.fv.SetupSimulationWithUserCodedRefinementCriterion(
                solver,
                volume_size=AdaptiveVolumeSize,
                min_volume_size=AdaptiveMinVolumeSize,
            )
        )
    else:
        solver.grid_construction = exahype2.api.fv.SetupSimulationOverRegularGrid(
            solver, volume_size=VolumeSize
        )
    solver.plotter = exahype2.api.fv.PlotAsVTKFile(
        solver, first_plot_time_stamp=0.0, time_in_between_plots=TimeInBetweenPlots
    )
    return solver
 
 
def build_and_run(directory, adaptive=False, dynamic=False, gradient=False):
    """! Generate, compile and run the tutorial setup within directory"""
    end_time = DynamicEndTime if dynamic or gradient else EndTime
    project = exahype2.api.Project(
        namespace=["tutorials", "exahype2", "euler"],
        project_name="Euler",
        directory=directory,
        executable="euler",
    )
    solver = create_solver(adaptive, dynamic, gradient)
    project.add_solver(solver)
    project.set_global_simulation_parameters(
        dimensions=2,
        offset=[-0.5, -0.5],
        size=[1.0, 1.0],
        min_end_time=end_time,
        max_end_time=end_time,
        periodic_BC=[True, True, True],
    )
    project.set_Peano4_installation(PeanoRoot, peano4.api.output.CompileMode.RELEASE)
    peano4_project = project.generate_Peano4_project()
    peano4_project.generate()
    subprocess.run(["make", "-j"], cwd=directory, check=True, stdout=subprocess.DEVNULL)
    subprocess.run(["./euler"], cwd=directory, check=True, stdout=subprocess.DEVNULL)
 
 
def build_and_run_fd4(directory, init_point_wisely):
    """!
 
    Generate, compile and run the fourth-order Finite Differences over the
    regular mesh within directory
 
    The run writes its output into directory/run.log, so we can parse the
    time step sizes afterwards. A run of the point explosion breaks down
    early, so we do not check its exit code.
 
    """
    project = exahype2.api.Project(
        namespace=["tutorials", "exahype2", "euler"],
        project_name="Euler",
        directory=directory,
        executable="euler",
    )
    solver = exahype2.api.Solver("EulerFD4", 4, 0)
    solver.add_user_solver_includes("#include <cmath>\n#include <algorithm>\n")
    solver.add_pde_term(exahype2.api.Flux(implementation=Flux, stateless=False))
    solver.add_pde_term(
        exahype2.api.MaxEigenvalue(implementation=MaxEigenvalue, stateless=False)
    )
    solver.add_pde_term(
        exahype2.api.InitPointWisely(implementation=init_point_wisely, stateless=False)
    )
    solver.numerical_scheme = exahype2.api.rkfd4.CentralDifferencesWithRungeKutta(
        solver, patch_size=9, rk_order=4
    )
    solver.grid_construction = exahype2.api.rkfd4.SetupSimulationOverRegularGrid(
        solver, volume_size=VolumeSize
    )
    solver.plotter = exahype2.api.rkfd4.PlotAsVTKFile(
        solver, first_plot_time_stamp=0.0, time_in_between_plots=TimeInBetweenPlots
    )
    project.add_solver(solver)
    project.set_global_simulation_parameters(
        dimensions=2,
        offset=[-0.5, -0.5],
        size=[1.0, 1.0],
        min_end_time=EndTime,
        max_end_time=EndTime,
        periodic_BC=[True, True, True],
    )
    project.set_Peano4_installation(PeanoRoot, peano4.api.output.CompileMode.RELEASE)
    peano4_project = project.generate_Peano4_project()
    peano4_project.generate()
    subprocess.run(["make", "-j"], cwd=directory, check=True, stdout=subprocess.DEVNULL)
    with open(os.path.join(directory, "run.log"), "w") as log:
        subprocess.run(["./euler"], cwd=directory, stdout=log, stderr=subprocess.STDOUT)
 
 
def read_snapshot(directory, snapshot, patch_file=PatchFile):
    """! Return the time stamp and all patches (offset, size, values) of a snapshot"""
    snapshots = tarch.plotter.api.notebook.read_snapshots(
        os.path.join(directory, patch_file)
    )
    time_stamp, patch_files = snapshots[snapshot]
    patches = [
        p for f in patch_files for p in tarch.plotter.api.notebook.read_patches(f)
    ]
    return time_stamp, patches
 
 
def number_of_snapshots(directory, patch_file=PatchFile):
    return len(
        tarch.plotter.api.notebook.read_snapshots(os.path.join(directory, patch_file))
    )
 
 
def density_along_x_axis(directory, snapshot, patch_file=PatchFile):
    """! Return time stamp plus volume centres and densities along y=0, x>=0"""
    time_stamp, patches = read_snapshot(directory, snapshot, patch_file)
    x, rho = [], []
    for offset, size, values in patches:
        if offset[1] <= 0.0 < offset[1] + size[1]:
            dof = values.shape[0]
            h = size[0] / dof
            row = int((0.0 - offset[1]) / (size[1] / dof))
            for i in range(dof):
                centre = offset[0] + (i + 0.5) * h
                if centre >= 0.0:
                    x.append(centre)
                    rho.append(values[row, i, 0])
    order = np.argsort(x)
    return time_stamp, np.array(x)[order], np.array(rho)[order]
 
 
def plot_mesh(directory):
    """! Initial density plus the patch boundaries"""
    file_name = os.path.join(directory, PatchFile)
    ax = tarch.plotter.api.notebook.plot_snapshot(file_name, 0, 0, label="density")
    _, patches = read_snapshot(directory, 0)
    for offset, size, values in patches:
        ax.add_patch(
            matplotlib.patches.Rectangle(
                offset, size[0], size[1], fill=False, edgecolor="white", linewidth=0.4
            )
        )
    ax.set_title(
        "Initial density\n{0}x{0} patches, {1}x{1} volumes per patch".format(
            round(1.0 / patches[0][1][0]), patches[0][2].shape[0]
        )
    )
    ax.figure.savefig(
        os.path.join(Here, "13_euler-with-amr-mesh.png"), dpi=100, bbox_inches="tight"
    )
    plt.close(ax.figure)
 
 
def plot_snapshots(directory):
    """! Density at three time stamps side by side"""
    file_name = os.path.join(directory, PatchFile)
    snapshots = number_of_snapshots(directory)
    fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
    for ax, snapshot in zip(
        axes, [(snapshots - 1) // 4, (snapshots - 1) // 2, snapshots - 1]
    ):
        tarch.plotter.api.notebook.plot_snapshot(
            file_name, snapshot, 0, label="density", ax=ax
        )
    fig.tight_layout()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-snapshots.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_profile(directory):
    """! Density along the positive x-axis at a few time stamps"""
    fig, ax = plt.subplots(figsize=(7, 4))
    snapshots = number_of_snapshots(directory)
    for snapshot, colour in zip(
        [0, (snapshots - 1) // 4, (snapshots - 1) // 2, snapshots - 1],
        ["C0", "C1", "C2", "C3"],
    ):
        t, x, rho = density_along_x_axis(directory, snapshot)
        ax.plot(
            x,
            rho,
            "o-",
            color=colour,
            markersize=2,
            linewidth=0.8,
            label="t={:.2f}".format(t),
        )
    ax.set_xlabel("distance r from the centre (along y=0)")
    ax.set_ylabel("density")
    ax.set_title("Density profile of the point explosion")
    ax.legend()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-profile.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_active_patches(directory):
    """!
 
    Highlight the patches which hold steep density gradients at the end
 
    A patch counts as active if the largest density jump between two
    neighbouring volumes within the patch exceeds 20% of the largest such
    jump in the whole domain. This is the kind of quantity a refinement
    criterion would look at.
 
    """
    file_name = os.path.join(directory, PatchFile)
    snapshot = number_of_snapshots(directory) - 1
    ax = tarch.plotter.api.notebook.plot_snapshot(
        file_name, snapshot, 0, label="density"
    )
    t, patches = read_snapshot(directory, snapshot)
 
    def max_jump(values):
        rho = values[..., 0]
        return max(
            np.abs(np.diff(rho, axis=0)).max(), np.abs(np.diff(rho, axis=1)).max()
        )
 
    jumps = [max_jump(values) for _, _, values in patches]
    threshold = 0.2 * max(jumps)
    active = 0
    for (offset, size, values), jump in zip(patches, jumps):
        if jump > threshold:
            active += 1
            ax.add_patch(
                matplotlib.patches.Rectangle(
                    offset,
                    size[0],
                    size[1],
                    facecolor=(1.0, 0.0, 0.0, 0.35),
                    edgecolor="red",
                    linewidth=0.5,
                )
            )
    ax.set_title(
        "t={:.2f}: {} of {} patches ({:.0f}%) hold steep gradients".format(
            t, active, len(patches), 100.0 * active / len(patches)
        )
    )
    ax.figure.savefig(
        os.path.join(Here, "13_euler-with-amr-active-patches.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(ax.figure)
    print("{} of {} patches are active".format(active, len(patches)))
 
 
def plot_adaptive_mesh(directory):
    """! Initial density on the adaptive mesh plus the patch boundaries"""
    file_name = os.path.join(directory, PatchFile)
    ax = tarch.plotter.api.notebook.plot_snapshot(file_name, 0, 0, label="density")
    _, patches = read_snapshot(directory, 0)
    for offset, size, values in patches:
        ax.add_patch(
            matplotlib.patches.Rectangle(
                offset, size[0], size[1], fill=False, edgecolor="white", linewidth=0.3
            )
        )
    sizes = sorted({round(size[0], 6) for _, size, _ in patches}, reverse=True)
    counts = [sum(1 for _, size, _ in patches if round(size[0], 6) == h) for h in sizes]
    ax.set_title(
        "Initial density on the adaptive mesh\n"
        + ", ".join(
            "{} patches of size 1/{}".format(count, round(1.0 / h))
            for h, count in zip(sizes, counts)
        ),
        fontsize=9,
    )
    ax.figure.savefig(
        os.path.join(Here, "13_euler-with-amr-adaptive-mesh.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(ax.figure)
    print("adaptive mesh: {}".format(dict(zip(sizes, counts))))
 
 
def plot_adaptive_snapshots(directory):
    """! Density at three time stamps on the adaptive mesh plus patch boundaries"""
    file_name = os.path.join(directory, PatchFile)
    snapshots = number_of_snapshots(directory)
    fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
    for ax, snapshot in zip(
        axes, [(snapshots - 1) // 4, (snapshots - 1) // 2, snapshots - 1]
    ):
        tarch.plotter.api.notebook.plot_snapshot(
            file_name, snapshot, 0, label="density", ax=ax
        )
        _, patches = read_snapshot(directory, snapshot)
        for offset, size, values in patches:
            ax.add_patch(
                matplotlib.patches.Rectangle(
                    offset,
                    size[0],
                    size[1],
                    fill=False,
                    edgecolor="white",
                    linewidth=0.2,
                )
            )
    fig.tight_layout()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-adaptive-snapshots.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_profile_comparison(regular_directory, adaptive_directory):
    """! Density along the positive x-axis: regular against adaptive mesh"""
    fig, ax = plt.subplots(figsize=(7, 4))
    snapshots = number_of_snapshots(regular_directory)
    assert snapshots == number_of_snapshots(adaptive_directory)
    for snapshot, colour in zip(
        [(snapshots - 1) // 4, (snapshots - 1) // 2, snapshots - 1],
        ["C1", "C2", "C3"],
    ):
        t, x, rho = density_along_x_axis(regular_directory, snapshot)
        ax.plot(
            x,
            rho,
            "-",
            color=colour,
            linewidth=1.0,
            label="regular, t={:.2f}".format(t),
        )
        t, x, rho = density_along_x_axis(adaptive_directory, snapshot)
        ax.plot(
            x,
            rho,
            "o",
            color=colour,
            markersize=1.5,
            label="adaptive, t={:.2f}".format(t),
        )
    ax.set_xlabel("distance r from the centre (along y=0)")
    ax.set_ylabel("density")
    ax.set_title("Regular mesh (lines) against adaptive mesh (dots)")
    ax.legend(fontsize=8)
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-profile-comparison.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def print_conservation(directory):
    """! Total mass and energy of the first and the last snapshot"""
    for snapshot in [0, number_of_snapshots(directory) - 1]:
        t, patches = read_snapshot(directory, snapshot)
        mass = sum(
            values[..., 0].sum() * (size[0] / values.shape[0]) ** 2
            for _, size, values in patches
        )
        energy = sum(
            values[..., 3].sum() * (size[0] / values.shape[0]) ** 2
            for _, size, values in patches
        )
        print("t={:.3f}: mass={:.8f}, energy={:.8f}".format(t, mass, energy))
 
 
def plot_setup_task_graph():
    """! Construct the setup task graph with the refinement criterion and draw it"""
    peano4.api.Logger().file_name = os.path.join(
        tempfile.gettempdir(), "13-euler-with-amr-log.md"
    )
    with tempfile.TemporaryDirectory() as directory:
        peano_project = peano4.api.Project(
            ["illustration"], "illustration", directory=directory
        )
        solver = create_solver(adaptive=True)
        exahype2_project = types.SimpleNamespace(periodic_BC=[True, True, True])
        setup_simulation = peano4.api.taskgraph.Graph()
        solver.grid_construction.construct_output(
            exahype2_project, peano_project, setup_simulation
        )
    setup_simulation.visualise(
        title="Setup task graph with point-wise refinement criterion"
    )
    output = os.path.join(Here, "13_euler-with-amr-setup-task-graph.png")
    plt.gcf().savefig(output, dpi=90)
    plt.close(plt.gcf())
    print("Wrote", output)
 
 
def patch_counts(directory):
    """! Time stamps and numbers of patches of all snapshots"""
    times, counts = [], []
    for snapshot in range(number_of_snapshots(directory)):
        t, patches = read_snapshot(directory, snapshot)
        times.append(t)
        counts.append(len(patches))
    return times, counts
 
 
def plot_dynamic_patch_count(directory):
    """! Number of patches per snapshot of the dynamically adaptive run"""
    times, counts = patch_counts(directory)
    fig, ax = plt.subplots(figsize=(7, 3.5))
    ax.step(times, counts, where="post", marker="o", markersize=3)
    ax.set_xlabel("time")
    ax.set_ylabel("number of patches")
    ax.set_title("Patches of the dynamically adaptive mesh")
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-dynamic-patch-count.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
    print("patch counts: {}".format(list(zip([round(t, 3) for t in times], counts))))
 
 
def plot_dynamic_mesh_change(directory):
    """! Density and patch boundaries right before the first mesh change and at the end"""
    times, counts = patch_counts(directory)
    changes = [i for i in range(1, len(counts)) if counts[i] != counts[i - 1]]
    if not changes:
        print("the dynamically adaptive mesh has not changed")
        return
    file_name = os.path.join(directory, PatchFile)
    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
    for ax, snapshot in zip(axes, [changes[0] - 1, len(counts) - 1]):
        tarch.plotter.api.notebook.plot_snapshot(
            file_name, snapshot, 0, label="density", ax=ax, clim=(0.0, 0.4)
        )
        _, patches = read_snapshot(directory, snapshot)
        for offset, size, values in patches:
            ax.add_patch(
                matplotlib.patches.Rectangle(
                    offset,
                    size[0],
                    size[1],
                    fill=False,
                    edgecolor="white",
                    linewidth=0.3,
                )
            )
        ax.set_xlim(-0.3, 0.3)
        ax.set_ylim(-0.3, 0.3)
        ax.set_title(
            "t={:.2f}: {} patches".format(times[snapshot], counts[snapshot]),
            fontsize=10,
        )
    fig.tight_layout()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-dynamic-mesh-change.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def gradient_snapshots(directory):
    """!
 
    Snapshots for the comparison of the point-wise and the gradient-based
    criterion: t=0.1, 0.2 and 0.25. We stop before t=0.3, as the shock
    then reaches the periodic boundary.
 
    """
    times, _ = patch_counts(directory)
    return [
        min(range(len(times)), key=lambda i: abs(times[i] - t))
        for t in [0.1, 0.2, 0.25]
    ]
 
 
def plot_gradient_profile_comparison(pointwise_directory, gradient_directory):
    """! Density along the positive x-axis: point-wise against gradient-based criterion"""
    fig, ax = plt.subplots(figsize=(7, 4))
    for snapshot, colour in zip(
        gradient_snapshots(gradient_directory), ["C1", "C2", "C3"]
    ):
        t, x, rho = density_along_x_axis(pointwise_directory, snapshot)
        ax.plot(
            x,
            rho,
            "-",
            color=colour,
            linewidth=1.0,
            label="point-wise, t={:.2f}".format(t),
        )
        t, x, rho = density_along_x_axis(gradient_directory, snapshot)
        ax.plot(
            x,
            rho,
            "o",
            color=colour,
            markersize=1.5,
            label="gradient, t={:.2f}".format(t),
        )
    ax.set_xlabel("distance r from the centre (along y=0)")
    ax.set_ylabel("density")
    ax.set_title("Point-wise criterion (lines) against gradient criterion (dots)")
    ax.legend(fontsize=8)
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-gradient-profile.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_gradient_patch_counts(pointwise_directory, gradient_directory):
    """! Number of patches per snapshot of both dynamically adaptive runs"""
    fig, ax = plt.subplots(figsize=(7, 3.5))
    for directory, label in [
        (pointwise_directory, "point-wise criterion"),
        (gradient_directory, "gradient criterion"),
    ]:
        times, counts = patch_counts(directory)
        ax.step(times, counts, where="post", marker="o", markersize=3, label=label)
        print("{}: {}".format(label, list(zip([round(t, 3) for t in times], counts))))
    ax.set_xlabel("time")
    ax.set_ylabel("number of patches")
    ax.set_title("Patches of the dynamically adaptive meshes")
    ax.legend()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-gradient-patch-count.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_gradient_mesh(pointwise_directory, gradient_directory):
    """!
 
    Density and patch boundaries of both runs at t=0.2
 
    We zoom into the upper right quadrant, as the boundaries of the fine
    patches would otherwise hide the solution.
 
    """
    snapshot = gradient_snapshots(gradient_directory)[1]
    fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
    for ax, directory, label in zip(
        axes,
        [pointwise_directory, gradient_directory],
        ["point-wise criterion", "gradient criterion"],
    ):
        tarch.plotter.api.notebook.plot_snapshot(
            os.path.join(directory, PatchFile),
            snapshot,
            0,
            label="density",
            ax=ax,
            clim=(0.0, 0.3),
        )
        t, patches = read_snapshot(directory, snapshot)
        for offset, size, values in patches:
            ax.add_patch(
                matplotlib.patches.Rectangle(
                    offset,
                    size[0],
                    size[1],
                    fill=False,
                    edgecolor="white",
                    linewidth=0.4,
                )
            )
        ax.set_xlim(0.0, 0.5)
        ax.set_ylim(0.0, 0.5)
        ax.set_title(
            "{}, t={:.2f}: {} patches".format(label, t, len(patches)), fontsize=10
        )
    fig.tight_layout()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-gradient-mesh.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def create_gradient_pictures(dynamic_directory, gradient_directory, skip_run):
    """!
 
    Run the gradient-based criterion and compare it to the point-wise one
 
    Expects the point-wise dynamically adaptive run within dynamic_directory.
 
    """
    if not skip_run:
        print("build and run with gradient criterion in {}".format(gradient_directory))
        build_and_run(gradient_directory, adaptive=True, gradient=True)
    plot_gradient_profile_comparison(dynamic_directory, gradient_directory)
    plot_gradient_patch_counts(dynamic_directory, gradient_directory)
    plot_gradient_mesh(dynamic_directory, gradient_directory)
    print_conservation(gradient_directory)
 
 
def plot_fd4_explosion_time_step_sizes(directory):
    """! Time step sizes of the fourth-order Finite Differences for the point explosion"""
    times, sizes = [], []
    with open(os.path.join(directory, "run.log")) as log:
        for line in log:
            match = re.search(r"\bt=([0-9.eE+-]+), new dt=([0-9.eE+-]+)", line)
            if match:
                times.append(float(match.group(1)))
                sizes.append(float(match.group(2)))
    fig, ax = plt.subplots(figsize=(7, 3.5))
    ax.semilogy(times, sizes, marker="o", markersize=2, linewidth=0.8)
    ax.set_xlabel("time")
    ax.set_ylabel("time step size")
    ax.set_title("Fourth-order Finite Differences for the point explosion")
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-fd4-explosion-time-step-size.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
    print(
        "FD4 explosion: {} time steps, last t={}, first dt={}, last dt={}".format(
            len(times), times[-1], sizes[0], sizes[-1]
        )
    )
 
 
def plot_fd4_snapshots(directory):
    """! Density of the Gaussian pulse at three time stamps side by side"""
    file_name = os.path.join(directory, FD4PatchFile)
    snapshots = number_of_snapshots(directory, FD4PatchFile)
    fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
    for ax, snapshot in zip(
        axes, [(snapshots - 1) // 4, (snapshots - 1) // 2, snapshots - 1]
    ):
        tarch.plotter.api.notebook.plot_snapshot(
            file_name, snapshot, 0, label="density", ax=ax
        )
    fig.tight_layout()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-fd4-snapshots.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_fd4_profile(directory):
    """! Density of the Gaussian pulse along the positive x-axis at a few time stamps"""
    fig, ax = plt.subplots(figsize=(7, 4))
    snapshots = number_of_snapshots(directory, FD4PatchFile)
    for snapshot, colour in zip(
        [0, (snapshots - 1) // 4, (snapshots - 1) // 2, snapshots - 1],
        ["C0", "C1", "C2", "C3"],
    ):
        t, x, rho = density_along_x_axis(directory, snapshot, FD4PatchFile)
        ax.plot(
            x,
            rho,
            "o-",
            color=colour,
            markersize=2,
            linewidth=0.8,
            label="t={:.2f}".format(t),
        )
        print(
            "FD4 pulse t={:.2f}: density in [{:.4f},{:.4f}], maximum at r={:.3f}".format(
                t, rho.min(), rho.max(), x[np.argmax(rho)]
            )
        )
    ax.set_xlabel("distance r from the centre (along y=0)")
    ax.set_ylabel("density")
    ax.set_title("Density profile of the Gaussian pulse (FD4)")
    ax.legend()
    fig.savefig(
        os.path.join(Here, "13_euler-with-amr-fd4-profile.png"),
        dpi=100,
        bbox_inches="tight",
    )
    plt.close(fig)
 
 
def plot_mesh_adaptation_task_graph():
    """! Construct the mesh adaptation task graph of the dynamically adaptive setup and draw it"""
    peano4.api.Logger().file_name = os.path.join(
        tempfile.gettempdir(), "13-euler-with-amr-log.md"
    )
    with tempfile.TemporaryDirectory() as directory:
        exahype2_project = exahype2.api.Project(
            ["illustration"], "illustration", "illustration", directory=directory
        )
        exahype2_project.set_global_simulation_parameters(
            dimensions=2,
            offset=[-0.5, -0.5],
            size=[1.0, 1.0],
            min_end_time=DynamicEndTime,
            periodic_BC=[True, True, True],
        )
        peano_project = peano4.api.Project(
            ["illustration"], "illustration", directory=directory
        )
        solver = create_solver(adaptive=True, dynamic=True)
        graphs = [peano4.api.taskgraph.Graph() for _ in range(4)]
        solver.generate_solver_realisation(exahype2_project, peano_project, *graphs)
    graphs[3].visualise(title="Mesh adaptation and load balancing task graph")
    output = os.path.join(Here, "13_euler-with-amr-mesh-adaptation-task-graph.png")
    plt.gcf().savefig(output, dpi=90)
    plt.close(plt.gcf())
    print("Wrote", output)
 
 
if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description="Create the PNGs of the tutorial 13_euler-with-amr.ipynb"
    )
    parser.add_argument(
        "--work-directory",
        default=None,
        help="Directory for builds and outputs (default: temporary directory)",
    )
    parser.add_argument(
        "--skip-run",
        action="store_true",
        help="Reuse the output of a previous run within the work directory",
    )
    parser.add_argument(
        "--fd4-only",
        action="store_true",
        help="Create only the pictures of the section on the fourth-order Finite Differences",
    )
    parser.add_argument(
        "--gradient-only",
        action="store_true",
        help="Create only the pictures of the section on the gradient-based refinement",
    )
    args = parser.parse_args()
 
    work_directory = args.work_directory or tempfile.mkdtemp(
        prefix="13_euler-with-amr-"
    )
    dynamic_directory = os.path.join(work_directory, "dynamic-mesh")
    gradient_directory = os.path.join(work_directory, "gradient-mesh")
    if args.gradient_only:
        os.makedirs(dynamic_directory, exist_ok=True)
        os.makedirs(gradient_directory, exist_ok=True)
        if not args.skip_run:
            print(
                "build and run over dynamically adaptive mesh in {}".format(
                    dynamic_directory
                )
            )
            build_and_run(dynamic_directory, adaptive=True, dynamic=True)
        create_gradient_pictures(dynamic_directory, gradient_directory, args.skip_run)
        print("wrote PNGs into {}".format(Here))
        sys.exit(0)
 
    fd4_explosion_directory = os.path.join(work_directory, "fd4-explosion")
    fd4_pulse_directory = os.path.join(work_directory, "fd4-pulse")
    os.makedirs(fd4_explosion_directory, exist_ok=True)
    os.makedirs(fd4_pulse_directory, exist_ok=True)
    if not args.skip_run:
        print(
            "build and run FD4 in {} and {}".format(
                fd4_explosion_directory, fd4_pulse_directory
            )
        )
        build_and_run_fd4(fd4_explosion_directory, InitPointWisely)
        build_and_run_fd4(fd4_pulse_directory, FD4InitPointWisely)
    plot_fd4_explosion_time_step_sizes(fd4_explosion_directory)
    plot_fd4_snapshots(fd4_pulse_directory)
    plot_fd4_profile(fd4_pulse_directory)
    if args.fd4_only:
        print("wrote PNGs into {}".format(Here))
        sys.exit(0)
 
    directory = os.path.join(work_directory, "regular-grid")
    os.makedirs(directory, exist_ok=True)
    if not args.skip_run:
        print("build and run in {}".format(directory))
        build_and_run(directory)
 
    plot_mesh(directory)
    plot_snapshots(directory)
    plot_profile(directory)
    plot_active_patches(directory)
 
    adaptive_directory = os.path.join(work_directory, "adaptive-mesh")
    os.makedirs(adaptive_directory, exist_ok=True)
    if not args.skip_run:
        print("build and run over adaptive mesh in {}".format(adaptive_directory))
        build_and_run(adaptive_directory, adaptive=True)
    plot_adaptive_mesh(adaptive_directory)
    plot_adaptive_snapshots(adaptive_directory)
    plot_profile_comparison(directory, adaptive_directory)
    print_conservation(adaptive_directory)
    plot_setup_task_graph()
 
    os.makedirs(dynamic_directory, exist_ok=True)
    if not args.skip_run:
        print(
            "build and run over dynamically adaptive mesh in {}".format(
                dynamic_directory
            )
        )
        build_and_run(dynamic_directory, adaptive=True, dynamic=True)
    plot_dynamic_patch_count(dynamic_directory)
    plot_dynamic_mesh_change(dynamic_directory)
    print_conservation(dynamic_directory)
    plot_mesh_adaptation_task_graph()
 
    os.makedirs(gradient_directory, exist_ok=True)
    create_gradient_pictures(dynamic_directory, gradient_directory, args.skip_run)
    print("wrote PNGs into {}".format(Here))