diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/README b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/README new file mode 100644 index 0000000000000000000000000000000000000000..8e537fc3e69f50cdf618cd34be105b97f4d601ac --- /dev/null +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/README @@ -0,0 +1,9 @@ +# Example to test the external NFW + MN external potential implementation. +# The gas disc is initialized so that the particles are in cicular orbit with that of the NFW + MN potential. +# The initial thickness of the gaseous disk is arbitrarius, so that it readjust while runnig the simulation. +# Parameters of the MN disk should be self-explanatory from the parameter file. + +tested with the following options: + +./configure --with-hydro=gadget2 --with-ext-potential=nfw_mn --disable-hand-vec + diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/generate_MN.py b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/generate_MN.py new file mode 100644 index 0000000000000000000000000000000000000000..a231b3a3b63f59e1438047cb383ecb7483f3b705 --- /dev/null +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/generate_MN.py @@ -0,0 +1,178 @@ +############################################################################### + # This file is part of SWIFT. + # Copyright (c) 2019 Alejandro Benitez-Llambay (alejandro.b.llambay@durham.ac.uk) + # + # 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 numpy as np +from scipy import special +import matplotlib.pyplot as plt +import h5py + +G = 4.299e-6 + +def mcum(r, c): + a = np.log(1+c*r)-c*r/(1.0+c*r) + b = np.log(1+c)-c/(1.0+c) + return a/b + +def MN_Vcirc(R, z, Rd=4.0, Md=3e10, Zd=0.4): + return np.sqrt(G*Md*R**2/(R**2+(Rd+np.sqrt(z**2+Zd**2))**2)**(3./2.)) + +def get_vcirc_expo(R, Mgas=3e10, Rd=4.0): + sigma0 = Mgas/(2.0*np.pi*Rd**2) + sigma = sigma0*np.exp(-R/Rd) + y = R/(2.0*Rd) + I0 = special.i0(y) + K0 = special.k0(y) + I1 = special.i1(y) + K1 = special.k1(y) + return np.sqrt(4.0*np.pi*G*sigma0*Rd*y**2*(I0*K0-I1*K1)) + +def dP1(x): + return x*np.exp(-x) + +def dP2(x): + return 1.0/np.cosh(x)**2 + +hlittle = 0.704 +rhoc = 2.775e2 * hlittle**2 + +#halo mass and concentration +#------------------- +m200 = 1.5e12 #Msun +c = 10.0 +r200 = (m200/(200*rhoc*4./3.*np.pi))**(1./3.) +#------------------- + +#Gaseous disk +#----------------- +npart = int(1e5) #This controls the resolution +md = 7e8 +rd = 0.02*r200 +zd = 0.1*rd +#----------------- + +#MN parameters +#----------------- +MN_Md = 3e10 #Msun +MN_Rd = 4.0 #kpc +MN_Zd = 0.4 #kpc +#----------------- + +#Outputfile +#------------------------ +outputfile = 'MN_ICs.hdf5' #Output file for SWIFT +#------------------------ + +x = [] +y = [] +z = [] +phi = [] + +while(len(x) <= npart): + n = int(4.5e5) + rmin = 0 + rmax = 8 + + zmin = -6 + zmax = 6 + + r0 = rmin+(rmax-rmin)*np.random.rand(n) + z0 = zmin+(zmax-zmin)*np.random.rand(n) + + ymin = 0 + ymax1 = np.max(dP1(r0)) + ymax2 = np.max(dP2(z0)) + + y1 = ymin+(ymax1-ymin)*np.random.rand(n) + y2 = ymin+(ymax2-ymin)*np.random.rand(n) + + k, = np.where( (y1 <= dP1(r0)) & (y2 <= dP2(z0) ) ) + + phi0 = 2*np.pi*np.random.rand(len(k)) + phi.extend(phi0) + x.extend(r0[k]*np.cos(phi0)) + y.extend(r0[k]*np.sin(phi0)) + z.extend(z0[k]) + +sample = np.random.randint(0,len(x), npart) +x = np.asarray(x)[:npart] +y = np.asarray(y)[:npart] +z = np.asarray(z)[:npart] +phi = np.asarray(phi)[:npart] + +r = np.sqrt(x**2+y**2+z**2) + +#Set positions +pgas = np.ones([npart,3]) +pgas[:,0] = x*rd +pgas[:,1] = y*rd +pgas[:,2] = z*zd +vgas = np.zeros([npart,3]) +rhogas = np.zeros(npart) +ugas = np.ones(npart)*0 +mgas = np.ones(npart)*md/npart*1e-10 + +#Set velocities +Rgas = np.sqrt(pgas[:,0]**2+pgas[:,1]**2) +ksort = np.argsort(Rgas) + +vcirc_gas = get_vcirc_expo(Rgas, Mgas=md, Rd=rd) #Approximate the velocity of the disk (better than spherical) +vcirc_dm = np.sqrt(G*mcum(Rgas/r200,c)*m200/(Rgas)) +vcirc_MN = MN_Vcirc(Rgas, pgas[:,2], Md=MN_Md, Rd=MN_Rd, Zd=MN_Zd) +vcirc_tot = np.sqrt(vcirc_dm**2+vcirc_gas**2+vcirc_MN**2) + +vgas[:,0] = vcirc_tot*np.sin(phi) +vgas[:,1] = -vcirc_tot*np.cos(phi) + + +boxsize = 400.0 #Boxsize for SWIFT + +for i in range(3): + pgas[:,i] += 200 + + +with h5py.File(outputfile,'w') as snap: + header = snap.create_group('Header') + + particles = snap.create_group('PartType0') + particles.create_dataset('Coordinates',data=pgas,dtype=np.float32) + particles.create_dataset('Velocities',data=vgas,dtype=np.float32) + particles.create_dataset('InternalEnergy',data=ugas,dtype=np.float32) + particles.create_dataset('Densities',data=rhogas,dtype=np.float32) + particles.create_dataset('Masses',data=mgas,dtype=np.float32) + particles.create_dataset('SmoothingLength',data=np.ones(len(mgas)),dtype=np.float32) + particles.create_dataset('ParticleIDs',data=np.arange(len(mgas)),dtype='uint') + + header.attrs['NumPart_ThisFile'] = np.array([len(mgas),0,0,0,0,0], dtype='uint') + header.attrs['NumPart_Total'] = np.array([len(mgas),0,0,0,0,0], dtype='uint') + header.attrs['MassTable'] = np.array([0,0,0,0,0,0], dtype='double') + header.attrs['Time'] = 0 + header.attrs['Redshift'] = 0 + header.attrs['Flag_Entropy_ICs'] = 0 + header.attrs['NumFilesPerSnapshot'] = 1 + header.attrs['NumPart_Total_HighWord'] = np.array([0,0,0,0,0,0], dtype='int') + header.attrs['BoxSize'] = boxsize + + snap.close() + + + + + + + diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/getIC.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/getIC.sh new file mode 100755 index 0000000000000000000000000000000000000000..7cfd203108fc04aa6a61e340b499e6d66ccc5c03 --- /dev/null +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/getIC.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/Isolated_NFW_MN.hdf5 diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/isolated_galaxy.yml new file mode 100644 index 0000000000000000000000000000000000000000..d2ca7aae80862e1bc1e89e645c6a1bb3adfcba1b --- /dev/null +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/isolated_galaxy.yml @@ -0,0 +1,61 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.9891E43 # 10^10 solar masses + UnitLength_in_cgs: 3.08567758E21 # 1 kpc + UnitVelocity_in_cgs: 1E5 # km/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters for the self-gravity scheme +Gravity: + eta: 0.025 # Constant dimensionless multiplier for time integration. + theta: 0.7 # Opening angle (Multipole acceptance criterion). + max_physical_baryon_softening: 0.100 # Physical softening length (in internal units). + +# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 2. # The end time of the simulation (in internal units). + dt_min: 1e-6 # 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: output # Common part of the name of output files + time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units) + delta_time: 0.001 # Time difference between consecutive outputs (in internal units) + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1e-2 # Time between statistics output + time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units) + +Scheduler: + max_top_level_cells: 16 + +# Parameters related to the initial conditions +InitialConditions: + file_name: Isolated_NFW_MN.hdf5 # The file to read + periodic: 0 # Are we running with periodic ICs? + +# Parameters for the hydrodynamics scheme +SPH: + minimal_temperature: 1000. # Minimum allowed gas temperature + 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.2 # Courant-Friedrich-Levy condition for time integration. + h_min_ratio: 0.1 # Minimal smoothing in units of softening. + h_max: 100. + + +#NFWPotential: +NFW_MNPotential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) + timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration + epsilon: 0.01 # Softening size (internal units) + concentration: 10.0 # concentration of the Halo + M_200: 150.0 # M200 of the galaxy disk + critical_density: 1.37E-8 # Critical density of the Universe in internal units + Mdisk: 3.0 # Disk mass (internal units) + Rdisk: 4.0 # Disk size (internal units) + Zdisk: 0.4704911 # Disk scale-height (internal units) diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/run.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..492eddc7b7fcf2e0a192a7408bdf888f1ef5a1a6 --- /dev/null +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_NFW_MN/run.sh @@ -0,0 +1,9 @@ +#!/bin/bash + +if [ ! -e Isolated_NFW_MN.hdf5 ] +then + echo "Fetching initial conditions for the isolated galaxy example..." + ./getIC.sh +fi + +../../swift -g -G -s --threads=16 isolated_galaxy.yml