[yt-svn] commit/yt: ngoldbaum: Merged in Astrodude87/yt (pull request #2346)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Sep 7 11:17:28 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/146a5e02f7fe/
Changeset: 146a5e02f7fe
Branch: yt
User: ngoldbaum
Date: 2016-09-07 18:17:01+00:00
Summary: Merged in Astrodude87/yt (pull request #2346)
Fixing Ramses time values
Affected #: 7 files
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 .hgignore
--- a/.hgignore
+++ b/.hgignore
@@ -63,6 +63,7 @@
yt/utilities/lib/quad_tree.c
yt/utilities/lib/ray_integrators.c
yt/utilities/lib/ragged_arrays.c
+yt/utilities/lib/cosmology_time.c
yt/utilities/lib/grid_traversal.c
yt/utilities/lib/marching_cubes.c
yt/utilities/lib/png_writer.h
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 setup.py
--- a/setup.py
+++ b/setup.py
@@ -155,6 +155,8 @@
Extension("yt.utilities.lib.primitives",
["yt/utilities/lib/primitives.pyx"],
libraries=std_libs),
+ Extension("yt.utilities.lib.cosmology_time",
+ ["yt/utilities/lib/cosmology_time.pyx"]),
Extension("yt.utilities.lib.origami",
["yt/utilities/lib/origami.pyx",
"yt/utilities/lib/origami_tags.c"],
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 yt/fields/particle_fields.py
--- a/yt/fields/particle_fields.py
+++ b/yt/fields/particle_fields.py
@@ -175,6 +175,14 @@
function=function, units=unit_system["velocity"], take_log=False,
validators=[ValidateSpatial(0)])
+ for method, name in zip(("cic", "sum"), ("cic", "nn")):
+ function = _get_density_weighted_deposit_field(
+ "age", "s", method)
+ registry.add_field(
+ ("deposit", ("%s_"+name+"_age") % (ptype)),
+ function=function, units=unit_system["time"], take_log=False,
+ validators=[ValidateSpatial(0)])
+
# Now some translation functions.
def particle_ones(field, data):
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 yt/frontends/ramses/data_structures.py
--- a/yt/frontends/ramses/data_structures.py
+++ b/yt/frontends/ramses/data_structures.py
@@ -43,6 +43,9 @@
RAMSESOctreeContainer
from yt.arraytypes import blankRecordArray
+from yt.utilities.lib.cosmology_time import \
+ friedman
+
class RAMSESDomainFile(object):
_last_mask = None
_last_selector_id = None
@@ -620,7 +623,6 @@
dom, mi, ma = f.readline().split()
self.hilbert_indices[int(dom)] = (float(mi), float(ma))
self.parameters.update(rheader)
- self.current_time = self.parameters['time'] * self.parameters['unit_t']
self.domain_left_edge = np.zeros(3, dtype='float64')
self.domain_dimensions = np.ones(3, dtype='int32') * \
2**(self.min_level+1)
@@ -643,6 +645,23 @@
self.max_level = rheader['levelmax'] - self.min_level - 1
f.close()
+
+ if self.cosmological_simulation == 0:
+ self.current_time = self.parameters['time'] * self.parameters['unit_t']
+ else :
+ self.tau_frw, self.t_frw, self.dtau, self.n_frw, self.time_tot = \
+ friedman( self.omega_matter, self.omega_lambda, 1. - self.omega_matter - self.omega_lambda )
+
+ age = self.parameters['time']
+ iage = 1 + int(10.*age/self.dtau)
+ iage = np.min([iage,self.n_frw/2 + (iage - self.n_frw/2)/10])
+
+ self.time_simu = self.t_frw[iage ]*(age-self.tau_frw[iage-1])/(self.tau_frw[iage]-self.tau_frw[iage-1])+ \
+ self.t_frw[iage-1]*(age-self.tau_frw[iage ])/(self.tau_frw[iage-1]-self.tau_frw[iage])
+
+ self.current_time = (self.time_tot + self.time_simu)/(self.hubble_constant*1e7/3.08e24)/self.parameters['unit_t']
+
+
@classmethod
def _is_valid(self, *args, **kwargs):
if not os.path.basename(args[0]).startswith("info_"): return False
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 yt/frontends/ramses/fields.py
--- a/yt/frontends/ramses/fields.py
+++ b/yt/frontends/ramses/fields.py
@@ -80,7 +80,7 @@
("particle_mass", ("code_mass", [], None)),
("particle_identifier", ("", ["particle_index"], None)),
("particle_refinement_level", ("", [], None)),
- ("particle_age", ("code_time", [], None)),
+ ("particle_age", ("code_time", ['age'], None)),
("particle_metallicity", ("", [], None)),
)
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 yt/frontends/ramses/io.py
--- a/yt/frontends/ramses/io.py
+++ b/yt/frontends/ramses/io.py
@@ -19,7 +19,10 @@
from yt.utilities.io_handler import \
BaseIOHandler
from yt.utilities.logger import ytLogger as mylog
+from yt.utilities.physical_ratios import cm_per_km, cm_per_mpc
import yt.utilities.fortran_utils as fpu
+from yt.utilities.lib.cosmology_time import \
+ get_ramses_ages
from yt.extern.six import PY3
if PY3:
@@ -101,4 +104,15 @@
tr[field] = fpu.read_vector(f, dt)
if field[1].startswith("particle_position"):
np.divide(tr[field], subset.domain.ds["boxlen"], tr[field])
+ cosmo = subset.domain.ds.cosmological_simulation
+ if cosmo == 1 and field[1] == "particle_age":
+ tf = subset.domain.ds.t_frw
+ dtau = subset.domain.ds.dtau
+ tauf = subset.domain.ds.tau_frw
+ tsim = subset.domain.ds.time_simu
+ h100 = subset.domain.ds.hubble_constant
+ nOver2 = subset.domain.ds.n_frw/2
+ t_scale = 1./(h100 * 100 * cm_per_km / cm_per_mpc)/subset.domain.ds['unit_t']
+ ages = tr[field]
+ tr[field] = get_ramses_ages(tf,tauf,dtau,tsim,t_scale,ages,nOver2,len(ages))
return tr
diff -r c144003a5fd281a25c7d283899141cca10cb340b -r 146a5e02f7fe72bf4496b89871244c5f795a2ff6 yt/utilities/lib/cosmology_time.pyx
--- /dev/null
+++ b/yt/utilities/lib/cosmology_time.pyx
@@ -0,0 +1,100 @@
+cimport numpy as np
+import numpy as np
+
+
+cdef double dadtau(double aexp_tau,double O_mat_0,double O_vac_0,double O_k_0):
+ return ( aexp_tau**3 * (O_mat_0 + O_vac_0*aexp_tau**3 + O_k_0*aexp_tau) )**0.5
+
+cdef double dadt(double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
+ return ( (1./aexp_t)*(O_mat_0 + O_vac_0*aexp_t**3 + O_k_0*aexp_t) )**0.5
+
+
+cdef step_cosmo(double alpha,double tau,double aexp_tau,double t,double aexp_t,double O_mat_0,double O_vac_0,double O_k_0):
+ cdef double dtau,aexp_tau_pre,dt,aexp_t_pre
+
+ dtau = alpha * aexp_tau / dadtau(aexp_tau,O_mat_0,O_vac_0,O_k_0)
+ aexp_tau_pre = aexp_tau - dadtau(aexp_tau,O_mat_0,O_vac_0,O_k_0)*dtau/2.0
+ aexp_tau = aexp_tau - dadtau(aexp_tau_pre,O_mat_0,O_vac_0,O_k_0)*dtau
+ tau = tau - dtau
+
+ dt = alpha * aexp_t / dadt(aexp_t,O_mat_0,O_vac_0,O_k_0)
+ aexp_t_pre = aexp_t - dadt(aexp_t,O_mat_0,O_vac_0,O_k_0)*dt/2.0
+ aexp_t = aexp_t - dadt(aexp_t_pre,O_mat_0,O_vac_0,O_k_0)*dt
+ t = t - dt
+
+ return tau,aexp_tau,t,aexp_t
+
+
+cpdef friedman(double O_mat_0,double O_vac_0,double O_k_0):
+ cdef double alpha=1.e-5,aexp_min=1.e-3,aexp_tau=1.,aexp_t=1.,tau=0.,t=0.
+ cdef int nstep=0,ntable=1000,n_out
+ cdef np.ndarray[double,mode='c'] t_out=np.zeros([ntable+1]),tau_out=np.zeros([ntable+1])
+ cdef double age_tot,delta_tau,next_tau
+
+ while aexp_tau >= aexp_min or aexp_t >= aexp_min:
+ nstep = nstep + 1
+ tau,aexp_tau,t,aexp_t = step_cosmo(alpha,tau,aexp_tau,t,aexp_t,O_mat_0,O_vac_0,O_k_0)
+
+ age_tot=-t
+ if nstep < ntable :
+ ntable = nstep
+ alpha = alpha / 2.
+
+ delta_tau = 20.*tau/ntable/11.
+
+ aexp_tau = 1.
+ aexp_t = 1.
+ tau = 0.
+ t = 0.
+
+ n_out = 0
+ t_out[n_out] = t
+ tau_out[n_out] = tau
+
+ next_tau = tau + delta_tau/10.
+
+ while n_out < ntable/2 :
+ tau,aexp_tau,t,aexp_t = step_cosmo(alpha,tau,aexp_tau,t,aexp_t,O_mat_0,O_vac_0,O_k_0)
+
+ if tau < next_tau:
+ n_out = n_out + 1
+ t_out[n_out] = t
+ tau_out[n_out] = tau
+ next_tau = next_tau + delta_tau/10.
+
+ while aexp_tau >= aexp_min or aexp_t >= aexp_min:
+ tau,aexp_tau,t,aexp_t = step_cosmo(alpha,tau,aexp_tau,t,aexp_t,O_mat_0,O_vac_0,O_k_0)
+
+ if tau < next_tau:
+ n_out = n_out + 1
+ t_out[n_out] = t
+ tau_out[n_out] = tau
+ next_tau = next_tau + delta_tau
+
+ n_out = ntable
+ t_out[n_out] = t
+ tau_out[n_out] = tau
+
+ return tau_out,t_out,delta_tau,ntable,age_tot
+
+cpdef get_ramses_ages(np.ndarray[double,mode='c'] tf,
+ np.ndarray[double,mode='c'] tauf,
+ double dtau,
+ double tsim,
+ double t_scale,
+ np.ndarray[double,mode='c'] ages,
+ int nOver2,
+ int ntot):
+
+ cdef np.ndarray[double,mode='c'] t
+ cdef np.ndarray[double,mode='c'] dage
+ cdef np.ndarray[int,mode='c'] iage
+
+ dage = 1 + (10*ages/dtau)
+ dage = np.minimum(dage, nOver2 + (dage - nOver2)/10.)
+ iage = np.array(dage,dtype=np.int32)
+
+ t = (tf[iage]*(ages - tauf[iage - 1]) / (tauf[iage] - tauf[iage - 1]))
+ t = t + (tf[iage-1]*(ages-tauf[iage]) / (tauf[iage-1]-tauf[iage]))
+ return (tsim - t)*t_scale
+
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list