Added 2D and 3D SineWavePotential tests. Need to run them at higher resolution...

Added 2D and 3D SineWavePotential tests. Need to run them at higher resolution and for longer, but they seem to produce reasonable results.

examples/SineWavePotential_1D/plotSolution.py

@@ -42,21 +42,30 @@ 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"])
 
 x = np.linspace(0., 1., 1000)
 rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
 
 P = cs2*rho
 
+ids_reverse = np.argsort(ids)
+
 gradP = np.zeros(P.shape)
 for i in range(len(P)):
-  im1 = i-1
+  iself = int(ids[i])
+  corr = 0.
+  im1 = iself-1
   if im1 < 0:
     im1 = len(P)-1
-  ip1 = i+1
+    corr = 1.
+  ip1 = iself+1
   if ip1 == len(P):
     ip1 = 0
-  gradP[i] = (P[ip1]-P[im1])/(coords[2,0]-coords[0,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)
 

examples/SineWavePotential_2D/makeIC.py

+###############################################################################
+# 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_1D = 50
+boxSize = [1., 1.]
+numPart = numPart_1D**2
+
+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_1D):
+  for j in range(numPart_1D):
+    coords[numPart_1D*i+j,0] = (i+0.5)/numPart_1D
+    coords[numPart_1D*i+j,1] = (j+0.5)/numPart_1D
+
+V = 1./numPart
+for i in range(numPart):
+  rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+  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"] = 2
+
+#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()

examples/SineWavePotential_2D/plotSolution.py

+###############################################################################
+# 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"])
+
+# ids_reverse gives the index original particle 0 now has in the particle arrays
+# and so on
+# note that this will only work if the particles do not move away too much from
+# there initial positions
+ids_reverse = np.argsort(ids)
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+n1D = int(np.sqrt(len(P)))
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+  iself = int(ids[i]/n1D)
+  jself = int(ids[i]-n1D*iself)
+  corr = 0.
+  im1 = iself-1
+  if im1 < 0:
+    im1 = n1D-1
+    corr = 1.
+  ip1 = iself+1
+  if ip1 == n1D:
+    ip1 = 0
+    corr = 1.
+  idxp1 = ids_reverse[ip1*n1D+jself]
+  idxm1 = ids_reverse[im1*n1D+jself]
+  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]+corr)
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+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]))

examples/SineWavePotential_2D/run.sh

+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+  echo "Generating initial conditions for the 1D SineWavePotential example..."
+  python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+  python plotSolution.py $f
+done

examples/SineWavePotential_2D/sineWavePotential.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.
+  UnitLength_in_cgs:   1.
+  UnitVelocity_in_cgs: 1.
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1.
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.   # 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:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sineWavePotential  # Common part of the name of output files
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          0.1          # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  sineWavePotential.hdf5       # The file to read
+ 
+# External potential parameters
+SineWavePotential:
+  amplitude: 10.
+  timestep_limit: 1.

examples/SineWavePotential_3D/makeIC.py

+###############################################################################
+# 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_1D = 20
+boxSize = [1., 1.]
+numPart = numPart_1D**3
+
+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_1D):
+#  coords[i,0] = (i+np.random.random())/numPart
+  for j in range(numPart_1D):
+    for k in range(numPart_1D):
+      coords[numPart_1D**2*i+numPart_1D*j+k,0] = (i+0.5)/numPart_1D
+      coords[numPart_1D**2*i+numPart_1D*j+k,1] = (j+0.5)/numPart_1D
+      coords[numPart_1D**2*i+numPart_1D*j+k,2] = (k+0.5)/numPart_1D
+
+V = 1./numPart
+for i in range(numPart):
+  # reasonable mass estimate (actually not really good, but better than assuming
+  # a constant volume)
+#  if i == 0:
+#    V = 0.5*(coords[1,0]-coords[-1,0]+1.)
+#  elif i == numPart-1:
+#    V = 0.5*(coords[0,0]+1.-coords[-2,0])
+#  else:
+#    V = 0.5*(coords[i+1,0] - coords[i-1,0])
+  rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+  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"] = 3
+
+#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()

examples/SineWavePotential_3D/plotSolution.py

+###############################################################################
+# 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"])
+
+# ids_reverse gives the index original particle 0 now has in the particle arrays
+# and so on
+# note that this will only work if the particles do not move away too much from
+# there initial positions
+ids_reverse = np.argsort(ids)
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+n1D = int(np.cbrt(len(P)))
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+  iself = int(ids[i]/n1D/n1D)
+  jself = int(int(ids[i]-n1D*iself)/n1D)
+  kself = int(ids[i]-n1D**2*iself-n1D*jself)
+  corr = 0.
+  im1 = iself-1
+  if im1 < 0:
+    im1 = n1D-1
+    corr = 1.
+  ip1 = iself+1
+  if ip1 == n1D:
+    ip1 = 0
+    corr = 1.
+  idxp1 = ids_reverse[ip1*n1D**2+jself*n1D+kself]
+  idxm1 = ids_reverse[im1*n1D**2+jself*n1D+kself]
+  gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]+corr)
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+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]))

examples/SineWavePotential_3D/run.sh

+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+  echo "Generating initial conditions for the 1D SineWavePotential example..."
+  python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+  python plotSolution.py $f
+done

examples/SineWavePotential_3D/sineWavePotential.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.
+  UnitLength_in_cgs:   1.
+  UnitVelocity_in_cgs: 1.
+  UnitCurrent_in_cgs:  1.
+  UnitTemp_in_cgs:     1.
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.1   # 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:     1.e-2  # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sineWavePotential  # Common part of the name of output files
+  time_first:          0.           # Time of the first output (in internal units)
+  delta_time:          0.01          # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  sineWavePotential.hdf5       # The file to read
+ 
+# External potential parameters
+SineWavePotential:
+  amplitude: 10.
+  timestep_limit: 1.