Capture progress on disk-patch

examples/IsothermalPotential/GravityOnly/test.pro

@@ -127,7 +127,7 @@ print,' relative Lz    change: (per cent) ',minmax(dl) * 100.
 if(iplot eq 1) then begin
 ; plot results on energy conservation for some particles
    nplot = min(10, nfollow)
-   wset,0
+   win,0
    xr = [min(tout), max(tout)]
    yr = [-2,2]*1d-2             ; in percent
    plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='time',ytitle='dE/E, dL/L (%)'
@@ -137,7 +137,7 @@ if(iplot eq 1) then begin
    screen_to_png,'e-time.png'
 
 ;  plot orbits of those particles
-   wset,2
+   win,2
    xr = [-100,100]
    yr = xr
    plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
@@ -146,7 +146,7 @@ if(iplot eq 1) then begin
    screen_to_png,'orbit.png'
 
 ; plot radial position of these particles
-   wset,4
+   win,4
    xr = [min(tout), max(tout)]
    yr = [0,80]
    plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='r'
@@ -155,7 +155,7 @@ for i=0,nplot-1 do begin dr = sqrt(reform(xout[i,*])^2 + reform(yout[i,*])^2) &
    screen_to_png,'r-time.png'
 
 ; make histogram of energy changes at end
-   wset,6
+   win,6
    ohist,de,x,y,-0.05,0.05,0.001
    plot,x,y,psym=10,xtitle='de (%)'
    screen_to_png,'de-hist.png'

src/const.h

@@ -65,8 +65,9 @@
 #define const_gravity_eta 0.025f
 
 /* External gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
+//#define EXTERNAL_POTENTIAL_POINTMASS
 //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+#define EXTERNAL_POTENTIAL_DISK_PATCH
 
 /* Are we debugging ? */
 //#define SWIFT_DEBUG_CHECKS

src/gravity/Default/gravity.h

@@ -46,7 +46,9 @@ gravity_compute_timestep_external(const struct external_potential* potential,
   dt = fminf(dt, external_gravity_isothermalpotential_timestep(potential,
                                                                phys_const, gp));
 #endif
-
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+  dt = fmin(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
+#endif
   return dt;
 }
 
@@ -110,7 +112,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
  * @param const_G Newton's constant
  */
 __attribute__((always_inline)) INLINE static void gravity_end_force(
-    struct gpart* gp, double const_G) {
+	struct gpart* gp, const double const_G) {
 
   /* Let's get physical... */
   gp->a_grav[0] *= const_G;
@@ -137,6 +139,9 @@ __attribute__((always_inline)) INLINE static void external_gravity(
 #ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
   external_gravity_isothermalpotential(potential, phys_const, gp);
 #endif
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+  external_gravity_disk_patch_potential(potential, phys_const, gp);
+#endif
 }
 
 /**

src/potentials.c

@@ -65,6 +65,17 @@ void potential_init(const struct swift_params* parameter_file,
       parameter_file, "IsothermalPotential:timestep_mult");
 
 #endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+  potential -> disk_patch_potential.surface_density = 
+	 parser_get_param_double(parameter_file, "Disk-PatchPotential:surface_density");
+  potential -> disk_patch_potential.scale_height = 
+	 parser_get_param_double(parameter_file, "Disk-PatchPotential:scale_height");
+  potential -> disk_patch_potential.z_disk = 
+	 parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
+  potential -> disk_patch_potential.timestep_mult = 
+	 parser_get_param_double(parameter_file, "Disk-PatchPotential:timestep_mult");
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
+
 }
 
 /**
@@ -94,4 +105,11 @@ void potential_print(const struct external_potential* potential) {
       potential->isothermal_potential.timestep_mult);
 
 #endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+  message("Disk-patch potential properties are surface_density = %e disk height= %e scale height= %e timestep multiplier= %e",
+			 potential->disk_patch_potential.surface_density,
+			 potential->disk_patch_potential.z_disk,
+			 potential->disk_patch_potential.scale_height,
+			 potential->disk_patch_potential.timestep_mult);
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
 }

src/potentials.h

@@ -34,6 +34,7 @@
 #include "part.h"
 #include "physical_constants.h"
 #include "units.h"
+#include "math.h"
 
 /* External Potential Properties */
 struct external_potential {
@@ -53,9 +54,66 @@ struct external_potential {
     float timestep_mult;
   } isothermal_potential;
 #endif
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+  struct {
+	 double surface_density;
+	 double scale_height;
+	 double z_disk;
+	 double timestep_mult;
+  } disk_patch_potential;
+#endif
 };
 
-/* Include external isothermal potential */
+/* external potential: disk_patch */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+/**
+ * @brief Computes the time-step due to the acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
+ *
+ * @param phys_cont The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline))
+    INLINE static float external_gravity_disk_patch_timestep(
+		  const struct external_potential* potential,
+        const struct phys_const* const phys_const,
+        const struct gpart* const g) {
+
+  const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
+  const float z_accel =  2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
+	 * tanh(dz / potential->disk_patch_potential.scale_height);
+  const float dz_accel_over_dt = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density /
+	 potential->disk_patch_potential.scale_height / cosh(dz / potential->disk_patch_potential.scale_height) / cosh(dz / potential->disk_patch_potential.scale_height) * fabs(g->v_full[2]);
+
+  /* demand that time step * derivative of acceleration < fraction of acceleration */
+  const float dt1 = potential->disk_patch_potential.timestep_mult * z_accel / dz_accel_over_dt;
+	 
+  /* demand that dt^2 * acceleration < fraction of scale height of disk */
+  const float dt2 = potential->disk_patch_potential.timestep_mult * potential->disk_patch_potential.scale_height / fabs(z_accel);
+  
+  float dt = dt1;
+  if(dt2 < dt) dt = dt2;
+  return dt2;
+}
+/**
+ * @brief Computes the gravitational acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
+ *
+ * @param phys_cont The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_disk_patch_potential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+
+  const float dz    = g->x[2] - potential->disk_patch_potential.z_disk;
+  g->a_grav[0]  = 0;
+  g->a_grav[1]  = 0;
+  const float z_accel =  2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
+	 * tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
+  if (dz > 0)
+	 g->a_grav[2] -= z_accel;
+  if (dz < 0)
+	 g->a_grav[2] += z_accel;
+}
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
+
 #ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 
 /**