Commit c0205edb by Bert Vandenbroucke

### Added scripts for a SineWavePotential test version with an ideal gas equation of state (in 1D).

parent 653a1eb8
 ############################################################################### # 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 . # ############################################################################## # 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 . # ############################################################################## # 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!