[yt-svn] commit/yt: drudd: Merged in drudd/yt (pull request #1172)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Sep 2 15:22:06 PDT 2014
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/5130e0892eb0/
Changeset: 5130e0892eb0
Branch: yt
User: drudd
Date: 2014-09-03 00:21:55
Summary: Merged in drudd/yt (pull request #1172)
Two approvers will have to do. ARTIO cosmological conversion improvements
Affected #: 5 files
diff -r 66661ac7b4208e71cd9b3ce97f36c1bdc06f5deb -r 5130e0892eb0cda34d5bff8a94e424d737e1518d yt/frontends/artio/_artio_caller.pyx
--- a/yt/frontends/artio/_artio_caller.pyx
+++ b/yt/frontends/artio/_artio_caller.pyx
@@ -25,6 +25,33 @@
cdef extern from "alloca.h":
void *alloca(int)
+cdef extern from "cosmology.h":
+ ctypedef struct CosmologyParameters "CosmologyParameters" :
+ pass
+ CosmologyParameters *cosmology_allocate()
+ void cosmology_free(CosmologyParameters *c)
+ void cosmology_set_fixed(CosmologyParameters *c)
+
+ void cosmology_set_OmegaM(CosmologyParameters *c, double value)
+ void cosmology_set_OmegaB(CosmologyParameters *c, double value)
+ void cosmology_set_OmegaL(CosmologyParameters *c, double value)
+ void cosmology_set_h(CosmologyParameters *c, double value)
+ void cosmology_set_DeltaDC(CosmologyParameters *c, double value)
+
+ double abox_from_auni(CosmologyParameters *c, double a) nogil
+ double tcode_from_auni(CosmologyParameters *c, double a) nogil
+ double tphys_from_auni(CosmologyParameters *c, double a) nogil
+
+ double auni_from_abox(CosmologyParameters *c, double v) nogil
+ double auni_from_tcode(CosmologyParameters *c, double v) nogil
+ double auni_from_tphys(CosmologyParameters *c, double v) nogil
+
+ double abox_from_tcode(CosmologyParameters *c, double tcode) nogil
+ double tcode_from_abox(CosmologyParameters *c, double abox) nogil
+
+ double tphys_from_abox(CosmologyParameters *c, double abox) nogil
+ double tphys_from_tcode(CosmologyParameters *c, double tcode) nogil
+
cdef extern from "artio.h":
ctypedef struct artio_fileset_handle "artio_fileset" :
pass
@@ -138,6 +165,10 @@
cdef public object parameters
cdef artio_fileset_handle *handle
+ # cosmology parameter for time unit conversion
+ cdef CosmologyParameters *cosmology
+ cdef float tcode_to_years
+
# common attributes
cdef public int num_grid
cdef int64_t num_root_cells
@@ -178,6 +209,19 @@
self.sfc_min = 0
self.sfc_max = self.num_root_cells-1
+ # initialize cosmology module
+ if self.parameters.has_key("abox"):
+ self.cosmology = cosmology_allocate()
+ cosmology_set_OmegaM(self.cosmology, self.parameters['OmegaM'][0])
+ cosmology_set_OmegaL(self.cosmology, self.parameters['OmegaL'][0])
+ cosmology_set_OmegaB(self.cosmology, self.parameters['OmegaB'][0])
+ cosmology_set_h(self.cosmology, self.parameters['hubble'][0])
+ cosmology_set_DeltaDC(self.cosmology, self.parameters['DeltaDC'][0])
+ cosmology_set_fixed(self.cosmology)
+ else:
+ self.cosmology = NULL
+ self.tcode_to_years = self.parameters['time_unit'][0]/(365.25*86400)
+
# grid detection
self.min_level = 0
self.max_level = self.parameters['grid_max_level'][0]
@@ -239,6 +283,8 @@
if self.primary_variables : free(self.primary_variables)
if self.secondary_variables : free(self.secondary_variables)
+ if self.cosmology : cosmology_free(self.cosmology)
+
if self.handle : artio_fileset_close(self.handle)
def read_parameters(self) :
@@ -288,6 +334,100 @@
self.parameters[key] = parameter
+ def abox_from_auni(self, np.float64_t a):
+ if self.cosmology:
+ return abox_from_auni(self.cosmology, a)
+ else:
+ raise RuntimeError("abox_from_auni called for non-cosmological ARTIO fileset!")
+
+ def tcode_from_auni(self, np.float64_t a):
+ if self.cosmology:
+ return tcode_from_auni(self.cosmology, a)
+ else:
+ raise RuntimeError("tcode_from_auni called for non-cosmological ARTIO fileset!")
+
+ def tphys_from_auni(self, np.float64_t a):
+ if self.cosmology:
+ return tphys_from_auni(self.cosmology, a)
+ else:
+ raise RuntimeError("tphys_from_auni called for non-cosmological ARTIO fileset!")
+
+ def auni_from_abox(self, np.float64_t v):
+ if self.cosmology:
+ return auni_from_abox(self.cosmology, v)
+ else:
+ raise RuntimeError("auni_from_abox called for non-cosmological ARTIO fileset!")
+
+ def auni_from_tcode(self, np.float64_t v):
+ if self.cosmology:
+ return auni_from_tcode(self.cosmology, v)
+ else:
+ raise RuntimeError("auni_from_tcode called for non-cosmological ARTIO fileset!")
+
+ @cython.wraparound(False)
+ @cython.boundscheck(False)
+ def auni_from_tcode_array(self, np.ndarray[np.float64_t] tcode):
+ cdef int i, nauni
+ cdef np.ndarray[np.float64_t, ndim=1] auni
+ cdef CosmologyParameters *cosmology = self.cosmology
+ if not cosmology:
+ raise RuntimeError("auni_from_tcode_array called for non-cosmological ARTIO fileset!")
+ auni = np.empty_like(tcode)
+ nauni = auni.shape[0]
+ with nogil:
+ for i in range(nauni):
+ auni[i] = auni_from_tcode(self.cosmology, tcode[i])
+ return auni
+
+ def auni_from_tphys(self, np.float64_t v):
+ if self.cosmology:
+ return auni_from_tphys(self.cosmology, v)
+ else:
+ raise RuntimeError("auni_from_tphys called for non-cosmological ARTIO fileset!")
+
+ def abox_from_tcode(self, np.float64_t abox):
+ if self.cosmology:
+ return abox_from_tcode(self.cosmology, abox)
+ else:
+ raise RuntimeError("abox_from_tcode called for non-cosmological ARTIO fileset!")
+
+ def tcode_from_abox(self, np.float64_t abox):
+ if self.cosmology:
+ return tcode_from_abox(self.cosmology, abox)
+ else:
+ raise RuntimeError("tcode_from_abox called for non-cosmological ARTIO fileset!")
+
+ def tphys_from_abox(self, np.float64_t abox):
+ if self.cosmology:
+ return tphys_from_abox(self.cosmology, abox)
+ else:
+ raise RuntimeError("tphys_from_abox called for non-cosmological ARTIO fileset!")
+
+ def tphys_from_tcode(self, np.float64_t tcode):
+ if self.cosmology:
+ return tphys_from_tcode(self.cosmology, tcode)
+ else:
+ return self.tcode_to_years*tcode
+
+ @cython.wraparound(False)
+ @cython.boundscheck(False)
+ def tphys_from_tcode_array(self, np.ndarray[np.float64_t] tcode):
+ cdef int i, ntphys
+ cdef np.ndarray[np.float64_t, ndim=1] tphys
+ cdef CosmologyParameters *cosmology = self.cosmology
+ tphys = np.empty_like(tcode)
+ ntphys = tcode.shape[0]
+
+ if cosmology:
+ tphys = np.empty_like(tcode)
+ ntphys = tcode.shape[0]
+ with nogil:
+ for i in range(ntphys):
+ tphys[i] = tphys_from_tcode(cosmology, tcode[i])
+ return tphys
+ else:
+ return tcode*self.tcode_to_years
+
# @cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
diff -r 66661ac7b4208e71cd9b3ce97f36c1bdc06f5deb -r 5130e0892eb0cda34d5bff8a94e424d737e1518d yt/frontends/artio/artio_headers/cosmology.c
--- /dev/null
+++ b/yt/frontends/artio/artio_headers/cosmology.c
@@ -0,0 +1,524 @@
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+
+#ifndef ERROR
+#include <stdio.h>
+#define ERROR(msg) { fprintf(stderr,"%s\n",msg); exit(1); }
+#endif
+
+#ifndef ASSERT
+#include <stdio.h>
+#define ASSERT(exp) { if(!(exp)) { fprintf(stderr,"Failed assertion %s, line: %d\n",#exp,__LINE__); exit(1); } }
+#endif
+
+#ifndef HEAPALLOC
+#include <stdlib.h>
+#define HEAPALLOC(type,size) (type *)malloc((size)*sizeof(type))
+#endif
+
+#ifndef NEWARR
+#include <stdlib.h>
+#define NEWARR(size) HEAPALLOC(double,size)
+#endif
+
+#ifndef DELETE
+#include <stdlib.h>
+#define DELETE(ptr) free(ptr)
+#endif
+
+
+#include "cosmology.h"
+
+
+struct CosmologyParametersStruct
+{
+ int set;
+ int ndex;
+ int size;
+ double *la;
+ double *aUni;
+ double *aBox;
+ double *tCode;
+ double *tPhys;
+ double *dPlus;
+ double *qPlus;
+ double aLow;
+ double tCodeOffset;
+
+ double OmegaM;
+ double OmegaD;
+ double OmegaB;
+ double OmegaL;
+ double OmegaK;
+ double OmegaR;
+ double h;
+ double DeltaDC;
+ int flat;
+ double Omh2;
+ double Obh2;
+};
+
+void cosmology_clear_table(CosmologyParameters *c);
+void cosmology_fill_table(CosmologyParameters *c, double amin, double amax);
+void cosmology_fill_table_abox(CosmologyParameters *c, int istart, int n);
+
+CosmologyParameters *cosmology_allocate() {
+ CosmologyParameters *c = HEAPALLOC(CosmologyParameters,1);
+ if ( c != NULL ) {
+ memset(c, 0, sizeof(CosmologyParameters));
+
+ c->ndex = 200;
+ c->aLow = 1.0e-2;
+ }
+ return c;
+}
+
+void cosmology_free(CosmologyParameters *c) {
+ cosmology_clear_table(c);
+ DELETE(c);
+}
+
+int cosmology_is_set(CosmologyParameters *c)
+{
+ return (c->OmegaM>0.0 && c->OmegaB>0.0 && c->h>0.0);
+}
+
+
+void cosmology_fail_on_reset(const char *name, double old_value, double new_value)
+{
+ char str[150];
+ sprintf(str,"Trying to change %s from %lg to %lg...\nCosmology has been fixed and cannot be changed.\n",name,old_value,new_value);
+ ERROR(str);
+}
+
+
+void cosmology_set_OmegaM(CosmologyParameters *c, double v)
+{
+ if(v < 1.0e-3) v = 1.0e-3;
+ if(fabs(c->OmegaM-v) > 1.0e-5)
+ {
+ if(c->set) cosmology_fail_on_reset("OmegaM",c->OmegaM,v);
+ c->OmegaM = v;
+ c->flat = (fabs(c->OmegaM+c->OmegaL-1.0) > 1.0e-5) ? 0 : 1;
+ cosmology_clear_table(c);
+ }
+}
+
+
+void cosmology_set_OmegaL(CosmologyParameters *c, double v)
+{
+ if(fabs(c->OmegaL-v) > 1.0e-5)
+ {
+ if(c->set) cosmology_fail_on_reset("OmegaL",c->OmegaL,v);
+ c->OmegaL = v;
+ c->flat = (fabs(c->OmegaM+c->OmegaL-1.0) > 1.0e-5) ? 0 : 1;
+ cosmology_clear_table(c);
+ }
+}
+
+
+void cosmology_set_OmegaB(CosmologyParameters *c, double v)
+{
+ if(v < 0.0) v = 0.0;
+ if(fabs(c->OmegaB-v) > 1.0e-5)
+ {
+ if(c->set) cosmology_fail_on_reset("OmegaB",c->OmegaB,v);
+ c->OmegaB = v;
+ cosmology_clear_table(c);
+ }
+}
+
+
+void cosmology_set_h(CosmologyParameters *c, double v)
+{
+ if(fabs(c->h-v) > 1.0e-5)
+ {
+ if(c->set) cosmology_fail_on_reset("h",c->h,v);
+ c->h = v;
+ cosmology_clear_table(c);
+ }
+}
+
+
+void cosmology_set_DeltaDC(CosmologyParameters *c, double v)
+{
+ if(fabs(c->DeltaDC-v) > 1.0e-3)
+ {
+ if(c->set) cosmology_fail_on_reset("DeltaDC",c->DeltaDC,v);
+ c->DeltaDC = v;
+ cosmology_clear_table(c);
+ }
+}
+
+
+void cosmology_init(CosmologyParameters *c)
+{
+ if(c->size == 0) /* reset only if the state is dirty */
+ {
+ if(!cosmology_is_set(c)) ERROR("Not all of the required cosmological parameters have been set; the minimum required set is (OmegaM,OmegaB,h).");
+
+ if(c->OmegaB > c->OmegaM) c->OmegaB = c->OmegaM;
+ c->OmegaD = c->OmegaM - c->OmegaB;
+ if(c->flat)
+ {
+ c->OmegaK = 0.0;
+ c->OmegaL = 1.0 - c->OmegaM;
+ }
+ else
+ {
+ c->OmegaK = 1.0 - (c->OmegaM+c->OmegaL);
+ }
+ c->OmegaR = 4.166e-5/(c->h*c->h);
+
+ c->Omh2 = c->OmegaM*c->h*c->h;
+ c->Obh2 = c->OmegaB*c->h*c->h;
+
+ cosmology_fill_table(c,c->aLow,1.0);
+
+ c->tCodeOffset = 0.0; /* Do need to set it to zero first */
+#ifndef NATIVE_TCODE_NORMALIZATION
+ c->tCodeOffset = 0.0 - tCode(c,inv_aBox(c,1.0));
+#endif
+ }
+}
+
+
+void cosmology_set_fixed(CosmologyParameters *c)
+{
+ cosmology_init(c);
+ c->set = 1;
+}
+
+
+double cosmology_mu(CosmologyParameters *c, double a)
+{
+ return sqrt(((a*a*c->OmegaL+c->OmegaK)*a+c->OmegaM)*a+c->OmegaR);
+}
+
+
+double cosmology_dc_factor(CosmologyParameters *c, double dPlus)
+{
+ double dc = 1.0 + dPlus*c->DeltaDC;
+ return 1.0/pow((dc>0.001)?dc:0.001,1.0/3.0);
+}
+
+
+void cosmology_fill_table_integrate(CosmologyParameters *c, double a, double y[], double f[])
+{
+ double mu = cosmology_mu(c, a);
+ double abox = a*cosmology_dc_factor(c, y[2]);
+
+ f[0] = a/(abox*abox*mu);
+ f[1] = a/mu;
+ f[2] = y[3]/(a*mu);
+ f[3] = 1.5*c->OmegaM*y[2]/mu;
+}
+
+
+void cosmology_fill_table_piece(CosmologyParameters *c, int istart, int n)
+{
+ int i, j;
+ double tPhysUnit = (3.0856775813e17/(365.25*86400))/c->h; /* 1/H0 in Julian years */
+
+ double x, aeq = c->OmegaR/c->OmegaM;
+ double tCodeFac = 1.0/sqrt(aeq);
+ double tPhysFac = tPhysUnit*aeq*sqrt(aeq)/sqrt(c->OmegaM);
+
+ double da, a0, y0[4], y1[4];
+ double f1[4], f2[4], f3[4], f4[4];
+
+ for(i=istart; i<n; i++)
+ {
+ c->aUni[i] = pow(10.0,c->la[i]);
+ }
+
+ /*
+ // Small a regime, use analytical formulae for matter + radiation model
+ */
+ for(i=istart; c->aUni[i]<(c->aLow+1.0e-9) && i<n; i++)
+ {
+ x = c->aUni[i]/aeq;
+
+ c->tPhys[i] = tPhysFac*2*x*x*(2+sqrt(x+1))/(3*pow(1+sqrt(x+1),2.0));
+ c->dPlus[i] = aeq*(x + 2.0/3.0 + (6*sqrt(1+x)+(2+3*x)*log(x)-2*(2+3*x)*log(1+sqrt(1+x)))/(log(64.0)-9)); /* long last term is the decaying mode generated after euality; it is very small for x > 10, I keep ot just for completeness; */
+ c->qPlus[i] = c->aUni[i]*cosmology_mu(c,c->aUni[i])*(1 + ((2+6*x)/(x*sqrt(1+x))+3*log(x)-6*log(1+sqrt(1+x)))/(log(64)-9)); /* this is a^2*dDPlus/dt/H0 */
+
+ c->aBox[i] = c->aUni[i]*cosmology_dc_factor(c,c->dPlus[i]);
+ c->tCode[i] = 1.0 - tCodeFac*asinh(sqrt(aeq/c->aBox[i]));
+ }
+
+ /*
+ // Large a regime, solve ODEs
+ */
+ ASSERT(i > 0);
+
+ tCodeFac = 0.5*sqrt(c->OmegaM);
+ tPhysFac = tPhysUnit;
+
+ y1[0] = c->tCode[i-1]/tCodeFac;
+ y1[1] = c->tPhys[i-1]/tPhysFac;
+ y1[2] = c->dPlus[i-1];
+ y1[3] = c->qPlus[i-1];
+
+ for(; i<n; i++)
+ {
+ a0 = c->aUni[i-1];
+ da = c->aUni[i] - a0;
+
+ /* RK4 integration */
+ for(j=0; j<4; j++) y0[j] = y1[j];
+ cosmology_fill_table_integrate(c, a0,y1,f1);
+
+ for(j=0; j<4; j++) y1[j] = y0[j] + 0.5*da*f1[j];
+ cosmology_fill_table_integrate(c, a0+0.5*da,y1,f2);
+
+ for(j=0; j<4; j++) y1[j] = y0[j] + 0.5*da*f2[j];
+ cosmology_fill_table_integrate(c, a0+0.5*da,y1,f3);
+
+ for(j=0; j<4; j++) y1[j] = y0[j] + da*f3[j];
+ cosmology_fill_table_integrate(c, a0+da,y1,f4);
+
+ for(j=0; j<4; j++) y1[j] = y0[j] + da*(f1[j]+2*f2[j]+2*f3[j]+f4[j])/6.0;
+
+ c->tCode[i] = tCodeFac*y1[0];
+ c->tPhys[i] = tPhysFac*y1[1];
+ c->dPlus[i] = y1[2];
+ c->qPlus[i] = y1[3];
+
+ c->aBox[i] = c->aUni[i]*cosmology_dc_factor(c,c->dPlus[i]);
+ }
+}
+
+
+void cosmology_fill_table(CosmologyParameters *c, double amin, double amax)
+{
+ int i, imin, imax, iold;
+ double dla = 1.0/c->ndex;
+ double lamin, lamax;
+ double *old_la = c->la;
+ double *old_aUni = c->aUni;
+ double *old_aBox = c->aBox;
+ double *old_tCode = c->tCode;
+ double *old_tPhys = c->tPhys;
+ double *old_dPlus = c->dPlus;
+ double *old_qPlus = c->qPlus;
+ int old_size = c->size;
+
+ if(amin > c->aLow) amin = c->aLow;
+ lamin = dla*floor(c->ndex*log10(amin));
+ lamax = dla*ceil(c->ndex*log10(amax));
+
+ c->size = 1 + (int)(0.5+c->ndex*(lamax-lamin));
+ ASSERT(fabs(lamax-lamin-dla*(c->size-1)) < 1.0e-14);
+
+ c->la = NEWARR(c->size); ASSERT(c->la != NULL);
+ c->aUni = NEWARR(c->size); ASSERT(c->aUni != NULL);
+ c->aBox = NEWARR(c->size); ASSERT(c->aBox != NULL);
+ c->tCode = NEWARR(c->size); ASSERT(c->tCode != NULL);
+ c->tPhys = NEWARR(c->size); ASSERT(c->tPhys != NULL);
+ c->dPlus = NEWARR(c->size); ASSERT(c->dPlus != NULL);
+ c->qPlus = NEWARR(c->size); ASSERT(c->qPlus != NULL);
+
+ /*
+ // New log10(aUni) table
+ */
+ for(i=0; i<c->size; i++)
+ {
+ c->la[i] = lamin + dla*i;
+ }
+
+ if(old_size == 0)
+ {
+ /*
+ // Filling the table for the first time
+ */
+ cosmology_fill_table_piece(c,0,c->size);
+ }
+ else
+ {
+ /*
+ // Find if we need to expand the lower end
+ */
+ if(lamin < old_la[0])
+ {
+ imin = (int)(0.5+c->ndex*(old_la[0]-lamin));
+ ASSERT(fabs(old_la[0]-lamin-dla*imin) < 1.0e-14);
+ }
+ else imin = 0;
+
+ /*
+ // Find if we need to expand the upper end
+ */
+ if(lamax > old_la[old_size-1])
+ {
+ imax = (int)(0.5+c->ndex*(old_la[old_size-1]-lamin));
+ ASSERT(fabs(old_la[old_size-1]-lamin-dla*imax) < 1.0e-14);
+ }
+ else imax = c->size - 1;
+
+ /*
+ // Re-use the rest
+ */
+ if(lamin > old_la[0])
+ {
+ iold = (int)(0.5+c->ndex*(lamin-old_la[0]));
+ ASSERT(fabs(lamin-old_la[0]-dla*iold) < 1.0e-14);
+ }
+ else iold = 0;
+
+ memcpy(c->aUni+imin,old_aUni+iold,sizeof(double)*(imax-imin+1));
+ memcpy(c->aBox+imin,old_aBox+iold,sizeof(double)*(imax-imin+1));
+ memcpy(c->tCode+imin,old_tCode+iold,sizeof(double)*(imax-imin+1));
+ memcpy(c->tPhys+imin,old_tPhys+iold,sizeof(double)*(imax-imin+1));
+ memcpy(c->dPlus+imin,old_dPlus+iold,sizeof(double)*(imax-imin+1));
+ memcpy(c->qPlus+imin,old_qPlus+iold,sizeof(double)*(imax-imin+1));
+
+ DELETE(old_la);
+ DELETE(old_aUni);
+ DELETE(old_aBox);
+ DELETE(old_tCode);
+ DELETE(old_tPhys);
+ DELETE(old_dPlus);
+ DELETE(old_qPlus);
+
+ /*
+ // Fill in additional pieces
+ */
+ if(imin > 0) cosmology_fill_table_piece(c,0,imin);
+ if(imax < c->size-1) cosmology_fill_table_piece(c,imax,c->size);
+ }
+}
+
+
+void cosmology_clear_table(CosmologyParameters *c)
+{
+ if(c->size > 0)
+ {
+ DELETE(c->la);
+ DELETE(c->aUni);
+ DELETE(c->aBox);
+ DELETE(c->tCode);
+ DELETE(c->tPhys);
+ DELETE(c->dPlus);
+ DELETE(c->qPlus);
+
+ c->size = 0;
+ c->la = NULL;
+ c->aUni = NULL;
+ c->aBox = NULL;
+ c->tCode = NULL;
+ c->tPhys = NULL;
+ c->dPlus = NULL;
+ c->qPlus = NULL;
+ }
+}
+
+
+void cosmology_check_range(CosmologyParameters *c, double a)
+{
+ ASSERT((a > 1.0e-9) && (a < 1.0e9));
+
+ if(c->size == 0) cosmology_init(c);
+
+ if(a < c->aUni[0])
+ {
+ cosmology_fill_table(c,a,c->aUni[c->size-1]);
+ }
+
+ if(a > c->aUni[c->size-1])
+ {
+ cosmology_fill_table(c,c->aUni[0],a);
+ }
+}
+
+
+void cosmology_set_thread_safe_range(CosmologyParameters *c, double amin, double amax)
+{
+ cosmology_check_range(c, amin);
+ cosmology_check_range(c, amax);
+}
+
+
+double cosmology_get_value_from_table(CosmologyParameters *c, double a, double table[])
+{
+ int idx = (int)(c->ndex*(log10(a)-c->la[0]));
+
+ ASSERT(idx>=0 && idx<c->size);
+
+ /*
+ // Do it as a function of aUni rather than la to ensure exact inversion
+ */
+ return table[idx] + (table[idx+1]-table[idx])/(c->aUni[idx+1]-c->aUni[idx])*(a-c->aUni[idx]);
+}
+
+
+int cosmology_find_index(CosmologyParameters *c, double v, double table[])
+{
+ int ic, il = 0;
+ int ih = c->size - 1;
+
+ if(v < table[0])
+ {
+ return -1;
+ }
+ if(v > table[c->size-1])
+ {
+ return c->size + 1;
+ }
+
+ while((ih-il) > 1)
+ {
+ ic = (il+ih)/2;
+ if(v > table[ic]) /* special, not fully optimal form to avoid checking that il < c->size-1 */
+ il = ic;
+ else
+ ih = ic;
+ }
+
+ ASSERT(il+1 < c->size);
+
+ return il;
+}
+
+
+/*
+// Direct and inverse functions
+*/
+#define DEFINE_FUN(name,offset) \
+double name(CosmologyParameters *c, double a) \
+{ \
+ cosmology_check_range(c,a); \
+ return cosmology_get_value_from_table(c,a,c->name) + offset; \
+} \
+double inv_##name(CosmologyParameters *c, double v) \
+{ \
+ int idx; \
+ double *table; \
+ if(c->size == 0) cosmology_init(c); \
+ v -= offset; \
+ table = c->name; \
+ idx = cosmology_find_index(c,v,table); \
+ while(idx < 0) \
+ { \
+ cosmology_check_range(c,0.5*c->aUni[0]); \
+ table = c->name; \
+ idx = cosmology_find_index(c,v,table); \
+ } \
+ while(idx > c->size) \
+ { \
+ cosmology_check_range(c,2.0*c->aUni[c->size-1]); \
+ table = c->name; \
+ idx = cosmology_find_index(c,v,table); \
+ } \
+ return c->aUni[idx] + (c->aUni[idx+1]-c->aUni[idx])/(table[idx+1]-table[idx])*(v-table[idx]); \
+}
+
+DEFINE_FUN(aBox,0.0);
+DEFINE_FUN(tCode,c->tCodeOffset);
+DEFINE_FUN(tPhys,0.0);
+DEFINE_FUN(dPlus,0.0);
+DEFINE_FUN(qPlus,0.0);
+
+#undef DEFINE_FUN
diff -r 66661ac7b4208e71cd9b3ce97f36c1bdc06f5deb -r 5130e0892eb0cda34d5bff8a94e424d737e1518d yt/frontends/artio/artio_headers/cosmology.h
--- /dev/null
+++ b/yt/frontends/artio/artio_headers/cosmology.h
@@ -0,0 +1,101 @@
+#ifndef __COSMOLOGY_H__
+#define __COSMOLOGY_H__
+
+typedef struct CosmologyParametersStruct CosmologyParameters;
+
+#define COSMOLOGY_DECLARE_PRIMARY_PARAMETER(name) \
+void cosmology_set_##name(CosmologyParameters *c, double value)
+
+#define cosmology_set(c,name,value) \
+cosmology_set_##name(c,value)
+
+COSMOLOGY_DECLARE_PRIMARY_PARAMETER(OmegaM);
+COSMOLOGY_DECLARE_PRIMARY_PARAMETER(OmegaB);
+COSMOLOGY_DECLARE_PRIMARY_PARAMETER(OmegaL);
+COSMOLOGY_DECLARE_PRIMARY_PARAMETER(h);
+COSMOLOGY_DECLARE_PRIMARY_PARAMETER(DeltaDC);
+
+#undef COSMOLOGY_DECLARE_PRIMARY_PARAMETER
+
+CosmologyParameters *cosmology_allocate();
+void cosmology_free(CosmologyParameters *c);
+
+/*
+// Check that all required cosmological parameters have been set.
+// The minimum set is OmegaM, OmegaB, and h. By default, zero OmegaL,
+// OmegaK, and the DC mode are assumed.
+*/
+int cosmology_is_set(CosmologyParameters *c);
+
+
+/*
+// Freeze the cosmology and forbid any further changes to it.
+// In codes that include user-customizable segments (like plugins),
+// this function van be used for insuring that a user does not
+// change the cosmology in mid-run.
+*/
+void cosmology_set_fixed(CosmologyParameters *c);
+
+
+/*
+// Manual initialization. This does not need to be called,
+// the initialization is done automatically on the first call
+// to a relevant function.
+*/
+void cosmology_init(CosmologyParameters *c);
+
+
+/*
+// Set the range of global scale factors for thread-safe
+// calls to direct functions until the argument leaves the range.
+*/
+void cosmology_set_thread_safe_range(CosmologyParameters *c, double amin, double amax);
+
+/*
+// Direct functions take the global cosmological scale factor as the argument.
+// These functionsare are thread-safe if called with the argument in the
+// range set by a prior call to cosmology_set_thread_safe_range(...).
+// Calling them with the argument outside that range is ok, but breaks
+// thread-safety assurance.
+*/
+
+#define DEFINE_FUN(name) \
+double name(CosmologyParameters *c, double a); \
+double inv_##name(CosmologyParameters *c, double v);
+
+DEFINE_FUN(aBox);
+DEFINE_FUN(tCode);
+DEFINE_FUN(tPhys);
+DEFINE_FUN(dPlus);
+DEFINE_FUN(qPlus); /* Q+ = a^2 dD+/(H0 dt) */
+
+#undef DEFINE_FUN
+
+/*
+// Conversion macros
+*/
+#define abox_from_auni(c,a) aBox(c,a)
+#define tcode_from_auni(c,a) tCode(c,a)
+#define tphys_from_auni(c,a) tPhys(c,a)
+#define dplus_from_auni(c,a) dPlus(c,a)
+
+#define auni_from_abox(c,v) inv_aBox(c,v)
+#define auni_from_tcode(c,v) inv_tCode(c,v)
+#define auni_from_tphys(c,v) inv_tPhys(c,v)
+#define auni_from_dplus(c,v) inv_dPlus(c,v)
+
+#define abox_from_tcode(c,tcode) aBox(c,inv_tCode(c,tcode))
+#define tcode_from_abox(c,abox) tCode(c,inv_aBox(c,abox))
+
+#define tphys_from_abox(c,abox) tPhys(c,inv_aBox(c,abox))
+#define tphys_from_tcode(c,tcode) tPhys(c,inv_tCode(c,tcode))
+#define dplus_from_tcode(c,tcode) dPlus(c,inv_tCode(c,tcode))
+
+/*
+// Hubble parameter in km/s/Mpc; defined as macro so that it can be
+// undefined if needed to avoid the name clash.
+*/
+double cosmology_mu(CosmologyParameters *c, double a);
+#define Hubble(c,a) (100*c->h*cosmology_mu(c,a)/(a*a))
+
+#endif /* __COSMOLOGY_H__ */
diff -r 66661ac7b4208e71cd9b3ce97f36c1bdc06f5deb -r 5130e0892eb0cda34d5bff8a94e424d737e1518d yt/frontends/artio/data_structures.py
--- a/yt/frontends/artio/data_structures.py
+++ b/yt/frontends/artio/data_structures.py
@@ -26,8 +26,7 @@
from yt.utilities.definitions import \
mpc_conversion, sec_conversion
from .fields import \
- ARTIOFieldInfo, \
- b2t
+ ARTIOFieldInfo
from yt.fields.particle_fields import \
standard_particle_fields
@@ -387,37 +386,40 @@
self.particle_types = ()
self.particle_types_raw = self.particle_types
- self.current_time = b2t(self.artio_parameters["tl"][0])
+ self.current_time = self.quan(self._handle.tphys_from_tcode(self.artio_parameters["tl"][0]),"yr")
# detect cosmology
if "abox" in self.artio_parameters:
+ self.cosmological_simulation = True
+
abox = self.artio_parameters["abox"][0]
- self.cosmological_simulation = True
self.omega_lambda = self.artio_parameters["OmegaL"][0]
self.omega_matter = self.artio_parameters["OmegaM"][0]
self.hubble_constant = self.artio_parameters["hubble"][0]
- self.current_redshift = 1.0/self.artio_parameters["abox"][0] - 1.0
+ self.current_redshift = 1.0/self.artio_parameters["auni"][0] - 1.0
+ self.current_redshift_box = 1.0/abox - 1.0
self.parameters["initial_redshift"] =\
1.0 / self.artio_parameters["auni_init"][0] - 1.0
self.parameters["CosmologyInitialRedshift"] =\
self.parameters["initial_redshift"]
- else:
- self.cosmological_simulation = False
- #units
- if self.cosmological_simulation:
self.parameters['unit_m'] = self.artio_parameters["mass_unit"][0]
self.parameters['unit_t'] =\
self.artio_parameters["time_unit"][0] * abox**2
self.parameters['unit_l'] =\
self.artio_parameters["length_unit"][0] * abox
+
+ if self.artio_parameters["DeltaDC"][0] != 0:
+ mylog.warn("DeltaDC != 0, which implies auni != abox. Be sure you understand which expansion parameter is appropriate for your use! (Gnedin, Kravtsov, & Rudd 2011)")
else:
+ self.cosmological_simulation = False
+
self.parameters['unit_l'] = self.artio_parameters["length_unit"][0]
self.parameters['unit_t'] = self.artio_parameters["time_unit"][0]
self.parameters['unit_m'] = self.artio_parameters["mass_unit"][0]
- # hard coded assumption of 3D periodicity (add to dataset)
+ # hard coded assumption of 3D periodicity
self.periodicity = (True, True, True)
@classmethod
diff -r 66661ac7b4208e71cd9b3ce97f36c1bdc06f5deb -r 5130e0892eb0cda34d5bff8a94e424d737e1518d yt/frontends/artio/fields.py
--- a/yt/frontends/artio/fields.py
+++ b/yt/frontends/artio/fields.py
@@ -19,6 +19,8 @@
from yt.funcs import mylog
from yt.fields.field_info_container import \
FieldInfoContainer
+from yt.fields.field_detector import \
+ FieldDetector
from yt.units.yt_array import \
YTArray
@@ -62,7 +64,8 @@
("MASS", ("code_mass", ["particle_mass"], None)),
("PID", ("", ["particle_index"], None)),
("SPECIES", ("", ["particle_type"], None)),
- ("BIRTH_TIME", ("code_time", ["creation_time"], None)),
+ ("BIRTH_TIME", ("", [], None)), # code-units defined as dimensionless to
+ # avoid incorrect conversion
("INITIAL_MASS", ("code_mass", ["initial_mass"], None)),
("METALLICITY_SNIa", ("", ["metallicity_snia"], None)),
("METALLICITY_SNII", ("", ["metallicity_snii"], None)),
@@ -111,123 +114,35 @@
take_log=True)
def setup_particle_fields(self, ptype):
+ if ptype == "STAR":
+ def _creation_time(field,data):
+ # this test is necessary to avoid passing invalid tcode values
+ # to the function tphys_from_tcode during field detection
+ # (1.0 is not a valid tcode value)
+ if isinstance(data, FieldDetector):
+ return data["STAR","BIRTH_TIME"]
+ return YTArray(data.ds._handle.tphys_from_tcode_array(data["STAR","BIRTH_TIME"]),"yr")
- def _particle_age(field, data):
- return b2t(data[ptype,"creation_time"])
- self.add_field((ptype, "particle_age"), function=_particle_age, units="s",
- particle_type=True)
+ def _age(field, data):
+ if isinstance(data, FieldDetector):
+ return data["STAR","creation_time"]
+ return data.ds.current_time - data["STAR","creation_time"]
+
+ self.add_field((ptype, "creation_time"), function=_creation_time, units="yr",
+ particle_type=True)
+ self.add_field((ptype, "age"), function=_age, units="yr",
+ particle_type=True)
+
+ if self.ds.cosmological_simulation:
+ def _creation_redshift(field,data):
+ # this test is necessary to avoid passing invalid tcode values
+ # to the function auni_from_tcode during field detection
+ # (1.0 is not a valid tcode value)
+ if isinstance(data, FieldDetector):
+ return data["STAR","BIRTH_TIME"]
+ return 1.0/data.ds._handle.auni_from_tcode_array(data["STAR","BIRTH_TIME"]) - 1.0
+
+ self.add_field((ptype, "creation_redshift"), function=_creation_redshift,
+ particle_type=True)
super(ARTIOFieldInfo, self).setup_particle_fields(ptype)
-
-#stolen from frontends/art/
-#All of these functions are to convert from hydro time var to
-#proper time
-sqrt = np.sqrt
-sign = np.sign
-
-
-def find_root(f, a, b, tol=1e-6):
- c = (a+b)/2.0
- last = -np.inf
- assert(sign(f(a)) != sign(f(b)))
- while np.abs(f(c)-last) > tol:
- last = f(c)
- if sign(last) == sign(f(b)):
- b = c
- else:
- a = c
- c = (a+b)/2.0
- return c
-
-
-def quad(fintegrand, xmin, xmax, n=1e4):
- spacings = np.logspace(np.log10(xmin), np.log10(xmax), n)
- integrand_arr = fintegrand(spacings)
- val = np.trapz(integrand_arr, dx=np.diff(spacings))
- return val
-
-
-def a2b(at, Om0=0.27, Oml0=0.73, h=0.700):
- def f_a2b(x):
- val = 0.5*sqrt(Om0) / x**3.0
- val /= sqrt(Om0/x**3.0 + Oml0 + (1.0-Om0-Oml0)/x**2.0)
- return val
- #val, err = si.quad(f_a2b,1,at)
- val = quad(f_a2b, 1, at)
- return val
-
-
-def b2a(bt, **kwargs):
- #converts code time into expansion factor
- #if Om0 ==1and OmL == 0 then b2a is (1 / (1-td))**2
- #if bt < -190.0 or bt > -.10: raise 'bt outside of range'
- f_b2a = lambda at: a2b(at, **kwargs)-bt
- return find_root(f_b2a, 1e-4, 1.1)
- #return so.brenth(f_b2a,1e-4,1.1)
- #return brent.brent(f_b2a)
-
-
-def a2t(at, Om0=0.27, Oml0=0.73, h=0.700):
- integrand = lambda x: 1./(x*sqrt(Oml0+Om0*x**-3.0))
- #current_time,err = si.quad(integrand,0.0,at,epsabs=1e-6,epsrel=1e-6)
- current_time = quad(integrand, 1e-4, at)
- #spacings = np.logspace(-5,np.log10(at),1e5)
- #integrand_arr = integrand(spacings)
- #current_time = np.trapz(integrand_arr,dx=np.diff(spacings))
- current_time *= 9.779/h
- return current_time
-
-
-def b2t(tb, n=1e2, logger=None, **kwargs):
- tb = np.array(tb)
- if len(np.atleast_1d(tb)) == 1:
- return a2t(b2a(tb))
- if tb.shape == ():
- return None
- if len(tb) < n:
- n = len(tb)
- age_min = a2t(b2a(tb.max(), **kwargs), **kwargs)
- age_max = a2t(b2a(tb.min(), **kwargs), **kwargs)
- tbs = -1.*np.logspace(np.log10(-tb.min()),
- np.log10(-tb.max()), n)
- ages = []
- for i, tbi in enumerate(tbs):
- ages += a2t(b2a(tbi)),
- if logger:
- logger(i)
- ages = np.array(ages)
- fb2t = np.interp(tb, tbs, ages)
- #fb2t = interp1d(tbs,ages)
- return fb2t*1e9*31556926
-
-
-def spread_ages(ages, logger=None, spread=.0e7*365*24*3600):
- #stars are formed in lumps; spread out the ages linearly
- da = np.diff(ages)
- assert np.all(da <= 0)
- #ages should always be decreasing, and ordered so
- agesd = np.zeros(ages.shape)
- idx, = np.where(da < 0)
- idx += 1
- #mark the right edges
- #spread this age evenly out to the next age
- lidx = 0
- lage = 0
- for i in idx:
- n = i-lidx
- #n stars affected
- rage = ages[i]
- lage = max(rage-spread, 0.0)
- agesd[lidx:i] = np.linspace(lage, rage, n)
- lidx = i
- #lage=rage
- if logger:
- logger(i)
- #we didn't get the last iter
- i = ages.shape[0]-1
- n = i-lidx
- #n stars affected
- rage = ages[i]
- lage = max(rage-spread, 0.0)
- agesd[lidx:i] = np.linspace(lage, rage, n)
- return agesd
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