diff --git a/examples/Disk-Patch/HydroStatic/disk-patch.yml b/examples/Disk-Patch/HydroStatic/disk-patch.yml
index d825a6c46e3b8607ef17810fffe20108aab73f60..f7df277501825a677c9b6c36730ee0a6a1d8a0ff 100644
--- a/examples/Disk-Patch/HydroStatic/disk-patch.yml
+++ b/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
diff --git a/examples/Disk-Patch/HydroStatic/makeIC.py b/examples/Disk-Patch/HydroStatic/makeIC.py
index 1dabfad77a4c73beee9cf8f8859dc85169e180ab..b7bf117ad56b68a5a792c1c16e63f2bc3ea4fe8d 100644
--- a/examples/Disk-Patch/HydroStatic/makeIC.py
+++ b/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
diff --git a/examples/Disk-Patch/HydroStatic/readme.text b/examples/Disk-Patch/HydroStatic/readme.text
index acb052a36a863d4b41a3162f4da9b45a7c66aa83..314334fd0383db0a74746da85cb3ab298aa3b060 100644
--- a/examples/Disk-Patch/HydroStatic/readme.text
+++ b/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
diff --git a/examples/Disk-Patch/HydroStatic/test.pro b/examples/Disk-Patch/HydroStatic/test.pro
index 90a62d3135a77a8dd1a3355fddb6755a14d32692..31e027e3a308b04f1c59222e1a339786857061ac 100644
--- a/examples/Disk-Patch/HydroStatic/test.pro
+++ b/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
diff --git a/examples/main.c b/examples/main.c
index 30e27cc55a6ec508dc352f1acc7f2eed71d71bdd..4aa7e7cd53c6df8b2ea5e611e902fce627a6dea6 100644
--- a/examples/main.c
+++ b/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 */
diff --git a/src/const.h b/src/const.h
index 08b56dbbaee0cb3334e5d4d04d45fcf78ee0d17d..f98d2ba28d0236448fe1c23e54072f8bd41156d0 100644
--- a/src/const.h
+++ b/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
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index d9f1a46f240a8bdf965749a5513ca3448d08aeb9..5103d6e475f38f2166447c30b87b9798373df8f9 100644
--- a/src/equation_of_state.h
+++ b/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
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 64664e77fed1430090386d1b48049a8e47092fe9..a415e725abdf6e9e2262200ee95e871ac751cb81 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/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 */
}
/**
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 1614a483ac3007d37770e04d40833f78d19a8cd2..17aae18ab07c3c6476c0e48ea1a657ca65df71dc 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/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; */
}
/**
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 0cd830c51a1836c0a56e7b4097990bf4d07a101f..3145390e3fb9e3b05dadd1cac56b98cec1d59e38 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/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;
diff --git a/src/potentials.c b/src/potentials.c
index 9a4480adadf8186e34174fa1c4d0b766a4a03d4c..26507622c3f93f86445181291172dd958f31f17a 100644
--- a/src/potentials.c
+++ b/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 */
}
diff --git a/src/potentials.h b/src/potentials.h
index 3c7cd0e750faf093a5b36e7ed8d08b16695c32cb..33bc64d57f981ef8e64c84e43852fd2f77dcc4ff 100644
--- a/src/potentials.h
+++ b/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);