#!/usr/bin/env python3 ############################################################################### # This file is part of SWIFT. # Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.com) # 2022 Tsang Keung Chan (chantsangkeung@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 . # ############################################################################## # ------------------------------------------------------------- # Add initial conditions for photon energies and fluxes # for 2D advection of photons. # First photon group: Top hat function with zero as the # baseline, advects in x direction # Second photon group: Top hat function with nonzero value # as the baseline, advcts in y direction. # Third photon group: Gaussian advecting diagonally # Fourth photon group: Circle moving radially from the center # ------------------------------------------------------------- import h5py import numpy as np import unyt from swiftsimio import Writer # define unit system to use unitsystem = unyt.unit_systems.cgs_unit_system # define box size boxsize = 1e10 * unitsystem["length"] # number of photon groups nPhotonGroups = 4 # filename of ICs to be generated outputfilename = "advection_2D.hdf5" def initial_condition(x): """ The initial conditions that will be advected x: particle position. 3D unyt array returns: E: photon energy for each photon group. List of scalars with size of nPhotonGroups F: photon flux for each photon group. List with size of nPhotonGroups of numpy arrays of shape (3,) """ # you can make the photon quantities unitless, the units will # already have been written down in the writer. E_list = [] F_list = [] c = unyt.c.to(unitsystem["length"] / unitsystem["time"]) # Group 1 Photons: # ------------------- in_x = 0.33 * boxsize < x[0] < 0.66 * boxsize in_y = 0.33 * boxsize < x[1] < 0.66 * boxsize if in_x and in_y: E = 1.0 else: E = 0.0 # Assuming all photons flow in only one direction # (optically thin regime, "free streaming limit"), # we have that |F| = c * E F = np.zeros(3, dtype=np.float64) F[0] = c * E E_list.append(E) F_list.append(F) # Group 2 Photons: # ------------------- in_x = 0.33 * boxsize < x[0] < 0.66 * boxsize in_y = 0.33 * boxsize < x[1] < 0.66 * boxsize if in_x and in_y: E = 2.0 else: E = 1.0 F = np.zeros(3, dtype=np.float64) F[1] = c * E E_list.append(E) F_list.append(F) # Group 3 Photons: # ------------------- sigma = 0.1 * boxsize mean = 0.5 * boxsize amplitude = 2.0 baseline = 1.0 E = ( amplitude * np.exp(-((x[0] - mean) ** 2 + (x[1] - mean) ** 2) / (2 * sigma ** 2)) + baseline ) F = np.zeros(3, dtype=np.float64) F[0] = c * E / 1.414213562 # sqrt(2) F[1] = c * E / 1.414213562 # sqrt(2) E_list.append(E) F_list.append(F) # Group 4 Photons: # ------------------- circle_radius = 0.15 * boxsize center = 0.5 * boxsize dx = x[0] - center dy = x[1] - center r = np.sqrt(dx ** 2 + dy ** 2) if r <= circle_radius: unit_vector = (dx / r, dy / r) E = 1.0 F = np.zeros(3, dtype=np.float64) F[0] = unit_vector[0] * c * E F[1] = unit_vector[1] * c * E else: E = 0.0 F = np.zeros(3, dtype=np.float64) E_list.append(E) F_list.append(F) return E_list, F_list if __name__ == "__main__": glass = h5py.File("glassPlane_128.hdf5", "r") # Read particle positions and h from the glass pos = glass["/PartType0/Coordinates"][:, :] h = glass["/PartType0/SmoothingLength"][:] glass.close() pos *= boxsize h *= boxsize numPart = np.size(h) w = Writer(unyt.unit_systems.cgs_unit_system, boxsize, dimension=2) w.gas.coordinates = pos w.gas.velocities = np.zeros((numPart, 3)) * (unyt.cm / unyt.s) w.gas.masses = np.ones(numPart, dtype=np.float64) * 1000 * unyt.g w.gas.internal_energy = ( np.ones(numPart, dtype=np.float64) * (300.0 * unyt.kb * unyt.K) / unyt.g ) # Generate initial guess for smoothing lengths based on MIPS w.gas.smoothing_length = h # If IDs are not present, this automatically generates w.write(outputfilename) # Now open file back up again and add photon groups # you can make them unitless, the units have already been # written down in the writer. In this case, it's in cgs. F = h5py.File(outputfilename, "r+") header = F["Header"] nparts = header.attrs["NumPart_ThisFile"][0] parts = F["/PartType0"] for grp in range(nPhotonGroups): dsetname = "PhotonEnergiesGroup{0:d}".format(grp + 1) energydata = np.zeros(nparts, dtype=np.float32) parts.create_dataset(dsetname, data=energydata) dsetname = "PhotonFluxesGroup{0:d}".format(grp + 1) # if dsetname not in parts.keys(): fluxdata = np.zeros((nparts, 3), dtype=np.float32) parts.create_dataset(dsetname, data=fluxdata) for p in range(nparts): E, Flux = initial_condition(pos[p]) for g in range(nPhotonGroups): Esetname = "PhotonEnergiesGroup{0:d}".format(g + 1) parts[Esetname][p] = E[g] Fsetname = "PhotonFluxesGroup{0:d}".format(g + 1) parts[Fsetname][p] = Flux[g] F.close()