updated readme.txt

examples/Disk-Patch/HydroStatic/disk-patch.yml

@@ -20,7 +20,7 @@ Statistics:
 # Parameters governing the snapshots
 Snapshots:
   basename:            Disk-Patch-dynamic # Common part of the name of output files
-  time_first:          968.         # Time of the first output (in internal units)
+  time_first:          968.       # Time of the first output (in internal units)
   delta_time:         24.         # Time difference between consecutive outputs (in internal units)
 
 # Parameters for the hydrodynamics scheme

examples/Disk-Patch/HydroStatic/makeIC.py

@@ -203,6 +203,8 @@ if (entropy_flag == 1):
 else:
     ds[()] = u    
 
+print u
+
 ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
 ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
 ds[()] = ids

examples/Disk-Patch/HydroStatic/readme.text

@@ -4,9 +4,13 @@ equations from Creasey, Theuns & Bower 2013.
 
 Run the swift using setting
 #define EXTERNAL_POTENTIAL_DISK_PATCH
-#define ISOTHERMAL_GLASS
-in const.h, geneate the initial conditions (a copy of the glass), and
-run this for 15 dynamical times
+#define EXTERNAL_POTENTIAL_DISK_PATCH_ICs
+#define ISOTHERMAL_GLASS (20.2615290634))
+undefine EOS_IDEAL_GAS
+in const.h.
+
+Generate ICs (pythin makeIC.py)
+run this using disk-patch-icc.py 
 (1) imposing the gas remains isothermal
 (2) particles suffer a viscous force
 (3) the gravitational force increases from 0 to its full value in 5
@@ -16,5 +20,5 @@ dynamica times
 Then, use the output as new initial conditions, run using
 #define EXTERNAL_POTENTIAL_DISK_PATCH
 only, and see whether it stays in hydrotatic equilibrium
-
+Use disk-patch.yml as parameter file
 

examples/Disk-Patch/HydroStatic/test.pro

@@ -8,7 +8,7 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f
 @physunits
 
 indir    = './'
-basefile = 'Disk-Patch-dynamic_'
+basefile = 'Disk-Patch_'
 
 ; set properties of potential
 uL   = phys.pc                  ; unit of length
@@ -121,72 +121,21 @@ x0 = reform(xout[0,*])
 y0 = reform(xout[1,*])
 z0 = reform(xout[2,*])
 
-; calculate relative energy change
-de    = 0.0 * eout
-dl    = 0.0 * lout
-nsave = isave
-for ifile=1, nsave-1 do de[*,ifile] = (eout[*,ifile]-eout[*,0])/eout[*,0]
-for ifile=1, nsave-1 do dl[*,ifile] = (lout[*,ifile] - lout[*,0])/lout[*,0]
-
-
-; calculate statistics of energy changes
-print,' relatve energy change: (per cent) ',minmax(de) * 100.
-print,' relative Lz    change: (per cent) ',minmax(dl) * 100.
-
-; plot enery and Lz conservation for some particles
-if(iplot eq 1) then begin
-   nplot = nfollow
-
-   ; plot density profile
-   wset,0
-   xr   = [0, 3*scale_height]
-   nbins = 100
-   zpos  = findgen(nbins)/float(nbins-1) * max(xr)
-   dens  = (surface_density/(2.d0*scale_height)) * 1./cosh(zpos/scale_height)^2
-   plot,[0],[0],xr=xr,/xs,yr=[0,max(dens)*1.4],/ys,/nodata,xtitle='|z|',ytitle='density'
-   oplot,zpos,dens,color=black,thick=3
-   oplot,abs(zout[*,1]),rout[*,1],psym=3
-   oplot,abs(zout[*,nsave-1]),rout[*,nsave-1],psym=3,color=red
-
-;; ; plot results on energy conservation for some particles
-;;    nplot = min(100, nfollow)
-;;    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 (%)'
-;;    for i=0,nplot-1 do oplot,tout,de[i,*]
-;;    for i=0,nplot-1 do oplot,tout,dl[i,*],color=red
-;;    legend,['dE/E','dL/L'],linestyle=[0,0],color=[black,red],box=0,/bottom,/left
-;;    screen_to_png,'e-time.png'
-
-;  plot vertical oscillation
-   wset,2
-   xr = [min(tout), max(tout)]
-   yr = [-3,3]*scale_height
-   plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='z(t)'
-   color = floor(findgen(nplot)*255/float(nplot))
-   for i=0,nplot-1,50 do oplot,tout,zout[i,*],color=color(i)
-   screen_to_png,'orbit.png'
-
-
-;   plot evolution of density for some particles close to disk
-   wset, 6
-   rr = (surface_density/(2.d0*scale_height)) * 1./cosh(abs(zout)/scale_height)^2 ; desired density
-   gd = where(abs(zout[*,0]) lt 50, ng)
-   plot,[0],[0],xr=[0, max(tout)], yr=10.^[-1,1], /xs, /ys, /nodata, /yl
-   nplot = min(40, ng)
-   color = floor(findgen(nplot)/float(nplot)*255)
-   for i=0, nplot-1 do oplot,tout[1:*],rout[gd[i],1:*]/rr[gd[i], 1:*], color=color[i],psym=-4
-
-
-;; ; make histogram of energy changes at end
-;;    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'
-
-
-endif
+
+; plot density profile and compare to analytic profile
+nplot = nfollow
+
+                                ; plot density profile
+wset,0
+xr   = [0, 3*scale_height]
+nbins = 100
+zpos  = findgen(nbins)/float(nbins-1) * max(xr)
+dens  = (surface_density/(2.d0*scale_height)) * 1./cosh(zpos/scale_height)^2
+plot,[0],[0],xr=xr,/xs,yr=[0,max(dens)*1.4],/ys,/nodata,xtitle='|z|',ytitle='density'
+oplot,zpos,dens,color=black,thick=3
+;oplot,abs(zout[*,1]),rout[*,1],psym=3 ; initial profile
+oplot,abs(zout[*,nsave-1]),rout[*,nsave-1],psym=3,color=red
+
 
 end
 

examples/main.c

@@ -329,10 +329,11 @@ int main(int argc, char *argv[]) {
   /* Initialise the hydro properties */
   struct hydro_props hydro_properties;
   hydro_props_init(&hydro_properties, params);
+  if (with_hydro && myrank == 0) eos_print(); 
 
   /* Initialise the external potential properties */
   struct external_potential potential;
-  if (with_external_gravity) potential_init(params, &us, &potential);
+  if (with_external_gravity) potential_init(params, &prog_const, &us, &potential);
   if (with_external_gravity && myrank == 0) potential_print(&potential);
 
   /* Read particles and space information from (GADGET) ICs */

src/const.h

@@ -43,7 +43,8 @@
 
 /* Equation of state choice */
 #define EOS_IDEAL_GAS
-//#define EOS_ISOTHERMAL_GAS
+/* if EOS_ISOTHERMAL_GAS is defined, keep thermal energy per unit mass equal to EOS_ISOTHERMAL_GAS in programme units */
+//#define EOS_ISOTHERMAL_GAS (20.2615290634)
 
 /* Kernel function to use */
 #define CUBIC_SPLINE_KERNEL
@@ -68,7 +69,9 @@
 //#define EXTERNAL_POTENTIAL_POINTMASS
 //#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
 #define EXTERNAL_POTENTIAL_DISK_PATCH
-//#define ISOTHERMAL_GLASS
+/* Add viscuous force to gas particles to speed-up glass making for disk-patch ICs */
+//#define EXTERNAL_POTENTIAL_DISK_PATCH_ICS
+
 
 /* Are we debugging ? */
 //#define SWIFT_DEBUG_CHECKS

src/equation_of_state.h

@@ -33,6 +33,9 @@
 
 /* ------------------------------------------------------------------------- */
 #if defined(EOS_IDEAL_GAS)
+#if defined(EOS_ISOTHERMAL_GAS)
+#error: cannot have ideal and isothermal gas at the same time!
+#endif
 
 /**
  * @brief Returns the internal energy given density and entropy
@@ -118,7 +121,17 @@ gas_soundspeed_from_internal_energy(float density, float u) {
 
   return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
 }
-
+__attribute__((always_inline)) INLINE static void eos_print(){
+  message(" assuming an equation of state of an ideal gas, with gamma = %e", hydro_gamma);
+}
+__attribute__((always_inline)) INLINE static float hydro_update_entropy(float density, float entropy)
+{
+  return entropy;
+}
+__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(float density, float entropy)
+{
+  return hydro_gamma_minus_one * pow_minus_gamma_minus_one(density);
+}
 /* ------------------------------------------------------------------------- */
 #elif defined(EOS_ISOTHERMAL_GAS)
 
@@ -131,10 +144,8 @@ gas_soundspeed_from_internal_energy(float density, float u) {
 __attribute__((always_inline)) INLINE static float
 gas_internal_energy_from_entropy(float density, float entropy) {
 
-  error("Missing definition !");
-  return 0.f;
+  return EOS_ISOTHERMAL_GAS;
 }
-
 /**
  * @brief Returns the pressure given density and entropy
  *
@@ -144,8 +155,7 @@ gas_internal_energy_from_entropy(float density, float entropy) {
 __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
     float density, float entropy) {
 
-  error("Missing definition !");
-  return 0.f;
+  return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS * density;
 }
 
 /**
@@ -157,8 +167,7 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
 __attribute__((always_inline)) INLINE static float
 gas_entropy_from_internal_energy(float density, float u) {
 
-  error("Missing definition !");
-  return 0.f;
+  return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS / pow(density, hydro_gamma-1);
 }
 
 /**
@@ -186,7 +195,18 @@ gas_soundspeed_from_internal_energy(float density, float u) {
   error("Missing definition !");
   return 0.f;
 }
-
+__attribute__((always_inline)) INLINE static void eos_print(){
+  message(" assuming an equation of state for an isothermal gas with utherm= %e and gamma = %e", EOS_ISOTHERMAL_GAS, hydro_gamma);
+}
+__attribute__((always_inline)) INLINE static float hydro_update_entropy(float density, float entropy)
+{
+  float thermal_energy = EOS_ISOTHERMAL_GAS;
+  return gas_entropy_from_internal_energy(density, thermal_energy);
+}
+__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(float density, float entropy)
+{
+  return 0;
+}
 /* ------------------------------------------------------------------------- */
 #else
 

src/hydro/Gadget2/hydro.h

@@ -308,8 +308,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part *restrict p) {
 
   p->force.h_dt *= p->h * 0.333333333f;
-
-  p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
+  p->entropy     = hydro_update_entropy(p->rho, p->entropy);
+  p->entropy_dt *= hydro_update_entropy_rate(p->rho, p->entropy);
+
+	 
+/* #ifdef EOS_ISOTHERMAL_GAS */
+/*   p->entropy_dt = 0; */
+/*   float dummy = 1; */
+/*   p->entropy    = gas_entropy_from_internal_energy(p->rho, dummy); */
+/* #else */
+/*   p->entropy_dt *= hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho); */
+/* #endif */
 }
 
 /**

src/hydro/Gadget2/hydro_io.h

@@ -20,6 +20,7 @@
 #include "adiabatic_index.h"
 #include "io_properties.h"
 #include "kernel_hydro.h"
+#include "equation_of_state.h"
 
 /**
  * @brief Specifies which particle fields to read from a dataset
@@ -54,8 +55,10 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
 
 float convert_u(struct engine* e, struct part* p) {
 
-  return p->entropy * pow_gamma_minus_one(p->rho) *
-         hydro_one_over_gamma_minus_one;
+  return gas_internal_energy_from_entropy(p->rho, p->entropy);
+
+  /* return p->entropy * pow_gamma_minus_one(p->rho) * */
+  /*        hydro_one_over_gamma_minus_one; */
 }
 
 /**

src/hydro/Gadget2/hydro_part.h

@@ -64,7 +64,7 @@ struct part {
   /* Derivative of the density with respect to smoothing length. */
   float rho_dh;
 
-  union {
+  //  union {
 
     struct {
 
@@ -100,7 +100,6 @@ struct part {
       float h_dt;
 
     } force;
-  };
 
   /* Particle ID. */
   long long id;

src/potentials.c

@@ -33,6 +33,7 @@
  * @param potential The external potential properties to initialize
  */
 void potential_init(const struct swift_params* parameter_file,
+						  const struct phys_const* const phys_const,
                     struct UnitSystem* us,
                     struct external_potential* potential) {
 
@@ -74,6 +75,7 @@ void potential_init(const struct swift_params* parameter_file,
 	 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");
+  potential -> disk_patch_potential.dynamical_time = sqrt(potential->disk_patch_potential.scale_height / (phys_const->const_newton_G * potential->disk_patch_potential.surface_density));
 #endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
 
 }
@@ -111,5 +113,8 @@ void potential_print(const struct external_potential* potential) {
 			 potential->disk_patch_potential.z_disk,
 			 potential->disk_patch_potential.scale_height,
 			 potential->disk_patch_potential.timestep_mult);
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
+  message("Disk-patch potential: imposing growth of gravity over time, and adding viscous force to gravity");
+#endif
 #endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
 }

src/potentials.h

@@ -36,8 +36,6 @@
 #include "physical_constants.h"
 #include "units.h"
 #include "math.h"
-#define FINLINE __attribute__((always_inline)) INLINE
-/* #define FINLINE */
 /* External Potential Properties */
 struct external_potential {
 
@@ -61,6 +59,7 @@ struct external_potential {
 	 double surface_density;
 	 double scale_height;
 	 double z_disk;
+	 double dynamical_time;
 	 double timestep_mult;
   } disk_patch_potential;
 #endif
@@ -74,17 +73,16 @@ struct external_potential {
  * @param phys_cont The physical constants in internal units.
  * @param gp Pointer to the g-particle data.
  */
-FINLINE static float external_gravity_disk_patch_timestep(
+__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) {
 
   
   /* initilize time step to disk dynamical time */
-  const float dt_dyn = sqrt(potential->disk_patch_potential.scale_height / (phys_const->const_newton_G * potential->disk_patch_potential.surface_density));
+  const float dt_dyn = potential->disk_patch_potential.dynamical_time;
   float dt = dt_dyn;
 
-
   /* absolute value of height above disk */
   const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
 
@@ -127,27 +125,23 @@ FINLINE static float external_gravity_disk_patch_timestep(
  * @param phys_cont The physical constants in internal units.
  * @param gp Pointer to the g-particle data.
  */
-FINLINE static void external_gravity_disk_patch_potential(const double time, const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+__attribute__((always_inline)) INLINE static void external_gravity_disk_patch_potential(const double time, const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
 
   const float G_newton = phys_const->const_newton_G;
-  const float surface_density= potential->disk_patch_potential.surface_density;
-  const float scale_height   = potential->disk_patch_potential.z_disk;
-  const double const_G       = phys_const->const_newton_G;
-  const float t_dyn          = sqrt(scale_height / (const_G * surface_density));
-
-
+  const float t_dyn = potential->disk_patch_potential.dynamical_time;
   const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
   g->a_grav[0]  += 0;
   g->a_grav[1]  += 0;
 
-#ifdef ISOTHERMAL_GLASS
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
   float reduction_factor = 1.;
   if(time < 5 * t_dyn)
 	 reduction_factor = time / (5 * t_dyn);
 #else
   const float reduction_factor = 1.;
 #endif
-  const float z_accel =  reduction_factor * 2 * phys_const->const_newton_G * M_PI * potential->disk_patch_potential.surface_density
+
+  const float z_accel =  reduction_factor * 2 * G_newton * M_PI * potential->disk_patch_potential.surface_density
   	 * tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
 
   if (dz > 0)
@@ -159,16 +153,12 @@ FINLINE static void external_gravity_disk_patch_potential(const double time, con
   if(abs(g->id_or_neg_offset) == 1)
 	 message(" time= %e, rf= %e, az= %e", time, reduction_factor, g->a_grav[2]);
 
-#ifdef ISOTHERMAL_GLASS
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
   /* TT: add viscous drag */
-  
-  g->a_grav[0]  -= g->v_full[0] / (2. * t_dyn) / const_G;
-  g->a_grav[1]  -= g->v_full[1] / (2. * t_dyn) / const_G;
-  g->a_grav[2]  -= g->v_full[2] / (2. * t_dyn) / const_G;
-
+  g->a_grav[0]  -= g->v_full[0] / (2 * t_dyn) / G_newton;
+  g->a_grav[1]  -= g->v_full[1] / (2 * t_dyn) / G_newton;
+  g->a_grav[2]  -= g->v_full[2] / (2 * t_dyn) / G_newton;
 #endif  
-
-
 }
 #endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
 
@@ -219,7 +209,7 @@ external_gravity_isothermalpotential_timestep(
  * @param phys_const The physical constants in internal units.
  * @param g Pointer to the g-particle data.
  */
-FINLINE static void external_gravity_isothermalpotential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+__attribute__((always_inline)) INLINE static void external_gravity_isothermalpotential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
 __attribute__((always_inline)) INLINE static void
 external_gravity_isothermalpotential(const struct external_potential* potential,
                                      const struct phys_const* const phys_const,
@@ -252,7 +242,7 @@ external_gravity_isothermalpotential(const struct external_potential* potential,
  * @param phys_const The physical constants in internal units.
  * @param g Pointer to the g-particle data.
  */
-FINLINE static float external_gravity_pointmass_timestep(const struct external_potential* potential,
+__attribute__((always_inline)) INLINE static float external_gravity_pointmass_timestep(const struct external_potential* potential,
 																					  const struct phys_const* const phys_const,
                                     const struct gpart* const g) {
 
@@ -287,7 +277,7 @@ FINLINE static float external_gravity_pointmass_timestep(const struct external_p
  * @param phys_const The physical constants in internal units.
  * @param g Pointer to the g-particle data.
  */
-FINLINE static void external_gravity_pointmass(
+__attribute__((always_inline)) INLINE static void external_gravity_pointmass(
     const struct external_potential* potential,
     const struct phys_const* const phys_const, struct gpart* g) {
 
@@ -306,6 +296,7 @@ FINLINE static void external_gravity_pointmass(
 
 /* Now, some generic functions, defined in the source file */
 void potential_init(const struct swift_params* parameter_file,
+						  const struct phys_const* const phys_const,
                     struct UnitSystem* us,
                     struct external_potential* potential);