diff --git a/configure.ac b/configure.ac index bc63e5143b691f355efd0b7cf074da8d40432e20..136d607645febdabe98c9ba1e485907ef0ba9ece 100644 --- a/configure.ac +++ b/configure.ac @@ -2624,7 +2624,7 @@ esac # External potential AC_ARG_WITH([ext-potential], [AS_HELP_STRING([--with-ext-potential=<pot>], - [external potential @<:@none, point-mass, point-mass-softened, isothermal, nfw, nfw-mn, hernquist, hernquist-sdmh05, disc-patch, sine-wave, constant, default: none@:>@] + [external potential @<:@none, point-mass, point-mass-softened, isothermal, nfw, nfw-mn, hernquist, hernquist-sdmh05, disc-patch, sine-wave, MWPotential2014, constant, default: none@:>@] )], [with_potential="$withval"], [with_potential="none"] @@ -2660,6 +2660,9 @@ case "$with_potential" in point-mass-softened) AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_SOFT], [1], [Softened point-mass potential with form 1/(r^2 + softening^2).]) ;; + MWPotential2014) + AC_DEFINE([EXTERNAL_POTENTIAL_MWPotential2014], [1], [Milky-Way like potential composed of a Navarro-Frenk-White + Miyamoto-Nagai disk + Power spherical cuttoff external potential.]) + ;; constant) AC_DEFINE([EXTERNAL_POTENTIAL_CONSTANT], [1], [Constant gravitational acceleration.]) ;; diff --git a/doc/RTD/source/ExternalPotentials/index.rst b/doc/RTD/source/ExternalPotentials/index.rst index 0454da53ab5a591452bb68bc1860859ecd83dc26..14537920850d8aa513f7269b62f8c09608c38dfa 100644 --- a/doc/RTD/source/ExternalPotentials/index.rst +++ b/doc/RTD/source/ExternalPotentials/index.rst @@ -344,7 +344,46 @@ follows the definitions of `Creasey, Theuns & Bower (2013) <https://adsabs.harvard.edu/abs/2013MNRAS.429.1922C>`_ equations (16) and (17). The potential is implemented along the x-axis. +12. MWPotential2014 (``MWPotential2014``) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This potential is based on ``galpy``'s ``MWPotential2014`` from `Jo Bovy (2015) <https://ui.adsabs.harvard.edu/abs/2015ApJS..216...29B>`_ and consists in a NFW potential for the halo, an axisymmetric Miyamoto-Nagai potential for the disk and a bulge modeled by a power spherical law with exponential cut-off. The bulge is given by the density: + +:math:`\rho(r) = A \left( \frac{r_1}{r} \right)^\alpha \exp \left( - \frac{r^2}{r_c^2} \right)`, + +where :math:`A` is an amplitude, :math:`r_1` is a reference radius for amplitude, :math:`\alpha` is the inner power and :math:`r_c` is the cut-off radius. + +The resulting potential is: + +:math:`\Phi_{\mathrm{MW}}(R, z) = f_1 \Phi_{\mathrm{NFW}} + f_2 \Phi_{\mathrm{MN}} + f_3 \Phi_{\text{bulge}}`, +where :math:`R^2 = x^2 + y^2` is the projected radius and :math:`f_1`, :math:`f_2` and :math:`f_3` are three coefficients that adjust the strength of each individual component. + +The parameters of the model are: + +.. code:: YAML + + MWPotential2014Potential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] # Location of centre of potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) + timestep_mult: 0.005 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration + epsilon: 0.001 # Softening size (internal units) + concentration: 9.823403437774843 # concentration of the Halo + M_200_Msun: 147.41031542774076e10 # M200 of the galaxy disk (in M_sun) + H: 1.2778254614201471 # Hubble constant in units of km/s/Mpc + Mdisk_Msun: 6.8e10 # Mass of the disk (in M_sun) + Rdisk_kpc: 3.0 # Effective radius of the disk (in kpc) + Zdisk_kpc: 0.280 # Scale-height of the disk (in kpc) + amplitude: 1.0 # Amplitude of the bulge + r_1_kpc: 1.0 # Reference radius for amplitude of the bulge (in kpc) + alpha: 1.8 # Exponent of the power law of the bulge + r_c_kpc: 1.9 # Cut-off radius of the bulge (in kpc) + potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] #Coefficients that adjust the strength of the halo (1st component), the disk (2nd component) and the bulge (3rd component) + +Note that the default value of the "Hubble constant" here seems odd. As it +enters multiplicatively with the :math:`f_1` term, the absolute normalisation is +actually not important. + How to implement your own potential ----------------------------------- diff --git a/examples/GravityTests/Hernquist_circularorbit/makeIC.py b/examples/GravityTests/Hernquist_circularorbit/makeIC.py index 474450f0e23704bfc43730872a978107f28704e9..4c69e7c966f95ca90bce20754c95fce47c45213e 100755 --- a/examples/GravityTests/Hernquist_circularorbit/makeIC.py +++ b/examples/GravityTests/Hernquist_circularorbit/makeIC.py @@ -16,9 +16,6 @@ # along with this program. If not, see <http://www.gnu.org/licenses/>. # ################################################################################ -from galpy.potential import NFWPotential -from galpy.orbit import Orbit -from galpy.util import bovy_conversion import numpy as np import matplotlib.pyplot as plt from astropy import units @@ -27,7 +24,7 @@ import h5py as h5 C = 8.0 M_200 = 2.0 N_PARTICLES = 3 -print("Initial conditions written to 'test_nfw.hdf5'") +print("Initial conditions written to 'circularorbitshernquist.hdf5'") pos = np.zeros((3, 3)) pos[0, 2] = 50.0 diff --git a/examples/GravityTests/MWPotential2014_circularorbit/README b/examples/GravityTests/MWPotential2014_circularorbit/README new file mode 100644 index 0000000000000000000000000000000000000000..056a242ad78337e3ea883351b265782261053232 --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/README @@ -0,0 +1,24 @@ +# Context + +This example tests the MWPotential2014, a reference Milky-Way potential implemented on galpy (https://docs.galpy.org). +Details of the parameters can be found in galpy's paper, p.7: +galpy: A Python Library for Galactic Dynamics, Jo Bovy (2015), Astrophys. J. Supp., 216, 29 (<arXiv/1412.3451>) + +# How to run this example + +In a terminal at the root of the "swiftsim" directory, type: + +`./autogen.sh` + +Then, configure swift to take into account external potentials. Type: + +`./configure --with-ext-potential=MWPotential2014` + +Feel free to adapt other configurations parameters depending on your system. Then, build the code by typing: + +`make` + +Finally, to run this example, open a terminal in the present directory and type: + +`./run.sh` + diff --git a/examples/GravityTests/MWPotential2014_circularorbit/makeIC.py b/examples/GravityTests/MWPotential2014_circularorbit/makeIC.py new file mode 100755 index 0000000000000000000000000000000000000000..02356b98227a89398fa63aa335d1cdd7d0c50769 --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/makeIC.py @@ -0,0 +1,75 @@ +############################################################################### +# This file is part of SWIFT. +# Copyright (c) 2023 Darwin Roduit (darwin.roduit@epfl.ch) +# +# 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/>. +# +################################################################################ +import numpy as np +import h5py as h5 + +box_size = 1000.0 +N_PARTICLES = 3 +print("Initial conditions written to 'circular_orbits_MW.hdf5'") + +pos = np.zeros((3, 3)) +pos[0, 2] = 5.0 +pos[1, 1] = 5.0 +pos[2, 0] = 30.0 +pos += np.array( + [box_size / 2, box_size / 2, box_size / 2] +) # shifts the particles to the center of the box +vel = np.zeros((3, 3)) +vel[0, 0] = 198.5424557586175 +vel[1, 0] = 225.55900735974072 +vel[2, 1] = 188.5272441374569 + +ids = np.array([1.0, 2.0, 3.0]) +mass = np.array([1.0, 1.0, 1.0]) * 1e-10 + +# File +file = h5.File("circular_orbits_MW.hdf5", "w") + +# Units +grp = file.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = 3.086e21 +grp.attrs["Unit mass in cgs (U_M)"] = 1.98848e43 +grp.attrs["Unit time in cgs (U_t)"] = 3.086e16 +grp.attrs["Unit current in cgs (U_I)"] = 1.0 +grp.attrs["Unit temperature in cgs (U_T)"] = 1.0 + +# Header +grp = file.create_group("/Header") +grp.attrs["BoxSize"] = box_size +grp.attrs["NumPart_Total"] = [0, N_PARTICLES, 0, 0, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [0, N_PARTICLES, 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, 0, 0, 0, 0, 0] +grp.attrs["Dimension"] = 3 + +# Runtime parameters +grp = file.create_group("/RuntimePars") +grp.attrs["PeriodicBoundariesOn"] = 1 + +# Particle group +grp1 = file.create_group("/PartType1") +ds = grp1.create_dataset("Velocities", (N_PARTICLES, 3), "f", data=vel) +ds = grp1.create_dataset("Masses", (N_PARTICLES,), "f", data=mass) +ds = grp1.create_dataset("ParticleIDs", (N_PARTICLES,), "L", data=ids) +ds = grp1.create_dataset("Coordinates", (N_PARTICLES, 3), "d", data=pos) + +file.close() diff --git a/examples/GravityTests/MWPotential2014_circularorbit/makePlots.py b/examples/GravityTests/MWPotential2014_circularorbit/makePlots.py new file mode 100755 index 0000000000000000000000000000000000000000..182c93900e9f70f51254fd099e5d242d2a574041 --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/makePlots.py @@ -0,0 +1,205 @@ +############################################################################### +# This file is part of SWIFT. +# Copyright (c) 2023 Darwin Roduit (darwin.roduit@epfl.ch) +# +# 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/>. +# +################################################################################ +import numpy as np +import h5py +import matplotlib.pyplot as plt + + +#%%Functions + + +def get_positions_and_time(N_snapshots, N_part, output_dir, boxsize): + xx = np.zeros((N_part, N_snapshots)) + yy = np.zeros((N_part, N_snapshots)) + zz = np.zeros((N_part, N_snapshots)) + time = np.zeros(N_snapshots) + pot = np.zeros((N_part, N_snapshots)) + for i in range(0, N_snapshots): + Data = h5py.File(output_dir + "/output_%04d.hdf5" % i, "r") + header = Data["Header"] + time[i] = header.attrs["Time"][0] + particles = Data["PartType1"] + positions = particles["Coordinates"] + xx[:, i] = positions[:, 0] - boxsize / 2.0 + yy[:, i] = positions[:, 1] - boxsize / 2.0 + zz[:, i] = positions[:, 2] - boxsize / 2.0 + pot[:, i] = particles["Potentials"][:] + return xx, yy, zz, time, pot + + +def plot_orbits(x, y, z, color, save_fig_name_suffix): + # Plots the orbits + fig, ax = plt.subplots(nrows=1, ncols=3, num=1, figsize=(12, 4.1)) + fig.suptitle("Orbits", fontsize=15) + ax[0].clear() + ax[1].clear() + + for i in range(0, N_part): + ax[0].plot(x[i, :], y[i, :], color[i]) + + ax[0].set_aspect("equal", "box") + ax[0].set_xlim([-35, 35]) + ax[0].set_ylim([-35, 35]) + ax[0].set_ylabel("y (kpc)") + ax[0].set_xlabel("x (kpc)") + + for i in range(0, N_part): + ax[1].plot(x[i, :], z[i, :], col[i]) + + ax[1].set_aspect("equal", "box") + ax[1].set_xlim([-35, 35]) + ax[1].set_ylim([-35, 35]) + ax[1].set_ylabel("z (kpc)") + ax[1].set_xlabel("x (kpc)") + ax[1].legend( + [ + "Particule 1, $R = 5$ kpc", + "Particule 2, $R = 5$ kpc", + "Particule 3, $R = 30$ kpc", + ] + ) + + for i in range(0, N_part): + ax[2].plot(y[i, :], z[i, :], col[i]) + + ax[2].set_aspect("equal", "box") + ax[2].set_xlim([-35, 35]) + ax[2].set_ylim([-35, 35]) + ax[2].set_ylabel("z (kpc)") + ax[2].set_xlabel("y (kpc)") + plt.tight_layout() + plt.savefig("circular_orbits" + save_fig_name_suffix + ".png", bbox_inches="tight") + plt.close() + + +def plot_deviation_from_circular_orbits(x, y, z, time, color, save_fig_name_suffix): + # Plot of the deviation from circular orbit + R_1_th = 5.0 # kpc + R_2_th = 5.0 # kpc + R_3_th = 30.0 # kpc + + fig2, ax2 = plt.subplots(nrows=1, ncols=2, num=2, figsize=(12, 4.5)) + fig2.suptitle("Deviation from circular orbit", fontsize=15) + ax2[0].clear() + ax2[1].clear() + + # Gather the x,y and z components of each particule into one array + pos_1 = np.array([x[0, :], y[0, :], z[0, :]]) + pos_2 = np.array([x[1, :], y[1, :], z[1, :]]) + pos_3 = np.array([x[2, :], y[2, :], z[2, :]]) + + # Compute the radii + r_1 = np.linalg.norm(pos_1, axis=0) + error_1 = np.abs(r_1 - R_1_th) / R_1_th * 100 + r_2 = np.linalg.norm(pos_2, axis=0) + error_2 = np.abs(r_2 - R_2_th) / R_2_th * 100 + r_3 = np.linalg.norm(pos_3, axis=0) + error_3 = np.abs(r_3 - R_3_th) / R_3_th * 100 + + ax2[0].plot(time, error_1, color[0]) + ax2[1].plot(time, error_2, color[1]) + ax2[1].plot(time, error_3, color[2]) + ax2[0].set_ylabel("Deviation (\%)") + ax2[0].set_xlabel("Time (Gyr)") + ax2[1].set_ylabel("Deviation (\%)") + ax2[1].set_xlabel("Time (Gyr)") + ax2[0].legend(["Particule 1, $R = 5$ kpc"]) + ax2[1].legend(["Particule 2, $R = 5$ kpc", "Particule 3, $R = 30$ kpc"]) + + plt.tight_layout() + plt.savefig("deviation" + save_fig_name_suffix + ".png", bbox_inches="tight") + plt.close() + return r_1, r_2, r_3 + + +def plot_deviation_from_original_data(r_1, r_2, r_3, time, color, save_fig_name_suffix): + """Make a comparison with the obtained data and ours to check nothing is broken.""" + filename = "original_radii.txt" + r_1_original, r_2_original, r_3_original = np.loadtxt(filename) + + # Plots the deviation wrt the original data + fig3, ax3 = plt.subplots(nrows=1, ncols=3, num=3, figsize=(12, 4.3)) + fig3.suptitle("Deviation from the original data", fontsize=15) + ax3[0].clear() + ax3[1].clear() + ax3[2].clear() + + error_1 = np.abs(r_1 - r_1_original) / r_1_original * 100 + error_2 = np.abs(r_2 - r_2_original) / r_2_original * 100 + error_3 = np.abs(r_3 - r_3_original) / r_3_original * 100 + + ax3[0].plot(time, error_1, col[0]) + ax3[1].plot(time, error_2, col[1]) + ax3[2].plot(time, error_3, col[2]) + ax3[0].set_ylabel("Deviation (\%)") + ax3[0].set_xlabel("Time (Gyr)") + ax3[1].set_ylabel("Deviation (\%)") + ax3[1].set_xlabel("Time (Gyr)") + ax3[2].set_ylabel("Deviation (\%)") + ax3[2].set_xlabel("Time (Gyr)") + ax3[0].legend(["Particule 1, $R = 5$ kpc"]) + ax3[1].legend(["Particule 2, $R = 5$ kpc"]) + ax3[2].legend(["Particule 3, $R = 30$ kpc"]) + plt.tight_layout() + plt.savefig( + "deviation_from_original_data" + save_fig_name_suffix + ".png", + bbox_inches="tight", + ) + plt.close() + + +#%%Plots the orbits, the deviation from the circular orbit and the deviation from the original precomputed data +# Notice that in this examples, the ouputs are set in suitable units in the parameters files. + +# General parameters +N_snapshots = 201 +N_part = 3 +boxsize = 1000.0 # kpc +col = ["b", "r", "c", "y", "k"] + +# First type of units (kpc) +output_dir = "output_1" +save_fig_name_suffix = "_simulation_kpc" +x_1, y_1, z_1, time_1, pot_1 = get_positions_and_time( + N_snapshots, N_part, output_dir, boxsize +) +plot_orbits(x_1, y_1, z_1, col, save_fig_name_suffix) +r_11, r_21, r_31 = plot_deviation_from_circular_orbits( + x_1, y_1, z_1, time_1, col, save_fig_name_suffix +) +plot_deviation_from_original_data(r_11, r_21, r_31, time_1, col, save_fig_name_suffix) + +# Second type of units (Mpc) (no need for units conversion to kpc, this is already done by swift in the snapshots) +output_dir = "output_2" +save_fig_name_suffix = "_simulation_Mpc" +x_2, y_2, z_2, time_2, pot_2 = get_positions_and_time( + N_snapshots, N_part, output_dir, boxsize +) +plot_orbits(x_2, y_2, z_2, col, save_fig_name_suffix) +r_12, r_22, r_32 = plot_deviation_from_circular_orbits( + x_2, y_2, z_2, time_2, col, save_fig_name_suffix +) +# plot_deviation_from_original_data(r_12, r_22, r_32, time_2, col, save_fig_name_suffix) #does not make sense since the original data are in kpc, not in Mpc + +#%%Saves our data to be the reference ones (precomputed) +# Uncomment only if corrections of the precomputed data must occur ! +# Original data : If some corrections occur in the potential default parameters, allows to correct +# the data +# filename = "original_radii.txt" +# np.savetxt(filename, (r_1, r_2, r_3)) diff --git a/examples/GravityTests/MWPotential2014_circularorbit/original_radii.txt b/examples/GravityTests/MWPotential2014_circularorbit/original_radii.txt new file mode 100644 index 0000000000000000000000000000000000000000..8f2e6a09e86727c1f64aa8cc0c3cb072073b387f --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/original_radii.txt @@ -0,0 +1,3 @@ +5.000000000000000000e+00 5.001166599459545559e+00 5.020579887129778207e+00 5.100189116397097600e+00 5.299839849332817820e+00 5.578710433997593476e+00 5.783030230255659987e+00 5.827598846292087131e+00 5.682565987909128147e+00 5.357447130085647657e+00 4.912578397139871988e+00 4.497394792994665380e+00 4.397021356282944105e+00 4.775431657567319910e+00 5.300932384346391579e+00 5.731042038011933570e+00 5.965091504370781728e+00 5.970489901937201971e+00 5.752272156326732500e+00 5.356957312236751534e+00 4.906750290982262008e+00 4.657121571444537089e+00 4.714457558888792477e+00 4.933301629890139317e+00 5.176061235968707486e+00 5.362349356536657119e+00 5.461122477067911873e+00 5.485087743207845534e+00 5.485339856546893600e+00 5.523077935189706800e+00 5.537107938749801228e+00 5.464804174270136095e+00 5.285227013585691580e+00 5.020384810243945672e+00 4.743263219255731578e+00 4.592936367545448206e+00 4.749715571294080618e+00 5.194699542807357240e+00 5.642069433457655769e+00 5.933223989683256150e+00 6.005807399351317244e+00 5.844582655179585196e+00 5.470114662333674715e+00 4.956958096079127962e+00 4.496399606097051027e+00 4.436328410208022710e+00 4.769513345769095025e+00 5.203928883967520846e+00 5.560443291869900051e+00 5.758575039733128342e+00 5.777233497154392161e+00 5.634900936928413984e+00 5.399734327667077949e+00 5.204348571913461008e+00 5.108383765639208818e+00 5.057693052077008034e+00 5.015035997691488667e+00 4.970823689765934361e+00 4.947419029227536846e+00 4.995207789090853723e+00 5.182270748184618192e+00 5.494607608495452489e+00 5.759878945928849525e+00 5.871757751107428369e+00 5.788023974259131066e+00 5.507664932242153810e+00 5.073360594415489366e+00 4.604596875155709590e+00 4.356042545916934294e+00 4.609781734172386791e+00 5.121533352837495556e+00 5.594559699714237055e+00 5.895503032313782477e+00 5.976985446160555604e+00 5.833832600072278218e+00 5.499788789299549840e+00 5.068750025511791435e+00 4.754624075584575671e+00 4.723949943832726817e+00 4.875445028387068724e+00 5.081882737105429371e+00 5.259925739770850761e+00 5.372637274737933843e+00 5.425242400507194418e+00 5.460118287161413342e+00 5.536220897651173090e+00 5.595059701109168948e+00 5.560241134615987235e+00 5.401506769780430872e+00 5.130567420218448582e+00 4.809371855699721365e+00 4.569652033714569406e+00 4.612229578661045437e+00 5.019851467780982901e+00 5.505814358798329700e+00 5.865468786224973918e+00 6.016174230579117932e+00 5.931623634176135695e+00 5.622177861123924814e+00 5.141870044592171674e+00 4.638547645705167533e+00 4.428307819955180413e+00 4.654500832829794099e+00 5.055593990419128048e+00 5.427492129376623176e+00 5.668492207620567491e+00 5.744061565586670248e+00 5.662423276901259683e+00 5.476861486048790084e+00 5.299864737421422056e+00 5.202470997364078364e+00 5.134308457492519295e+00 5.058899583759844276e+00 4.970512274614167225e+00 4.896917303931519250e+00 4.898286265130379569e+00 5.054476498797095374e+00 5.384729626924356083e+00 5.707038692608169761e+00 5.888268525975762557e+00 5.871970152607735471e+00 5.647189263247745394e+00 5.241550304616231060e+00 4.747405841925699477e+00 4.378927045831982667e+00 4.474186584819125123e+00 4.940610149752312275e+00 5.439852850044650800e+00 5.799438267533770563e+00 5.953178490303065118e+00 5.885318802304978725e+00 5.618141372733455263e+00 5.224255162563557597e+00 4.875392411458151720e+00 4.768956166135601471e+00 4.850718648821461976e+00 5.009868149351704325e+00 5.165493314050841889e+00 5.278258318859305298e+00 5.348100423409463566e+00 5.410955922599304913e+00 5.522902505871598144e+00 5.630663735706390227e+00 5.642294107886619337e+00 5.517118052537815842e+00 5.256074072268789088e+00 4.908246395925647043e+00 4.590870106450615218e+00 4.510269017110958956e+00 4.841676564404478356e+00 5.347528554683028013e+00 5.768526167155091144e+00 5.995637396520894669e+00 5.989832219840335092e+00 5.751871797307715362e+00 5.320044178608156926e+00 4.806704151149128634e+00 4.474934801273977136e+00 4.576629640621881379e+00 4.921538717217767811e+00 5.290768217179389055e+00 5.562191457746648915e+00 5.686534351951245014e+00 5.661795715468069368e+00 5.528388603164072101e+00 5.380407491001560238e+00 5.294789296099782661e+00 5.223565412511718797e+00 5.127422326009080322e+00 5.001198016563644266e+00 4.875299081804920220e+00 4.818639828265833813e+00 4.925709984592579360e+00 5.252798132458202041e+00 5.624850406729384744e+00 5.875379086292169539e+00 5.930742494670786513e+00 5.769927451901470938e+00 5.407923458724494203e+00 4.913709661908932880e+00 4.459744560941122593e+00 4.382798389554843688e+00 4.768537003581841738e+00 5.272907504565521819e+00 5.680248961594461576e+00 5.900602226243958803e+00 5.906378074467261818e+00 5.709143330590744192e+00 5.364799048397715708e+00 5.008887298284076195e+00 4.843398709646222322e+00 4.859874698091035228e+00 4.965584565534517658e+00 5.086695693063617796e+00 5.185431020385737000e+00 5.259341369321015591e+00 5.340599814904859244e+00 5.482731832086808588e+00 5.640753798791949158e+00 5.705643537489685002e+00 5.624938857189243357e+00 +5.000000000000000000e+00 4.999888359284162753e+00 5.000641684297063350e+00 5.000994447130620024e+00 5.001675693602675388e+00 5.001624794631814197e+00 5.001558065482559989e+00 5.000888414115478575e+00 5.000403424123103235e+00 4.999857363290737489e+00 4.999815117935768605e+00 5.000030078010763823e+00 5.000590428170696988e+00 5.001185604320395939e+00 5.001578829707598928e+00 5.001700786476619420e+00 5.001329335974549650e+00 5.000875735587428927e+00 5.000172435947884964e+00 4.999938404239019363e+00 4.999745387639388028e+00 5.000283699222444866e+00 5.000643388144556667e+00 5.001470334121859374e+00 5.001565799279720004e+00 5.001856736927319247e+00 5.001184047933382892e+00 5.000858585776699172e+00 5.000045524839611844e+00 4.999938201977301766e+00 4.999789233051740389e+00 5.000353429170261599e+00 5.000792119848084027e+00 5.001452348355899247e+00 5.001632422736415506e+00 5.001608535971686642e+00 5.001127158676252904e+00 5.000553905860360970e+00 5.000029329624915420e+00 4.999762162118537745e+00 4.999946677627926661e+00 5.000335136347326070e+00 5.001035783553454017e+00 5.001422330001844330e+00 5.001779305392438424e+00 5.001442290734069651e+00 5.001159754497918719e+00 5.000348863737741567e+00 5.000122347527026889e+00 4.999700132208437786e+00 5.000215353627103454e+00 5.000418901608895794e+00 5.001274128363529137e+00 5.001475223091539135e+00 5.001811379870214402e+00 5.001359505875552891e+00 5.000995546297854588e+00 5.000250370590015869e+00 4.999942233798068791e+00 4.999755726773103426e+00 5.000113010258659507e+00 5.000609015959638803e+00 5.001219673802144605e+00 5.001614354202936319e+00 5.001636066331053954e+00 5.001364334692166125e+00 5.000733606126373409e+00 5.000258441419644484e+00 4.999780701257404480e+00 4.999929042417343261e+00 5.000122602854585097e+00 5.000898089680752889e+00 5.001250786342685828e+00 5.001830340890685456e+00 5.001536553836007926e+00 5.001450960015414893e+00 5.000564462723501968e+00 5.000273841477640424e+00 4.999729946197625452e+00 5.000019708256767004e+00 5.000238197444725863e+00 5.000993609020794928e+00 5.001370698158823735e+00 5.001741598978786563e+00 5.001520071329959016e+00 5.001141309344722785e+00 5.000492961311684859e+00 5.000007097885235119e+00 4.999789471734320756e+00 4.999926978056772775e+00 5.000452501234716074e+00 5.000982517444656494e+00 5.001568938123119423e+00 5.001630494678826366e+00 5.001581127018102535e+00 5.000918624733289519e+00 5.000520523906761383e+00 4.999849635612912735e+00 4.999963116585512957e+00 4.999945908278064621e+00 5.000768878036669918e+00 5.001060851264059082e+00 5.001749243021976099e+00 5.001594604383085674e+00 5.001528201965431464e+00 5.000787308820795474e+00 5.000362783584799686e+00 4.999809889944675234e+00 4.999870251224309392e+00 5.000081273895164991e+00 5.000706925098050704e+00 5.001234607158520085e+00 5.001628277065568895e+00 5.001642902217743547e+00 5.001272369050412614e+00 5.000752316183323387e+00 5.000115144777662834e+00 4.999877063130510280e+00 4.999785278748652750e+00 5.000314975961055453e+00 5.000732923535549546e+00 5.001486668371692978e+00 5.001580864064437471e+00 5.001767526649143925e+00 5.001102524724919896e+00 5.000815405699214899e+00 4.999977098680566279e+00 4.999965430706085456e+00 4.999820211656420987e+00 5.000468706648866224e+00 5.000862674405269992e+00 5.001540997553139967e+00 5.001619721607925229e+00 5.001589058085179396e+00 5.001022878428104335e+00 5.000494347648318794e+00 4.999950437180141094e+00 4.999782011985283603e+00 4.999968771211793950e+00 5.000436397880152484e+00 5.001086192651295725e+00 5.001485312783913173e+00 5.001737933662662172e+00 5.001395316031677751e+00 5.001032550774532126e+00 5.000272753974341455e+00 5.000032632149568279e+00 4.999711660403833413e+00 5.000227368341616518e+00 5.000504073344837153e+00 5.001396382624039738e+00 5.001505342442887248e+00 5.001831185827574799e+00 5.001279357006295001e+00 5.000939059512067075e+00 5.000156498511878489e+00 4.999939419656115547e+00 4.999760777523306388e+00 5.000213242684194981e+00 5.000677422077679068e+00 5.001316266444979952e+00 5.001612758519516255e+00 5.001623006669117188e+00 5.001255205432448392e+00 5.000655856805036770e+00 5.000151079135440213e+00 4.999768706765014059e+00 4.999923501590752828e+00 5.000205665288950385e+00 5.000941154591838078e+00 5.001315826587424240e+00 5.001794596942993110e+00 5.001491317612160259e+00 5.001316595758519057e+00 5.000470387348008749e+00 5.000254580857065534e+00 4.999712075477516393e+00 5.000102976947414568e+00 5.000306391555606744e+00 5.001111012410059509e+00 5.001407461804255661e+00 5.001773556648744012e+00 5.001449382725613013e+00 5.001085296333972252e+00 5.000387772118549456e+00 4.999982482345314061e+00 4.999767868097886314e+00 5.000003634425123522e+00 5.000507888369177145e+00 5.001079803140623170e+00 5.001579362536846318e+00 5.001633342805087423e+00 5.001483416394396642e+00 5.000841451375783286e+00 5.000401977832478195e+00 4.999818075316714072e+00 4.999936405692012364e+00 5.000014185319412441e+00 5.000809312551328212e+00 5.001135437752519231e+00 5.001827481510867202e+00 +3.000000000000000000e+01 2.999779458340697147e+01 2.999913108927551875e+01 2.999822513999371054e+01 2.999869101241313629e+01 2.999908157608862069e+01 2.999867255062164162e+01 3.000035432115671696e+01 2.999906436031323764e+01 3.000130726868057707e+01 2.999985124246509116e+01 3.000120011581939750e+01 3.000101355345713472e+01 3.000145740135981853e+01 3.000252851377030439e+01 3.000205476548366335e+01 3.000437014149219905e+01 3.000296575532702192e+01 3.000506656450175669e+01 3.000416345420759257e+01 3.000531256555992954e+01 3.000561911665013071e+01 3.000580235611645818e+01 3.000730436832045811e+01 3.000650795996002174e+01 3.000919171647245420e+01 3.000740253616102748e+01 3.000908722559980291e+01 3.000846124470268705e+01 3.000913539394623442e+01 3.000966175811804604e+01 3.000931550053518393e+01 3.001098523842369659e+01 3.000961070226721716e+01 3.001169411489794570e+01 3.001000798292042404e+01 3.001105493839978067e+01 3.001049866425505996e+01 3.001050661892987037e+01 3.001107862320722219e+01 3.001004719415356803e+01 3.001174860804977840e+01 3.000967995094207907e+01 3.001106890886319079e+01 3.000941297065721614e+01 3.000977242153797064e+01 3.000925845687551785e+01 3.000859589580250741e+01 3.000923251594434049e+01 3.000755748209316565e+01 3.000935520671618661e+01 3.000667892784183266e+01 3.000748164415950114e+01 3.000598473874775962e+01 3.000580528773574329e+01 3.000550129646441633e+01 3.000435340448430566e+01 3.000525642811016169e+01 3.000315424345881254e+01 3.000455578094075548e+01 3.000223634049040911e+01 3.000270501649730548e+01 3.000162796403852283e+01 3.000117720180458747e+01 3.000135582863928008e+01 2.999999786615864394e+01 3.000144392723029085e+01 2.999919002924884381e+01 3.000046804618323293e+01 2.999877353580113848e+01 2.999916917874825728e+01 2.999876482042754233e+01 2.999828485654238719e+01 2.999917646486512979e+01 2.999782544062704659e+01 3.000001627767003498e+01 2.999779633743827034e+01 2.999911844063658251e+01 2.999819830488571526e+01 2.999865030911109187e+01 2.999902751315119787e+01 2.999860581120435299e+01 3.000027557069040540e+01 2.999897424318053041e+01 3.000120647142659891e+01 2.999974052798087598e+01 3.000108036609949735e+01 3.000088565466134227e+01 3.000132237211635555e+01 3.000238750303375213e+01 3.000190931347448497e+01 3.000422185867911296e+01 3.000281580647389745e+01 3.000491608599069338e+01 3.000401361861365146e+01 3.000516478912791385e+01 3.000547462841453950e+01 3.000566249461711266e+01 3.000717035022314505e+01 3.000638099269628611e+01 3.000907299236889969e+01 3.000729313447100211e+01 3.000898818458806971e+01 3.000837351407666986e+01 3.000905988756514375e+01 3.000959933442729266e+01 3.000926695882656148e+01 3.001095125287651655e+01 3.000959180150701044e+01 3.001169059343623147e+01 3.001002003559164066e+01 3.001108278198048041e+01 3.001054248304392047e+01 3.001056653531087548e+01 3.001115458351709364e+01 3.001013901847990439e+01 3.001185589737429282e+01 3.000980246935861828e+01 3.001120595492282561e+01 3.000956393020903334e+01 3.000993653574142783e+01 3.000943488286652894e+01 3.000878378812352665e+01 3.000943101649924571e+01 3.000776555887522434e+01 3.000957175115255993e+01 3.000690282240585205e+01 3.000771163940339292e+01 3.000621963347996513e+01 3.000604381912500784e+01 3.000574221323055824e+01 3.000459546464529481e+01 3.000549825809542170e+01 3.000339446894059847e+01 3.000479306306118232e+01 3.000246913202072463e+01 3.000293187427373454e+01 3.000184761572847592e+01 3.000138830507564691e+01 3.000155700520485880e+01 3.000018799611447662e+01 3.000162189224056064e+01 2.999935478803817901e+01 3.000061862764304976e+01 2.999890904656334101e+01 2.999928887416733758e+01 2.999886807861885174e+01 2.999837114075570099e+01 2.999924530142823542e+01 2.999787658246655297e+01 3.000004954612104058e+01 2.999781168046782653e+01 2.999911595907233774e+01 2.999817825072545574e+01 2.999861297335620591e+01 2.999897334516337466e+01 2.999853524782260195e+01 3.000018922740121141e+01 2.999887294236322077e+01 3.000109111179866161e+01 2.999961196637686811e+01 3.000093950505366180e+01 3.000073363009094862e+01 3.000116029442814991e+01 3.000221660273533431e+01 3.000173085754470392e+01 3.000403703825239532e+01 3.000262587834090766e+01 3.000472258970311401e+01 3.000381811956868461e+01 3.000496875334294344e+01 3.000527961996388271e+01 3.000546996417365975e+01 3.000698169725789910e+01 3.000619758041846197e+01 3.000889621429398701e+01 3.000712430259497765e+01 3.000882861345186825e+01 3.000822447317817065e+01 3.000892257912635230e+01 3.000947488798068008e+01 3.000915634131698795e+01 3.001085531632530490e+01 3.000951134895363381e+01 3.001162660663587900e+01 3.000997326126640274e+01 3.001105384355043526e+01 3.001053154124043232e+01 3.001057371847517885e+01 3.001118006937515403e+01 3.001018279679330192e+01 3.001191804078371206e+01 3.000988258949455556e+01 3.001130365826893609e+01 3.000967875120074169e+01 3.001006775551087458e+01 3.000958182273279817e+01 3.000894558530716338e+01 3.000960665058126864e+01 3.000795400716476991e+01 3.000977178013717150e+01 diff --git a/examples/GravityTests/MWPotential2014_circularorbit/params_unit_1.yml b/examples/GravityTests/MWPotential2014_circularorbit/params_unit_1.yml new file mode 100755 index 0000000000000000000000000000000000000000..9e7ff5e2bfebb0e19e7231d72bbf8efe3162278c --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/params_unit_1.yml @@ -0,0 +1,53 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e+43 # 10^10 Solar masses + UnitLength_in_cgs: 3.086e+21 # kpc + UnitVelocity_in_cgs: 1e5 # 1 km / s in cm/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 2.0 # The end time of the simulation (in internal units). + dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units). + dt_max: 1e0 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: output_1/output # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) + UnitMass_in_cgs: 1.98848e+43 + UnitLength_in_cgs: 3.086e+21 + UnitVelocity_in_cgs: 1e5 + UnitCurrent_in_cgs: 1 + UnitTemp_in_cgs: 1 + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1.0 # Time between statistics output + +# Parameters related to the initial conditions +InitialConditions: + file_name: circular_orbits_MW.hdf5 # The file to read + periodic: 1 + +# NFW_MN_PSC potential parameters +MWPotential2014Potential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] #Centre of the potential with respect to centre of the box + timestep_mult: 0.005 # Dimensionless pre-factor for the time-step condition + epsilon: 0.001 # Softening size (internal units) + #The following parameters are the default ones. + # concentration: 9.823403437774843 + # M_200_Msun: 147.41031542774076e10 + # H: 127.78254614201471e-2 + # Mdisk_Msun: 6.8e10 + # Rdisk_kpc: 3.0 + # Zdisk_kpc: 0.280 + # amplitude: 1.0 + # r_1_kpc: 1.0 + # alpha: 1.8 + # r_c_kpc: 1.9 + # potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] diff --git a/examples/GravityTests/MWPotential2014_circularorbit/params_unit_2.yml b/examples/GravityTests/MWPotential2014_circularorbit/params_unit_2.yml new file mode 100755 index 0000000000000000000000000000000000000000..8d6b232b11f301d20f1a41c4872b056cfe408c4f --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/params_unit_2.yml @@ -0,0 +1,52 @@ +# Define the system of units to use internally. +InternalUnitSystem: + UnitMass_in_cgs: 1.98848e+43 # 10^10 Solar masses + UnitLength_in_cgs: 3.086e+24 # Mpc + UnitVelocity_in_cgs: 1e5 # 1 km / s in cm/s + UnitCurrent_in_cgs: 1 # Amperes + UnitTemp_in_cgs: 1 # Kelvin + +# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) +TimeIntegration: + time_begin: 0. # The starting time of the simulation (in internal units). + time_end: 2.0e-3 # The end time of the simulation (in internal units). + dt_min: 1e-13 # The minimal time-step size of the simulation (in internal units). + dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units). + +# Parameters governing the snapshots +Snapshots: + basename: output_2/output # Common part of the name of output files + time_first: 0. # Time of the first output (in internal units) + delta_time: 1e-5 # Time difference between consecutive outputs (in internal units) + UnitMass_in_cgs: 1.98848e+43 + UnitLength_in_cgs: 3.086e+21 + UnitVelocity_in_cgs: 1e5 + UnitCurrent_in_cgs: 1 + UnitTemp_in_cgs: 1 +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 1.0 # Time between statistics output + +# Parameters related to the initial conditions +InitialConditions: + file_name: circular_orbits_MW.hdf5 # The file to read + periodic: 1 + +# NFW_MN_PSC potential parameters +MWPotential2014Potential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] #Centre of the potential with respect to centre of the box + timestep_mult: 0.0005 # Dimensionless pre-factor for the time-step condition + epsilon: 0.001e-3 # Softening size (internal units) + #The following parameters are the default ones. + # concentration: 9.823403437774843 + # M_200_Msun: 147.41031542774076e10 + # H: 127.78254614201471e-2 + # Mdisk_Msun: 6.8e10 + # Rdisk_kpc: 3.0 + # Zdisk_kpc: 0.280 + # amplitude: 1.0 + # r_1_kpc: 1.0 + # alpha: 1.8 + # r_c_kpc: 1.9 + # potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] diff --git a/examples/GravityTests/MWPotential2014_circularorbit/run.sh b/examples/GravityTests/MWPotential2014_circularorbit/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..a2821996a9cfde87bfdaaff36df9bd576a083241 --- /dev/null +++ b/examples/GravityTests/MWPotential2014_circularorbit/run.sh @@ -0,0 +1,86 @@ +#!/bin/bash + +#Creates a directory for the outputs +DIR=output_1 #First test of units conversion +if [ -d "$DIR" ]; +then + echo "$DIR directory exists. Its content will be removed." + rm $DIR/output_* +else + echo "$DIR directory does not exists. It will be created." + mkdir $DIR +fi + +DIR=output_2 #Second test of units conversion +if [ -d "$DIR" ]; +then + echo "$DIR directory exists. Its content will be removed." + rm $DIR/output_* +else + echo "$DIR directory does not exists. It will be created." + mkdir $DIR +fi + +#Clears the previous figures +echo "Clearing existing figures." +if [ -f "circular_orbits_simulation_kpc.png" ]; +then + rm circular_orbits_simulation_kpc.png +fi + +if [ -f "circular_orbits_simulation_Mpc.png" ]; +then + rm circular_orbits_simulation_Mpc.png +fi + +if [ -f "deviation_simulation_kpc.png" ]; +then + rm deviation_simulation_kpc.png +fi +if [ -f "deviation_simulation_Mpc.png" ]; +then + rm deviation_simulation_Mpc.png +fi + +if [ -f "deviation_from_original_data_simulation_kpc.png" ]; +then + rm deviation_from_original_data_simulation_kpc.png +fi + +if [ -f "deviation_from_original_data_simulation_Mpc.png" ]; +then + rm deviation_from_original_data_simulation_Mpc.png +fi + +#Clears the IC file +if [ -f "circular_orbits_MW.hdf5" ]; +then + rm circular_orbits_MW.hdf5 +fi + + +#Generates the initial conditions +echo "Generate initial conditions for circular orbits" +if command -v python3 &>/dev/null; then + python3 makeIC.py +else + python3 makeIC.py +fi + +#Runs the simulation +# self gravity G, external potential g, hydro s, threads t and high verbosity v +echo "Starting the simulation with the first type of units (kpc)... Look at output_1.log for the simulation details." +../../../swift --external-gravity --threads=8 params_unit_1.yml 2>&1 > output_1.log +echo "Simulation ended." + +echo "Starting the simulation with the second type of units (Mpc)... Look at output_2.log for the simulation details." +../../../swift --external-gravity --threads=8 params_unit_2.yml 2>&1 > output_2.log +echo "Simulation ended." + +#Saves the plots +echo "Save plots of the circular orbits and of the errors" +if command -v python3 &>/dev/null; then + python3 makePlots.py +else + python3 makePlots.py +fi diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index b10e6d82e79947a4b568ce2bbe68229edcf52700..0e5cac61fcce905f60ac8cbc7ed986081dfa4ed0 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -420,6 +420,23 @@ SineWavePotential: timestep_limit: 1. # Time-step dimensionless pre-factor. growth_time: 0. # (Optional) Time for the potential to grow to its final size. +# MWPotential2014 potential +MWPotential2014Potential: + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + position: [0.,0.,0.] # Location of centre of potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) + timestep_mult: 0.005 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration + epsilon: 0.001 # Softening size (internal units) + concentration: 9.823403437774843 # concentration of the Halo + M_200_Msun: 147.41031542774076e10 # M200 of the galaxy disk (in M_sun) + H: 1.2778254614201471 # Hubble constant in units of km/s/Mpc + Mdisk_Msun: 6.8e10 # Mass of the disk (in M_sun) + Rdisk_kpc: 3.0 # Effective radius of the disk (in kpc) + Zdisk_kpc: 0.280 # Scale-height of the disk (in kpc) + amplitude: 1.0 # Amplitude of the bulge + r_1_kpc: 1.0 # Reference radius for amplitude of the bulge (in kpc) + alpha: 1.8 # Exponent of the power law of the bulge + r_c_kpc: 1.9 # Cut-off radius of the bulge (in kpc) + potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] # Coefficients that adjust the strength of the halo (1st component), the disk (2nd component) and the bulge (3rd component) # Parameters related to entropy floors ---------------------------------------------- diff --git a/src/potential.h b/src/potential.h index 9011eee11123f0d04883a69119ac6ee86aeb354d..311c3b3b4c824424939cfa218e96ea8e7ae81581 100644 --- a/src/potential.h +++ b/src/potential.h @@ -42,6 +42,8 @@ #include "./potential/nfw/potential.h" #elif defined(EXTERNAL_POTENTIAL_NFW_MN) #include "./potential/nfw_mn/potential.h" +#elif defined(EXTERNAL_POTENTIAL_MWPotential2014) +#include "./potential/MWPotential2014/potential.h" #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH) #include "./potential/disc_patch/potential.h" #elif defined(EXTERNAL_POTENTIAL_SINE_WAVE) diff --git a/src/potential/MWPotential2014/potential.h b/src/potential/MWPotential2014/potential.h new file mode 100644 index 0000000000000000000000000000000000000000..09b80014dde76b0d77b542582465080e4289fe41 --- /dev/null +++ b/src/potential/MWPotential2014/potential.h @@ -0,0 +1,497 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2023 Darwin Roduit + * + * 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_POTENTIAL_MWPotential2014_H +#define SWIFT_POTENTIAL_MWPotential2014_H + +/* Config parameters. */ +#include <config.h> + +/* Some standard headers. */ +#include <float.h> +#include <math.h> + +/* Local includes. */ +#include "error.h" +#include "gravity.h" +#include "parser.h" +#include "part.h" +#include "physical_constants.h" +#include "space.h" +#include "units.h" + +#ifdef HAVE_LIBGSL +#include <gsl/gsl_sf_gamma.h> +#endif + +/** + * @brief External Potential Properties - MWPotential2014 composed by + * NFW + Miyamoto-Nagai + Power Spherical cut-off potentials + * + * halo --> rho_NFW(r) = rho_0 / ( (r/R_s)*(1+r/R_s)^2 ) + * disk --> Phi_MN(R,z) = -G * Mdisk / (R^2 + (Rdisk + + * (z^2+Zdisk^2)^1/2)^2)^(1/2) bulge --> rho_PSC(r) = + * amplitude*(r_1/r)^alpha*exp(-(r/r_c)^2) + * + * We however parametrise this in terms of c and virial_mass, Mdisk, Rdisk + * and Zdisk. Also, each potential is given a contribution amplitude such that + * the resulting potential is: + * Phi_tot = f_1 * Phi_NFW + f_2 * Phi_MN + f_3 * Phi_PSC, + * with f_1, f_2 and f_3 contained in the array f. + * + * This potential is inspired by the following article: + * galpy: A Python Library for Galactic Dynamics, Jo Bovy (2015), + * Astrophys. J. Supp., 216, 29 (arXiv/1412.3451). + */ +struct external_potential { + + /*! Position of the centre of potential */ + double x[3]; + + /*! The scale radius of the NFW potential */ + double r_s; + + /*! The pre-factor \f$ 4 \pi G \rho_0 \r_s^3 \f$ */ + double pre_factor; + + /*! Hubble parameter */ + double H; + + /*! The concentration parameter */ + double c_200; + + /*! The virial mass */ + double M_200; + + /*! Disk Size */ + double Rdisk; + + /*! Disk height */ + double Zdisk; + + /*! Disk Mass */ + double Mdisk; + + /*! Amplitude for the PSC potential */ + double amplitude; + + /*! Reference radius for amplitude */ + double r_1; + + /*! Inner power */ + double alpha; + + /*! Cut-off radius */ + double r_c; + + /*! Contribution of each potential : f[0]*NFW + f[1]*MN + f[2]*PSP */ + double f[3]; + + /*! Prefactor \f$ 2 \pi amplitude r_1^\alpha r_c^(3-\alpha) \f$ */ + double prefactor_psc_1; + + /*! Prefactor \f$ 2 \pi amplitude r_1^\alpha r_c^(2-\alpha) \f$ */ + double prefactor_psc_2; + + /*! Gamma function evaluation \f$ \Gamma((3-\alpha)/2 \f$ */ + double gamma_psc; + + /*! Time-step condition pre_factor, this factor is used to multiply times the + * orbital time, so in the case of 0.01 we take 1% of the orbital time as + * the time integration steps */ + double timestep_mult; + + /*! Minimum time step based on the orbital time at the softening times + * the timestep_mult */ + double mintime; + + /*! Common log term \f$ \ln(1+c_{200}) - \frac{c_{200}}{1 + c_{200}} \f$ */ + double log_c200_term; + + /*! Softening length */ + double eps; +}; + +/** + * @brief Computes the time-step due to the acceleration from the NFW + MN + PSC + * potential as a fraction (timestep_mult) of the circular orbital time of that + * particle. + * + * @param time The current time. + * @param potential The #external_potential used in the run. + * @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) { + +#ifdef HAVE_LIBGSL + + const float dx = g->x[0] - potential->x[0]; + const float dy = g->x[1] - potential->x[1]; + const float dz = g->x[2] - potential->x[2]; + + const float R2 = dx * dx + dy * dy; + const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps); + + /* Vcirc for NFW */ + const float M_NFW = potential->pre_factor * (logf(1.0f + r / potential->r_s) - + r / (r + potential->r_s)); + const float Vcirc_NFW = sqrtf((phys_const->const_newton_G * M_NFW) / r); + + /* Now for MN */ + const float R = sqrtf(R2); + const float sqrt_term = sqrtf(dz * dz + potential->Zdisk * potential->Zdisk); + const float MN_denominator = + powf(R2 + powf(potential->Rdisk + sqrt_term, 2.0f), 1.5f); + const float dPhi_dR_MN = potential->Mdisk * R / MN_denominator; + const float dPhi_dz_MN = potential->Mdisk * dz * + (potential->Rdisk + sqrt_term) / + (sqrt_term * MN_denominator); + const float Vcirc_MN = sqrtf(phys_const->const_newton_G * R * dPhi_dR_MN + + phys_const->const_newton_G * dz * dPhi_dz_MN); + + /* Now for PSC */ + const float r2 = r * r; + const float M_psc = + potential->prefactor_psc_1 * + (potential->gamma_psc - + gsl_sf_gamma_inc(1.5f - 0.5f * potential->alpha, + r2 / (potential->r_c * potential->r_c))); + const float Vcirc_PSC = sqrtf(phys_const->const_newton_G * M_psc / r); + + /* Total circular velocity */ + const float Vcirc = sqrtf(potential->f[0] * Vcirc_NFW * Vcirc_NFW + + potential->f[1] * Vcirc_MN * Vcirc_MN + + potential->f[2] * Vcirc_PSC * Vcirc_PSC); + + const float period = 2.0f * M_PI * r / Vcirc; + + /* Time-step as a fraction of the circular period */ + const float time_step = potential->timestep_mult * period; + + return max(time_step, potential->mintime); + +#else + error("Code not compiled with GSL. Can't compute MWPotential2014."); + return 0.0; +#endif +} + +/** + * @brief Computes the gravitational acceleration from an NFW Halo potential + + * MN disk + PSC bulge. + * + * Note that the accelerations are multiplied by Newton's G constant + * later on. + * + * a_x = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * x + * - dphi_PSC/dr*x/r + * a_y = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * y + * - dphi_PSC/dr*y/r + * a_z = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * z + * - dphi_PSC/dr*z/r + * + * + * @param time The current time. + * @param potential The #external_potential used in the run. + * @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) { + +#ifdef HAVE_LIBGSL + + const float dx = g->x[0] - potential->x[0]; + const float dy = g->x[1] - potential->x[1]; + const float dz = g->x[2] - potential->x[2]; + + /* First for the NFW part */ + const float R2 = dx * dx + dy * dy; + const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps); + + const float r_inv = 1.0f / r; + const float M_NFW = potential->pre_factor * (logf(1.0f + r / potential->r_s) - + r / (r + potential->r_s)); + const float dpot_dr_NFW = M_NFW * r_inv * r_inv; + const float pot_nfw = + -potential->pre_factor * logf(1.0f + r / potential->r_s) * r_inv; + g->a_grav[0] -= potential->f[0] * dpot_dr_NFW * dx * r_inv; + g->a_grav[1] -= potential->f[0] * dpot_dr_NFW * dy * r_inv; + g->a_grav[2] -= potential->f[0] * dpot_dr_NFW * dz * r_inv; + gravity_add_comoving_potential(g, potential->f[0] * pot_nfw); + + /* Now the the MN disk */ + const float f1 = sqrtf(potential->Zdisk * potential->Zdisk + dz * dz); + const float f2 = potential->Rdisk + f1; + const float f3 = powf(R2 + f2 * f2, -1.5f); + const float mn_term = potential->Rdisk + sqrtf(potential->Zdisk + dz * dz); + const float pot_mn = -potential->Mdisk / sqrtf(R2 + mn_term * mn_term); + + g->a_grav[0] -= potential->f[1] * potential->Mdisk * f3 * dx; + g->a_grav[1] -= potential->f[1] * potential->Mdisk * f3 * dy; + g->a_grav[2] -= potential->f[1] * potential->Mdisk * f3 * (f2 / f1) * dz; + gravity_add_comoving_potential(g, potential->f[1] * pot_mn); + + /* Now the the PSC bulge */ + const float r2 = r * r; + const float M_psc = + potential->prefactor_psc_1 * + (potential->gamma_psc - + gsl_sf_gamma_inc(1.5f - 0.5f * potential->alpha, + r2 / (potential->r_c * potential->r_c))); + const float dpot_dr = M_psc / r2; + const float pot_psc = + -M_psc / r - potential->prefactor_psc_2 * + gsl_sf_gamma_inc(1.0f - 0.5f * potential->alpha, + r2 / (potential->r_c * potential->r_c)); + + g->a_grav[0] -= potential->f[2] * dpot_dr * dx * r_inv; + g->a_grav[1] -= potential->f[2] * dpot_dr * dy * r_inv; + g->a_grav[2] -= potential->f[2] * dpot_dr * dz * r_inv; + gravity_add_comoving_potential(g, potential->f[2] * pot_psc); + +#else + error("Code not compiled with GSL. Can't compute MWPotential2014."); +#endif +} + +/** + * @brief Computes the gravitational potential energy of a particle in an + * NFW potential + MN potential. + * + * phi = f[0] * (-4 * pi * G * rho_0 * r_s^3 * ln(1+r/r_s)) - f[1] * (G * Mdisk + * / sqrt(R^2 + (Rdisk + sqrt(z^2 + Zdisk^2))^2)) + f[2] * [- G / r * (2 * pi * + * amplitude * r_1^alpha r_c^(3 - alpha) * gamma_inf((3 - alpha)/2, r^2 / r_c^2) + * ) - 2 * pi * G * amplitude * r_1^alpha * r_c^(2 - alpha) + * Gamma_sup((2-alpha)/2, r^2 / r_c^2 ) ] + * + * @param time The current time (unused here). + * @param potential The #external_potential used in the run. + * @param phys_const Physical constants in internal units. + * @param g 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* g) { + +#ifdef HAVE_LIBGSL + + const float dx = g->x[0] - potential->x[0]; + const float dy = g->x[1] - potential->x[1]; + const float dz = g->x[2] - potential->x[2]; + + /* First for the NFW profile */ + const float R2 = dx * dx + dy * dy; + const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps); + const float r_inv = 1.0f / r; + const float pot_nfw = + -potential->pre_factor * logf(1.0f + r / potential->r_s) * r_inv; + + /* Now for the MN disk */ + const float sqrt_term = sqrtf(dz * dz + potential->Zdisk * potential->Zdisk); + const float MN_denominator = + sqrtf(R2 + powf(potential->Rdisk + sqrt_term, 2.0f)); + const float mn_pot = -potential->Mdisk / MN_denominator; + + /* Now for PSC bulge */ + const float r2 = r * r; + const float M_psc = + potential->prefactor_psc_1 * + (potential->gamma_psc - + gsl_sf_gamma_inc(1.5f - 0.5f * potential->alpha, + r2 / (potential->r_c * potential->r_c))); + const float psc_pot = + -M_psc / r - potential->prefactor_psc_2 * + gsl_sf_gamma_inc(1.0f - 0.5f * potential->alpha, + r2 / (potential->r_c * potential->r_c)); + + return phys_const->const_newton_G * + (potential->f[0] * pot_nfw + potential->f[1] * mn_pot + + potential->f[2] * psc_pot); + +#else + error("Code not compiled with GSL. Can't compute MWPotential2014."); + return 0.0; +#endif +} + +/** + * @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( + struct swift_params* parameter_file, const struct phys_const* phys_const, + const struct unit_system* us, const struct space* s, + struct external_potential* potential) { + +#ifdef HAVE_LIBGSL + + /* Read in the position of the centre of potential */ + parser_get_param_double_array( + parameter_file, "MWPotential2014Potential:position", 3, potential->x); + + /* Is the position absolute or relative to the centre of the box? */ + const int useabspos = parser_get_param_int( + parameter_file, "MWPotential2014Potential:useabspos"); + + /* Define the default value in the above system of units*/ + const double c_200_default = 9.823403437774843; + const double M_200_Msun_default = 147.41031542774076e10; /* M_sun */ + const double H_default = 127.78254614201471e-2; /* no unit */ + const double Mdisk_Msun_default = 6.8e10; /* M_sun */ + const double Rdisk_kpc_default = 3.0; /* kpc */ + const double Zdisk_kpc_default = 0.280; /* kpc */ + const double amplitude_default = 1.0; /* no unit */ + const double r_1_kpc_default = 1.0; /* kpc */ + const double alpha_default = 1.8; /* no unit */ + const double r_c_kpc_default = 1.9; /* kpc */ + potential->f[0] = 0.4367419745056084; /* no unit */ + potential->f[1] = 1.002641971008805; /* no unit */ + potential->f[2] = 0.022264787598364262; /* no unit */ + + if (!useabspos) { + potential->x[0] += s->dim[0] / 2.; + potential->x[1] += s->dim[1] / 2.; + potential->x[2] += s->dim[2] / 2.; + } + + /* Read the other parameters of the model */ + potential->timestep_mult = parser_get_param_double( + parameter_file, "MWPotential2014Potential:timestep_mult"); + + /* Bug fix : Read the softening length from the params file */ + potential->eps = parser_get_param_double(parameter_file, + "MWPotential2014Potential:epsilon"); + + potential->c_200 = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:concentration", c_200_default); + potential->M_200 = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:M_200_Msun", + M_200_Msun_default); + potential->H = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:H", H_default); + potential->Mdisk = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:Mdisk_kpc", Mdisk_Msun_default); + potential->Rdisk = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:Rdisk_kpc", Rdisk_kpc_default); + potential->Zdisk = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:Zdisk_kpc", Zdisk_kpc_default); + potential->amplitude = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:amplitude", amplitude_default); + potential->r_1 = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:r_1_kpc", r_1_kpc_default); + potential->alpha = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:alpha", alpha_default); + potential->r_c = parser_get_opt_param_double( + parameter_file, "MWPotential2014Potential:r_c_kpc", r_c_kpc_default); + parser_get_opt_param_double_array( + parameter_file, "MWPotential2014Potential:potential_factors", 3, + potential->f); + + /* Convert to internal system of units by using the + * physical constants defined in this system */ + potential->M_200 *= phys_const->const_solar_mass; + potential->H *= phys_const->const_reduced_hubble; + potential->Mdisk *= phys_const->const_solar_mass; + potential->Rdisk *= 1000. * phys_const->const_parsec; + potential->Zdisk *= 1000. * phys_const->const_parsec; + potential->r_1 *= 1000. * phys_const->const_parsec; + potential->r_c *= 1000. * phys_const->const_parsec; + + /* Compute rho_c */ + const double rho_c = 3.0 * potential->H * potential->H / + (8.0 * M_PI * phys_const->const_newton_G); + + /* Compute R_200 */ + const double R_200 = + cbrtf(3.0 * potential->M_200 / (4. * M_PI * 200.0 * rho_c)); + + /* NFW scale-radius */ + potential->r_s = R_200 / potential->c_200; + const double r_s3 = potential->r_s * potential->r_s * potential->r_s; + + /* Log(c_200) term appearing in many expressions */ + potential->log_c200_term = + log(1. + potential->c_200) - potential->c_200 / (1. + potential->c_200); + + const double rho_0 = + potential->M_200 / (4.f * M_PI * r_s3 * potential->log_c200_term); + + /* Pre-factor for the accelerations (note G is multiplied in later on) */ + potential->pre_factor = 4.0f * M_PI * rho_0 * r_s3; + + /* Prefactor for the mass of the PSC profile */ + potential->prefactor_psc_1 = 2.0 * M_PI * potential->amplitude * + pow(potential->r_1, potential->alpha) * + pow(potential->r_c, 3.0 - potential->alpha); + + /* Gamma function value for the mass of the PSC profile */ + potential->gamma_psc = gsl_sf_gamma(1.5 - 0.5 * potential->alpha); + + /* Prefactor for the potential of the PSC profile */ + potential->prefactor_psc_2 = 2.0 * M_PI * potential->amplitude * + pow(potential->r_1, potential->alpha) * + pow(potential->r_c, 2.0 - potential->alpha); + + /* Compute the orbital time at the softening radius */ + const double sqrtgm = sqrt(phys_const->const_newton_G * potential->M_200); + const double epslnthing = log(1.0 + potential->eps / potential->r_s) - + potential->eps / (potential->eps + potential->r_s); + + potential->mintime = 2. * M_PI * potential->eps * sqrtf(potential->eps) * + sqrtf(potential->log_c200_term / epslnthing) / sqrtgm * + potential->timestep_mult; +#else + error("Code not compiled with GSL. Can't compute MWPotential2014."); +#endif +} + +/** + * @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 'MWPotential2014' " + "with properties (in internal units) are " + "(x,y,z) = " + "(%e, %e, %e), c_200 = %e, M_200 = %e, H = %e, M_disk = %e, R_disk = %e, " + "z_disk = %e, amplitude = %e, r_1 = %e, alpha = %e, r_c = %e, timestep " + "multiplier = %e mintime = %e", + potential->x[0], potential->x[1], potential->x[2], potential->c_200, + potential->M_200, potential->H, potential->Mdisk, potential->Rdisk, + potential->Zdisk, potential->amplitude, potential->r_1, potential->alpha, + potential->r_c, potential->timestep_mult, potential->mintime); +} + +#endif /* SWIFT_POTENTIAL_MWPotential2014_H */