Minimal resolution we want to look at for the MHD cloud collapse

examples/MHDTests/MinimalResMagnetisedCloudCollapse/README.md

+Magnetised Cloud Collapse
+
+---------------------------------------------------------------------
+
+Test case for investigation of coupling of MHD to gravity.
+
+The test has been ivestigated with the Direct Induction (Price 2018)
+formulation of MHD, and default parameters in this directory have been
+chosen to best couple with a Wendland-C2 kernel. Note the need for
+a barotropic equation of state and an adiabatic index of 4/3. We
+thus suggest you compile SWIFT as :
+
+./configure --with-spmhd=direct-induction --with-kernel=wendland-C2 --with-equation-of-state=barotropic-gas --with-adiabatic-index=4/3
+
+to run this test.
+
+Generate ICs with makeIC.py. You can provide as a command line
+argument the particle number on a side of the cube out of which
+we carve the collapsing spherical cloud (all other attributes
+of the set up are computed automatically) e.g. :
+
+python3 makeIC.py -n 128 
+
+for a high resolution setup. The script also prints an
+indicative value of the gravitational softening length necessary to
+resolve the densities of physical interest; addapt accordingly in
+the parameter file.
+
+The simulation can then be run as :
+
+../../../swift -s -G --threads=12 magnetised_cloud.yml
+
+Results for snapshot snap.hdf5 can then be plotted and storred in
+im.png as :
+
+python3 plotSolution.py snap.hdf5 im.png
+
+The directory also contains a bash script that does all this for you.

examples/MHDTests/MinimalResMagnetisedCloudCollapse/getGlass.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_16.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_64.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_128.hdf5
+

examples/MHDTests/MinimalResMagnetisedCloudCollapse/magnetised_cloud.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841e33       # Grams
+  UnitLength_in_cgs:   3.08567758149e18 # Centimeters
+  UnitVelocity_in_cgs: 1e5              # Centimeters per second
+  UnitCurrent_in_cgs:  1e10             # 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:   0.05       # The end time of the simulation (in internal units).
+  dt_min:     1e-14      # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3       # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            magnetised_cloud # Common part of the name of output files
+  time_first:          0.               # Time of the first output (in internal units)
+  delta_time:          0.0025           # Time difference between consecutive outputs (in internal units)
+  compression:         1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          0.00001 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target h in units of the mean inter-particle separation
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  h_max:		 0.1
+
+# Parameters for MHD scheme
+MHD:
+  monopole_subtraction:   1.0
+  artificial_diffusion:   1.0
+  hyperbolic_dedner:      1.0
+  hyperbolic_dedner_divv: 0.5
+  parabolic_dedner:       1.0
+  resistive_eta:          0.0
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                0.025
+  MAC:                adaptive
+  theta_cr:           0.7
+  epsilon_fmm:        0.001
+  max_physical_baryon_softening: 0.00001    # Physical softening length (in internal units).
+  mesh_side_length:   64
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./magnetised_cloud.hdf5       # The file to read
+  periodic:   1
+
+# Parameters related to the barotropic equation of state
+EoS:
+  barotropic_vacuum_sound_speed:   0.2             # Internal system of units (aka. km/s)
+  barotropic_core_density:         1.47756194451e8 # Internal system of units (aka. Msun / pc**3)

examples/MHDTests/MinimalResMagnetisedCloudCollapse/makeIC.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+import h5py
+import argparse as ap
+from numpy import *
+from scipy.spatial.transform import Rotation
+
+# Constants
+G = 6.67430e-8
+mu0 = 1.25663706127e-1
+c1 = 0.53
+
+# Parameters (taken from Hopkins 2016)
+Rcloud = 4.628516371e16
+Lbox = 10.0 * Rcloud
+
+gamma = 4.0 / 3.0  # Gas adiabatic index
+
+M = 1.99e33  # total mass of the sphere
+T = 4.7e5  # initial orbital period in years
+Omega = 2 * pi / (T * 3.1536e7)  # initial angular frequency of cloud
+
+mu = (
+    10
+)  # mass to magnetic field flux through sphere, normalised to a critical value for collapse. Refer to e.g. Henebelle & Fromang 2008 for details.
+Bini = 3.0 / c1 * sqrt(mu0 * G / 5.0) * M / (Rcloud * Rcloud) * 1 / mu
+
+# Barotropic EoS parameters
+cs0 = 2e4
+inv_rho_c = 1e14
+
+# Attributes of cloud and ambient medium
+volume_cloud = (4 / 3) * pi * Rcloud ** 3
+volume_cloud_box = (2 * Rcloud) ** 3
+volume_sim_box = Lbox ** 3
+
+rho_in = M / volume_cloud
+rho_out_to_rho_in = 1 / 360
+rho_out = rho_out_to_rho_in * rho_in
+
+P_in = rho_in * cs0 * cs0 * sqrt(1.0 + (rho_in * inv_rho_c) ** gamma)
+P_out = rho_out * cs0 * cs0 * sqrt(1.0 + (rho_out * inv_rho_c) ** gamma)
+
+# Read glass files
+fileName = "magnetised_cloud.hdf5"
+
+glass = h5py.File("glassCube_16.hdf5", "r")
+pos_gf = glass["/PartType0/Coordinates"][:, :]
+h_gf = glass["/PartType0/SmoothingLength"][:]
+
+# Position cloud and ambient medium particles
+cloud_box_side = 2.0 * Rcloud
+atmosphere_box_side = (1.0 / cbrt(rho_out_to_rho_in)) * cloud_box_side
+
+pos_in = cloud_box_side * pos_gf
+h_in = cloud_box_side * h_gf
+
+pos_in -= 0.5 * cloud_box_side
+
+mask_in = pos_in[:, 0] ** 2 + pos_in[:, 1] ** 2 + pos_in[:, 2] ** 2 < Rcloud * Rcloud
+
+pos_in = pos_in[mask_in]
+h_in = h_in[mask_in]
+
+numPart_in = int(len(h_in))
+
+pos_out = atmosphere_box_side * pos_gf
+h_out = atmosphere_box_side * h_gf
+
+pos_out -= 0.5 * atmosphere_box_side
+
+mask_out = (
+    (pos_out[:, 0] ** 2 + pos_out[:, 1] ** 2 + pos_out[:, 2] ** 2 > Rcloud * Rcloud)
+    & (abs(pos_out[:, 0]) < 0.5 * Lbox)
+    & (abs(pos_out[:, 1]) < 0.5 * Lbox)
+    & (abs(pos_out[:, 2]) < 0.5 * Lbox)
+)
+
+pos_out = pos_out[mask_out]
+h_out = h_out[mask_out]
+
+pos = concatenate((pos_in, pos_out), axis=0)
+h = concatenate((h_in, h_out), axis=0)
+
+numPart = int(len(h))
+
+# Solid body rotation for cloud particles
+mask = pos[:, 0] ** 2 + pos[:, 1] ** 2 + pos[:, 2] ** 2 < Rcloud * Rcloud
+
+x = pos[:, 0]
+y = pos[:, 1]
+
+R = sqrt(x * x + y * y)
+
+phi = arctan2(y, x)
+cos_phi = cos(phi)
+sin_phi = sin(phi)
+
+v = zeros((numPart, 3))
+v[mask][:, 0] = -Omega * R[mask] * sin_phi[mask]
+v[mask][:, 1] = Omega * R[mask] * cos_phi[mask]
+
+pos += 0.5 * Lbox
+
+# Other attributes
+ids = linspace(1, numPart, numPart)
+m = ones(numPart) * M / numPart_in
+
+u = ones(numPart)
+u[mask] *= P_in / ((gamma - 1) * rho_in)
+u[~mask] *= P_out / ((gamma - 1) * rho_out)
+
+B = zeros((numPart, 3))
+B[:, 2] = Bini
+
+epsilon_lim = cbrt(M / (numPart_in * 1e-11)) / 3.086e18
+print(epsilon_lim)
+print(
+    "The softening length you need to correctly resolve densities up to 1e-11 g cm^-3 is %f pc"
+    % epsilon_lim
+)
+
+# File
+file = h5py.File(fileName, "w")
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = [Lbox, Lbox, Lbox]
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+grp.attrs["Dimension"] = 3
+
+# Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.0
+grp.attrs["Unit mass in cgs (U_M)"] = 1.0
+grp.attrs["Unit time in cgs (U_t)"] = 1.0
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+# Particle group
+grp = file.create_group("/PartType0")
+grp.create_dataset("Coordinates", data=pos, dtype="d")
+grp.create_dataset("Velocities", data=v, dtype="f")
+grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("SmoothingLength", data=h, dtype="f")
+grp.create_dataset("InternalEnergy", data=u, dtype="f")
+grp.create_dataset("ParticleIDs", data=ids, dtype="L")
+grp.create_dataset("MagneticFluxDensities", data=B, dtype="f")
+
+file.close()

examples/MHDTests/MinimalResMagnetisedCloudCollapse/plotSolution.py

+import numpy as np
+import h5py
+import sys
+import unyt
+import matplotlib.pyplot as plt
+
+from swiftsimio import load
+from swiftsimio.visualisation.rotation import rotation_matrix_from_vector
+from swiftsimio.visualisation.slice import slice_gas
+
+filename = sys.argv[1]
+data = load(filename)
+center = 0.5 * data.metadata.boxsize
+
+R0 = 0.015 * 3.086e18 * unyt.cm
+
+tff = 3e4 * 3.156e7 * unyt.s
+tsim = data.metadata.time
+tplot = tsim / tff
+print("Showing results at %2f free fall times" % tplot)
+
+# First create a mass-weighted temperature dataset
+r = data.gas.coordinates - center
+v = data.gas.velocities
+norm_r = np.sqrt(r[:, 0] ** 2 + r[:, 1] ** 2 + r[:, 2] ** 2)
+vr = np.sum(v * r, axis=1) / norm_r
+
+B = data.gas.magnetic_flux_densities
+normB = np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2 + B[:, 2] ** 2)
+divB = data.gas.magnetic_divergences
+h = data.gas.smoothing_lengths
+
+data.gas.mass_weighted_densities = data.gas.masses * data.gas.densities
+data.gas.mass_weighted_error = data.gas.masses * np.log10(
+    np.maximum(h * abs(divB) / normB, 1e-6)
+)
+data.gas.mass_weighted_vr = data.gas.masses * vr
+data.gas.mass_weighted_Btheta = data.gas.masses * np.sqrt(B[:, 0] ** 2 + B[:, 1] ** 2)
+data.gas.mass_weighted_Bz = data.gas.masses * abs(B[:, 2])
+"""
+rotation_matrix = [[1,0,0],
+                   [0,1,0],
+                   [0,0,1]]
+"""
+rotation_matrix = [[0, 0, 1], [0, 1, 0], [-1, 0, 0]]
+
+visualise_region = [
+    center[0] - 0.15 * R0,
+    center[0] + 0.15 * R0,
+    center[2] - 0.15 * R0,
+    center[2] + 0.15 * R0,
+]
+
+visualise_region_zoom = [
+    center[0] - 0.015 * R0,
+    center[0] + 0.015 * R0,
+    center[2] - 0.015 * R0,
+    center[2] + 0.015 * R0,
+]
+
+common_arguments = dict(
+    data=data,
+    resolution=512,
+    parallel=True,
+    rotation_center=unyt.unyt_array(center),
+    rotation_matrix=rotation_matrix,
+    region=visualise_region,
+    z_slice=0.0 * unyt.cm,
+)
+common_arguments_zoom = dict(
+    data=data,
+    resolution=512,
+    parallel=True,
+    rotation_center=unyt.unyt_array(center),
+    rotation_matrix=rotation_matrix,
+    region=visualise_region_zoom,
+    z_slice=0.0 * unyt.cm,
+)
+
+# Map in msun / mpc^3
+mass_map = slice_gas(**common_arguments, project="masses")
+
+mass_zoom_map = slice_gas(**common_arguments_zoom, project="masses")
+
+mass_weighted_density_map = slice_gas(
+    **common_arguments, project="mass_weighted_densities"
+)
+
+mass_weighted_error_map = slice_gas(**common_arguments, project="mass_weighted_error")
+
+mass_weighted_vr_map = slice_gas(**common_arguments, project="mass_weighted_vr")
+
+mass_weighted_vr_zoom_map = slice_gas(
+    **common_arguments_zoom, project="mass_weighted_vr"
+)
+
+mass_weighted_Btheta_map = slice_gas(**common_arguments, project="mass_weighted_Btheta")
+
+mass_weighted_Bz_map = slice_gas(**common_arguments, project="mass_weighted_Bz")
+
+density_map = mass_weighted_density_map / mass_map
+error_map = mass_weighted_error_map / mass_map
+vr_map = mass_weighted_vr_map / mass_map
+
+# vr_map[vr_map.value == -np.inf] = 0.0 * unyt.cm / unyt.s
+
+vr_zoom_map = mass_weighted_vr_zoom_map / mass_zoom_map
+Btheta_map = mass_weighted_Btheta_map / mass_map
+Bz_map = mass_weighted_Bz_map / mass_map
+
+density_map.convert_to_units(unyt.g * unyt.cm ** (-3))
+vr_map.convert_to_units(unyt.km / unyt.s)
+vr_zoom_map.convert_to_units(unyt.km / unyt.s)
+Btheta_map.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+Bz_map.convert_to_units(1e-7 * unyt.g / (unyt.statA * unyt.s * unyt.s))
+
+fig, ax = plt.subplots(3, 2, sharey=True, figsize=(10, 15))
+
+fig.suptitle("Plotting at %.2f free fall times" % tplot)
+
+a00 = ax[0, 0].contourf(
+    np.log10(density_map.value),
+    cmap="jet",
+    extend="both",
+    levels=np.linspace(-17.0, -11.0, 100),
+)
+a01 = ax[0, 1].contourf(
+    error_map.value, cmap="jet", extend="both", levels=np.arange(-6.0, 3.0, 1.0)
+)
+a10 = ax[1, 0].contourf(
+    vr_map.value, cmap="jet", extend="both", levels=np.linspace(-1.5, 2.5, 100)
+)
+a11 = ax[1, 1].contourf(
+    vr_zoom_map.value, cmap="jet", extend="both", levels=np.linspace(-1.5, 2.5, 100)
+)
+a20 = ax[2, 0].contourf(
+    np.log10(np.maximum(Btheta_map.value, 1.0)),
+    cmap="jet",
+    extend="both",
+    levels=np.linspace(0.0, 5.0, 100),
+)
+a21 = ax[2, 1].contourf(
+    np.log10(np.maximum(Bz_map.value, 1.0)),
+    cmap="jet",
+    extend="both",
+    levels=np.linspace(1.0, 4.5, 100),
+)
+
+# locs = [10.67, 32, 53.33]
+# locs = [21.33, 64, 106.67]
+locs = [85.33, 256.00, 426.67]
+# locs   = [170.67, 512.00, 853.33]
+# locs   = [341.33, 1024, 1706.67]
+labels = [-0.1, 0.0, 0.1]
+
+# locs   = [128, 256, 384]
+# labels = [-1.0, 0.0, 1.0]
+
+for ii in range(3):
+    ax[ii, 0].set_ylabel(r"$z/R_0$")
+    ax[ii, 0].set_yticks(locs, labels)
+
+for ii in range(3):
+    for jj in range(2):
+        # R = 0.001/15 * 1707
+        # print(R)
+        # circle = plt.Circle((256,256),R, color='k', fill=False)
+        # ax[ii,jj].add_patch(circle)
+        ax[ii, jj].set_xlabel(r"$x/R_0$")
+        ax[ii, jj].set_xticks(locs, labels)
+        ax[ii, jj].set_aspect("equal")
+
+cbar1 = plt.colorbar(
+    a00, ax=ax[0, 0], fraction=0.046, pad=0.04, ticks=np.linspace(-17, -11, 7)
+)
+cbar1.set_label(r"$\mathrm{log}_{10} \rho \; [g {cm}^{-3}]$")
+cbar2 = plt.colorbar(a01, ax=ax[0, 1], fraction=0.046, pad=0.04)
+cbar2.set_label(r"$\mathrm{log}_{10} \frac{h \: \nabla \cdot B}{|B|}$")
+cbar3 = plt.colorbar(
+    a10, ax=ax[1, 0], fraction=0.046, pad=0.04, ticks=np.linspace(-1.5, 2.5, 9)
+)
+cbar3.set_label(r"$v_r \; [km s^{-1}$")
+cbar4 = plt.colorbar(
+    a11, ax=ax[1, 1], fraction=0.046, pad=0.04, ticks=np.linspace(-1.5, 2.5, 9)
+)
+cbar4.set_label(r"$v_r \; [km s^{-1}$")
+cbar5 = plt.colorbar(
+    a20, ax=ax[2, 0], fraction=0.046, pad=0.04, ticks=np.linspace(0.0, 5.0, 6)
+)
+cbar5.set_label(r"$\mathrm{log}_{10} B_\theta \; [\mu G]$")
+cbar6 = plt.colorbar(
+    a21, ax=ax[2, 1], fraction=0.046, pad=0.04, ticks=np.linspace(1.0, 4.5, 8)
+)
+cbar6.set_label(r"$\mathrm{log}_{10} B_z \; [\mu G]$")
+
+plt.tight_layout()
+
+plt.savefig(sys.argv[2])

examples/MHDTests/MinimalResMagnetisedCloudCollapse/run.sh

+#!/bin/bash
+
+if [ ! -e glassCube_16.hdf5 ]
+then
+    echo "Fetching glass cubes to be used in IC generation for the magnetised cloud collapse example ..."
+    ./getGlass.sh
+fi
+
+# Generate the initial conditions if they are not present.
+if [ ! -e magnetised_cloud.hdf5 ]
+then
+    echo "Generating initial conditions for the magnetised cloud collapse example ..."
+    python3 makeIC.py
+fi
+
+# Run SWIFT
+../../../swift --hydro --self-gravity --limiter --threads=16 magnetised_cloud.yml 2>&1 | tee output.log