diff --git a/configure.ac b/configure.ac index 386a18c0dfa7cb160046f748865cf19990479892..d5a08c7c0863bbb6dfcdbde74d68f1fe5d909907 100644 --- a/configure.ac +++ b/configure.ac @@ -751,7 +751,7 @@ esac # External potential AC_ARG_WITH([ext-potential], [AS_HELP_STRING([--with-ext-potential=<pot>], - [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch default: none@:>@] + [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch, sine-wave default: none@:>@] )], [with_potential="$withval"], [with_potential="none"] @@ -769,6 +769,9 @@ case "$with_potential" in disc-patch) AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential]) ;; + sine-wave) + AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D]) + ;; *) AC_MSG_ERROR([Unknown external potential: $with_potential]) ;; diff --git a/examples/SineWavePotential_1D/makeIC.py b/examples/SineWavePotential_1D/makeIC.py new file mode 100644 index 0000000000000000000000000000000000000000..feab6d095d28abfa2d44291f1c7dd71aa9c7f0f4 --- /dev/null +++ b/examples/SineWavePotential_1D/makeIC.py @@ -0,0 +1,99 @@ +############################################################################### +# 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+np.random.random())/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"] = 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.py b/examples/SineWavePotential_1D/plotSolution.py new file mode 100644 index 0000000000000000000000000000000000000000..62cbd814dc9035aaca1164fdf4e4060c73b665bb --- /dev/null +++ b/examples/SineWavePotential_1D/plotSolution.py @@ -0,0 +1,68 @@ +############################################################################### +# 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"]) + +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 + +gradP = np.zeros(P.shape) +for i in range(len(P)): + im1 = i-1 + if im1 < 0: + im1 = len(P)-1 + ip1 = i+1 + if ip1 == len(P): + ip1 = 0 + gradP[i] = (P[ip1]-P[im1])/(coords[2,0]-coords[0,0]) + +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_1D/run.sh b/examples/SineWavePotential_1D/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37 --- /dev/null +++ b/examples/SineWavePotential_1D/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_1D/sineWavePotential.yml b/examples/SineWavePotential_1D/sineWavePotential.yml new file mode 100644 index 0000000000000000000000000000000000000000..1f98b47a0e3c72e77a2024aebf4dcd863ed02e9b --- /dev/null +++ b/examples/SineWavePotential_1D/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: 10. # 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: 1e-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: 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/src/potential.h b/src/potential.h index 4c20f4c7b4eaa7ca134b3d90ae790a4710a1f217..e6ab9a5bd6bd91801eae0b3f1e3d8f65778f5065 100644 --- a/src/potential.h +++ b/src/potential.h @@ -36,6 +36,8 @@ #include "./potential/isothermal/potential.h" #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH) #include "./potential/disc_patch/potential.h" +#elif defined(EXTERNAL_POTENTIAL_SINE_WAVE) +#include "./potential/sine_wave/potential.h" #else #error "Invalid choice of external potential" #endif diff --git a/src/potential/sine_wave/potential.h b/src/potential/sine_wave/potential.h new file mode 100644 index 0000000000000000000000000000000000000000..e2e2b8ffcc170c28a5facc8373a81746811a9991 --- /dev/null +++ b/src/potential/sine_wave/potential.h @@ -0,0 +1,133 @@ +/******************************************************************************* + * 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/>. + * + ******************************************************************************/ +#ifndef SWIFT_SINE_WAVE_H +#define SWIFT_SINE_WAVE_H + +/* Config parameters. */ +#include "../config.h" + +/* Some standard headers. */ +#include <float.h> +#include <math.h> + +/* Local includes. */ +#include "const.h" +#include "error.h" +#include "parser.h" +#include "part.h" +#include "physical_constants.h" +#include "space.h" +#include "units.h" + +/** + * @brief External Potential Properties - Sine wave case + */ +struct external_potential { + + /*! Amplitude of the sine wave. */ + double amplitude; + + /*! Time-step limiting factor. */ + double timestep_limit; +}; + +/** + * @brief Computes the time-step from the acceleration due to a sine wave. + * + * @param time The current time. + * @param potential The properties of the potential. + * @param phys_const The physical constants in internal units. + * @param g Pointer to the g-particle data. + */ +__attribute__((always_inline)) INLINE static float external_gravity_timestep( + double time, const struct external_potential* restrict potential, + const struct phys_const* restrict phys_const, + const struct gpart* restrict g) { + + return potential->timestep_limit; +} + +/** + * @brief Computes the gravitational acceleration along x given by the sine + * wave. + * + * @param time The current time in internal units. + * @param potential The properties of the potential. + * @param phys_const The physical constants in internal units. + * @param g Pointer to the g-particle data. + */ +__attribute__((always_inline)) INLINE static void external_gravity_acceleration( + double time, const struct external_potential* restrict potential, + const struct phys_const* restrict phys_const, struct gpart* restrict g) { + + g->a_grav[0] = potential->amplitude * sin(2. * M_PI * g->x[0]) / + phys_const->const_newton_G; +} + +/** + * @brief Computes the gravitational potential energy of a particle in the + * sine wave. + * + * @param time The current time. + * @param potential The #external_potential used in the run. + * @param phys_const Physical constants in internal units. + * @param gp Pointer to the particle data. + */ +__attribute__((always_inline)) INLINE static float +external_gravity_get_potential_energy( + double time, const struct external_potential* potential, + const struct phys_const* const phys_const, const struct gpart* gp) { + + /* this potential does not really have a potential energy */ + return 0.; +} + +/** + * @brief Initialises the external potential properties in the internal system + * of units. + * + * @param parameter_file The parsed parameter file + * @param phys_const Physical constants in internal units + * @param us The current internal system of units + * @param potential The external potential properties to initialize + */ +static INLINE void potential_init_backend( + const struct swift_params* parameter_file, + const struct phys_const* phys_const, const struct unit_system* us, + const struct space* s, struct external_potential* potential) { + + potential->amplitude = + parser_get_param_double(parameter_file, "SineWavePotential:amplitude"); + potential->timestep_limit = parser_get_param_double( + parameter_file, "SineWavePotential:timestep_limit"); +} + +/** + * @brief Prints the properties of the external potential to stdout. + * + * @param potential The external potential properties. + */ +static INLINE void potential_print_backend( + const struct external_potential* potential) { + + message("External potential is a sine wave with amplitude %g", + potential->amplitude); +} + +#endif /* SWIFT_SINE_WAVE_H */