Commit aca9d796 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Removed ideal gas specific SineWavePotential scripts and incorporated them in...

Removed ideal gas specific SineWavePotential scripts and incorporated them in the isothermal script. Ideal gas now seems to work.
parent c0205edb
......@@ -38,7 +38,7 @@ 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)
u = np.zeros(numPart) + uconst
ids = np.arange(numPart, dtype = 'L')
rho = np.zeros(numPart)
......
###############################################################################
# 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()
###############################################################################
# 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]))
#!/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
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