diff --git a/examples/Disk-Patch/HydroStatic/disk-patch.yml b/examples/Disk-Patch/HydroStatic/disk-patch.yml
index 46cef05130eb337de25d15d9f34e2c346a9e314e..d825a6c46e3b8607ef17810fffe20108aab73f60 100644
--- a/examples/Disk-Patch/HydroStatic/disk-patch.yml
+++ b/examples/Disk-Patch/HydroStatic/disk-patch.yml
@@ -8,10 +8,10 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 480. # The end time of the simulation (in internal units).
- dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units).
- dt_max: 4. # The maximal time-step size of the simulation (in internal units).
+ time_begin: 968 # The starting time of the simulation (in internal units).
+ time_end: 12000. # The end time of the simulation (in internal units).
+ dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1. # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
@@ -19,9 +19,9 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
- basename: Disk-Patch # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
- delta_time: 8. # Time difference between consecutive outputs (in internal units)
+ basename: Disk-Patch-dynamic # Common part of the name of output files
+ 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
SPH:
@@ -33,7 +33,7 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
- file_name: Disk-Patch.hdf5 # The file to read
+ file_name: Disk-Patch-dynamic.hdf5 # The file to read
h_scaling: 1. # A scaling factor to apply to all smoothing lengths in the ICs.
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
diff --git a/examples/Disk-Patch/HydroStatic/dynamic.pro b/examples/Disk-Patch/HydroStatic/dynamic.pro
new file mode 100644
index 0000000000000000000000000000000000000000..b598530d766d5c092802b9b0acc0aae9d067d8c6
--- /dev/null
+++ b/examples/Disk-Patch/HydroStatic/dynamic.pro
@@ -0,0 +1,104 @@
+;
+; test energy / angular momentum conservation of test problem
+;
+
+iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare final and initial energy
+
+; set physical constants
+@physunits
+
+indir = './'
+basefile = 'Disk-Patch-dynamic_'
+
+; set properties of potential
+uL = phys.pc ; unit of length
+uM = phys.msun ; unit of mass
+uV = 1d5 ; unit of velocity
+
+; properties of patch
+surface_density = 10.
+scale_height = 100.
+z_disk = 200.;
+; derived units
+constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
+pcentre = [0.,0.,z_disk] * pc / uL
+
+;
+infile = indir + basefile + '*'
+spawn,'ls -1 '+infile,res
+nfiles = n_elements(res)
+
+
+; choose: calculate change of energy and Lz, comparing first and last
+; snapshots for all particles, or do so for a subset
+
+; compare all
+ifile = 0
+inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+id = h5rd(inf,'PartType0/ParticleIDs')
+nfollow = n_elements(id)
+
+
+; compute anlytic profile
+nbins = 100
+zbins = findgen(nbins)/float(nbins-1) * 2 * scale_height
+rbins = (surface_density/(2.*scale_height)) / cosh(abs(zbins)/scale_height)^2
+
+
+; plot analytic profile
+wset,0
+plot,[0],[0],xr=[0,2*scale_height],yr=[0,max(rbins)],/nodata,xtitle='|z|',ytitle=textoidl('\rho')
+oplot,zbins,rbins
+
+ifile = 0
+nskip = nfiles - 1
+isave = 0
+nplot = 8192 ; randomly plot particles
+color = floor(findgen(nfiles)/float(nfiles-1)*255)
+;for ifile=0,nfiles-1,nskip do begin
+tsave = [0.]
+for ifile=0,nfiles-1,nskip do begin
+ inf = indir + basefile + strtrim(string(ifile,'(i3.3)'),1) + '.hdf5'
+ time = h5ra(inf, 'Header','Time')
+ tsave = [tsave, time]
+ print,' time= ',time
+ p = h5rd(inf,'PartType0/Coordinates')
+ v = h5rd(inf,'PartType0/Velocities')
+ id = h5rd(inf,'PartType0/ParticleIDs')
+ rho = h5rd(inf,'PartType0/Density')
+ h = h5rd(inf,'PartType0/SmoothingLength')
+ utherm = h5rd(inf,'PartType0/InternalEnergy')
+ indx = sort(id)
+
+; substract disk centre
+ for ic=0,2 do p[ic,*]=p[ic,*] - pcentre[ic]
+
+
+;; ; if you want to sort particles by ID
+;; id = id[indx]
+;; rho = rho[indx]
+;; utherm = utherm[indx]
+;; h = h[indx]
+;; for ic=0,2 do begin
+;; tmp = reform(p[ic,*]) & p[ic,*] = tmp[indx]
+;; tmp = reform(v[ic,*]) & v[ic,*] = tmp[indx]
+;; endfor
+
+ ip = floor(randomu(ifile+1,nplot)*n_elements(rho))
+ color = red
+ if(ifile eq 0) then color=black
+ oplot,abs(p[2,ip]), rho[ip], psym=3, color=color
+
+ isave = isave + 1
+
+endfor
+tsave = tsave[1:*]
+
+label = ['']
+for i=0,n_elements(tsave)-1 do label=[label,'t='+string(tsave[i],format='(f8.0)')]
+label = label[1:*]
+legend,[label[0],label[1]],linestyle=[0,0],color=[black,red],box=0,/top,/right
+
+end
+
+
diff --git a/examples/Disk-Patch/HydroStatic/makeIC.py b/examples/Disk-Patch/HydroStatic/makeIC.py
index cfc3ac638870535b7fb273f1e99bced2353e0f3f..1dabfad77a4c73beee9cf8f8859dc85169e180ab 100644
--- a/examples/Disk-Patch/HydroStatic/makeIC.py
+++ b/examples/Disk-Patch/HydroStatic/makeIC.py
@@ -57,7 +57,7 @@ print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
# parameters of potential
surface_density = 10.
-scale_height = 400.
+scale_height = 100.
gamma = 5./3.
# derived units
@@ -69,7 +69,7 @@ v_disp = numpy.sqrt(2 * utherm)
soundspeed = numpy.sqrt(utherm / (gamma * (gamma-1.)))
t_dyn = numpy.sqrt(scale_height / (const_G * surface_density))
t_cross = scale_height / soundspeed
-print 'dynamical time = ',t_dyn,' sound crossing time = ',t_cross,' sound speed= ',soundspeed,' 3D velocity dispersion = ',v_disp
+print 'dynamical time = ',t_dyn,' sound crossing time = ',t_cross,' sound speed= ',soundspeed,' 3D velocity dispersion = ',v_disp,' thermal_energy= ',utherm
# Parameters
@@ -117,22 +117,25 @@ for ix in range(0,ntile):
glass_p = numpy.concatenate(glass_p, axis=0)
glass_h = numpy.concatenate(glass_h, axis=0)
+# random shuffle of glas ICs
+numpy.random.seed(12345)
+indx = numpy.random.rand(numpy.shape(glass_h)[0])
+indx = numpy.argsort(indx)
+glass_p = glass_p[indx, :]
+glass_h = glass_h[indx]
-# use some of these
-numGas = 5000
-
-
-# use these as ICs
+# select numGas of them
+numGas = 8192
pos = glass_p[0:numGas,:]
h = glass_h[0:numGas]
numGas = numpy.shape(pos)[0]
-print ' numGas = ', numGas
-
# compute furthe properties of ICs
column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
enclosed_mass = column_density * boxSize * boxSize
pmass = enclosed_mass / numGas
+meanrho = enclosed_mass / boxSize**3
+print 'pmass= ',pmass,' mean(rho) = ', meanrho,' entropy= ', (gamma-1) * utherm / meanrho**(gamma-1)
# desired density
rho = surface_density / (2. * scale_height) / numpy.cosh(abs(pos[:,2])/scale_height)**2
diff --git a/examples/Disk-Patch/HydroStatic/test.pro b/examples/Disk-Patch/HydroStatic/test.pro
index 675c30f6aabb1accabecb50f4314a1644c35d577..90a62d3135a77a8dd1a3355fddb6755a14d32692 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_'
+basefile = 'Disk-Patch-dynamic_'
; set properties of potential
uL = phys.pc ; unit of length
@@ -21,7 +21,7 @@ scale_height = 100.
; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
-pcentre = [0.,0.,300.] * pc / uL
+pcentre = [0.,0.,200.] * pc / uL
;
infile = indir + basefile + '*'
@@ -39,7 +39,7 @@ id = h5rd(inf,'PartType0/ParticleIDs')
nfollow = n_elements(id)
; follow a subset
-nfollow = min(4000, nfollow) ; number of particles to follow
+; nfollow = min(4000, nfollow) ; number of particles to follow
;
if (iplot eq 1) then begin
@@ -171,7 +171,7 @@ if(iplot eq 1) then begin
; plot evolution of density for some particles close to disk
wset, 6
- rr = (surface_density/(2.d0*scale_height)) * 1./cosh(abs(zout)/100.)^2 ; desired density
+ 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)
diff --git a/src/const.h b/src/const.h
index c14419ac58f0130857643c76929cd56d5373a81a..08b56dbbaee0cb3334e5d4d04d45fcf78ee0d17d 100644
--- a/src/const.h
+++ b/src/const.h
@@ -68,7 +68,7 @@
//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
#define EXTERNAL_POTENTIAL_DISK_PATCH
-#define ISOTHERMAL_GLASS
+//#define ISOTHERMAL_GLASS
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 74a68fdc335b95e393c1cfad6cd8d07c32ce5458..1fb67220977d029ba821ef66e93ee61cae1394e2 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -62,17 +62,19 @@ __attribute__((always_inline)) INLINE static float
gravity_compute_timestep_self(const struct phys_const* const phys_const,
const struct gpart* const gp) {
- const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
- gp->a_grav[1] * gp->a_grav[1] +
- gp->a_grav[2] * gp->a_grav[2];
+ /* const float ac2 = gp->a_grav[0] * gp->a_grav[0] + */
+ /* gp->a_grav[1] * gp->a_grav[1] + */
+ /* gp->a_grav[2] * gp->a_grav[2]; */
- const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
+ /* const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN; */
- const float dt = sqrt(2.f * const_gravity_eta * gp->epsilon / ac);
+ /* const float dt = sqrt(2.f * const_gravity_eta * gp->epsilon / ac); */
- return dt;
+ /* return dt; */
+ return FLT_MAX;
}
+
/**
* @brief Initialises the g-particles for the first time
*
@@ -130,6 +132,7 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
* @param gp The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void external_gravity(
+ const double time,
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* gp) {
@@ -140,7 +143,7 @@ __attribute__((always_inline)) INLINE static void external_gravity(
external_gravity_isothermalpotential(potential, phys_const, gp);
#endif
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- external_gravity_disk_patch_potential(potential, phys_const, gp);
+ external_gravity_disk_patch_potential(time, potential, phys_const, gp);
#endif
}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index a1448d8686a30de88e24629558f01661daca223a..64664e77fed1430090386d1b48049a8e47092fe9 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -225,7 +225,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
-
const float irho = 1.f / p->rho;
/* Divide the pressure by the density and density gradient */
diff --git a/src/kick.h b/src/kick.h
index 7c6dbd16fbf613b63ca3225a03ebb73f2a486b36..b57e13d4ebf27d3a366d571e7fd4cd819653f726 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -46,31 +46,11 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
gp->ti_begin = gp->ti_end;
gp->ti_end = gp->ti_begin + new_dti;
-
-#ifdef ISOTHERMAL_GLASS
- /* TT: add viscous drag */
- float a_tot[3] = {gp->a_grav[0], gp->a_grav[1], gp->a_grav[2]};
- const float surface_density= 10.;
- const float scale_height = 100;
- const float t_dyn = sqrt(scale_height / (const_G * surface_density));
- const double time = ti_start * timeBase;
- a_tot[0] -= gp->v_full[0] / (0.05 * t_dyn);
- a_tot[1] -= gp->v_full[1] / (0.05 * t_dyn);
- a_tot[2] -= gp->v_full[2] / (0.05 * t_dyn);
-
- /* Kick particles in momentum space */
- gp->v_full[0] += a_tot[0] * dt;
- gp->v_full[1] += a_tot[1] * dt;
- gp->v_full[2] += a_tot[2] * dt;
-
-#else
-
/* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt;
gp->v_full[1] += gp->a_grav[1] * dt;
gp->v_full[2] += gp->a_grav[2] * dt;
-#endif
/* Extra kick work */
gravity_kick_extra(gp, dt, half_dt);
}
@@ -109,17 +89,6 @@ __attribute__((always_inline)) INLINE static void kick_part(
a_tot[2] += p->gpart->a_grav[2];
}
-#ifdef ISOTHERMAL_GLASS
- /* TT: add viscous drag */
- const float surface_density= 10.;
- const float scale_height = 100;
- const float t_dyn = sqrt(scale_height / (const_G * surface_density));
- const double time = ti_start * timeBase;
- a_tot[0] -= p->gpart->v_full[0] / (0.05 * t_dyn);
- a_tot[1] -= p->gpart->v_full[1] / (0.05 * t_dyn);
- a_tot[2] -= p->gpart->v_full[2] / (0.05 * t_dyn);
-#endif
-
/* Kick particles in momentum space */
xp->v_full[0] += a_tot[0] * dt;
xp->v_full[1] += a_tot[1] * dt;
diff --git a/src/potentials.h b/src/potentials.h
index cc49545e95b9a8859277028849bd86b73ce03391..2a9999c72068295f8af93deb6a5bddadd2f92f98 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -130,21 +130,48 @@ 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 struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+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) {
- const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
+ 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 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); */
- /* impose *no* graitational acceleration */
- const float z_accel = 0.;
+#ifdef ISOTHERMAL_GLASS
+ 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
+ * tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
if (dz > 0)
- g->a_grav[2] -= z_accel;
+ g->a_grav[2] -= z_accel / G_newton; /* returned acceleraton is multiplied by G later on */
if (dz < 0)
- g->a_grav[2] += z_accel;
+ g->a_grav[2] += z_accel / G_newton;
+
+
+ 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
+ /* 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;
+
+#endif
+
+
}
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
diff --git a/src/runner.c b/src/runner.c
index 6e99d8cece2cf086b1209ac2f8314301b9779e20..7212f2baeb02b3ac00dfa10c1e7b78617fb4d9dc 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -105,7 +105,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
const int ti_current = r->e->ti_current;
const struct external_potential *potential = r->e->external_potential;
const struct phys_const *constants = r->e->physical_constants;
-
+ const double time = r->e->time;
TIMER_TIC;
/* Recurse? */
@@ -128,7 +128,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
/* Is this part within the time step? */
if (g->ti_end <= ti_current) {
- external_gravity(potential, constants, g);
+ external_gravity(time, potential, constants, g);
}
}
diff --git a/src/timestep.h b/src/timestep.h
index 569120cf9cf989b633da35529f7693c8f62a1910..fc1869156fc356f51b924379b5f4666233ab0491 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -111,6 +111,13 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
const float new_dt_self = FLT_MAX; // MATTHIEU
new_dt_grav = fminf(new_dt_external, new_dt_self);
+
+ if(p->id == -1)
+ {
+ printParticle_single(p, xp);
+ message(" dt_hydro= %e, dt_grav_external= %e dt_grav_self=%e ",new_dt_hydro,new_dt_external,new_dt_self);
+ }
+
}
/* Final time-step is minimum of hydro and gravity */
@@ -123,6 +130,9 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
: FLT_MAX;
new_dt = fminf(new_dt, dt_h_change);
+ if(p->id == -1)
+ message(" new_dt= %e", new_dt);
+
/* Limit timestep within the allowed range */
new_dt = fminf(new_dt, e->dt_max);