Commit 602ad119 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'add-blob-test' into 'master'

Adds the blob test to examples/HydroTests/BlobTest_3D

See merge request !847
parents 84ac95e5 9d80be7d
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`.
# 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
"""
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)
"""
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")
#!/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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment