diff --git a/examples/SineWavePotential_1D/makeIC_ideal.py b/examples/SineWavePotential_1D/makeIC_ideal.py new file mode 100644 index 0000000000000000000000000000000000000000..5f41b6b1ceabf48f64af3927ef09a38ee4785db2 --- /dev/null +++ b/examples/SineWavePotential_1D/makeIC_ideal.py @@ -0,0 +1,93 @@ +############################################################################### +# 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/>. +# +############################################################################## + +# Generates a distorted 1D grid with a density profile that balances out the +# external sine wave potential if run with an isothermal equation of state. + +import numpy as np +import h5py + +# constant thermal energy +# the definition below assumes the same thermal energy is defined in const.h, +# and also that the code was configured with an adiabatic index of 5./3. +uconst = 20.2615290634 +cs2 = 2.*uconst/3. +A = 10. + +fileName = "sineWavePotential.hdf5" +numPart = 100 +boxSize = 1. + +coords = np.zeros((numPart, 3)) +v = np.zeros((numPart, 3)) +m = np.zeros(numPart) + 1. +h = np.zeros(numPart) + 2./numPart +u = np.zeros(numPart) +ids = np.arange(numPart, dtype = 'L') +rho = np.zeros(numPart) + +# first set the positions, as we try to do a reasonable volume estimate to +# set the masses +for i in range(numPart): + coords[i,0] = (i+0.5)/numPart + +V = 1./numPart +for i in range(numPart): + rho[i] = 1000. + u[i] = 1.5*(-0.5*A/np.pi*np.cos(2.*np.pi*coords[i,0]) + 2.*A) + m[i] = rho[i]*V + +#File +file = h5py.File(fileName, 'w') + +# Header +grp = file.create_group("/Header") +grp.attrs["BoxSize"] = boxSize +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"] = 1 + +#Runtime parameters +grp = file.create_group("/RuntimePars") +grp.attrs["PeriodicBoundariesOn"] = 1 + +#Units +grp = file.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = 1. +grp.attrs["Unit mass in cgs (U_M)"] = 1. +grp.attrs["Unit time in cgs (U_t)"] = 1. +grp.attrs["Unit current in cgs (U_I)"] = 1. +grp.attrs["Unit temperature in cgs (U_T)"] = 1. + +#Particle group +grp = file.create_group("/PartType0") +grp.create_dataset('Coordinates', data=coords, 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('Density', data=rho, dtype='f') + +file.close() diff --git a/examples/SineWavePotential_1D/plotSolution_ideal.py b/examples/SineWavePotential_1D/plotSolution_ideal.py new file mode 100644 index 0000000000000000000000000000000000000000..14039200fca787f9f99d2268c105931e4390322f --- /dev/null +++ b/examples/SineWavePotential_1D/plotSolution_ideal.py @@ -0,0 +1,73 @@ +############################################################################### +# 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/>. +# +############################################################################## + +# Plots some quantities for the snapshot file which is passed on as a command +# line argument (full name) + +import numpy as np +import h5py +import sys +import pylab as pl + +# these should be the same as in makeIC.py +uconst = 20.2615290634 +cs2 = 2.*uconst/3. +A = 10. + +if len(sys.argv) < 2: + print "Need to provide a filename argument!" + exit() + +fileName = sys.argv[1] + +file = h5py.File(fileName, 'r') +coords = np.array(file["/PartType0/Coordinates"]) +rho = np.array(file["/PartType0/Density"]) +u = np.array(file["/PartType0/InternalEnergy"]) +agrav = np.array(file["/PartType0/GravAcceleration"]) +m = np.array(file["/PartType0/Masses"]) +ids = np.array(file["/PartType0/ParticleIDs"]) + +P = 2.*rho*u/3. + +ids_reverse = np.argsort(ids) + +gradP = np.zeros(P.shape) +for i in range(len(P)): + iself = int(ids[i]) + corr = 0. + im1 = iself-1 + if im1 < 0: + im1 = len(P)-1 + corr = 1. + ip1 = iself+1 + if ip1 == len(P): + ip1 = 0 + corr = 1. + idxp1 = ids_reverse[ip1] + idxm1 = ids_reverse[im1] + gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]) + +fig, ax = pl.subplots(2, 2) + +ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5) +ax[0][1].plot(coords[:,0], gradP/rho, "b.") +ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5) +ax[1][1].plot(coords[:,0], m, "y.") +pl.savefig("{fileName}.png".format(fileName = fileName[:-5])) diff --git a/examples/SineWavePotential_1D/run_ideal.sh b/examples/SineWavePotential_1D/run_ideal.sh new file mode 100755 index 0000000000000000000000000000000000000000..a335a8a35ba5874fe00cae0c4ecf10901bd38c33 --- /dev/null +++ b/examples/SineWavePotential_1D/run_ideal.sh @@ -0,0 +1,14 @@ +#!/bin/bash + +if [ ! -e sineWavePotential.hdf5 ] +then + echo "Generating initial conditions for the 1D SineWavePotential example..." + python makeIC_ideal.py +fi + +../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log + +for f in sineWavePotential_*.hdf5 +do + python plotSolution_ideal.py $f +done