From e960c7f6a8bbc8fadd57a83369aa010aff41d1bf Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Mon, 16 Dec 2013 14:20:47 +0100 Subject: Initial commit. --- cactus_utils.py | 1156 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1156 insertions(+) create mode 100644 cactus_utils.py diff --git a/cactus_utils.py b/cactus_utils.py new file mode 100644 index 0000000..100f31f --- /dev/null +++ b/cactus_utils.py @@ -0,0 +1,1156 @@ +import h5py +import numpy as np +import pylab +import scipy.interpolate +import scipy.integrate +import scipy.optimize + +import collections +import os +import os.path +import re + +axis_prop = { 'x' : 0, 'y' : 1, 'z' : 2 } +axis_data = { 'x' : 2, 'y' : 1, 'z' : 0 } + +EPSILON = 1e-3 + + +def get_default_gridfunc(datafile): + """ + Get the name of the default grid function from a given data file + """ + funcs = datafile['Parameters and Global Attributes']['Datasets'] + + if funcs.dtype.type != np.string_: + funcs_str = ''.join(map(chr, funcs)) + else: + funcs_str = funcs.value + + return funcs_str.strip().split('\n')[0] + + +def _get_iteration_from_time(datafile, time): + """ + Get the simulation iteration with the given time. + """ + for item in datafile.values(): + try: + if item.attrs['time'] == time: + return item.attrs['timestep'] + except KeyError: + pass + + raise ValueError("No timelevel with time %d" % (time)) + + +def _get_last_iteration(datafile, name, rl = 0): + """ + Get the last iteration for which the data for the given gridfunction + on the given refinement level exists in the data file. + """ + last_it = -1 + + for item in datafile.values(): + try: + if (item.attrs['name'] == name and + item.attrs['level'] == rl): + last_it = max(last_it, item.attrs['timestep']) + except KeyError: + pass + + if last_it < 0: + raise KeyError("No compatible datafile in the file") + + return last_it + + +def _extract_paren(string): + """ + From a string starting with a parentheses-like pairwise construct -- + '(,[,{,<', extract the part inside it, taking into account nesting. + """ + if string[0] == '[': + parens = ('[', ']') + elif string[0] == '(': + parens = ('(', ')') + elif string[0] == '{': + parens = ('{', '}') + elif string[0] == '<': + parens = ('<', '>') + else: + raise ValueError("String '%s' does not start with a recognized pairwise delimiter" % string) + + cnt = 0 + for i in xrange(len(string)): + if string[i] == parens[0]: + cnt += 1 + elif string[i] == parens[1]: + cnt -= 1 + + if cnt == 0: + return string[:i + 1] + + raise ValueError("Unbalanced parentheses") + + +def _parse_grid_struct(string): + """ + Parse the grid structure v5 string into a dict. + """ + ret = {} + s = string + while len(s) > 0: + try: + key, tail = s.split(':', 1) + except ValueError: + if len(ret) > 0: + return ret + else: + return string + + val = _extract_paren(tail) + ret[key] = val + + s = tail[len(val) + 1:] + + return ret + + +def _get_timestep(datafile): + """ + Get the time step on the coarsest grid. + """ + gridstruct = datafile['Parameters and Global Attributes']['Grid Structure v5'] + + if gridstruct.dtype.type != np.string_: + s = ''.join(map(chr, gridstruct)) + else: + s = gridstruct.value + + match = re.search('grid_delta_times:\[(.*?)\],', s) + if match is None: + raise ValueError("Invalid grid structure entry") + + deltas = eval('list(%s)' % (match.group(1))) + return deltas[0] + + +def _get_grid_spacing(datafile): + """ + Get the numbers of grid points (not counting ghosts) in each direction. + """ + gridstruct = datafile['Parameters and Global Attributes']['Grid Structure v5'] + + if gridstruct.dtype.type != np.string_: + s = ''.join(map(chr, gridstruct)) + else: + s = gridstruct.value + + struct = _parse_grid_struct(s) + match = re.search('\[([0-9]*),([0-9]*),([0-9]*)\]/[0-9]*\)', struct['grid_structure']) + gridpoints = np.array(map(int, (match.group(1), match.group(2), match.group(3)))) + + match = re.search('\[([0-9]*),([0-9]*),([0-9]*)\]', struct['grid_ghosts']) + ghosts = np.array(map(int, (match.group(1), match.group(2), match.group(3)))) + + gridpoints -= 2 * ghosts + + return gridpoints + + +def _get_data_axis(dataset, axis, axvalue): + """ + Get a tuple of (coordinates, data) along the specified axis ('x', 'y', or 'z'). + axvalue is the value of the other two axes, it may be a number for the same + value on both or an iterable with at least two elements for distinct values. + """ + slicelist = [None] * 3 + try: + if len(axvalue) > 2: + raise ValueError("Need at most two axis values") + elif len(axvalue) == 1: + axvalue = [axvalue] * 2 + except TypeError: + axvalue = [axvalue] * 2 + + # data is stored in [z, y, x] order, properties are in normal [x, y, z] order + axes = collections.OrderedDict([["x", {"data_idx" : 2, "idx" : 0 }], + ["y", {"data_idx" : 1, "idx" : 1 }], + ["z", {"data_idx" : 0, "idx" : 2 }]]) + + # for the requested axis, just remove the ghost zones + ax = axes.pop(axis) + ghosts = dataset.attrs["cctk_nghostzones"][ax["idx"]] + delta = dataset.attrs["delta"][ax["idx"]] + origin = dataset.attrs["origin"][ax["idx"]] + slicelist[ax["data_idx"]] = slice(ghosts, -ghosts) + + # for the other two axes, locate the proper values + while len(axes): + axname, ax = axes.popitem(last = False) + + o = dataset.attrs["origin"][ax["idx"]] + d = dataset.attrs["delta"][ax["idx"]] + idx = int((axvalue.pop(0) - o) / d) + + if idx >= dataset.shape[ax["data_idx"]]: + raise ValueError("Axis %s is outside the domain [%f, %f]" % (axname, o, o + d * dataset.shape[ax["data_idx"]])) + + slicelist[ax["data_idx"]] = idx + + data = dataset[tuple(slicelist)] + + # now compute the corresponding distances from coordinate origin + r = np.array([(origin + ghosts * delta + n * delta) for n in xrange(data.shape[0])]) + + return r, data + + +def plot(datafile, axis = "x", axisvalue = 0.0, iteration = "last", time = -1, gridfunc = None, rl = 0, c = -1, mod = None, label = None): + """ + Draw a 2D plot. + """ + # guess dataset name + if gridfunc is None: + gridfunc = get_default_gridfunc(datafile) + + if time >= 0: + iteration = _get_iteration_from_time(datafile, time) + elif iteration == "last": + iteration = _get_last_iteration(datafile, gridfunc, rl) + + dataset = datafile['%s it=%d tl=0 rl=%d%s' % (gridfunc, iteration, rl, " c=%d" % c if c >= 0 else "")] + if axis == "x" or axis == "y" or axis == "z": + r, data = _get_data_axis(dataset, axis, axisvalue) +# elif len(axis) == 2 and re.match("[xyz]{2}", axis): +# r, data = get_data_diagonal(dataset, axis) + else: + raise ValueError("Invalid axis specifier %s" % axis) + + if mod: + data = mod(data) + + if not label: + label = "%s time=%f level=%d" % (gridfunc, dataset.attrs['time'], rl) + + pylab.plot(r, data, label = label) + pylab.legend() + + +def plot_diff(datafile, axis = "x", axisvalue = 0.0, it1 = 0, it2 = "last", time1 = -1, time2 = -1, gridfunc = None, rl = 0, absdiff = 0): + """ + Plot a difference between two times or iterations. + """ + # guess dataset name + if gridfunc is None: + gridfunc = get_default_gridfunc(datafile) + + if time1 >= 0: + it1 = _get_iteration_from_time(datafile, time1) + elif it1 == "last": + it1 = _get_last_iteration(datafile, gridfunc, rl) + + if time2 >= 0: + it2 = _get_iteration_from_time(datafile, time2) + elif it2 == "last": + it2 = _get_last_iteration(datafile, gridfunc, rl) + + dataset1 = datafile['%s it=%d tl=0 rl=%d' % (gridfunc, it1, rl)] + dataset2 = datafile['%s it=%d tl=0 rl=%d' % (gridfunc, it2, rl)] + if axis == "x" or axis == "y" or axis == "z": + r, data1 = _get_data_axis(dataset1, axis, axisvalue) + r, data2 = _get_data_axis(dataset2, axis, axisvalue) +# elif len(axis) == 2 and re.match("[xyz]{2}", axis): +# r, data = get_data_diagonal(dataset, axis) + else: + raise ValueError("Invalid axis specifier %s" % axis) + + data = data1 - data2 + if absdiff: + data = np.abs(data) + + pylab.plot(r, data, label = "%s" % (gridfunc)) + pylab.legend() + + +def plot_convergence_time(datafiles, axis, gridfunc, rl): + pass + + +def plot_convergence_space(datafiles, axis, gridfunc, rl): + pass + + +def plot_convergence(datafiles, axis = "x", gridfunc = None, rl = 0, time = "last", mode = "selfconsistent"): + if mode == "selfconsistent" and len(datafiles) < 3: + raise ValueError("At least three data sets needed for self-consistent convergence plot") + elif len(datafiles) < 2: + raise ValueError("At least two data sets needed for convergence plot") + + if gridfunc is None: + gridfunc = get_default_gridfunc(datafiles[0]) + + # check if the datafiles have changing time steps + prev_timestep = 0 + timestep_changing = 0 + for df in datafiles: + timestep = _get_timestep(df) + if prev_timestep == 0: + prev_timestep = timestep + elif prev_timestep != timestep: + timestep_changing = 1 + break + + datasets = [f['%s it=0 tl=0 rl=%d' % (gridfunc, rl)] for f in datafiles] + + # check if our chosen datasets have different grid spacings + prev_spacing = None + spacing_changing = 0 + for ds in datasets: + spacing = ds.attrs['delta'] + if prev_spacing is None: + prev_spacing = spacing + elif not np.all(spacing == prev_spacing): + spacing_changing = 1 + break + del datasets + + if timestep_changing and spacing_changing: + raise ValueError("Both the time step and grid spacing change between files.") + if not timestep_changing and not spacing_changing: + raise ValueError("Neither the time step nor grid spacing change between files.") + + if timestep_changing: + plot_convergence_time(datafiles, axis, gridfunc, rl) + else: + plot_convergence_space(datafiles, axis, gridfunc, rl) + + +def dminus(arr, dir, dx): + """ + Fourth order centered derivative in the given direction. + Second order centered derivatives are used for the points next to the + border ones and first order for the border points. + """ + ret = np.empty_like(arr) + slicelist_ret = [slice(None)] * len(arr.shape) + slicelist_ret[dir] = slice(1, None) + slicelist_m1 = [slice(None)] * len(arr.shape) + slicelist_m1[dir] = slice(0, -1) + slicelist_0 = [slice(None)] * len(arr.shape) + slicelist_0[dir] = slice(1, None) + + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_0)] - arr[tuple(slicelist_m1)]) / dx + + slicelist_ret[dir] = slice(0, 1) + slicelist_m1[dir] = slice(1, 2) + slicelist_0[dir] = slice(0, 1) + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_m1)] - arr[tuple(slicelist_0)]) / dx + + return ret + + +def dplus(arr, dir, dx): + """ + Fourth order centered derivative in the given direction. + Second order centered derivatives are used for the points next to the + border ones and first order for the border points. + """ + ret = np.empty_like(arr) + slicelist_ret = [slice(None)] * len(arr.shape) + slicelist_ret[dir] = slice(0, -1) + slicelist_p1 = [slice(None)] * len(arr.shape) + slicelist_p1[dir] = slice(1, None) + slicelist_0 = [slice(None)] * len(arr.shape) + slicelist_0[dir] = slice(0, -1) + + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_p1)] - arr[tuple(slicelist_0)]) / dx + + slicelist_ret[dir] = slice(-1, None) + slicelist_p1[dir] = slice(-2, -1) + slicelist_0[dir] = slice(-1, None) + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_0)] - arr[tuple(slicelist_p1)]) / dx + + return ret + + +def diff4(arr, dir, dx, coeff = None, div = 12.): + """ + Fourth order centered derivative in the given direction. + Second order centered derivatives are used for the points next to the + border ones and first order for the border points. + """ + ret = np.zeros_like(arr) + slicelist_ret = [slice(None)] * len(arr.shape) + slicelist_ret[dir] = slice(2, -2) + slicelist_0 = [slice(None)] * len(arr.shape) + slicelist_0[dir] = slice(2, -2) + slicelist_p1 = [slice(None)] * len(arr.shape) + slicelist_p1[dir] = slice(3, -1) + slicelist_p2 = [slice(None)] * len(arr.shape) + slicelist_p2[dir] = slice(4, None) + slicelist_m1 = [slice(None)] * len(arr.shape) + slicelist_m1[dir] = slice(1, -3) + slicelist_m2 = [slice(None)] * len(arr.shape) + slicelist_m2[dir] = slice(0, -4) + + sl = (slicelist_p2, slicelist_p1, slicelist_0, slicelist_m1, slicelist_m2) + if coeff is None: + coeff = (-1, 8, 0, -8, 1) + + for s, c in zip(sl, coeff): + ret[tuple(slicelist_ret)] += c * arr[tuple(s)] + ret[tuple(slicelist_ret)] /= (div * dx) + + slicelist_ret[dir] = slice(0, 1) + slicelist_p1[dir] = slice(1, 2) + slicelist_m1[dir] = slice(0, 1) + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_p1)] - arr[tuple(slicelist_m1)]) / dx + + slicelist_ret[dir] = slice(-1, None) + slicelist_p1[dir] = slice(-1, None) + slicelist_m1[dir] = slice(-2, -1) + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_p1)] - arr[tuple(slicelist_m1)]) / dx + + slicelist_ret[dir] = slice(1, 2) + slicelist_p1[dir] = slice(2, 3) + slicelist_m1[dir] = slice(0, 1) + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_p1)] - arr[tuple(slicelist_m1)]) / (2 * dx) + + slicelist_ret[dir] = slice(-2, -1) + slicelist_p1[dir] = slice(-1, None) + slicelist_m1[dir] = slice(-3, -2) + ret[tuple(slicelist_ret)] = (arr[tuple(slicelist_p1)] - arr[tuple(slicelist_m1)]) / (2 * dx) + return ret + + +def get_meshgrid(datafile, gridfunc = None, rl = 0, c = -1, axis = "z"): + """ + Get a tuple of coordinate mesh grids for the given data file. + """ + if not gridfunc: + gridfunc = get_default_gridfunc(datafile) + dataset = datafile['%s it=0 tl=0 rl=%d%s' % (gridfunc, rl, " c=%d" % c if c >= 0 else "")] + o = dataset.attrs['origin'] + d = dataset.attrs['delta'] + if axis == 'z': + idx = (0, 1) + elif axis == 'y': + idx = (0, 2) + elif axis == 'x': + idx = (1, 2) + else: + raise ValueError('Invalid axis %s' % axis) + + shape = (dataset.shape[2], dataset.shape[1], dataset.shape[0]) + return np.meshgrid(np.array([o[idx[0]] + n * d[idx[0]] for n in xrange(shape[idx[0]])]), np.array([o[idx[1]] + n * d[idx[1]] for n in xrange(shape[idx[1]])])) + + +def _get_extrema(ar, rad = 3): + """ + Find local extrema of a 1-D array. + Returned is a tuple (maxima, minima) of the (coordinate, value) of the + extrema. rad determines the number of points in each direction over which + the data must be monotonous for the center to be classified as an extremum. + """ + maxima = [] + minima = [] + dar = diff4(ar, 0, 1) + spl = scipy.interpolate.UnivariateSpline(range(ar.shape[0]), ar) + + for i in xrange(rad, len(ar) - rad): + maximum = minimum = True + for j in xrange(1, rad): + for dir in (1, -1): + idx = i + j * dir + if ar[idx] > ar[idx - dir]: + maximum = False + elif ar[idx] < ar[idx - dir]: + minimum = False + + if maximum or minimum: + x = np.array(range(2 * rad + 1)) + y = dar[i - rad : i + rad + 1] + roots = scipy.interpolate.UnivariateSpline(x, y).roots() + if len(roots) > 0: + root = roots[0] + i - rad + (maxima if maximum else minima).append((root, spl(root))) + + return maxima, minima + + +def _get_movement(coord_old, val_old, ex_new): + """ + Return the number of points an extremum i in val_old moves in val_new + or None if no movement could be detected. + """ + dist = None + for coord_new, val_new in ex_new: + ratio = val_new / val_old + dist_new = coord_new - coord_old + if ratio < 1.1 and ratio > 0.9 and (dist is None or abs(dist_new) < abs(dist)): + dist = dist_new + return dist + + +def calc_movement(ar1, ar2): + """ + Find all local extrema in ar1 and calculate how much the move in ar2. + Return the list of matched movements, in points. + """ + maxima1, minima1 = _get_extrema(ar1) + maxima2, minima2 = _get_extrema(ar2) + + ret = [] + for coord, val in maxima1: + diff = _get_movement(coord, val, maxima2) + if diff is not None: + ret.append(diff) + for coord, val in minima1: + diff = _get_movement(coord, val, minima2) + if diff is not None: + ret.append(diff) + + return ret + + +def gf_itertime(datafile, rl = 0): + """ + Find all iterations that are stored in the given datafile. + Return an ordered list of (iteration, time) tuples. + """ + ret = set() + + rl = ' rl=%d' % rl + for key in datafile.keys(): + if not rl in key: + continue + + match = re.search('it=([0-9]*)', key) + if match: + ret.add((int(match.group(1)), datafile[key].attrs['time'])) + + return sorted(ret) + + +def gf_iter(gridfunc, rl = 0): + """ + Find all iterations that are stored in the given datafile. + Return an ordered list of iterations. + """ + return zip(*gf_itertime(gridfunc, rl))[0] + + +def gf_time(gridfunc, rl = 0): + """ + Find all iterations that are stored in the given datafile. + Return an ordered list of times. + """ + return zip(*gf_itertime(gridfunc, rl))[1] + + +def trumpet_r(R, M = 1.0): + fac1 = ((2 * R + M + np.sqrt(4 * (R ** 2) + 4 * M * R + 3 * (M ** 2))) / 4) + fac2 = ((4 + 3 * np.sqrt(2)) * (2 * R - 3 * M)) / (8 * R + 6 * M + 3 * np.sqrt(8 * (R ** 2) + 8 * M * R + 6 * (M ** 2))) + return fac1 * (fac2 ** (1. / np.sqrt(2))) + + +def trumpet_log1_C2(n): + s = np.sqrt(4 + 9 * (n ** 2)) + a_c = np.sqrt((s - 3 * n) / (s + 3 * n)) + C2 = ((3 * n + s) ** 3) / (128 * (n ** 3)) * np.exp(-2 * a_c / n) + return C2 + + +def _trumpet_log1_alphabar(alpha, R, M, n): + num = 3 * M - 2 * R + 2 * R * (alpha ** 2) + den = (R - 2 * M + n * R * alpha - R * (alpha ** 2)) + ret = - n * num / (R * den) + return ret + + +def _trumpet_log1_alpha_implicit(alpha, R, M, n): + C2 = trumpet_log1_C2(n) + return alpha ** 2 - 1 + 2 * M / R - C2 * np.exp(2 * alpha / n) / (R ** 4) + + +def trumpet_log1_lapse(R, M = 1.0, n = 2.0): + R_c = (3 * (n ** 2) * M + np.sqrt((2 * n * M) ** 2 + (3 * (n ** 2) * M) ** 2)) / (4 * (n ** 2)) + s = np.sqrt(4 + 9 * (n ** 2)) + a_c = np.sqrt((s - 3 * n) / (s + 3 * n)) + + for i in xrange(R.shape[0]): + if R[i] > R_c: + break + + tmp = scipy.integrate.odeint(_trumpet_log1_alphabar, a_c, [R_c, R[i - 1]], args = (M, n))[1] + ret1 = scipy.integrate.odeint(_trumpet_log1_alphabar, tmp, R[:i][::-1], args = (M, n))[::-1] + ret1 = ret1.reshape(ret1.shape[0]) + + tmp = scipy.integrate.odeint(_trumpet_log1_alphabar, a_c, [R_c, R[i]], args = (M, n))[1] + ret2 = scipy.integrate.odeint(_trumpet_log1_alphabar, tmp, R[i:], args = (M, n)) + ret2 = ret2.reshape(ret2.shape[0]) + + return np.concatenate((ret1, ret2)) + + +def trumpet_log1_trk(R, M = 1.0, n = 2.0): + a = trumpet_log1_lapse(R, M, n) + b = np.sqrt((a ** 2) - (1 - 2 * M / R)) + bbar = diff4(b, 0, R[1] - R[0]) + + K = 2 * b / R + bbar + return K + + +def trumpet_lapse(R, M = 1.0): + C2 = 27. * (M ** 4) / 16 + return np.sqrt(1 - 2 * M / R + C2 / (R ** 4)) + + +def trumpet_log1_r(R, M = 1.0, n = 2.0): + R_aux = np.linspace(1.32, 100, 10000) + alpha_aux = trumpet_log1_lapse(R_aux) + + R_alpha_aux = scipy.interpolate.UnivariateSpline(alpha_aux, R_aux, s = 0) + alpha_s = alpha_aux[20] + coeff = R_alpha_aux(alpha_s) ** (1. / alpha_s) + + def integrand_aux(a, R_alpha = R_alpha_aux, alpha = alpha_aux): + if a > alpha[-1]: + RR = 2 / (1 - a ** 2) + else: + RR = R_alpha(a) + ret = np.log(RR) / (a ** 2) + return ret + C0, err = scipy.integrate.quad(integrand_aux, alpha_s, 1) + + alpha = trumpet_log1_lapse(R, M, n) + R_alpha = scipy.interpolate.UnivariateSpline(alpha, R, s = 0) + Rbar_alpha = scipy.interpolate.UnivariateSpline(alpha[1:-1], (R[2:] - R[:-2]) / (2 * (R[1] - R[0])), s = 0) + + def integrand(a, R_alpha = R_alpha, alpha = alpha): + if a > alpha[-1]: + RR = 2 / (1 - a ** 2) + else: + RR = R_alpha(a) + ret = np.log(RR) / (a ** 2) + return ret + + def integrand2(a, R_alpha = R_alpha, alpha = alpha): + if a > alpha[-1]: + RR = 2 / (1 - a) + Rbar = - 4 * a / ((1 - a ** 2) ** 2) + else: + RR = R_alpha(a) + Rbar = Rbar_alpha(a) + return Rbar / (RR * a) + + r1 = [] + r2 = [] + + if alpha[0] >= alpha_s: + i = 0 + else: + for i in xrange(alpha.shape[0]): + if alpha[i] >= alpha_s: + break + + for j in xrange(i): + a = alpha[j] + integral, err, = scipy.integrate.quad(integrand2, a, alpha_s) + r1.append(np.exp(-integral - C0) * coeff) + + int_prev = 0 + alpha_prev = alpha_s + for j in xrange(i, alpha.shape[0]): + a = alpha[j] + integral, err, = scipy.integrate.quad(integrand, alpha_prev, a) + integral += int_prev + alpha_prev = a + int_prev = integral + + r2.append(np.exp(integral - C0) * R_alpha(a) ** (1. / a)) + + return np.array(r1 + r2) + + +def calc_kretschmann(s, axis = 2, xoffset = 0): + rl = s.rl + + sliceidx = [slice(0, 7), slice(0, 7), slice(xoffset, xoffset + 7)] + sliceidx[axis] = slice(None) + sliceidx = tuple(sliceidx) + + offsetidx = [3] * 3 + offsetidx[axis] = slice(None) + offsetidx = tuple(offsetidx) + + #sliceidx = [slice(None)]*3 + #sliceidx = tuple(sliceidx) + + #offsetidx = [slice(None)]*3 + #offsetidx = tuple(offsetidx) + + # load grid functions + phi = s['phi'][sliceidx] + trK = s['trK'][sliceidx] + gt = np.array([[s['gt33'][sliceidx], s['gt23'][sliceidx], s['gt13'][sliceidx]], + [s['gt23'][sliceidx], s['gt22'][sliceidx], s['gt12'][sliceidx]], + [s['gt13'][sliceidx], s['gt12'][sliceidx], s['gt11'][sliceidx]]]) + At = np.array([[s['At33'][sliceidx], s['At23'][sliceidx], s['At13'][sliceidx]], + [s['At23'][sliceidx], s['At22'][sliceidx], s['At12'][sliceidx]], + [s['At13'][sliceidx], s['At12'][sliceidx], s['At11'][sliceidx]]]) + + # compute derivatives needed later + dphi = np.array([diff4(phi, i, rl.dx[0]) for i in 0, 1, 2]) + d2phi = np.array([diff4(dphi, i, rl.dx[0]) for i in 1, 2, 3]) + dgt = np.array([diff4(gt, i, rl.dx[0]) for i in 2, 3, 4]) + d2gt = np.array([diff4(dgt, i, rl.dx[0]) for i in 3, 4, 5]) + dAt = np.array([diff4(At, i, rl.dx[0]) for i in 2, 3, 4]) + dtrK = np.array([diff4(trK, i, rl.dx[0]) for i in 0, 1, 2]) + + phi = phi[offsetidx] + trK = trK[offsetidx] + gt = gt[(slice(None),) * 2 + offsetidx] + At = At[(slice(None),) * 2 + offsetidx] + + dphi = dphi [(slice(None),) + offsetidx] + d2phi = d2phi[(slice(None),) * 2 + offsetidx] + dgt = dgt [(slice(None),) * 3 + offsetidx] + d2gt = d2gt [(slice(None),) * 4 + offsetidx] + dAt = dAt [(slice(None),) * 3 + offsetidx] + dtrK = dtrK [(slice(None),) + offsetidx] + + e4phi = 1. / phi ** 2 + gtu = np.zeros_like(gt) + for i in xrange(gt.shape[2]): + gtu[:, :, i] = np.linalg.inv(gt[:, :, i]) + #for i in xrange(gt.shape[2]): + # for j in xrange(gt.shape[3]): + # for k in xrange(gt.shape[4]): + # gtu[:, :, i, j, k] = np.linalg.inv(gt[:, :, i, j, k]) + + g = e4phi * gt + gu = gtu / e4phi + + K = e4phi * At + (1. / 3.) * trK * g + + Dphi = -0.5 * dphi / phi + D2phi = - 0.5 * (d2phi - np.einsum('i...,j...->ij...', dphi, dphi) / phi) / phi + + Dg = e4phi * (4 * np.einsum('i...,jk...->ijk...', Dphi, gt) + dgt) + D2g = e4phi * (4 * np.einsum('ij...,kl...->ijkl...', D2phi + 4 * np.einsum('i...,j...->ij...', Dphi, Dphi), gt) + 4 * np.einsum('i...,jkl...->ijkl...', Dphi, dgt) + 4 * np.einsum('i...,jkl...->jikl...', Dphi, dgt) + d2gt) + + Gl = 0.5 * (np.einsum('cab...->abc...', Dg) + np.einsum('bac...->abc...', Dg) - Dg) + G = np.einsum('ij...,jkl...->ikl...', gu, Gl) + + R = (0.5 * (np.einsum('adbc...->abcd...', D2g) + np.einsum('bcad...->abcd...', D2g) - np.einsum('bdac...->abcd...', D2g) - np.einsum('acbd...->abcd...', D2g)) + + np.einsum('ead...,ebc...->abcd...', Gl, G) - np.einsum('eac...,ebd...->abcd...', Gl, G)) + + A = R + np.einsum('ik...,jl...->ijkl...', K, K) - np.einsum('il...,jk...->ijkl...', K, K) + trA = np.einsum('kl...,kilj...->ij...', gu, A) + term1 = (np.einsum('ai...,bj...,ck...,dl...,abcd...,ijkl...', gu, gu, gu, gu, A, A)) + term3 = 4 * np.einsum('ai...,bj...,ab...,ij...', gu, gu, trA, trA) + + t1 = e4phi * (4 * np.einsum('i...,jk...->ijk...', Dphi, At) + dAt) + t2 = (1. / 3.) * (np.einsum('i...,jk...->ijk...', dtrK, g) + trK * Dg) + t3 = - np.einsum('iab...,ic...->bac...', G, K) - np.einsum('icb...,ia...->bac...', G, K) + CDK = t1 + t2 + t3 + CDKu = np.einsum('ai...,bj...,ck...,ijk...->abc...', gu, gu, gu, CDK) + + term2 = 8 * (np.einsum('abc...,bac...', CDK, CDKu) - np.einsum('abc...,abc...', CDK, CDKu)) + + return term1 + term2 + term3 + + +def calc_kretschmann2(s): + rl = s.rl + + sliceidx = [slice(0, 7), slice(None), slice(None)] + sliceidx = tuple(sliceidx) + + offsetidx = [3, slice(None), slice(None)] + offsetidx = tuple(offsetidx) + + #sliceidx = [slice(None)]*3 + #sliceidx = tuple(sliceidx) + + #offsetidx = [slice(None)]*3 + #offsetidx = tuple(offsetidx) + + # load grid functions + phi = s['phi'][sliceidx] + trK = s['trK'][sliceidx] + gt = np.array([[s['gt33'][sliceidx], s['gt23'][sliceidx], s['gt13'][sliceidx]], + [s['gt23'][sliceidx], s['gt22'][sliceidx], s['gt12'][sliceidx]], + [s['gt13'][sliceidx], s['gt12'][sliceidx], s['gt11'][sliceidx]]]) + At = np.array([[s['At33'][sliceidx], s['At23'][sliceidx], s['At13'][sliceidx]], + [s['At23'][sliceidx], s['At22'][sliceidx], s['At12'][sliceidx]], + [s['At13'][sliceidx], s['At12'][sliceidx], s['At11'][sliceidx]]]) + + # compute derivatives needed later + dphi = np.array([diff4(phi, i, rl.dx[0]) for i in 0, 1, 2]) + d2phi = np.array([diff4(dphi, i, rl.dx[0]) for i in 1, 2, 3]) + dgt = np.array([diff4(gt, i, rl.dx[0]) for i in 2, 3, 4]) + d2gt = np.array([diff4(dgt, i, rl.dx[0]) for i in 3, 4, 5]) + dAt = np.array([diff4(At, i, rl.dx[0]) for i in 2, 3, 4]) + dtrK = np.array([diff4(trK, i, rl.dx[0]) for i in 0, 1, 2]) + + phi = phi[offsetidx] + trK = trK[offsetidx] + gt = gt[(slice(None),) * 2 + offsetidx] + At = At[(slice(None),) * 2 + offsetidx] + + dphi = dphi [(slice(None),) + offsetidx] + d2phi = d2phi[(slice(None),) * 2 + offsetidx] + dgt = dgt [(slice(None),) * 3 + offsetidx] + d2gt = d2gt [(slice(None),) * 4 + offsetidx] + dAt = dAt [(slice(None),) * 3 + offsetidx] + dtrK = dtrK [(slice(None),) + offsetidx] + + e4phi = 1. / phi ** 2 + gtu = np.zeros_like(gt) + for i in xrange(gt.shape[2]): + for j in xrange(gt.shape[3]): + gtu[:, :, i, j] = np.linalg.inv(gt[:, :, i, j]) + #for i in xrange(gt.shape[2]): + # for j in xrange(gt.shape[3]): + # for k in xrange(gt.shape[4]): + # gtu[:, :, i, j, k] = np.linalg.inv(gt[:, :, i, j, k]) + + g = e4phi * gt + gu = gtu / e4phi + + K = e4phi * At + (1. / 3.) * trK * g + + Dphi = -0.5 * dphi / phi + D2phi = - 0.5 * (d2phi - np.einsum('i...,j...->ij...', dphi, dphi) / phi) / phi + + Dg = e4phi * (4 * np.einsum('i...,jk...->ijk...', Dphi, gt) + dgt) + D2g = e4phi * (4 * np.einsum('ij...,kl...->ijkl...', D2phi + 4 * np.einsum('i...,j...->ij...', Dphi, Dphi), gt) + 4 * np.einsum('i...,jkl...->ijkl...', Dphi, dgt) + 4 * np.einsum('i...,jkl...->jikl...', Dphi, dgt) + d2gt) + + Gl = 0.5 * (np.einsum('cab...->abc...', Dg) + np.einsum('bac...->abc...', Dg) - Dg) + G = np.einsum('ij...,jkl...->ikl...', gu, Gl) + + R = (0.5 * (np.einsum('adbc...->abcd...', D2g) + np.einsum('bcad...->abcd...', D2g) - np.einsum('bdac...->abcd...', D2g) - np.einsum('acbd...->abcd...', D2g)) + + np.einsum('ead...,ebc...->abcd...', Gl, G) - np.einsum('eac...,ebd...->abcd...', Gl, G)) + + A = R + np.einsum('ik...,jl...->ijkl...', K, K) - np.einsum('il...,jk...->ijkl...', K, K) + trA = np.einsum('kl...,kilj...->ij...', gu, A) + term1 = (np.einsum('ai...,bj...,ck...,dl...,abcd...,ijkl...', gu, gu, gu, gu, A, A)) + term3 = 4 * np.einsum('ai...,bj...,ab...,ij...', gu, gu, trA, trA) + + t1 = e4phi * (4 * np.einsum('i...,jk...->ijk...', Dphi, At) + dAt) + t2 = (1. / 3.) * (np.einsum('i...,jk...->ijk...', dtrK, g) + trK * Dg) + t3 = - np.einsum('iab...,ic...->bac...', G, K) - np.einsum('icb...,ia...->bac...', G, K) + CDK = t1 + t2 + t3 + CDKu = np.einsum('ai...,bj...,ck...,ijk...->abc...', gu, gu, gu, CDK) + + term2 = 8 * (np.einsum('abc...,bac...', CDK, CDKu) - np.einsum('abc...,abc...', CDK, CDKu)) + + return term1 + term2 + term3 + + +def l2norm(sd, it, gf, maxrl = None, minrl = None): + results = [] + + if maxrl is None: + maxrl = len(sd.rl) - 1 + if minrl is None: + minrl = 0 + + prev_bounds_low = [] + prev_bounds_high = [] + for rl in sd.rl[maxrl:minrl:-1]: + s = rl.slice(it) + trim = tuple([slice(rl.ghosts[i], -rl.ghosts[i]) for i in xrange(3)]) + data = s[gf][trim] ** 2 + x = s.x[trim[0]] + y = s.y[trim[1]] + z = s.z[trim[2]] + + if not prev_bounds_low: + while len(data.shape) > 0: + data = scipy.integrate.simps(data, dx = rl.dx[0]) + results.append(data) + else: + res = 0.0 + if prev_bounds_low[0] > x[0]: + for i, xx in enumerate(x): + if prev_bounds_low[0] <= xx: + break + + tmp = scipy.integrate.simps(data[:, :, :i], dx = rl.dx[0]) + while len(tmp.shape) > 0: + tmp = scipy.integrate.simps(tmp, dx = rl.dx[0]) + res += tmp + if prev_bounds_low[1] > y[0]: + for i, yy in enumerate(y): + if prev_bounds_low[1] <= yy: + break + + tmp = scipy.integrate.simps(data[:, :i], axis = 1, dx = rl.dx[0]) + while len(tmp.shape) > 0: + tmp = scipy.integrate.simps(tmp, dx = rl.dx[0]) + res += tmp + if prev_bounds_low[2] > z[0]: + for i, zz in enumerate(z): + if prev_bounds_low[2] <= zz: + break + + tmp = scipy.integrate.simps(data[:i], axis = 0, dx = rl.dx[0]) + while len(tmp.shape) > 0: + tmp = scipy.integrate.simps(tmp, dx = rl.dx[0]) + res += tmp + if prev_bounds_high[0] < x[-1]: + for i, xx in enumerate(x): + if prev_bounds_high[0] >= xx: + i -= 1 + break + + tmp = scipy.integrate.simps(data[:, :, i:], dx = rl.dx[0]) + while len(tmp.shape) > 0: + tmp = scipy.integrate.simps(tmp, dx = rl.dx[0]) + res += tmp + if prev_bounds_high[1] < y[-1]: + for i, yy in enumerate(y): + if prev_bounds_high[1] >= yy: + i -= 1 + break + + tmp = scipy.integrate.simps(data[:, i:], axis = 1, dx = rl.dx[0]) + while len(tmp.shape) > 0: + tmp = scipy.integrate.simps(tmp, dx = rl.dx[0]) + res += tmp + if prev_bounds_high[2] < z[-1]: + for i, zz in enumerate(z): + if prev_bounds_high[2] >= zz: + i -= 1 + break + + tmp = scipy.integrate.simps(data[i:], axis = 0, dx = rl.dx[0]) + while len(tmp.shape) > 0: + tmp = scipy.integrate.simps(tmp, dx = rl.dx[0]) + res += tmp + + results.append(res) + + prev_bounds_low = [x[0], y[0], z[0]] + prev_bounds_high = [x[-1], y[-1], z[-1]] + + return sum(results) + + +class RefinedSlice: + # simulation time on this slice + time = -1 + # simulation iteration on this slice + it = -1 + + # [x, y, z] coordinates of the origin of the grid on this slice + o = None + + # a list of arrays of coordinate points along each direction [x, y, z] + c = None + + # arrays of coordinate points along each direction + x = None + y = None + z = None + + rl = None + + # private + _sd = None + _rl = None + + def __init__(self, sd, rl, time = -1, iteration = -1): + self._sd = sd + self._rl = rl + self.rl = rl + + if iteration == 'last': + iteration = rl.itertimes[-1][0] + + if time >= 0: + for it, t in rl.itertimes: + if abs(t - time) < EPSILON: + self.time = time + self.it = it + break + if self.time < 0: + raise IndexError('No such time: ' + str(time)) + elif iteration >= 0: + for it, t in rl.itertimes: + if it == iteration: + self.time = t + self.it = it + break + if self.it < 0: + raise IndexError('No such iteration: ' + str(iteration)) + else: + raise TypeError('Neither time nor iteration provided') + + # pick a representative datafile and get the grid properties from it + if 'H' in self._sd.df: + df = self._sd.df['H'] + else: + df = self._sd.df.values()[0] + + gf = get_default_gridfunc(df) + ds = df['%s it=%d tl=0 rl=%d' % (gf, self.it, rl.n)] + + self.o = ds.attrs['origin'] + + self.c = [np.linspace(self.o[axis_prop[i]], self.o[axis_prop[i]] + rl.dx[axis_prop[i]] * ds.shape[axis_data[i]], ds.shape[axis_data[i]], endpoint = False) for i in ('x', 'y', 'z')] + self.x, self.y, self.z = self.c + + def __getitem__(self, key): + if key in self._sd.df: + gf = get_default_gridfunc(self._sd.df[key]) + return self._sd.df[key]['%s it=%d tl=0 rl=%d' % (gf, self.it, self._rl.n)] + elif key in self._sd.gf: + df = self._sd.gf[key] + return df['%s it=%d tl=0 rl=%d' % (key, self.it, self._rl.n)] + raise IndexError + + +class RefinementLevel: + # refinement level number + n = -1 + + # [dx, dy, dz] grid spacings + dx = None + # number of ghost gridpoints in [x, y, z] dimensions + ghosts = None + + # the time step + dt = 0.0 + + # a sorted tuple of all the (iteration, time) pairs for which the data + # is present on this refinement level + itertimes = None + + # private + _sd = None + + def __init__(self, sd, n): + self._sd = sd + self.n = n + + def slice(self, iteration = -1, time = -1): + return RefinedSlice(self._sd, self, time, iteration) + + def __str__(self): + return 'Refinement level %d; extent [%g, %g, %g] -- [%g, %g, %g] ([%g, %g, %g] -- [%g, %g, %g]); dx = [%g, %g, %g]; dt = %g' % \ + (self.n, self.x[self.ghosts[0]], self.y[self.ghosts[1]], self.z[self.ghosts[2]], + self.x[-self.ghosts[0] - 1], self.y[-self.ghosts[1] - 1], self.z[-self.ghosts[2] - 1], + self.x[0], self.y[0], self.z[0], self.x[-1], self.y[-1], self.z[-1], + self.dx[0], self.dx[1], self.dx[2], self.dt) + + +class SimulationData: + # directory name + dirname = None + + # datafiles + df = None + + # a dictionary of all the { gridfunction : datafile it is located in } + # pairs for this set of data + gf = None + + # courant factor of the time integration (dx / dt) + courant = 0 + + # per-refinement level parameters + rl = None + + def __init__(self, dirname): + self.dirname = os.path.abspath(dirname) + self.df = {} + self.gf = {} + + # open all the hdf5 files in the dir + for f in os.listdir(dirname): + if not f.endswith('.h5') or f.startswith('checkpoint'): + continue + + self.df[f[:-3]] = h5py.File('%s/%s' % (dirname, f), 'r') + if len(self.df) == 0: + raise ValueError('No HDF5 data files in the directory.') + + for df in self.df.values(): + funcs = df['Parameters and Global Attributes']['Datasets'] + + if funcs.dtype.type != np.string_: + funcs_str = ''.join(map(chr, funcs)).strip('\x00') + else: + funcs_str = funcs.value + + for ds in funcs_str.strip().split(): + if ds in self.gf: + raise ValueError('Gridfunction %s present in more than one datafile: %s and %s' % (ds, self.gf[ds].filename, df.filename)) + self.gf[ds] = df + + # pick a representative datafile and get the grid properties from it + if 'H' in self.df: + df = self.df['H'] + else: + df = self.df.values()[0] + + gf = get_default_gridfunc(df) + + # get the refinement levels, iterations and times + self.rl = [] + while True: + cur_rl = len(self.rl) + try: + ds = df['%s it=0 tl=0 rl=%d' % (gf, cur_rl)] + except KeyError: + break + + self.rl.append(RefinementLevel(self, len(self.rl))) + rl = self.rl[-1] + + rl.itertimes = gf_itertime(df, cur_rl) + rl.dx = ds.attrs['delta'] + rl.ghosts = ds.attrs["cctk_nghostzones"] + rl.dt = _get_timestep(df) / (1 << (len(self.rl) - 1)) + + self.courant = self.rl[0].dx[0] / self.rl[0].dt + + def __del__(self): + if self.df: + map(h5py.File.close, self.df.values()) + + def calc_velocities(self, get_data, rl = 0, t_start = 0, t_end = float('inf'), offsets = None): + rl = self.rl[rl] + dt = rl.itertimes[1][1] - rl.itertimes[0][1] + + ret = [] + for idx, (iter, t) in enumerate(rl.itertimes): + if t < t_start: + continue + if t > t_end: + break + + data1 = get_data(self, iter) + + if offsets is None: + step = int(self.courant) + offsets = xrange(step, 10 * step, step) + + for offset in offsets: + diff_t = offset * dt + try: + data2 = get_data(self, rl.itertimes[idx + offset][0]) + except IndexError: + break + mov = calc_movement(data1, data2) + + for it in mov: + ret.append(it * rl.dx[0] / diff_t) + return ret -- cgit v1.2.3