added iso_thermal flag, vicous force, updated make_icc.py

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.

examples/Disk-Patch/HydroStatic/dynamic.pro

+;
+;  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
+
+

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

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)

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

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
 }
 

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 */

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;

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 */
 

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);
     }
   }
 

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);