From a89dc2f3cdb4f50a2f47b80bd0d6544bb1999d6d Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 27 Jun 2023 16:46:59 +0200 Subject: simdata: add API for loading time evolution of data along an axis --- simdata.py | 216 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 215 insertions(+), 1 deletion(-) 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() -- cgit v1.2.3