Adds the blob test to examples/HydroTests/BlobTest_3D

examples/HydroTests/BlobTest_3D/README.md

+Blob Test
+=========
+
+This test is very similar to the one presented in Agertz+ 2007
+(https://ui.adsabs.harvard.edu/abs/2007MNRAS.380..963A/abstract) that tests
+un-seeded fluid mixing.
+
+You can use the makeIC.py script to generate the initial conditions, and
+then the makeMovie.py script to look at the results. Note that this is a very
+dynamic test, and similarly to the kelvin-helmholtz it does _not_ have an explicit
+answer and is best viewed as a movie with more than a grain of salt.
+
+The expected results are:
+
++ For schemes with surface tension terms, the blob stays together
++ For schemes without surface tension terms, the blob breaks up.
+
+This is at much lower resolution than the original test in Agertz+ 2007.
+To change the resolution, you will need to edit the `makeIC.py` file.
+A similar resolution to the GAD\_10M presented in that paper can be gained
+by using `num_on_side=128`.

examples/HydroTests/BlobTest_3D/blob.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   10   # The end time of the simulation (in internal units).
+  dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            blob # Common part of the name of output files
+  time_first:          0.       # Time of the first output (in internal units)
+  delta_time:          0.1      # Time difference between consecutive outputs (in internal units)
+  compression:         1
+  
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./blob.hdf5       # The file to read
+  periodic:   1
+  cleanup_smoothing_lengths: 1

examples/HydroTests/BlobTest_3D/makeIC.py

+"""
+Creates the BCC ICs for the blob test.
+"""
+
+import numpy as np
+import h5py
+
+from unyt import cm, g, s, erg
+from unyt.unit_systems import cgs_unit_system
+
+from swiftsimio import Writer
+
+
+def generate_cube(num_on_side, side_length=1.0):
+    """
+    Generates a cube
+    """
+
+    values = np.linspace(0.0, side_length, num_on_side + 1)[:-1]
+
+    positions = np.empty((num_on_side ** 3, 3), dtype=float)
+
+    for x in range(num_on_side):
+        for y in range(num_on_side):
+            for z in range(num_on_side):
+                index = x * num_on_side + y * num_on_side ** 2 + z
+
+                positions[index, 0] = values[x]
+                positions[index, 1] = values[y]
+                positions[index, 2] = values[z]
+
+    return positions
+
+
+def generate_bcc_lattice(num_on_side, side_length=1.0):
+    cube = generate_cube(num_on_side // 2, side_length)
+
+    mips = side_length / num_on_side
+
+    positions = np.concatenate([cube, cube + mips * 0.5])
+
+    return positions
+
+
+def generate_outside_of_blob(
+    num_on_side, num_copies_x, side_length, blob_centre_x, blob_radius, initial_velocity
+):
+    """
+    Generates the area outside of the blob, including a
+    cut out for where the blob will fit in.
+    """
+
+    bcc_lattice = generate_bcc_lattice(num_on_side)
+
+    # We now need to duplicate.
+    bcc_lattice_full = np.concatenate(
+        [bcc_lattice + x * np.array([1.0, 0.0, 0.0]) for x in range(num_copies_x)]
+    )
+
+    # Now size it appropriately
+    bcc_lattice_full *= side_length
+
+    # Now we need to chop out the region that our blob will live in
+    dx = bcc_lattice_full - np.array(
+        [blob_centre_x, 0.5 * side_length, 0.5 * side_length]
+    )
+    squared_radius = np.sum(dx * dx, axis=1)
+
+    cut_out_mask = squared_radius > blob_radius ** 2
+
+    positions = bcc_lattice_full[cut_out_mask]
+
+    # Now we can generate the velocities
+    velocities = np.zeros_like(positions)
+    velocities[:, 0] = initial_velocity
+
+    return positions, velocities
+
+
+def generate_blob(num_on_side, side_length, blob_centre_x, blob_radius):
+    """
+    Generate the positions and velocities for the blob.
+    """
+
+    bcc_lattice = generate_bcc_lattice(num_on_side)
+
+    # Update to respect side length
+    bcc_lattice *= side_length
+
+    # Find the radii
+    squared_radius = np.sum((bcc_lattice - 0.5 * side_length) ** 2, axis=1)
+    inside_sphere = squared_radius < blob_radius * blob_radius
+
+    # Now select out particles
+    bcc_lattice = bcc_lattice[inside_sphere]
+
+    # Move to the correct x_position
+    bcc_lattice[:, 0] = bcc_lattice[:, 0] + blob_centre_x - 0.5 * side_length
+    positions = bcc_lattice
+
+    # Generate velocities
+    velocities = np.zeros_like(positions)
+
+    return positions, velocities
+
+
+def write_out_ics(
+    filename,
+    num_on_side,
+    side_length=1.0,
+    duplications=4,
+    blob_radius=0.1,
+    blob_location_x=0.5,
+    inside_density=10,
+    velocity_outside=2.7,  # Actually the mach number
+):
+    x = Writer(
+        cgs_unit_system, [side_length * duplications, side_length, side_length] * cm
+    )
+
+    positions_outside, velocities_outside = generate_outside_of_blob(
+        num_on_side,
+        duplications,
+        side_length,
+        blob_location_x,
+        blob_radius,
+        initial_velocity=1.0,
+    )
+    positions_inside, velocities_inside = generate_blob(
+        int(num_on_side * np.cbrt(inside_density)),
+        side_length,
+        blob_location_x,
+        blob_radius,
+    )
+
+    coordinates = np.concatenate([positions_inside, positions_outside])
+
+    velocities = np.concatenate([velocities_inside, velocities_outside])
+
+    x.gas.coordinates = coordinates * cm
+
+    x.gas.velocities = velocities * cm / s
+
+    x.gas.masses = (
+        np.ones(coordinates.shape[0], dtype=float) * (side_length / num_on_side) * g
+    )
+
+    # Set to be (1 / n) c_s
+    x.gas.internal_energy = (
+        np.ones(coordinates.shape[0], dtype=float)
+        * ((1.0 / velocity_outside) * x.gas.velocities[-1][0]) ** 2
+        / (5.0 / 3.0 - 1)
+    )
+
+    # Need to correct the particles inside the sphere so that we are in pressure equilibrium everywhere
+    x.gas.internal_energy[: len(positions_inside)] /= inside_density
+
+    x.gas.generate_smoothing_lengths(
+        boxsize=(duplications / 2) * side_length * cm, dimension=3
+    )
+
+    x.write(filename)
+
+    return
+
+
+if __name__ == "__main__":
+    write_out_ics(filename="blob.hdf5", num_on_side=64)
+

examples/HydroTests/BlobTest_3D/makeMovie.py

+"""
+Makes a movie of the output of the blob test.
+
+Josh Borrow (joshua.borrow@durham.ac.uk) 2019
+
+LGPLv3
+"""
+
+from swiftsimio import load
+from swiftsimio.visualisation import scatter
+
+from p_tqdm import p_map
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from matplotlib.colors import LogNorm
+from matplotlib.animation import FuncAnimation
+
+info_frames = 15
+start_frame = 0
+end_frame = 101
+resolution = 1024
+snapshot_name = "blob"
+cmap = "Spectral_r"
+text_args = dict(color="black")
+# plot = "pressure"
+# name = "Pressure $P$"
+plot = "density"
+name = "Fluid Density $\\rho$"
+
+
+def get_image(n):
+    """
+    Gets the image for snapshot n, and also returns the associated
+    SWIFT metadata object.
+    """
+    filename = f"{snapshot_name}_{n:04d}.hdf5"
+
+    data = load(filename)
+    boxsize = data.metadata.boxsize[0].value
+
+    output = np.zeros((resolution, resolution * 4), dtype=float)
+
+    x, y, _ = data.gas.coordinates.value.T
+
+    # This is an oblong box but we can only make squares!
+    for box, box_edges in enumerate([[0.0, 1.1], [0.9, 2.1], [1.9, 3.1], [2.9, 4.0]]):
+        mask = np.logical_and(x >= box_edges[0], x <= box_edges[1])
+        masked_x = x[mask] - np.float64(box)
+        masked_y = y[mask]
+
+        hsml = data.gas.smoothing_length.value[mask]
+
+        if plot == "density":
+            mass = data.gas.masses.value[mask]
+            image = scatter(x=masked_y, y=masked_x, m=mass, h=hsml, res=resolution)
+        else:
+            quantity = getattr(data.gas, plot).value[mask]
+            # Need to divide out the particle density for non-projected density quantities
+            image = scatter(
+                x=masked_y, y=masked_x, m=quantity, h=hsml, res=resolution
+            ) / scatter(
+                x=masked_y, y=masked_x, m=np.ones_like(quantity), h=hsml, res=resolution
+            )
+
+        output[:, box * resolution : (box + 1) * resolution] = image
+
+    return output, data.metadata
+
+
+def get_data_dump(metadata):
+    """
+    Gets a big data dump from the SWIFT metadata
+    """
+
+    try:
+        viscosity = metadata.viscosity_info
+    except:
+        viscosity = "No info"
+
+    try:
+        diffusion = metadata.diffusion_info
+    except:
+        diffusion = "No info"
+
+    output = (
+        "$\\bf{Blob}$ $\\bf{Test}$\n\n"
+        "$\\bf{SWIFT}$\n"
+        + metadata.code_info
+        + "\n\n"
+        + "$\\bf{Compiler}$\n"
+        + metadata.compiler_info
+        + "\n\n"
+        + "$\\bf{Hydrodynamics}$\n"
+        + metadata.hydro_info
+        + "\n\n"
+        + "$\\bf{Viscosity}$\n"
+        + viscosity
+        + "\n\n"
+        + "$\\bf{Diffusion}$\n"
+        + diffusion
+    )
+
+    return output
+
+
+def time_formatter(metadata):
+    return f"$t = {metadata.t:2.2f}$"
+
+
+# Generate the frames and unpack our variables
+images, metadata = zip(*p_map(get_image, list(range(start_frame, end_frame))))
+
+# The edges are funny because of the non-periodicity.
+central_region = images[0][
+    resolution // 10 : resolution - resolution // 10,
+    resolution // 10 : resolution - resolution // 10,
+]
+norm = LogNorm(vmin=np.min(central_region), vmax=np.max(central_region), clip="black")
+
+fig, ax = plt.subplots(figsize=(8 * 4, 8), dpi=resolution // 8)
+
+fig.subplots_adjust(0, 0, 1, 1)
+ax.axis("off")
+
+# Set up the initial state
+image = ax.imshow(np.zeros_like(images[0]), norm=norm, cmap=cmap, origin="lower")
+
+description_text = ax.text(
+    0.5,
+    0.5,
+    get_data_dump(metadata[0]),
+    va="center",
+    ha="center",
+    **text_args,
+    transform=ax.transAxes,
+)
+
+time_text = ax.text(
+    (1 - 0.025 * 0.25),
+    0.975,
+    time_formatter(metadata[0]),
+    **text_args,
+    va="top",
+    ha="right",
+    transform=ax.transAxes,
+)
+
+info_text = ax.text(
+    0.025 * 0.25, 0.975, name, **text_args, va="top", ha="left", transform=ax.transAxes
+)
+
+
+def animate(n):
+    # Display just our original frames at t < 0
+    if n >= 0:
+        image.set_array(images[n])
+        description_text.set_text("")
+        time_text.set_text(time_formatter(metadata[n]))
+
+    return (image,)
+
+
+animation = FuncAnimation(
+    fig, animate, range(start_frame - info_frames, end_frame), interval=40
+)
+
+animation.save("blob.mp4")

examples/HydroTests/BlobTest_3D/run.sh

+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e blob.hdf5 ]
+then
+    echo "Generating initial conditions for the Sod shock example..."
+    python makeIC.py
+fi
+
+# Run SWIFT
+../../swift --hydro --threads=4 blob.yml
+
+python makeMovie.py