added Hernquist halo as an example for an external gravitational potential

Makefile

@@ -288,7 +288,7 @@ INCL    += time_integration/timestep.h time_integration/driftfac.h
 
 SUBDIRS += gravity
 OBJS    += gravity/gravity.o gravity/ewald.o gravity/ewald_test.o \
-           gravity/grav_forcetest.o \
+           gravity/grav_forcetest.o gravity/grav_external.o \
            gravity/grav_direct.o gravity/second_order_ics.o
 INCL    += gravity/ewald.h gravity/ewaldtensors.h gravity/grav_forcetest.h
 

Template-Config.sh

@@ -31,7 +31,7 @@ SELFGRAVITY                                   # switch to enable self-gravity of
 #EXTRA_HIGH_EWALD_ACCURACY                    # this uses third-order instead of second-order Taylor expansion to interpolate Ewald corrections from table 
 #ALLOW_DIRECT_SUMMATION                       # allows calculation of direct summation gravity force if only a tiny number of particles as active 
 #EXTERNALGRAVITY                              # switches on inclusion of external gravitational potential
- 
+#EXTERNALGRAVITY_STATICHQ                     # example for a simple external potential due to a Hernquist halo
 
 #--------------------------------------- TreePM Options
 

documentation/04_config-options.md

@@ -350,6 +350,13 @@ and/or parameters what this external field is.
 
 -------
 
+**EXTERNALGRAVITY_STATICHQ**
+
+Activates a simple external potential due to a Hernquist dark matter
+halo, with parameters specified in the parameterfile.
+
+-------
+
 
 TreePM Options                                             {#treepm}
 ==============

documentation/05_parameterfile.md

@@ -1165,3 +1165,25 @@ temperature onto the gas. This is however done in terms of a minimum
 energy per unit mass `u` , i.e. if `MinEgySpec` is set to some finite
 value, `u` will not be allowed to drop below this value.
 
+-------
+
+Special features                                      {#special}
+================
+
+
+**A_StaticHQHalo**        5.0
+
+In case the `EXTERNALGRAVITY_STATICHQ` option is activated, this
+specifies the scale length of a static Hernquist halo whose
+gravitational potential is added to the force calculation. The
+halo is centered at the origin, and the scale length is given
+in internal length units.
+
+-------
+
+**Mass_StaticHQHalo**     100.0
+
+This parameter is only active when `EXTERNALGRAVITY_STATICHQ` is
+enabled, and then gives the total mass (in internal units) of the
+halos that is added as a static potential to the force computation.
+

src/data/allvars.cc

@@ -194,6 +194,11 @@ void global_data_all_processes::register_parameters(void)
 #ifdef CREATE_GRID
   add_param("GridSize", &GridSize, PARAM_INT, PARAM_FIXED);
 #endif
+
+#ifdef EXTERNALGRAVITY_STATICHQ
+  add_param("A_StaticHQHalo", &A_StaticHQHalo, PARAM_DOUBLE, PARAM_FIXED);
+  add_param("Mass_StaticHQHalo", &Mass_StaticHQHalo, PARAM_DOUBLE, PARAM_FIXED);
+#endif
 }
 
 /*! \brief This function reads a table with a list of desired output times.

src/data/allvars.h

@@ -348,6 +348,11 @@ struct global_data_all_processes : public parameters
   int GridSize;
 #endif
 
+#ifdef EXTERNALGRAVITY_STATICHQ
+  double A_StaticHQHalo;
+  double Mass_StaticHQHalo;
+#endif
+
   void set_cosmo_factors_for_current_time(void);
   void register_parameters(void);
   void read_outputlist(char *fname);

src/data/constants.h

@@ -270,6 +270,10 @@
 #error "LIGHTCONE_OUTPUT_ACCELERATIONS only works with PMGRID if TREEPM_NOTIMESPLIT is used and HIERARCHICAL_GRAVITY is not used"
 #endif
 
+#if defined(EXTERNALGRAVITY_STATICHQ) && !defined(EXTERNALGRAVITY)
+#error "EXTERNALGRAVITY_STATICHQ only works when EXTERNALGRAVITY is activated"
+#endif
+
 #ifndef ASMTH
 /** ASMTH gives the scale of the short-range/long-range force split in units of FFT-mesh cells */
 #define ASMTH 1.25

src/gravity/grav_external.cc

+/*******************************************************************************
+ * \copyright   This file is part of the GADGET4 N-body/SPH code developed
+ * \copyright   by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
+ * \copyright   (vspringel@mpa-garching.mpg.de) and all contributing authors.
+ *******************************************************************************/
+
+/*! \file gravity_external.cc
+ *
+ * \brief can add an optional external gravity field
+ */
+
+#include "gadgetconfig.h"
+
+#ifdef EXTERNALGRAVITY
+
+#include <math.h>
+#include <mpi.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "../data/allvars.h"
+#include "../data/dtypes.h"
+#include "../data/intposconvert.h"
+#include "../data/mymalloc.h"
+#include "../domain/domain.h"
+#include "../fmm/fmm.h"
+#include "../gravtree/gravtree.h"
+#include "../logs/logs.h"
+#include "../logs/timer.h"
+#include "../main/simulation.h"
+#include "../mpi_utils/mpi_utils.h"
+#include "../pm/pm.h"
+#include "../system/system.h"
+#include "../time_integration/timestep.h"
+
+void sim::gravity_external(void)
+{
+  for(int i = 0; i < Sp.TimeBinsGravity.NActiveParticles; i++)
+    {
+      int target = Sp.TimeBinsGravity.ActiveParticleList[i];
+
+#if defined(EVALPOTENTIAL) || defined(OUTPUT_POTENTIAL)
+      Sp.P[target].ExtPotential = 0;
+#endif
+
+#ifdef EXTERNALGRAVITY_STATICHQ
+      {
+        vector<double> pos;
+        Sp.intpos_to_pos(Sp.P[target].IntPos, pos.da); /* converts the integer coordinate to floating point */
+
+        double r = sqrt(pos.r2());
+
+        double m = All.Mass_StaticHQHalo * pow(r / (r + All.A_StaticHQHalo), 2);
+
+        if(r > 0)
+          Sp.P[target].GravAccel += (-All.G * m / (r * r * r)) * pos;
+
+#if defined(EVALPOTENTIAL) || defined(OUTPUT_POTENTIAL)
+        Sp.P[target].ExtPotential += (-All.G * m / (r + All.A_StaticHQHalo));
+#endif
+      }
+#endif
+    }
+}
+
+#endif

src/main/simulation.h

@@ -176,7 +176,7 @@ class sim : public pinning, public test_io_bandwidth
 #endif
 
 #ifdef EXTERNALGRAVITY
-  void gravity_external(void) { Terminate("Currently, no external gravity field is implemented.\n"); }
+  void gravity_external(void);
 #endif
 };