summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--simdata.py216
1 files changed, 215 insertions, 1 deletions
diff --git a/simdata.py b/simdata.py
index 2464cb2..77fb5a8 100644
--- a/simdata.py
+++ b/simdata.py
@@ -4,6 +4,11 @@ This module provides a wrapper around a Cactus simulation output, i.e. a directo
import os
from collections.abc import MutableMapping
+from functools import cached_property
+
+import numpy as np
+
+from math_utils.array_utils import array_reflect
from . import datafile
@@ -59,6 +64,212 @@ class _RefinedSlice(object):
return self._sd.df[key].slice(self.it, rl = self.rl.n)
raise IndexError
+class _AxisData:
+ """
+ Wrapper providing convenient access to quantities along a given axis
+ as functions of time. Data is exported as (t, x) 2D arrays.
+
+
+ """
+
+ sd = None
+ "parent simulation data"
+
+ axis = None
+ "axis along which the data is loaded, one of: x/y/z"
+ axis_idx = None
+ "axis index: 0/1/2 for x/y/z"
+
+ rl = None
+ "refinement level number"
+
+ it = None
+ "1D array of simulation iterations on which data is available"
+
+ t = None
+ "1D array of coordinate times, same shape as it"
+
+ x = None
+ "1D array of spatial coordinates on the axis"
+
+ T = None
+ "2D array of times corresponding to returned data"
+
+ X = None
+ "2D array of spatial coordinates corresponding to returned data"
+
+ _reflect = None
+ _zero_shift = None
+ _data_sl = None
+
+ def __init__(self, sd, rl, *,
+ axis = 'x', maxtime = None, reflect = False,
+ zero_shift = None):
+ self.sd = sd
+ self.rl = sd.rl[rl].n
+
+ if axis == 'x':
+ self.axis_idx = 0
+ elif axis == 'y':
+ self.axis_idx = 1
+ elif axis == 'z':
+ self.axis_idx = 2
+ else:
+ raise ValueError('axis must be one of: x y z')
+ self.axis = axis
+
+ self._reflect = reflect
+ self._zero_shift = zero_shift
+
+ df_name = os.path.splitext(os.path.basename(sd._df.path))[0].split('.')[0]
+ df = sd.df['%s.%s' % (df_name, axis)]
+
+ itertimes = df.itertimes[rl]
+ if maxtime is not None:
+ itertimes = list(filter(lambda it: it[1] <= maxtime, itertimes))
+ if len(itertimes) < 1:
+ raise ValueError('No data for maxtime: %s' % maxtime)
+
+ self.it, self.t = map(np.array, zip(*itertimes))
+
+ s = df.slice(it = self.it[0], rl = self.rl)
+ x = s.layout.coords[0]
+
+ ghosts = np.nonzero(np.abs(x) < 1e-10)[0][0]
+ self._data_sl = slice(ghosts, -ghosts)
+
+ x = x[self._data_sl]
+ if reflect:
+ x = array_reflect(x, parity = -1.0)
+
+ self.x = x
+
+ self.T, self.X = np.meshgrid(self.t, self.x, indexing = 'ij')
+
+ def _load_data(self, gf, parity):
+ df = self.sd.df['%s.%s' % (gf, self.axis)]
+
+ ret = []
+ for it in self.it:
+ d = df.slice(it, rl = self.rl).data[self._data_sl]
+ if self._reflect:
+ d = array_reflect(d, parity = parity)
+
+ ret.append(d)
+ return np.array(ret)
+
+ # scalars
+ @cached_property
+ def alpha(self): return self._load_data('alpha', 1.0)
+ @cached_property
+ def W(self): return self._load_data('W', 1.0)
+ @cached_property
+ def phi(self): return self._load_data('phi', 1.0)
+ @cached_property
+ def trK(self): return self._load_data('trK', 1.0)
+ @cached_property
+ def zeta(self): return self._load_data('zeta', 1.0)
+ @cached_property
+ def Kretsch(self): return self._load_data('Kretsch', 1.0)
+ @cached_property
+ def H(self): return self._load_data('H', 1.0)
+
+ # vectors
+ @cached_property
+ def betax(self):
+ if self._zero_shift:
+ return np.zeros_like(self.T)
+ return self._load_data('beta1', -1.0 if self.axis == 'x' else 1.0)
+ @cached_property
+ def betay(self):
+ if self._zero_shift:
+ return np.zeros_like(self.T)
+ return self._load_data('beta2', -1.0 if self.axis == 'y' else 1.0)
+ @cached_property
+ def betaz(self):
+ if self._zero_shift:
+ return np.zeros_like(self.T)
+ return self._load_data('beta3', -1.0 if self.axis == 'z' else 1.0)
+
+ @cached_property
+ def conf_fact(self): return 1.0 / (self.phi * self.phi)
+
+ # 2-tensors
+ @cached_property
+ def gtxx(self): return self._load_data('gt11', 1.0)
+ @cached_property
+ def gtyy(self): return self._load_data('gt22', 1.0)
+ @cached_property
+ def gtzz(self): return self._load_data('gt33', 1.0)
+ @cached_property
+ def Atxx(self): return self._load_data('At11', 1.0)
+ @cached_property
+ def Atyy(self): return self._load_data('At22', 1.0)
+ @cached_property
+ def Atzz(self): return self._load_data('At33', 1.0)
+
+ @cached_property
+ def gtxy(self):
+ return self._load_data('gt12', 1.0 if self.axis == 'z' else -1.0)
+ @cached_property
+ def gtxz(self):
+ return self._load_data('gt13', 1.0 if self.axis == 'y' else -1.0)
+ @cached_property
+ def gtyz(self):
+ return self._load_data('gt23', 1.0 if self.axis == 'x' else -1.0)
+ @cached_property
+ def Atxy(self):
+ return self._load_data('At12', 1.0 if self.axis == 'z' else -1.0)
+ @cached_property
+ def Atxz(self):
+ return self._load_data('At13', 1.0 if self.axis == 'y' else -1.0)
+ @cached_property
+ def Atyz(self):
+ return self._load_data('At23', 1.0 if self.axis == 'x' else -1.0)
+
+ @cached_property
+ def gxx(self): return self.gtxx * self.conf_fact
+ @cached_property
+ def gyy(self): return self.gtyy * self.conf_fact
+ @cached_property
+ def gzz(self): return self.gtzz * self.conf_fact
+ @cached_property
+ def Kxx(self): return self.Atxx * self.conf_fact + self.gxx * self.trK / 3.
+ @cached_property
+ def Kyy(self): return self.Atyy * self.conf_fact + self.gyy * self.trK / 3.
+ @cached_property
+ def Kzz(self): return self.Atzz * self.conf_fact + self.gzz * self.trK / 3.
+
+ @cached_property
+ def gxy(self): return self.gtxy * self.conf_fact
+ @cached_property
+ def gxz(self): return self.gtxz * self.conf_fact
+ @cached_property
+ def gyz(self): return self.gtyz * self.conf_fact
+ @cached_property
+ def Kxy(self): return self.Atxy * self.conf_fact + self.gxy * self.trK / 3.
+ @cached_property
+ def Kxz(self): return self.Atxz * self.conf_fact + self.gxz * self.trK / 3.
+ @cached_property
+ def Kyz(self): return self.Atyz * self.conf_fact + self.gyz * self.trK / 3.
+
+ @cached_property
+ def gamma(self):
+ """
+ spatial metric
+ """
+ return np.array([[self.gxx, self.gxy, self.gxz],
+ [self.gxz, self.gyy, self.gyz],
+ [self.gxz, self.gyz, self.gzz]])
+
+ @cached_property
+ def K(self):
+ """
+ spatial extrinsic curvature
+ """
+ ret = np.array([[self.Axx, self.Axy, self.Axz],
+ [self.Axz, self.Ayy, self.Ayz],
+ [self.Axz, self.Ayz, self.Azz]])
class _RefinementLevel(object):
# refinement level number
@@ -150,7 +361,7 @@ class SimulationData(object):
if 'H' in self.df:
df = self.df['H']
else:
- df = self.df.values()[0]
+ df = next(iter(self.df.values()))
self._df = df
@@ -159,5 +370,8 @@ class SimulationData(object):
for rl in df.rl:
self.rl.append(_RefinementLevel(self, len(self.rl)))
+ def axis_data(self, rl, **kwargs):
+ return _AxisData(self, rl, **kwargs)
+
def close(self):
self.df.close()