From 9d80be7dba9c4b0098cecabaa32bae54e9ae16db Mon Sep 17 00:00:00 2001
From: Josh Borrow <joshua.borrow@durham.ac.uk>
Date: Thu, 20 Jun 2019 10:14:55 +0100
Subject: [PATCH] Adds the blob test to examples/HydroTests/BlobTest_3D
---
examples/HydroTests/BlobTest_3D/README.md | 21 +++
examples/HydroTests/BlobTest_3D/blob.yml | 36 ++++
examples/HydroTests/BlobTest_3D/makeIC.py | 169 +++++++++++++++++++
examples/HydroTests/BlobTest_3D/makeMovie.py | 169 +++++++++++++++++++
examples/HydroTests/BlobTest_3D/run.sh | 13 ++
5 files changed, 408 insertions(+)
create mode 100644 examples/HydroTests/BlobTest_3D/README.md
create mode 100644 examples/HydroTests/BlobTest_3D/blob.yml
create mode 100644 examples/HydroTests/BlobTest_3D/makeIC.py
create mode 100644 examples/HydroTests/BlobTest_3D/makeMovie.py
create mode 100755 examples/HydroTests/BlobTest_3D/run.sh
diff --git a/examples/HydroTests/BlobTest_3D/README.md b/examples/HydroTests/BlobTest_3D/README.md
new file mode 100644
index 0000000000..2bd59abd7a
--- /dev/null
+++ b/examples/HydroTests/BlobTest_3D/README.md
@@ -0,0 +1,21 @@
+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`.
diff --git a/examples/HydroTests/BlobTest_3D/blob.yml b/examples/HydroTests/BlobTest_3D/blob.yml
new file mode 100644
index 0000000000..b319618282
--- /dev/null
+++ b/examples/HydroTests/BlobTest_3D/blob.yml
@@ -0,0 +1,36 @@
+# 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
diff --git a/examples/HydroTests/BlobTest_3D/makeIC.py b/examples/HydroTests/BlobTest_3D/makeIC.py
new file mode 100644
index 0000000000..859612f8a6
--- /dev/null
+++ b/examples/HydroTests/BlobTest_3D/makeIC.py
@@ -0,0 +1,169 @@
+"""
+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)
+
diff --git a/examples/HydroTests/BlobTest_3D/makeMovie.py b/examples/HydroTests/BlobTest_3D/makeMovie.py
new file mode 100644
index 0000000000..9ae4a538e0
--- /dev/null
+++ b/examples/HydroTests/BlobTest_3D/makeMovie.py
@@ -0,0 +1,169 @@
+"""
+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")
diff --git a/examples/HydroTests/BlobTest_3D/run.sh b/examples/HydroTests/BlobTest_3D/run.sh
new file mode 100755
index 0000000000..b59594227c
--- /dev/null
+++ b/examples/HydroTests/BlobTest_3D/run.sh
@@ -0,0 +1,13 @@
+#!/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
--
GitLab