From 653a1eb86491a65d7fbaef45ea75a29ad89e22f0 Mon Sep 17 00:00:00 2001 From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be> Date: Fri, 3 Mar 2017 18:15:02 +0000 Subject: [PATCH] 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 | 15 ++- examples/SineWavePotential_2D/makeIC.py | 95 ++++++++++++++++ examples/SineWavePotential_2D/plotSolution.py | 83 ++++++++++++++ examples/SineWavePotential_2D/run.sh | 14 +++ .../sineWavePotential.yml | 39 +++++++ examples/SineWavePotential_3D/makeIC.py | 106 ++++++++++++++++++ examples/SineWavePotential_3D/plotSolution.py | 84 ++++++++++++++ examples/SineWavePotential_3D/run.sh | 14 +++ .../sineWavePotential.yml | 39 +++++++ 9 files changed, 486 insertions(+), 3 deletions(-) create mode 100644 examples/SineWavePotential_2D/makeIC.py create mode 100644 examples/SineWavePotential_2D/plotSolution.py create mode 100755 examples/SineWavePotential_2D/run.sh create mode 100644 examples/SineWavePotential_2D/sineWavePotential.yml create mode 100644 examples/SineWavePotential_3D/makeIC.py create mode 100644 examples/SineWavePotential_3D/plotSolution.py create mode 100755 examples/SineWavePotential_3D/run.sh create mode 100644 examples/SineWavePotential_3D/sineWavePotential.yml diff --git a/examples/SineWavePotential_1D/plotSolution.py b/examples/SineWavePotential_1D/plotSolution.py index 62cbd814dc..65e981e464 100644 --- a/examples/SineWavePotential_1D/plotSolution.py +++ b/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) diff --git a/examples/SineWavePotential_2D/makeIC.py b/examples/SineWavePotential_2D/makeIC.py new file mode 100644 index 0000000000..950dad41d9 --- /dev/null +++ b/examples/SineWavePotential_2D/makeIC.py @@ -0,0 +1,95 @@ +############################################################################### +# 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() diff --git a/examples/SineWavePotential_2D/plotSolution.py b/examples/SineWavePotential_2D/plotSolution.py new file mode 100644 index 0000000000..ee02f59c40 --- /dev/null +++ b/examples/SineWavePotential_2D/plotSolution.py @@ -0,0 +1,83 @@ +############################################################################### +# 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])) diff --git a/examples/SineWavePotential_2D/run.sh b/examples/SineWavePotential_2D/run.sh new file mode 100755 index 0000000000..077cf1c0cc --- /dev/null +++ b/examples/SineWavePotential_2D/run.sh @@ -0,0 +1,14 @@ +#!/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 diff --git a/examples/SineWavePotential_2D/sineWavePotential.yml b/examples/SineWavePotential_2D/sineWavePotential.yml new file mode 100644 index 0000000000..f5a565ce4d --- /dev/null +++ b/examples/SineWavePotential_2D/sineWavePotential.yml @@ -0,0 +1,39 @@ +# 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. diff --git a/examples/SineWavePotential_3D/makeIC.py b/examples/SineWavePotential_3D/makeIC.py new file mode 100644 index 0000000000..c205884650 --- /dev/null +++ b/examples/SineWavePotential_3D/makeIC.py @@ -0,0 +1,106 @@ +############################################################################### +# 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() diff --git a/examples/SineWavePotential_3D/plotSolution.py b/examples/SineWavePotential_3D/plotSolution.py new file mode 100644 index 0000000000..7ae5dcd2a5 --- /dev/null +++ b/examples/SineWavePotential_3D/plotSolution.py @@ -0,0 +1,84 @@ +############################################################################### +# 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])) diff --git a/examples/SineWavePotential_3D/run.sh b/examples/SineWavePotential_3D/run.sh new file mode 100755 index 0000000000..077cf1c0cc --- /dev/null +++ b/examples/SineWavePotential_3D/run.sh @@ -0,0 +1,14 @@ +#!/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 diff --git a/examples/SineWavePotential_3D/sineWavePotential.yml b/examples/SineWavePotential_3D/sineWavePotential.yml new file mode 100644 index 0000000000..391383568a --- /dev/null +++ b/examples/SineWavePotential_3D/sineWavePotential.yml @@ -0,0 +1,39 @@ +# 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. -- GitLab