From 63f30b2bf0de3e811610210531b3de79e7601182 Mon Sep 17 00:00:00 2001 From: Hengrui Zhu Date: Mon, 13 Jan 2025 10:09:02 -0500 Subject: [PATCH 01/72] analysis script --- scripts/cce/athk_to_specter.py | 684 +++++++++++++++++++++++++++++++++ 1 file changed, 684 insertions(+) create mode 100755 scripts/cce/athk_to_specter.py diff --git a/scripts/cce/athk_to_specter.py b/scripts/cce/athk_to_specter.py new file mode 100755 index 000000000..84f5462ea --- /dev/null +++ b/scripts/cce/athk_to_specter.py @@ -0,0 +1,684 @@ +#!/usr/bin/env python3 +## Alireza Rashti - Oct 2024 (C) +## usage: +## $ ./me -h +## + +import sys +import os +import numpy as np +from scipy import special +import math +import argparse +import h5py +import struct +import glob +from scipy.interpolate import interp1d + +# from itertools import product +# import matplotlib.pyplot as plt +# import glob +# import sympy +# import re + +## ---------------------------------------------------------------------- ## + +## field names +g_field_names = [ + "gxx", + "gxy", + "gxz", + "gyy", + "gyz", + "gzz", + "betax", + "betay", + "betaz", + "alp", +] + +g_name_map = { + "gxx": "gxx", + "gxy": "gxy", + "gxz": "gxz", + "gyy": "gyy", + "gyz": "gyz", + "gzz": "gzz", + "betax": "Shiftx", + "betay": "Shifty", + "betaz": "Shiftz", + "alp": "Lapse", +} + +g_name_index = { + "gxx": 4, + "gxy": 5, + "gxz": 6, + "gyy": 7, + "gyz": 8, + "gzz": 9, + "betax": 1, + "betay": 2, + "betaz": 3, + "alp": 0, +} + +## real/imag +g_re = 0 +g_im = 1 + +## args +g_args = None +## various attrs +g_attrs = None + +## debug +g_debug_max_l = 2 + + +def parse_cli(): + """ + arg parser + """ + p = argparse.ArgumentParser( + description="convert Athenak CCE dumps to Spectre CCE") + p.add_argument("-fpath", type=str, required=True, help="/path/to/cce/h5/dumps") + p.add_argument("-ftype", type=str, required=True, help="input file type:{h5,bin}") + p.add_argument("-d_out", type=str, required=True, help="/path/to/output/dir") + p.add_argument("-debug", type=str, default="n", help="debug=[y,n]") + p.add_argument( + "-radius", + type=float, + required=True, + help="interpolate all fields and their derivatives at this radius.", + ) + p.add_argument( + "-t_deriv", + type=str, + default="Fourier", + help="method to take the time derivative of fields:{Fourier,fin_diff}", + ) + p.add_argument( + "-r_deriv", + type=str, + default="ChebU", + help= + "method to take the radial derivative of fields:{ChebU:Chebyshev of second kind}", + ) + p.add_argument( + "-interpolation", + type=str, + default="ChebU", + help= + "method to interpolate fields at a given r:{ChebU:Chebyshev of second kind}", + ) + + args = p.parse_args() + return args + + +def load(fpath: str, field_name: str, attrs: dict) -> list: + """ + read the field accroding to attrs. + return convention: + ret[real/imag, time_level, n, lm], eg: + ret[g_re,3,2,:] = Re(C_2lm(t=3)) for all lm + ret[g_im,3,2,:] = Im(C_2lm(t=3)) for all lm + """ + + if attrs["file_type"] == "h5": + lev_t = attrs["lev_t"] + max_n = attrs["max_n"] + max_lm = attrs["max_lm"] + shape = (len([g_re, g_im]), lev_t, max_n, max_lm) + ret = np.empty(shape=shape, dtype=float) + with h5py.File(fpath, "r") as h5f: + # read & save + for i in range(0, lev_t): + key = f"{i}" + h5_re = h5f[f"{key}/{field_name}/re"] + h5_im = h5f[f"{key}/{field_name}/im"] + ret[g_re, i, :] = h5_re + ret[g_im, i, :] = h5_im + elif attrs["file_type"] == "bin": + lev_t = attrs["lev_t"] + max_n = attrs["max_n"] + max_lm = attrs["max_lm"] + shape = (len([g_re, g_im]), lev_t, max_n, max_lm) + ret = np.empty(shape=shape, dtype=float) + # Load the list of files + flist = sorted(glob.glob(fpath+'/cce_*.bin')) + data = [] + dat_real = [] + dat_imag = [] + t = [] + for f in flist: + nr, num_l_modes, time, rin, rout, data_real, data_imag, index_to_lm = read_cce_file(f) + t.append(time) + dat_real.append(data_real[:,g_name_index[field_name],:]) + dat_imag.append(data_imag[:,g_name_index[field_name],:]) + dat_real = np.array(dat_real) + dat_imag = np.array(dat_imag) + t = np.array(t) + f_real = interp1d(t, dat_real,axis=0) + f_imag = interp1d(t, dat_imag,axis=0) + + ret[g_re, :] = f_real(attrs["time"]) + ret[g_im, :] = f_imag(attrs["time"]) + else: + raise ValueError("no such option") + + # print(ret) + return ret + + +def lm_mode(l, m): + """ + l and m mode convention + """ + return l * l + l + m + + +def read_cce_file(filename): + """ + Reads binary data from a CCE output file generated by the provided C++ function. + + Parameters: + filename (str): The path to the binary file. + + Returns: + nr (int): Number of radial points. + num_l_modes (int): Number of l modes. + time (float): Time value from the simulation. + rin (float): Inner radial boundary. + rout (float): Outer radial boundary. + data_real (numpy.ndarray): Real part of the data, shape (nr, 10, num_angular_modes). + data_imag (numpy.ndarray): Imaginary part of the data, shape (nr, 10, num_angular_modes). + index_to_lm (dict): Mapping from data index to (l, m) values. + """ + with open(filename, 'rb') as f: + # Read number of radial points (nr) + nr_bytes = f.read(4) + nr = struct.unpack(' dict: + """ + find attributes such as num. of time level, and n, lm in C_nlm + also saves the time value at each slice. + """ + attrs = {} + if type == "h5": + attrs["file_type"] = "h5" + with h5py.File(fpath, "r") as h5f: + # find attribute about num. of time level, and n,l,m in C_nlm + attrs["lev_t"] = len(h5f.keys()) - 1 + attrs["max_n"], attrs["max_lm"] = h5f[f"1/{field_name}/re"].shape + attrs["max_l"] = int(math.sqrt(attrs["max_lm"])) - 1 + attrs["r_in"] = h5f["metadata"].attrs["Rin"] + attrs["r_out"] = h5f["metadata"].attrs["Rout"] + # read & save time + time = [] + for i in range(0, attrs["lev_t"]): + key = f"{i}" + t = h5f[key].attrs["Time"][0] + time.append(t) + attrs["time"] = np.array(time) + elif type == "bin": + attrs["file_type"] = "bin" + # Load the list of files + flist = sorted(glob.glob(fpath+'/cce_*.bin')) + t = [] + data = [] + for f in flist: + nr, num_l_modes, time, rin, rout, data_real, data_imag, index_to_lm = read_cce_file(f) + t.append(time) + #data.append(data_real + 1j * data_imag) + t = np.array(t) + #data = np.array(data) + attrs["lev_t"] = len(flist) + attrs["max_n"] = nr + attrs["max_l"] = num_l_modes + attrs["max_lm"] = (num_l_modes + 1) ** 2 + attrs["r_in"] = np.array([rin]) + attrs["r_out"] = np.array([rout]) + attrs["time"] = np.linspace(t.min(),t.max(),t.shape[0]) + else: + raise ValueError("no such option") + + # print(attrs) + return attrs + + +def time_derivative_findiff(field: np.array, field_name: str, attrs: dict, + args) -> np.array: + """ + return the time derivative of the given field using finite diff. 2nd order + field(t,rel/img,n,lm) + """ + + print(f"finite difference time derivative: {field_name}", flush=True) + _, len_t, len_n, len_lm = field.shape + time = attrs["time"] + dt = np.gradient(time, 2) + dfield = np.empty_like(field) + + for n in range(len_n): + for lm in range(len_lm): + dfield[g_re, :, n, lm] = np.gradient(field[g_re, :, n, lm], 2) / dt + dfield[g_im, :, n, lm] = np.gradient(field[g_im, :, n, lm], 2) / dt + + return dfield + + +def time_derivative_fourier(field: np.array, field_name: str, attrs: dict, + args) -> np.array: + """ + return the time derivative of the given field using Fourier method + field(t,rel/img,n,lm) + """ + + print(f"Fourier time derivative: {field_name}", flush=True) + _, len_t, len_n, len_lm = field.shape + dt = attrs["time"][2] - attrs["time"][1] + wm = math.pi * 2.0 / (len_t * dt) + + dfield = np.empty_like(field) + for n in range(len_n): + for lm in range(len_lm): + coeff = field[g_re, :, n, lm] + 1j * field[g_im, :, n, lm] + # F. transform + fft_coeff = np.fft.fft(coeff) + # if args["debug"] == 'y': + # print("debug: normalization?",round(coeff[1],6) == round(np.fft.ifft(fft_coeff)[1],6)) + + # time derivative + half = len_t // 2 + 1 + omega = np.empty(shape=half) + for i in range(0, half): + omega[i] = i * wm + + dfft_coeff = np.empty_like(fft_coeff) + dfft_coeff[0] = 0 + dfft_coeff[1:half] = (-np.imag(fft_coeff[1:half]) + + 1j * np.real(fft_coeff[1:half])) * omega[1:] + dfft_coeff[half:] = (np.imag(fft_coeff[half:]) - + 1j * np.real(fft_coeff[half:])) * omega[::-1][1:] + + # not optimized version + """ + dfft_coeff[0] = 0 + for i in range(1, half): + omega = i * wm + re = np.real(fft_coeff[i]) + im = np.imag(fft_coeff[i]) + re2 = np.real(fft_coeff[-i]) + im2 = np.imag(fft_coeff[-i]) + + dfft_coeff[i] = omega*complex(-im, re) + dfft_coeff[-i] = omega*complex(im2, -re2) + + """ + # F. inverse + coeff = np.fft.ifft(dfft_coeff) + dfield[g_re, :, n, lm] = np.real(coeff) + dfield[g_im, :, n, lm] = np.imag(coeff) + + if args["debug"] == "y": + for n in range(len_n): + for l in range(2, g_debug_max_l + 1): + for m in range(-l, l + 1): + hfile = (f"{args['d_out']}/debug_{field_name}_n{n}l{l}m{m}.txt") + write_data = np.column_stack(( + attrs["time"], + dfield[g_re, :, n, lm_mode(l, m)], + dfield[g_im, :, n, lm_mode(l, m)], + field[g_re, :, n, lm_mode(l, m)], + field[g_im, :, n, lm_mode(l, m)], + )) + np.savetxt(hfile, write_data, header="t dre/dt dim/dt re im") + + return dfield + + +def dUk_dx(order: int, x: float) -> float: + """ + d(Chebyshev of second kind)/dx + """ + assert x != 1 and x != -1 + t = special.chebyt(order + 1)(x) + u = special.chebyu(order)(x) + duk_dx = (order + 1) * t - x * u + duk_dx /= x**2 - 1 + + return duk_dx + + +def radial_derivative_at_r_chebu(field: np.array, + field_name: str, + attrs: dict, + args) -> np.array: + """ + return the radial derivative of the given field using Chebyshev of + 2nd kind method at the radius of interest. + + f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind + collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) + x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} + + field(t,rel/img,n,lm) + """ + + print(f"ChebyU radial derivative: {field_name}", flush=True) + _, len_t, len_n, len_lm = field.shape + + r_1 = attrs["r_in"][0] + r_2 = attrs["r_out"][0] + assert r_1 != r_2 + dx_dr = 2 / (r_2 - r_1) + + if args["debug"] == "y": + # populate collocation points, roots of U_i + x_i = np.empty(shape=len_n, dtype=float) + for i in range(len_n): + x_i[i] = math.cos(math.pi * (i + 1) / (len_n + 1)) + + # dU_k/dx|x=x_i + duk_dx = np.empty(shape=(len_n, len_n), dtype=float) + for k in range(len_n): + for i in range(len_n): + t = special.chebyt(k + 1)(x_i[i]) + u = special.chebyu(k)(x_i[i]) + duk_dx[k, i] = (k + 1) * t - x_i[i] * u + + duk_dx /= np.square(x_i) - 1 + + uk = np.empty(shape=len_n, dtype=float) + tk = np.empty(shape=len_n, dtype=float) + for k in range(len_n): + hfile = f"{args['d_out']}/cheb_k{k}.txt" + for i in range(len_n): + tk[i] = special.chebyt(k)(x_i[i]) + uk[i] = special.chebyu(k)(x_i[i]) + + write_data = np.column_stack((x_i, uk, tk, duk_dx[k, :])) + np.savetxt( + hfile, + write_data, + header=f"x_i uk{k}(x_i) tk{k}(x_i) duk{k}(x_i)/dx", + ) + + dfield = np.zeros(shape=(len([g_re, g_im]), len_t, len_lm)) + r = args["radius"] + x = (2 * r - r_1 - r_2) / (r_2 - r_1) + for k in range(len_n): + dfield[:, :, :] += field[:, :, k, :] * dUk_dx(k, x) + + return dfield * dx_dr + + +def time_derivative(field: np.array, field_name: str, attrs: dict, args): + """ + return the time derivative of the given field + """ + + if args["t_deriv"] == "Fourier": + return time_derivative_fourier(field, field_name, attrs, args) + elif args["t_deriv"] == "fin_diff": + return time_derivative_findiff(field, field_name, attrs, args) + else: + raise ValueError("no such option") + + +def radial_derivative_at_r(field: np.array, field_name: str, attrs: dict, args): + """ + return the radial derivative of the given field at R=r + """ + + if args["r_deriv"] == "ChebU": + return radial_derivative_at_r_chebu(field, field_name, attrs, args) + else: + raise ValueError("no such option") + + +class Interpolate_at_r: + + def __init__(self, attrs: dict, args: dict): + """ + interpolate the given field at R=r + """ + self.attrs = attrs + self.args = args + self.len_t = attrs["lev_t"] + self.len_n = attrs["max_n"] + self.len_lm = attrs["max_lm"] + r_1 = attrs["r_in"][0] + r_2 = attrs["r_out"][0] + self.r = r = args["radius"] + self.x = (2 * r - r_1 - r_2) / (r_2 - r_1) + + if args["interpolation"] == "ChebU": + self.Uk = np.empty(shape=self.len_n) + for k in range(self.len_n): + self.Uk[k] = special.chebyu(k)(self.x) + self.interp = self.interpolate_at_r_chebu + else: + raise ValueError("no such option") + + def interpolate_at_r_chebu(self, field: np.array, field_name: str): + """ + interpolate at R=r using Cheb U. + """ + print(f"Interpolating at R={self.r}: {field_name}", flush=True) + + field_r = np.zeros(shape=(len([g_re, g_im]), self.len_t, self.len_lm)) + for k in range(self.len_n): + field_r[:, :, :] += field[:, :, k, :] * self.Uk[k] + + return field_r + + def interpolate(self, field: np.array, field_name: str): + return self.interp(field, field_name) + + +def process_field(field_name: str) -> dict: + """ + - read data + - find time derives + - find radial derives + - interpolate at R=r + """ + + # return + attrs = g_attrs + args = g_args + db = {} + + # load data + field = load(args["fpath"], field_name, attrs) + # db[f"{field_name}"] = field + + # time derivative + dfield_dt = time_derivative(field, field_name, attrs, args) + + # interpolate at a specific radii + interpolate = Interpolate_at_r(attrs, args) + field_at_r = interpolate.interpolate(field, field_name) + db[f"{field_name}|r"] = field_at_r + + dfield_dt_at_r = interpolate.interpolate(dfield_dt, f"d{field_name}/dt") + db[f"d{field_name}/dt|r"] = dfield_dt_at_r + + # radial derivative at R=r + dfield_dr_at_r = radial_derivative_at_r(field, field_name, attrs, args) + db[f"d{field_name}/dr|r"] = dfield_dr_at_r + + return db + + +def h5_create_group(h5file, group_name: str): + """ + create a group for h5 + """ + h5group = None + + # create group if not exists + if h5file.get(group_name, default=None) == None: + h5group = h5file.create_group(group_name) + else: + raise ValueError("this group {group_name} is already exists.") + + return h5group + + +def h5_write_data(h5file, + data: np.array, + data_name: str, + attrs: dict, + args: dict): + """ + reminder: + data[real/imag, time_level, lm] + + write syntax, eg: + + h5["gxx.dat"] = + [time_level, ['time', 'gxx_Re(0,0)', 'gxx_Im(0,0)', 'gxx_Re(1,1)', 'gxx_Im(1,1)', ...] ] + + h5["gxx.dat"].attrs['Legend'] = the associated column = + array(['time', 'gxx_Re(0,0)', 'gxx_Im(0,0)', 'gxx_Re(1,1)', 'gxx_Im(1,1)', ...]) + + # => h5["gxx.dat"][3,0] = value of time at the dump level 3 + # => h5["gxx.dat"][4,1] = value of gxx_Re(0,0) at the dump level 4 + + """ + + dataset_conf = dict( + name=f"{data_name}", + shape=(attrs["lev_t"], len([g_re, g_im]) * (attrs["max_l"]**2) + 1), + dtype=float, # chunks=True, + # compression="gzip", + # shuffle=True, + ) + + data_attrs = ["time"] + + if args["debug"] == "y": + print(dataset_conf, flush=True) + + h5file.create_dataset(**dataset_conf) + + flat = 0 + h5file[f"{data_name}"][:, flat] = attrs["time"] + flat += 1 + for l in range(0, attrs["max_l"]): + for m in range(l, -l - 1,-1): + data_attrs.append(f"{data_name[:-4]}_Re({l},{m})") + data_attrs.append(f"{data_name[:-4]}_Im({l},{m})") + h5file[f"{data_name}"][:, flat] = data[g_re, :, lm_mode(l, m)] + h5file[f"{data_name}"][:, flat + 1] = data[g_im, :, lm_mode(l, m)] + flat += 2 + Legend = [s.encode('ascii', 'ignore') for s in data_attrs] + h5file[f"{data_name}"].attrs["Legend"] = Legend + + +def write(f: str, db: dict, attrs: dict, args: dict): + """ + write data on disk + """ + print(f"writing: {f}", flush=True) + + field_name = g_name_map[f"{f}"] + field_name_key = f"{field_name}.dat" + dfield_name_dr_key = f"Dr{field_name}.dat" + dfield_name_dt_key = f"Dt{field_name}.dat" + + r = args["radius"] + file_name = os.path.join(args["d_out"], f"CceR{r:07.2f}.h5") + with h5py.File(file_name, "a") as h5file: + + name = field_name_key + data = db[f"{f}|r"] + h5_write_data(h5file, data, name, attrs, args) + + name = dfield_name_dr_key + data = db[f"d{f}/dr|r"] + h5_write_data(h5file, data, name, attrs, args) + + name = dfield_name_dt_key + data = db[f"d{f}/dt|r"] + h5_write_data(h5file, data, name, attrs, args) + + +def main(args): + """ + create output required by Specter code + ref: https://spectre-code.org/tutorial_cce.html + """ + + global g_attrs + global g_args + + g_args = args + + # check if output dir exist, if not, mkdir + if not os.path.exists(g_args["d_out"]): + os.makedirs(g_args["d_out"]) + # find attribute for an arbitrary field + g_attrs = get_attribute(args["fpath"],type=args["ftype"]) + + # for each field + # I'm afraid this method takes too much memory + # from multiprocessing import Pool + # with Pool(processes=len(g_field_names)) as p: + # db = p.map(process_field, g_field_names) + + for f in g_field_names: + db = process_field(f) + # write on disk + write(f, db, g_attrs, g_args) + + +if __name__ == "__main__": + args = parse_cli() + main(args.__dict__) From 74b341705c8e89b1988f1a8c6be6fd4b1e8c3f97 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 10:35:46 -0500 Subject: [PATCH 02/72] feat: nan checker. --- scripts/cce/athk_to_specter.py | 4 ++++ scripts/cce/cce_run.sh | 13 +++++++++++++ scripts/cce/pre_run.sh | 7 +++++++ 3 files changed, 24 insertions(+) create mode 100755 scripts/cce/cce_run.sh create mode 100755 scripts/cce/pre_run.sh diff --git a/scripts/cce/athk_to_specter.py b/scripts/cce/athk_to_specter.py index 84f5462ea..57933554c 100755 --- a/scripts/cce/athk_to_specter.py +++ b/scripts/cce/athk_to_specter.py @@ -613,6 +613,10 @@ def h5_write_data(h5file, flat += 1 for l in range(0, attrs["max_l"]): for m in range(l, -l - 1,-1): + if args["debug"] == "y": + assert all(data[g_re, :, lm_mode(l, m)] != np.nan) + assert all(data[g_im, :, lm_mode(l, m)] != np.nan) + data_attrs.append(f"{data_name[:-4]}_Re({l},{m})") data_attrs.append(f"{data_name[:-4]}_Im({l},{m})") h5file[f"{data_name}"][:, flat] = data[g_re, :, lm_mode(l, m)] diff --git a/scripts/cce/cce_run.sh b/scripts/cce/cce_run.sh new file mode 100755 index 000000000..5037849f2 --- /dev/null +++ b/scripts/cce/cce_run.sh @@ -0,0 +1,13 @@ +#!/usr/bin/env bash + +cd runs + +echo "running preproc...\n" +./PreprocessCceWorldtube --input-file PreprocessCceWorldtube_test.yaml + +echo "-------------------------" + + +echo "running CCE..." +./CharacteristicExtract --input-file cce_par.yaml + diff --git a/scripts/cce/pre_run.sh b/scripts/cce/pre_run.sh new file mode 100755 index 000000000..beae51279 --- /dev/null +++ b/scripts/cce/pre_run.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env bash + +./athk_to_specter.py -ftype bin -d_out ./debug -debug y \ +-t_deriv fin_diff -r_deriv ChebU -interpolation ChebU \ +-radius 30 \ +-fpath ./dat/bin/cce_bbh/ + From d97f67752b66b5d5bd273e5d41188fe89af0bd3c Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 10:39:08 -0500 Subject: [PATCH 03/72] fix: name --- scripts/cce/{athk_to_specter.py => athk_to_spectre.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename scripts/cce/{athk_to_specter.py => athk_to_spectre.py} (100%) diff --git a/scripts/cce/athk_to_specter.py b/scripts/cce/athk_to_spectre.py similarity index 100% rename from scripts/cce/athk_to_specter.py rename to scripts/cce/athk_to_spectre.py From 0c8e36297f82f3dee0c4c65ae3d29e79642e18f5 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 11:32:19 -0500 Subject: [PATCH 04/72] feat: add debug script. --- scripts/cce/debug_athk_to_spectre.py | 120 +++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100755 scripts/cce/debug_athk_to_spectre.py diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py new file mode 100755 index 000000000..229e815bc --- /dev/null +++ b/scripts/cce/debug_athk_to_spectre.py @@ -0,0 +1,120 @@ +#!/usr/bin/env python3 +## Alireza Rashti - Jan 2025 (C) +## usage: +## $ ./me -h +## + +import sys +import os +import numpy as np +import math as m +import argparse +import h5py + +# import matplotlib.pyplot as plt +# import glob +# import sympy +## ---------------------------------------------------------------------- ## + +g_name_map = { + "gxx": "gxx", + "gxy": "gxy", + "gxz": "gxz", + "gyy": "gyy", + "gyz": "gyz", + "gzz": "gzz", + "betax": "Shiftx", + "betay": "Shifty", + "betaz": "Shiftz", + "alp": "Lapse", +} + + +def parse_cli(): + """ + arg parser + """ + p = argparse.ArgumentParser(description="debugging athk_to_spectre.py") + p.add_argument( + "-debug", + type=str, + default="plot_simple", + help="debug type=[plot_simple]", + ) + p.add_argument( + "-fpath", + type=str, + required=True, + help="path/to/output/athk_to_spectre.py/h5", + ) + p.add_argument( + "-field_name", + type=str, + default="gxx", + help="plot for this field [gxx,gxy,...]", + ) + p.add_argument( + "-field_mode", + type=str, + default="Re_(2,2)", + help="plot this mode[Re_(l,m),Im_(l,m)]", + ) + + args = p.parse_args() + return args + + +def find_h5_mode(h5f, field_name, mode_name, args): + mode = 0 + flag = False + for m in h5f[field_name].attrs["Legend"]: + if m.find(mode_name) != -1: + flag = True + break + mode += 1 + + assert flag == True + return mode + + +def read_h5_mode_and_derivs(args): + + field_name = g_name_map[f"{field_name}"] + field_name_key = f"{field_name}.dat" + dfield_name_dr_key = f"Dr{field_name}.dat" + dfield_name_dt_key = f"Dt{field_name}.dat" + + with h5py.File(args["fpath"], "r") as h5f: + h5file[f"{data_name}"].attrs["Legend"] = Legend + + mode = find_h5_mode(h5f, f"{field_name_key}", args["field_mode"], args) + f = h5f[f"{field_name_key}"][:, mode] + + mode = find_h5_mode(h5f, f"{dfield_name_dr_key}", args["field_mode"], args) + drf = h5f[f"{dfield_name_dr_key}"][:, mode] + + mode = find_h5_mode(h5f, f"{dfield_name_dt_key}", args["field_mode"], args) + dtf = h5f[f"{dfield_name_dt_key}"][:, mode] + + return f, drf, dtf + + +def debug_plot_simple(args): + dat = read_h5_mode_and_derivs(args) + plot_simple_live(dat, args) + + +def main(args): + """ + debug + """ + + if args["debug"] == "plot_simple": + debug_plot_simple(args) + else: + raise ValueError("no such option") + + +if __name__ == "__main__": + args = parse_cli() + main(args.__dict__) From be07710fee6eb78335f4306435c2546d5055058f Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 11:36:25 -0500 Subject: [PATCH 05/72] fix: bug. --- scripts/cce/debug_athk_to_spectre.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py index 229e815bc..e9904f764 100755 --- a/scripts/cce/debug_athk_to_spectre.py +++ b/scripts/cce/debug_athk_to_spectre.py @@ -56,8 +56,8 @@ def parse_cli(): p.add_argument( "-field_mode", type=str, - default="Re_(2,2)", - help="plot this mode[Re_(l,m),Im_(l,m)]", + default="Re(2,2)", + help="plot this mode[Re(l,m),Im(l,m)]", ) args = p.parse_args() @@ -69,6 +69,7 @@ def find_h5_mode(h5f, field_name, mode_name, args): flag = False for m in h5f[field_name].attrs["Legend"]: if m.find(mode_name) != -1: + print("found mode for", field_name, m,mode_name) flag = True break mode += 1 @@ -79,14 +80,12 @@ def find_h5_mode(h5f, field_name, mode_name, args): def read_h5_mode_and_derivs(args): - field_name = g_name_map[f"{field_name}"] + field_name = g_name_map[args["field_name"]] field_name_key = f"{field_name}.dat" dfield_name_dr_key = f"Dr{field_name}.dat" dfield_name_dt_key = f"Dt{field_name}.dat" with h5py.File(args["fpath"], "r") as h5f: - h5file[f"{data_name}"].attrs["Legend"] = Legend - mode = find_h5_mode(h5f, f"{field_name_key}", args["field_mode"], args) f = h5f[f"{field_name_key}"][:, mode] From 21680c224f1214cc8ff394e0d7866a32b1bda0ec Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 12:12:13 -0500 Subject: [PATCH 06/72] feat: plot --- scripts/cce/debug_athk_to_spectre.py | 56 ++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 4 deletions(-) diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py index e9904f764..d828fa992 100755 --- a/scripts/cce/debug_athk_to_spectre.py +++ b/scripts/cce/debug_athk_to_spectre.py @@ -10,8 +10,8 @@ import math as m import argparse import h5py +import matplotlib.pyplot as plt -# import matplotlib.pyplot as plt # import glob # import sympy ## ---------------------------------------------------------------------- ## @@ -87,6 +87,7 @@ def read_h5_mode_and_derivs(args): with h5py.File(args["fpath"], "r") as h5f: mode = find_h5_mode(h5f, f"{field_name_key}", args["field_mode"], args) + t = h5f[f"{field_name_key}"][:, 0] f = h5f[f"{field_name_key}"][:, mode] mode = find_h5_mode(h5f, f"{dfield_name_dr_key}", args["field_mode"], args) @@ -95,13 +96,60 @@ def read_h5_mode_and_derivs(args): mode = find_h5_mode(h5f, f"{dfield_name_dt_key}", args["field_mode"], args) dtf = h5f[f"{dfield_name_dt_key}"][:, mode] - return f, drf, dtf + return (f, drf, dtf, t) +def plot_simple_v_t(dat, args): + """ + plot value vs time + """ + fig, axes = plt.subplots(3, 1, sharex=True) + + # f + ax = axes[0] + label=args["field_name"] + conf = dict(ls="-", label=label, color="k") + ax.plot(dat[-1], dat[0], **conf) + ax.set_ylabel(label) + ax.grid(True) + + # drf + ax = axes[1] + label="d"+args["field_name"] + "/dr" + conf = dict(ls="-", label=label, color="k") + ax.plot(dat[-1], dat[1], **conf) + ax.set_ylabel(label) + ax.grid(True) + + # dtf + ax = axes[2] + label="d"+args["field_name"] + "/dt" + conf = dict(ls="-", label=label, color="k") + ax.plot(dat[-1], dat[2], **conf) + ax.grid(True) + ax.set_ylabel(label) + ax.set_xlabel("t/M") + + plt.tight_layout() + plt.show() + + """ + file_out = args["o"] + if args["debug"] == "y": + plt.show() + else: + plt.savefig(file_out) + os.system("pdfcrop " + file_out + " " + file_out) + """ + + + + + + def debug_plot_simple(args): dat = read_h5_mode_and_derivs(args) - plot_simple_live(dat, args) - + plot_simple_v_t(dat, args) def main(args): """ From 4641d85bff17b54aba5e23a909b1277f0d18169d Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 12:27:22 -0500 Subject: [PATCH 07/72] fix. --- scripts/cce/debug_athk_to_spectre.py | 42 +++++++++++++++------------- scripts/cce/pre_run.sh | 7 ++++- 2 files changed, 28 insertions(+), 21 deletions(-) diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py index d828fa992..fc463bf86 100755 --- a/scripts/cce/debug_athk_to_spectre.py +++ b/scripts/cce/debug_athk_to_spectre.py @@ -47,6 +47,12 @@ def parse_cli(): required=True, help="path/to/output/athk_to_spectre.py/h5", ) + p.add_argument( + "-dout", + type=str, + default="./", + help="path/to/output/dir", + ) p.add_argument( "-field_name", type=str, @@ -69,7 +75,7 @@ def find_h5_mode(h5f, field_name, mode_name, args): flag = False for m in h5f[field_name].attrs["Legend"]: if m.find(mode_name) != -1: - print("found mode for", field_name, m,mode_name) + print("found mode for", field_name, m, mode_name) flag = True break mode += 1 @@ -102,12 +108,12 @@ def read_h5_mode_and_derivs(args): def plot_simple_v_t(dat, args): """ plot value vs time - """ + """ fig, axes = plt.subplots(3, 1, sharex=True) - + # f ax = axes[0] - label=args["field_name"] + label = args["field_name"] conf = dict(ls="-", label=label, color="k") ax.plot(dat[-1], dat[0], **conf) ax.set_ylabel(label) @@ -115,7 +121,7 @@ def plot_simple_v_t(dat, args): # drf ax = axes[1] - label="d"+args["field_name"] + "/dr" + label = "d" + args["field_name"] + "/dr" conf = dict(ls="-", label=label, color="k") ax.plot(dat[-1], dat[1], **conf) ax.set_ylabel(label) @@ -123,34 +129,30 @@ def plot_simple_v_t(dat, args): # dtf ax = axes[2] - label="d"+args["field_name"] + "/dt" + label = "d" + args["field_name"] + "/dt" conf = dict(ls="-", label=label, color="k") ax.plot(dat[-1], dat[2], **conf) ax.grid(True) ax.set_ylabel(label) ax.set_xlabel("t/M") - - plt.tight_layout() - plt.show() - - """ - file_out = args["o"] - if args["debug"] == "y": - plt.show() - else: - plt.savefig(file_out) - os.system("pdfcrop " + file_out + " " + file_out) - """ + plt.tight_layout() + # plt.show() - + mode = args["field_mode"][3:-1] + lm = mode.split(",") + file_out = os.path.join( + args["dout"], + args["field_name"] + "_" + f"l{lm[0]}m{lm[1]}" + "_vs_time.png", + ) + plt.savefig(file_out, dpi=200) - def debug_plot_simple(args): dat = read_h5_mode_and_derivs(args) plot_simple_v_t(dat, args) + def main(args): """ debug diff --git a/scripts/cce/pre_run.sh b/scripts/cce/pre_run.sh index beae51279..fb0a23594 100755 --- a/scripts/cce/pre_run.sh +++ b/scripts/cce/pre_run.sh @@ -1,7 +1,12 @@ #!/usr/bin/env bash -./athk_to_specter.py -ftype bin -d_out ./debug -debug y \ +if false +then +./athk_to_spectre.py -ftype bin -d_out ./debug -debug y \ -t_deriv fin_diff -r_deriv ChebU -interpolation ChebU \ -radius 30 \ -fpath ./dat/bin/cce_bbh/ +fi +# plot fig +./debug_athk_to_spectre.py -debug plot_simple -fpath ./debug/CceR0030.00.h5 -dout ./debug From c618d78fdbf7045a516d6405dd31865e6944f558 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 16 Jan 2025 12:32:00 -0500 Subject: [PATCH 08/72] feat: add plot script. --- scripts/cce/pre_run.sh | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/scripts/cce/pre_run.sh b/scripts/cce/pre_run.sh index fb0a23594..098e7d536 100755 --- a/scripts/cce/pre_run.sh +++ b/scripts/cce/pre_run.sh @@ -8,5 +8,34 @@ then -fpath ./dat/bin/cce_bbh/ fi +file="./debug/CceR0030.00.h5" # plot fig -./debug_athk_to_spectre.py -debug plot_simple -fpath ./debug/CceR0030.00.h5 -dout ./debug +fname="gxx" +mode="Re(0,0)" +./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ + -fpath $file -field_name $fname -field_mode $mode + +mode="Re(2,2)" +./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ + -fpath $file -field_name $fname -field_mode $mode + + +fname="alp" +mode="Re(0,0)" +./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ + -fpath $file -field_name $fname -field_mode $mode + +mode="Re(2,2)" +./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ + -fpath $file -field_name $fname -field_mode $mode + + +fname="betax" +mode="Re(0,0)" +./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ + -fpath $file -field_name $fname -field_mode $mode + +mode="Re(2,2)" +./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ + -fpath $file -field_name $fname -field_mode $mode + From a7d66a10a8521b0177e3dbbb3e17509d28627a91 Mon Sep 17 00:00:00 2001 From: Alireza Date: Fri, 17 Jan 2025 08:26:04 -0500 Subject: [PATCH 09/72] feat: add mode convergence. --- scripts/cce/debug_athk_to_spectre.py | 104 +++++++++++++++++++++++++-- 1 file changed, 97 insertions(+), 7 deletions(-) diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py index fc463bf86..4fd110154 100755 --- a/scripts/cce/debug_athk_to_spectre.py +++ b/scripts/cce/debug_athk_to_spectre.py @@ -66,11 +66,18 @@ def parse_cli(): help="plot this mode[Re(l,m),Im(l,m)]", ) + p.add_argument( + "-time_dump", + type=int, + default=100, + help="how often dump for mode convergence", + ) + args = p.parse_args() return args -def find_h5_mode(h5f, field_name, mode_name, args): +def find_h5_1mode(h5f, field_name, mode_name, args): mode = 0 flag = False for m in h5f[field_name].attrs["Legend"]: @@ -84,7 +91,32 @@ def find_h5_mode(h5f, field_name, mode_name, args): return mode -def read_h5_mode_and_derivs(args): +def find_h5_all_modes(h5f, field_name, mode_name, args): + modes = [] + names = [] + re_or_im = mode_name[0:2] ## Re(...) + print(re_or_im) + assert re_or_im == "Re" or re_or_im == "Im" + + legends = h5f[field_name].attrs["Legend"] + for i in range(len(legends)): + m = legends[i] + if m.find(re_or_im) != -1: + names.append(m) + modes.append(i) + + assert len(modes) + + # print(modes) + # print(names) + + return modes + + +def read_h5_1mode(args): + """ + return t, f(t)|mode, df(t)/dr|mod, df(t)/dt|mode + """ field_name = g_name_map[args["field_name"]] field_name_key = f"{field_name}.dat" @@ -92,19 +124,72 @@ def read_h5_mode_and_derivs(args): dfield_name_dt_key = f"Dt{field_name}.dat" with h5py.File(args["fpath"], "r") as h5f: - mode = find_h5_mode(h5f, f"{field_name_key}", args["field_mode"], args) + mode = find_h5_1mode(h5f, f"{field_name_key}", args["field_mode"], args) t = h5f[f"{field_name_key}"][:, 0] f = h5f[f"{field_name_key}"][:, mode] - mode = find_h5_mode(h5f, f"{dfield_name_dr_key}", args["field_mode"], args) + mode = find_h5_1mode(h5f, f"{dfield_name_dr_key}", args["field_mode"], args) drf = h5f[f"{dfield_name_dr_key}"][:, mode] - mode = find_h5_mode(h5f, f"{dfield_name_dt_key}", args["field_mode"], args) + mode = find_h5_1mode(h5f, f"{dfield_name_dt_key}", args["field_mode"], args) dtf = h5f[f"{dfield_name_dt_key}"][:, mode] return (f, drf, dtf, t) +def read_h5_all_modes(args): + """ + return t, f(mode)|t, df(mode)/dr|t, df(mode)/dt|t + + f[t_i,modes] = at t_i what are the modes + """ + + field_name = g_name_map[args["field_name"]] + field_name_key = f"{field_name}.dat" + dfield_name_dr_key = f"Dr{field_name}.dat" + dfield_name_dt_key = f"Dt{field_name}.dat" + + with h5py.File(args["fpath"], "r") as h5f: + t = h5f[f"{field_name_key}"][:, 0] + n_dumps = len(t) // args["time_dump"] + n_dumps = n_dumps + 1 if len(t) % args["time_dump"] else n_dumps + + modes = find_h5_all_modes(h5f, + f"{field_name_key}", + args["field_mode"], + args) + f = np.empty(shape=(n_dumps, len(modes))) + k = 0 # note: f[k,.] is 0,1,... but i jumps + for i in range(0, len(t), args["time_dump"]): + for j in range(len(modes)): + f[k, j] = h5f[f"{field_name_key}"][i, modes[j]] + k += 1 + + modes = find_h5_all_modes(h5f, + f"{dfield_name_dr_key}", + args["field_mode"], + args) + drf = np.empty(shape=(n_dumps, len(modes))) + k = 0 + for i in range(0, len(t), args["time_dump"]): + for j in range(len(modes)): + drf[k, j] = h5f[f"{dfield_name_dr_key}"][i, modes[j]] + k += 1 + + modes = find_h5_all_modes(h5f, + f"{dfield_name_dt_key}", + args["field_mode"], + args) + dtf = np.empty(shape=(n_dumps, len(modes))) + k = 0 + for i in range(0, len(t), args["time_dump"]): + for j in range(len(modes)): + dtf[k, j] = h5f[f"{dfield_name_dt_key}"][i, modes[j]] + k += 1 + + return (f, drf, dtf, t) + + def plot_simple_v_t(dat, args): """ plot value vs time @@ -149,8 +234,13 @@ def plot_simple_v_t(dat, args): def debug_plot_simple(args): - dat = read_h5_mode_and_derivs(args) - plot_simple_v_t(dat, args) + # value vs time + # dat = read_h5_1mode(args) + # plot_simple_v_t(dat, args) + + # conv test + dat = read_h5_all_modes(args) + # plot_simple_mode(dat, args) def main(args): From 54791157c9696cb05752f39ec89a92d508c85199 Mon Sep 17 00:00:00 2001 From: Alireza Date: Fri, 17 Jan 2025 08:52:48 -0500 Subject: [PATCH 10/72] fix: plots. --- scripts/cce/debug_athk_to_spectre.py | 85 +++++++++++++++++++++++++--- 1 file changed, 77 insertions(+), 8 deletions(-) diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py index 4fd110154..5f3cc09f3 100755 --- a/scripts/cce/debug_athk_to_spectre.py +++ b/scripts/cce/debug_athk_to_spectre.py @@ -93,16 +93,16 @@ def find_h5_1mode(h5f, field_name, mode_name, args): def find_h5_all_modes(h5f, field_name, mode_name, args): modes = [] - names = [] + # names = [] re_or_im = mode_name[0:2] ## Re(...) - print(re_or_im) + # print(re_or_im) assert re_or_im == "Re" or re_or_im == "Im" legends = h5f[field_name].attrs["Legend"] for i in range(len(legends)): m = legends[i] if m.find(re_or_im) != -1: - names.append(m) + # names.append(m) modes.append(i) assert len(modes) @@ -137,13 +137,14 @@ def read_h5_1mode(args): return (f, drf, dtf, t) -def read_h5_all_modes(args): +def read_h5_all_modes(args) -> dict: """ return t, f(mode)|t, df(mode)/dr|t, df(mode)/dt|t f[t_i,modes] = at t_i what are the modes """ + dat = {} field_name = g_name_map[args["field_name"]] field_name_key = f"{field_name}.dat" dfield_name_dr_key = f"Dr{field_name}.dat" @@ -187,13 +188,21 @@ def read_h5_all_modes(args): dtf[k, j] = h5f[f"{dfield_name_dt_key}"][i, modes[j]] k += 1 - return (f, drf, dtf, t) + dat["t"] = t + dat["f"] = f + dat["drf"] = drf + dat["dtf"] = dtf + dat["n_dumps"] = n_dumps + + return dat def plot_simple_v_t(dat, args): """ plot value vs time """ + + print("plot value vs time ...") fig, axes = plt.subplots(3, 1, sharex=True) # f @@ -233,14 +242,74 @@ def plot_simple_v_t(dat, args): plt.savefig(file_out, dpi=200) +def plot_simple_modes(dat, args): + """ + plot modes for different times + """ + + print("plot modes for different times ...") + + t = dat["t"] + f = dat["f"] + drf = dat["drf"] + dtf = dat["dtf"] + n_dumps = dat["n_dumps"] + x = np.arange(0, len(f[0,:]), dtype=int) + p = 0 + for i in range(0, len(t), args["time_dump"]): + title = f"t_{t[i]}" + + fig, axes = plt.subplots(3, 1, sharex=True) + + # f + ax = axes[0] + ax.set_title(title) + label = args["field_name"] + conf = dict(ls="-", label=label, color="k") + ax.plot(x, np.abs(f[p, :]), **conf) + ax.set_ylabel(label) + ax.set_yscale('log') + ax.grid(True) + + # drf + ax = axes[1] + label = "d" + args["field_name"] + "/dr" + conf = dict(ls="-", label=label, color="k") + ax.plot(x, np.abs(drf[p,:]), **conf) + ax.set_ylabel(label) + ax.set_yscale('log') + ax.grid(True) + + # dtf + ax = axes[2] + label = "d" + args["field_name"] + "/dt" + conf = dict(ls="-", label=label, color="k") + ax.plot(x, np.abs(dtf[p,:]), **conf) + ax.grid(True) + ax.set_ylabel(label) + ax.set_yscale('log') + ax.set_xlabel("modes") + + plt.tight_layout() + # plt.show() + + file_out = os.path.join( + args["dout"], + args["field_name"] + "_" + "modes_" + f"{title}.png", + ) + plt.savefig(file_out, dpi=200) + + p += 1 + + def debug_plot_simple(args): # value vs time - # dat = read_h5_1mode(args) - # plot_simple_v_t(dat, args) + dat = read_h5_1mode(args) + plot_simple_v_t(dat, args) # conv test dat = read_h5_all_modes(args) - # plot_simple_mode(dat, args) + plot_simple_modes(dat, args) def main(args): From 193da9b17b640887cb0778af909ff3274bcf1dc5 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 07:24:20 -0500 Subject: [PATCH 11/72] fix: minor. --- scripts/cce/debug_athk_to_spectre.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/cce/debug_athk_to_spectre.py b/scripts/cce/debug_athk_to_spectre.py index 5f3cc09f3..b56cf1800 100755 --- a/scripts/cce/debug_athk_to_spectre.py +++ b/scripts/cce/debug_athk_to_spectre.py @@ -298,7 +298,7 @@ def plot_simple_modes(dat, args): args["field_name"] + "_" + "modes_" + f"{title}.png", ) plt.savefig(file_out, dpi=200) - + plt.close() p += 1 From 47636c01edcb8c3358fd36d372966d4f85e5f8d3 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 09:57:50 -0500 Subject: [PATCH 12/72] fix: minor bug. --- scripts/cce/athk_to_spectre.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 57933554c..d27c0ba80 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -95,7 +95,7 @@ def parse_cli(): p.add_argument( "-t_deriv", type=str, - default="Fourier", + default="fin_diff", help="method to take the time derivative of fields:{Fourier,fin_diff}", ) p.add_argument( From 2d01df15868d8568133e1e6e3716d83fbd3f3006 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 09:58:13 -0500 Subject: [PATCH 13/72] update debug scripts. --- scripts/cce/cce_run.sh | 28 ++++++++++++++++++++++++---- scripts/cce/pre_run.sh | 37 ++++++++++++++++++++++++------------- 2 files changed, 48 insertions(+), 17 deletions(-) diff --git a/scripts/cce/cce_run.sh b/scripts/cce/cce_run.sh index 5037849f2..919bf90e7 100755 --- a/scripts/cce/cce_run.sh +++ b/scripts/cce/cce_run.sh @@ -1,13 +1,33 @@ #!/usr/bin/env bash +set -x +if false +then cd runs - echo "running preproc...\n" -./PreprocessCceWorldtube --input-file PreprocessCceWorldtube_test.yaml - +./PreprocessCceWorldtube --input-file PreprocessCceWorldtube_bin.yaml echo "-------------------------" +echo "running CCE..." +./CharacteristicExtract --input-file cce_par_bin.yaml +# h5 test +elif false +then +cd runs +echo "running preproc...\n" +./PreprocessCceWorldtube --input-file PreprocessCceWorldtube_h5.yaml +echo "-------------------------" +echo "running CCE..." +./CharacteristicExtract --input-file cce_par_h5.yaml +# q2 +elif true +then +cd runs +echo "running preproc...\n" +./PreprocessCceWorldtube --input-file PreprocessCceWorldtube_q2_h5.yaml +echo "-------------------------" echo "running CCE..." -./CharacteristicExtract --input-file cce_par.yaml +./CharacteristicExtract --input-file cce_par_q2_h5.yaml +fi diff --git a/scripts/cce/pre_run.sh b/scripts/cce/pre_run.sh index 098e7d536..33a6c0cd8 100755 --- a/scripts/cce/pre_run.sh +++ b/scripts/cce/pre_run.sh @@ -9,33 +9,44 @@ then fi file="./debug/CceR0030.00.h5" +dout='./debug/bin' + +if false +then +./athk_to_spectre.py -fpath dat/bbh_q2.0_chizp0.0_chizm0.0_d10.0_lev13_n128_fixed_cce_decomp_shell_3.h5 \ +-ftype h5 -d_out ./debug/q2/ -radius 100 +fi + +file="./debug/q2/CceR0100.00.h5" +dout='./debug/q2' + # plot fig fname="gxx" -mode="Re(0,0)" -./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ - -fpath $file -field_name $fname -field_mode $mode - mode="Re(2,2)" -./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ +./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ -fpath $file -field_name $fname -field_mode $mode +#mode="Re(0,0)" +#./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ +# -fpath $file -field_name $fname -field_mode $mode + fname="alp" -mode="Re(0,0)" -./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ - -fpath $file -field_name $fname -field_mode $mode +#mode="Re(0,0)" +#./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ +# -fpath $file -field_name $fname -field_mode $mode mode="Re(2,2)" -./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ +./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ -fpath $file -field_name $fname -field_mode $mode fname="betax" -mode="Re(0,0)" -./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ - -fpath $file -field_name $fname -field_mode $mode +#mode="Re(0,0)" +#./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ +# -fpath $file -field_name $fname -field_mode $mode mode="Re(2,2)" -./debug_athk_to_spectre.py -debug plot_simple -dout ./debug \ +./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ -fpath $file -field_name $fname -field_mode $mode From 7482abdde26bca28c62a6982706848c6a52c9cf1 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 09:58:52 -0500 Subject: [PATCH 14/72] feat: add cce spectre plot. --- scripts/cce/cce_plot.py | 97 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 97 insertions(+) create mode 100755 scripts/cce/cce_plot.py diff --git a/scripts/cce/cce_plot.py b/scripts/cce/cce_plot.py new file mode 100755 index 000000000..6006847a3 --- /dev/null +++ b/scripts/cce/cce_plot.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +## Alireza Rashti - Jan 2025 (C) +## usage: +## $ ./me -h +## + +import sys +import os +import numpy as np +import math as m +import argparse +import re +import h5py + +# import matplotlib.pyplot as plt +# import glob +# import sympy +## ---------------------------------------------------------------------- ## + + +def parse_cli(): + """ + arg parser + """ + p = argparse.ArgumentParser(description="plotting cce spectre output") + p.add_argument("-f", + type=str, + required=True, + help="/path/to/specter/cce/output.h5") + p.add_argument("-dout", + type=str, + required=True, + help="/path/to/output/file.txt") + p.add_argument( + "-field", + type=str, + default="Strain", + help="field names:[News,Psi0,Psi1,Psi2,Psi3,Psi4,Strain]", + ) + p.add_argument( + "-mode", + type=str, + default="Real Y_2,2", + help='field modes:["Real Y_2,-2","Imag Y_2,-2", ...]', + ) + args = p.parse_args() + return args + + +def find_h5_1mode(h5f, field_name, mode_name, args): + mode = 0 + flag = False + for m in h5f[field_name].attrs["Legend"]: + if m.find(mode_name) != -1: + print("found mode for", field_name, m, mode_name) + flag = True + break + mode += 1 + + assert flag == True + return mode + + +def read_save_to_txt(args): + + with h5py.File(args["f"], "r") as h5f: + # find group + keys = h5f.keys() + group = [k for k in keys if k.find("Spectre") != -1][0] + group = "/" + group + "/" + dataset = group+args["field"] + mode = find_h5_1mode(h5f, dataset, args["mode"], args) + t = h5f[dataset][:, 0] + f = h5f[dataset][:, mode] + + mode_txt = args["mode"].replace(",", "_") + mode_txt = mode_txt.replace(" ", "") + out = os.path.join(args["dout"], args["field"] + "_" + mode_txt + ".txt") + #print(out) + stack = np.column_stack((t, f)) + np.savetxt(out, + stack, + header=f'# t {args["mode"]}', + fmt="%.16e %.16e") + + +def main(args): + """ + read cce spectre file + """ + + read_save_to_txt(args) + + +if __name__ == "__main__": + args = parse_cli() + main(args.__dict__) From af2703473910399c790e53bc6451590fdf792a51 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 10:42:56 -0500 Subject: [PATCH 15/72] fix: bugs. --- scripts/cce/cce_plot.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/scripts/cce/cce_plot.py b/scripts/cce/cce_plot.py index 6006847a3..2afcd65cf 100755 --- a/scripts/cce/cce_plot.py +++ b/scripts/cce/cce_plot.py @@ -34,14 +34,14 @@ def parse_cli(): p.add_argument( "-field", type=str, - default="Strain", + default="Psi4", help="field names:[News,Psi0,Psi1,Psi2,Psi3,Psi4,Strain]", ) p.add_argument( "-mode", type=str, - default="Real Y_2,2", - help='field modes:["Real Y_2,-2","Imag Y_2,-2", ...]', + default="Y_2,2", + help='field modes:["Y_2,-2","Y_2,-2", ...]', ) args = p.parse_args() return args @@ -69,20 +69,21 @@ def read_save_to_txt(args): group = [k for k in keys if k.find("Spectre") != -1][0] group = "/" + group + "/" dataset = group+args["field"] - mode = find_h5_1mode(h5f, dataset, args["mode"], args) + mode_re = find_h5_1mode(h5f, dataset, "Real "+args["mode"], args) + mode_im = find_h5_1mode(h5f, dataset, "Imag "+args["mode"], args) t = h5f[dataset][:, 0] - f = h5f[dataset][:, mode] + re = h5f[dataset][:, mode_re] + im = h5f[dataset][:, mode_im] mode_txt = args["mode"].replace(",", "_") mode_txt = mode_txt.replace(" ", "") out = os.path.join(args["dout"], args["field"] + "_" + mode_txt + ".txt") #print(out) - stack = np.column_stack((t, f)) + stack = np.column_stack((t, re, im)) np.savetxt(out, stack, - header=f'# t {args["mode"]}', - fmt="%.16e %.16e") - + header=f'# t re{args["mode"]} im{args["mode"]}', + fmt="%.16e %.16e %.16e") def main(args): """ From 11c04dca5286a9a2c7155e110d27da1d845efb53 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 10:43:19 -0500 Subject: [PATCH 16/72] renamed: cce_plot.py -> cce_txt.py --- scripts/cce/{cce_plot.py => cce_txt.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename scripts/cce/{cce_plot.py => cce_txt.py} (100%) diff --git a/scripts/cce/cce_plot.py b/scripts/cce/cce_txt.py similarity index 100% rename from scripts/cce/cce_plot.py rename to scripts/cce/cce_txt.py From a97cc9a7d5bddbb52ac0b3d568c8a8e7881984ea Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 11:03:29 -0500 Subject: [PATCH 17/72] renamed: cce_txt.py -> cce_extract_txt.py --- scripts/cce/{cce_txt.py => cce_extract_txt.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename scripts/cce/{cce_txt.py => cce_extract_txt.py} (100%) diff --git a/scripts/cce/cce_txt.py b/scripts/cce/cce_extract_txt.py similarity index 100% rename from scripts/cce/cce_txt.py rename to scripts/cce/cce_extract_txt.py From 436c93448bdc8ddce603fc2d59c689cc7c663802 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 11:05:14 -0500 Subject: [PATCH 18/72] fix: minors. --- scripts/cce/cce_extract_txt.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/scripts/cce/cce_extract_txt.py b/scripts/cce/cce_extract_txt.py index 2afcd65cf..93c13fdc2 100755 --- a/scripts/cce/cce_extract_txt.py +++ b/scripts/cce/cce_extract_txt.py @@ -22,7 +22,8 @@ def parse_cli(): """ arg parser """ - p = argparse.ArgumentParser(description="plotting cce spectre output") + p = argparse.ArgumentParser( + description="extract cce spectre output in a txt file.") p.add_argument("-f", type=str, required=True, @@ -68,9 +69,9 @@ def read_save_to_txt(args): keys = h5f.keys() group = [k for k in keys if k.find("Spectre") != -1][0] group = "/" + group + "/" - dataset = group+args["field"] - mode_re = find_h5_1mode(h5f, dataset, "Real "+args["mode"], args) - mode_im = find_h5_1mode(h5f, dataset, "Imag "+args["mode"], args) + dataset = group + args["field"] + mode_re = find_h5_1mode(h5f, dataset, "Real " + args["mode"], args) + mode_im = find_h5_1mode(h5f, dataset, "Imag " + args["mode"], args) t = h5f[dataset][:, 0] re = h5f[dataset][:, mode_re] im = h5f[dataset][:, mode_im] @@ -78,12 +79,15 @@ def read_save_to_txt(args): mode_txt = args["mode"].replace(",", "_") mode_txt = mode_txt.replace(" ", "") out = os.path.join(args["dout"], args["field"] + "_" + mode_txt + ".txt") - #print(out) + # print(out) stack = np.column_stack((t, re, im)) - np.savetxt(out, - stack, - header=f'# t re{args["mode"]} im{args["mode"]}', - fmt="%.16e %.16e %.16e") + np.savetxt( + out, + stack, + header=f'# t re{args["mode"]} im{args["mode"]}', + fmt="%.16e %.16e %.16e", + ) + def main(args): """ From 2de7fa8b938443ad5091d7ad4be4feaca239f698 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 22 Jan 2025 11:09:35 -0500 Subject: [PATCH 19/72] add: memo. --- scripts/cce/readme.txt | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 scripts/cce/readme.txt diff --git a/scripts/cce/readme.txt b/scripts/cce/readme.txt new file mode 100644 index 000000000..b7f202df5 --- /dev/null +++ b/scripts/cce/readme.txt @@ -0,0 +1,7 @@ +# h5 +./athk_to_spectre.py -fpath dat/bbh_q2.0_chizp0.0_chizm0.0_d10.0_lev13_n128_fixed_cce_decomp_shell_3.h5 -ftype h5 -d_out ./debug/q2/ -radius 100 + +# extract txt +./cce_extract_txt.py -f ./runs/CharacteristicExtractReduction.h5 -dout ./debug/q2/ -field Strain +./cce_extract_txt.py -f ./runs/CharacteristicExtractReduction.h5 -dout ./debug/q2/ -field Psi4 + From 49e9dcaa59ec608aa8fbf4882fbb9ad2a34eee86 Mon Sep 17 00:00:00 2001 From: Alireza Date: Mon, 27 Jan 2025 10:44:39 -0500 Subject: [PATCH 20/72] fix & feat: fix some bugs and improve bin read speed. --- scripts/cce/athk_to_spectre.py | 145 +++++++++++++++++---------------- scripts/cce/pre_run.sh | 43 ++++++---- 2 files changed, 103 insertions(+), 85 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index d27c0ba80..70efaa27b 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -3,6 +3,9 @@ ## usage: ## $ ./me -h ## +## BUG: +## -[] see TODOs +## -[] see NOTEs import sys import os @@ -82,8 +85,14 @@ def parse_cli(): """ p = argparse.ArgumentParser( description="convert Athenak CCE dumps to Spectre CCE") - p.add_argument("-fpath", type=str, required=True, help="/path/to/cce/h5/dumps") - p.add_argument("-ftype", type=str, required=True, help="input file type:{h5,bin}") + p.add_argument("-fpath", + type=str, + required=True, + help="/path/to/cce/dir/or/file/dumps") + p.add_argument("-ftype", + type=str, + required=True, + help="input file type:{h5,bin}") p.add_argument("-d_out", type=str, required=True, help="/path/to/output/dir") p.add_argument("-debug", type=str, default="n", help="debug=[y,n]") p.add_argument( @@ -141,28 +150,44 @@ def load(fpath: str, field_name: str, attrs: dict) -> list: ret[g_re, i, :] = h5_re ret[g_im, i, :] = h5_im elif attrs["file_type"] == "bin": - lev_t = attrs["lev_t"] - max_n = attrs["max_n"] - max_lm = attrs["max_lm"] - shape = (len([g_re, g_im]), lev_t, max_n, max_lm) - ret = np.empty(shape=shape, dtype=float) # Load the list of files - flist = sorted(glob.glob(fpath+'/cce_*.bin')) + ##TODO: this depends on file name + flist = sorted(glob.glob(fpath + "/cce_*.bin")) data = [] dat_real = [] dat_imag = [] t = [] for f in flist: - nr, num_l_modes, time, rin, rout, data_real, data_imag, index_to_lm = read_cce_file(f) + ( + nr, + num_l_modes, + time, + rin, + rout, + data_real, + data_imag, + index_to_lm, + ) = read_cce_file(f) t.append(time) - dat_real.append(data_real[:,g_name_index[field_name],:]) - dat_imag.append(data_imag[:,g_name_index[field_name],:]) + dat_real.append(data_real[:, g_name_index[field_name], :]) + dat_imag.append(data_imag[:, g_name_index[field_name], :]) dat_real = np.array(dat_real) dat_imag = np.array(dat_imag) t = np.array(t) - f_real = interp1d(t, dat_real,axis=0) - f_imag = interp1d(t, dat_imag,axis=0) + f_real = interp1d(t, dat_real, axis=0) + f_imag = interp1d(t, dat_imag, axis=0) + attrs["lev_t"] = len(flist) + attrs["max_n"] = nr + attrs["max_l"] = num_l_modes + attrs["max_lm"] = (num_l_modes + 1)**2 + attrs["r_in"] = np.array([rin]) + attrs["r_out"] = np.array([rout]) + attrs["time"] = np.linspace(t.min(), t.max(), t.shape[0]) + + shape = (len([g_re, g_im]), attrs["lev_t"], attrs["max_n"], attrs["max_lm"]) + ret = np.empty(shape=shape, dtype=float) + ## NOTE: is it fine that we do interpolation? ret[g_re, :] = f_real(attrs["time"]) ret[g_im, :] = f_imag(attrs["time"]) else: @@ -180,7 +205,7 @@ def lm_mode(l, m): def read_cce_file(filename): - """ + """ Reads binary data from a CCE output file generated by the provided C++ function. Parameters: @@ -196,48 +221,49 @@ def read_cce_file(filename): data_imag (numpy.ndarray): Imaginary part of the data, shape (nr, 10, num_angular_modes). index_to_lm (dict): Mapping from data index to (l, m) values. """ - with open(filename, 'rb') as f: - # Read number of radial points (nr) - nr_bytes = f.read(4) - nr = struct.unpack(' Date: Mon, 27 Jan 2025 11:48:46 -0500 Subject: [PATCH 21/72] minor. --- scripts/cce/pre_run.sh | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/scripts/cce/pre_run.sh b/scripts/cce/pre_run.sh index 72e4d8b4f..5daa3406c 100755 --- a/scripts/cce/pre_run.sh +++ b/scripts/cce/pre_run.sh @@ -1,7 +1,7 @@ #!/usr/bin/env bash # bin -if true +if false then # create dout='./debug/dev' @@ -57,5 +57,23 @@ mode="Re(2,2)" #mode="Re(2,2)" #./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ # -fpath $file -field_name $fname -field_mode $mode + +# h5 bbh, q=2, fourier time derivatives +elif true +then +file="./debug/q2_t_fourier/CceR0100.00.h5" +dout='./debug/q2_t_fourier' +mkdir -p ${dout} + +./athk_to_spectre.py -fpath dat/bbh_q2.0_chizp0.0_chizm0.0_d10.0_lev13_n128_fixed_cce_decomp_shell_3.h5 \ +-ftype h5 -d_out "$dout" -radius 100 -t_deriv "Fourier" + + +# plot fig: +fname="gxx" +mode="Re(2,2)" +./debug_athk_to_spectre.py -debug plot_simple -dout "${dout}" \ + -fpath $file -field_name $fname -field_mode $mode + fi From a560a82076daf833a98410f2eb779facad63e051 Mon Sep 17 00:00:00 2001 From: Alireza Date: Mon, 27 Jan 2025 11:55:12 -0500 Subject: [PATCH 22/72] fix: fourier time derives. --- scripts/cce/athk_to_spectre.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 70efaa27b..30ecb6345 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -185,7 +185,12 @@ def load(fpath: str, field_name: str, attrs: dict) -> list: attrs["r_out"] = np.array([rout]) attrs["time"] = np.linspace(t.min(), t.max(), t.shape[0]) - shape = (len([g_re, g_im]), attrs["lev_t"], attrs["max_n"], attrs["max_lm"]) + shape = ( + len([g_re, g_im]), + attrs["lev_t"], + attrs["max_n"], + attrs["max_lm"], + ) ret = np.empty(shape=shape, dtype=float) ## NOTE: is it fine that we do interpolation? ret[g_re, :] = f_real(attrs["time"]) @@ -279,6 +284,9 @@ def get_attribute(fpath: str, with h5py.File(fpath, "r") as h5f: # find attribute about num. of time level, and n,l,m in C_nlm attrs["lev_t"] = len(h5f.keys()) - 1 + if attrs["lev_t"] % 2 == 0: # for fourier transformation + attrs["lev_t"] -= 1 + attrs["max_n"], attrs["max_lm"] = h5f[f"1/{field_name}/re"].shape attrs["max_l"] = int(math.sqrt(attrs["max_lm"])) - 1 attrs["r_in"] = h5f["metadata"].attrs["Rin"] @@ -349,10 +357,11 @@ def time_derivative_fourier(field: np.array, field_name: str, attrs: dict, dfft_coeff = np.empty_like(fft_coeff) dfft_coeff[0] = 0 + dfft_coeff[1:half] = (-np.imag(fft_coeff[1:half]) + 1j * np.real(fft_coeff[1:half])) * omega[1:] dfft_coeff[half:] = (np.imag(fft_coeff[half:]) - - 1j * np.real(fft_coeff[half:])) * omega[::-1][1:] + 1j * np.real(fft_coeff[half:])) * (omega[::-1][1:]) # not optimized version """ From ab602d7b0cbcf28c700f77029a6d0654080b98a3 Mon Sep 17 00:00:00 2001 From: Alireza Date: Mon, 27 Jan 2025 12:03:38 -0500 Subject: [PATCH 23/72] add note. --- scripts/cce/athk_to_spectre.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 30ecb6345..783c2ef1d 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -334,7 +334,8 @@ def time_derivative_fourier(field: np.array, field_name: str, attrs: dict, return the time derivative of the given field using Fourier method field(t,rel/img,n,lm) """ - + + ## TODO: Fourier time derives not tested! print(f"Fourier time derivative: {field_name}", flush=True) _, len_t, len_n, len_lm = field.shape dt = attrs["time"][2] - attrs["time"][1] From 1ff1ea02defb116fc1ddb5fe477513a7100b23d5 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 29 Jan 2025 10:12:28 -0500 Subject: [PATCH 24/72] feat: add chebu expansion. --- scripts/cce/athk_to_spectre.py | 113 ++++++++++++++++++++++++++++++--- 1 file changed, 105 insertions(+), 8 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 783c2ef1d..25fca9a5d 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -89,10 +89,12 @@ def parse_cli(): type=str, required=True, help="/path/to/cce/dir/or/file/dumps") - p.add_argument("-ftype", - type=str, - required=True, - help="input file type:{h5,bin}") + p.add_argument( + "-ftype", + type=str, + required=True, + help="input file type:{h5:for pitnull,bin: for athenak}", + ) p.add_argument("-d_out", type=str, required=True, help="/path/to/output/dir") p.add_argument("-debug", type=str, default="n", help="debug=[y,n]") p.add_argument( @@ -158,6 +160,7 @@ def load(fpath: str, field_name: str, attrs: dict) -> list: dat_imag = [] t = [] for f in flist: + ##TODO: it reads all field and only use one of them. it's inefficient ( nr, num_l_modes, @@ -311,7 +314,7 @@ def time_derivative_findiff(field: np.array, field_name: str, attrs: dict, args) -> np.array: """ return the time derivative of the given field using finite diff. 2nd order - field(t,rel/img,n,lm) + field(rel/img,t,n,lm) """ print(f"finite difference time derivative: {field_name}", flush=True) @@ -332,9 +335,9 @@ def time_derivative_fourier(field: np.array, field_name: str, attrs: dict, args) -> np.array: """ return the time derivative of the given field using Fourier method - field(t,rel/img,n,lm) + field(rel/img,t, n,lm) """ - + ## TODO: Fourier time derives not tested! print(f"Fourier time derivative: {field_name}", flush=True) _, len_t, len_n, len_lm = field.shape @@ -425,7 +428,7 @@ def radial_derivative_at_r_chebu(field: np.array, collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} - field(t,rel/img,n,lm) + field(rel/img,t,n,lm) """ print(f"ChebyU radial derivative: {field_name}", flush=True) @@ -489,6 +492,97 @@ def time_derivative(field: np.array, field_name: str, attrs: dict, args): raise ValueError("no such option") +class ChebUExpansion: + """ + f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind + collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) + """ + + def __init__(self, attrs: dict, args: dict): + self.N = N = attrs["max_n"] # num. of coeffs or num of collocation pnts + + # roots: + x_i = np.empty(shape=(N), dtype=float) + for i in range(N): + x_i[i] = cos(pi * (i + 1) / (N + 1)) + self.x_i = x_i + + # weights: + w_i = np.empty(shape=(N), dtype=float) + for i in range(N): + w_i[i] = sin(pi * (i + 1) / (N + 1)) + w_i = np.square(w_i) + w_i *= pi / (N + 1) + self.w_i = w_i + + # ChebU_j(x_i): + Uj_xi = np.empty(shape=(N, N), dtype=float) + for j in range(N): + for i in range(N): + x_i = cos(pi * (i + 1) / (N + 1)) + Uj_xi[j, i] = special.chebyu(j)(x_i) + self.Uj_xi = Uj_xi + + def Coefficients(self, field: np.array) -> np.array: + """ + use Gauss Quadrature to expand the given field in ChebU + f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind + collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) + x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} + + """ + + coeffs = np.zeros(shape=(N), dtype=float) + + w_i = self.w_i + x_i = self.w_i + Uji = self.Uj_xi + + for j in range(self.N): + tmp = 0.0 + for i in range(self.N): + tmp += w_i * field[i] * Uji[j, i] + coeffs[j] = tmp + + return coeffs + + +def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, args): + """ + expands field in the radial direction using chebyshev of second kind + and returns the expansion coefficients. we have + + f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind + collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) + x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} + + field(rel/img,t,n,lm) + + """ + + print(f"ChebU expansion {field_name}", flush=True) + expand = ChebUExpansion(attrs, args) + _, len_t, _, len_lm = field.shape + + for t in range(len_t): + for lm in range(len_lm): + field[g_re, t, :, lm] = expand.Coefficients(field[g_re, t, :, lm]) + field[g_im, t, :, lm] = expand.Coefficients(field[g_im, t, :, lm]) + + +def radial_expansion(field: np.array, field_name: str, attrs: dict, args): + """ + expands field in the radial direction(if needed) and returns + """ + + if args["ftype"] == "bin": + return radial_expansion_chebu(field, field_name, attrs, args) + elif args["ftype"] == "h5": + return field + else: + raise ValueError("no such option") + + def radial_derivative_at_r(field: np.array, field_name: str, attrs: dict, args): """ return the radial derivative of the given field at R=r @@ -557,6 +651,9 @@ def process_field(field_name: str) -> dict: field = load(args["fpath"], field_name, attrs) # db[f"{field_name}"] = field + # radial expansion + field = radial_expansion(field, field_name, attrs, args) + # time derivative dfield_dt = time_derivative(field, field_name, attrs, args) From c30296344082db86e968cf829a12bd95980f3651 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 29 Jan 2025 10:12:45 -0500 Subject: [PATCH 25/72] minor. --- scripts/cce/pre_run.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/cce/pre_run.sh b/scripts/cce/pre_run.sh index 5daa3406c..1941f639b 100755 --- a/scripts/cce/pre_run.sh +++ b/scripts/cce/pre_run.sh @@ -1,7 +1,7 @@ #!/usr/bin/env bash # bin -if false +if true then # create dout='./debug/dev' @@ -59,7 +59,7 @@ mode="Re(2,2)" # -fpath $file -field_name $fname -field_mode $mode # h5 bbh, q=2, fourier time derivatives -elif true +elif false then file="./debug/q2_t_fourier/CceR0100.00.h5" dout='./debug/q2_t_fourier' From 997c08ce624fdb73dbe59577823ef52b67585cef Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 29 Jan 2025 10:18:17 -0500 Subject: [PATCH 26/72] fix: bugs. --- scripts/cce/athk_to_spectre.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 25fca9a5d..6798690ca 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -500,17 +500,18 @@ class ChebUExpansion: def __init__(self, attrs: dict, args: dict): self.N = N = attrs["max_n"] # num. of coeffs or num of collocation pnts + pi = math.pi # roots: x_i = np.empty(shape=(N), dtype=float) for i in range(N): - x_i[i] = cos(pi * (i + 1) / (N + 1)) + x_i[i] = math.cos(pi * (i + 1) / (N + 1)) self.x_i = x_i # weights: w_i = np.empty(shape=(N), dtype=float) for i in range(N): - w_i[i] = sin(pi * (i + 1) / (N + 1)) + w_i[i] = math.sin(pi * (i + 1) / (N + 1)) w_i = np.square(w_i) w_i *= pi / (N + 1) self.w_i = w_i @@ -519,7 +520,7 @@ def __init__(self, attrs: dict, args: dict): Uj_xi = np.empty(shape=(N, N), dtype=float) for j in range(N): for i in range(N): - x_i = cos(pi * (i + 1) / (N + 1)) + x_i = math.cos(pi * (i + 1) / (N + 1)) Uj_xi[j, i] = special.chebyu(j)(x_i) self.Uj_xi = Uj_xi @@ -532,7 +533,7 @@ def Coefficients(self, field: np.array) -> np.array: """ - coeffs = np.zeros(shape=(N), dtype=float) + coeffs = np.zeros(shape=(self.N), dtype=float) w_i = self.w_i x_i = self.w_i @@ -541,13 +542,14 @@ def Coefficients(self, field: np.array) -> np.array: for j in range(self.N): tmp = 0.0 for i in range(self.N): - tmp += w_i * field[i] * Uji[j, i] + tmp += w_i[i] * field[i] * Uji[j, i] coeffs[j] = tmp - return coeffs + return coeffs * (2.0 / math.pi) -def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, args): +def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, + args) -> np.array: """ expands field in the radial direction using chebyshev of second kind and returns the expansion coefficients. we have @@ -569,6 +571,8 @@ def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, args): field[g_re, t, :, lm] = expand.Coefficients(field[g_re, t, :, lm]) field[g_im, t, :, lm] = expand.Coefficients(field[g_im, t, :, lm]) + return field + def radial_expansion(field: np.array, field_name: str, attrs: dict, args): """ From 4f88e0563462bef5e080b5ec717f8731fdcc9701 Mon Sep 17 00:00:00 2001 From: Alireza Date: Wed, 29 Jan 2025 11:03:07 -0500 Subject: [PATCH 27/72] chores. --- scripts/cce/athk_to_spectre.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 6798690ca..16ba56386 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -431,11 +431,17 @@ def radial_derivative_at_r_chebu(field: np.array, field(rel/img,t,n,lm) """ - print(f"ChebyU radial derivative: {field_name}", flush=True) - _, len_t, len_n, len_lm = field.shape - r_1 = attrs["r_in"][0] r_2 = attrs["r_out"][0] + r = args["radius"] + + print( + f"ChebyU radial derivative: {field_name}, at r={r} in [{r_1},{r_2}]", + flush=True, + ) + + _, len_t, len_n, len_lm = field.shape + assert r_1 != r_2 dx_dr = 2 / (r_2 - r_1) From 60fd8d3e5f3dc4167259c79e0daddd5cd12cb3a9 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 08:46:04 -0500 Subject: [PATCH 28/72] fix: bugs add: tests. --- scripts/cce/athk_to_spectre.py | 45 +++++++++++++++++++++++++--------- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 16ba56386..ea62b563c 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -526,10 +526,12 @@ def __init__(self, attrs: dict, args: dict): Uj_xi = np.empty(shape=(N, N), dtype=float) for j in range(N): for i in range(N): - x_i = math.cos(pi * (i + 1) / (N + 1)) - Uj_xi[j, i] = special.chebyu(j)(x_i) + Uj_xi[j, i] = special.chebyu(j)(x_i[i]) self.Uj_xi = Uj_xi + if True: + self.__Debug() + def Coefficients(self, field: np.array) -> np.array: """ use Gauss Quadrature to expand the given field in ChebU @@ -539,19 +541,40 @@ def Coefficients(self, field: np.array) -> np.array: """ - coeffs = np.zeros(shape=(self.N), dtype=float) - w_i = self.w_i - x_i = self.w_i + x_i = self.x_i Uji = self.Uj_xi - for j in range(self.N): - tmp = 0.0 - for i in range(self.N): - tmp += w_i[i] * field[i] * Uji[j, i] - coeffs[j] = tmp + coeffs = Uji @ (w_i * field) + coeffs *= 2.0 / math.pi + + return coeffs + + def __Debug(self): + """ + test if the expansion works for different known functions + """ + + N = self.N + x_i = self.x_i + + # populate funcs: + fs = [] + for n in range(N): + f = special.chebyu(n)(x_i) + fs.append((f, n)) + + # find coeffs + cs = [] + for n in range(N): + c = self.Coefficients(fs[n][0]) + cs.append((c, fs[n][1])) + + for n in range(N): + order = cs[n][1] + print(f"chebyu_{order}, coeffs = {cs[n][0]}\n", flush=True) - return coeffs * (2.0 / math.pi) + exit(1) def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, From 1b09aac405c9af2124add7dd47a2791ecf5ded49 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 10:16:04 -0500 Subject: [PATCH 29/72] fix: nasty bugs in cheb. exp. --- scripts/cce/athk_to_spectre.py | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index ea62b563c..f39534484 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -508,28 +508,19 @@ def __init__(self, attrs: dict, args: dict): self.N = N = attrs["max_n"] # num. of coeffs or num of collocation pnts pi = math.pi - # roots: - x_i = np.empty(shape=(N), dtype=float) - for i in range(N): - x_i[i] = math.cos(pi * (i + 1) / (N + 1)) - self.x_i = x_i - - # weights: - w_i = np.empty(shape=(N), dtype=float) - for i in range(N): - w_i[i] = math.sin(pi * (i + 1) / (N + 1)) - w_i = np.square(w_i) - w_i *= pi / (N + 1) - self.w_i = w_i + # U_N roots: + a_i = np.arange(0, N, dtype=int) + 1 + self.x_i = x_i = np.cos(math.pi * a_i / (N + 1)) + # quadrature weights + self.w_i = (1 - np.square(x_i)) * pi / (N + 1) # ChebU_j(x_i): Uj_xi = np.empty(shape=(N, N), dtype=float) for j in range(N): - for i in range(N): - Uj_xi[j, i] = special.chebyu(j)(x_i[i]) + Uj_xi[j, :] = special.chebyu(j)(x_i) self.Uj_xi = Uj_xi - if True: + if args["debug"] == "y": self.__Debug() def Coefficients(self, field: np.array) -> np.array: @@ -558,7 +549,7 @@ def __Debug(self): N = self.N x_i = self.x_i - # populate funcs: + # populate funcs using bases themselves fs = [] for n in range(N): f = special.chebyu(n)(x_i) @@ -570,11 +561,10 @@ def __Debug(self): c = self.Coefficients(fs[n][0]) cs.append((c, fs[n][1])) + # expect to see only the n-th entry is 1 for the chebyu(n) of order n for n in range(N): order = cs[n][1] - print(f"chebyu_{order}, coeffs = {cs[n][0]}\n", flush=True) - - exit(1) + print(f"chebyu{order}, coeffs = {cs[n][0]}\n", flush=True) def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, From 72a0e3d3adeac9fe014a5d0b155ccc0279ad6ba4 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 10:18:21 -0500 Subject: [PATCH 30/72] minor. --- scripts/cce/athk_to_spectre.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index f39534484..59ae79399 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -546,7 +546,7 @@ def __Debug(self): test if the expansion works for different known functions """ - N = self.N + N = self.N+1 # +1 to see the roots of U_N x_i = self.x_i # populate funcs using bases themselves From ea74e4af64fc01e276321f0b1c7840c0af738ff6 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 10:37:24 -0500 Subject: [PATCH 31/72] fix: collocation points for interpolation. --- src/utils/chebyshev.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/chebyshev.hpp b/src/utils/chebyshev.hpp index 4b7426879..ecc4c1f27 100644 --- a/src/utils/chebyshev.hpp +++ b/src/utils/chebyshev.hpp @@ -15,7 +15,7 @@ KOKKOS_INLINE_FUNCTION Real ChebyshevSecondKindCollocationPoints(Real xmin, Real xmax, int N, int k) { // over the interval of -1 to 1 - Real x_k = std::cos( M_PI * (k + 1) / (N + 2)); + Real x_k = std::cos( M_PI * (k + 1) / (N + 1)); return 0.5 * ( (xmin - xmax) * x_k + (xmin + xmax)); } From 30385b688b8cc9dbcea02287ac99e796953f2645 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 11:59:56 -0500 Subject: [PATCH 32/72] fix: sign bug due to x=[-1,1] convention. --- scripts/cce/athk_to_spectre.py | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 59ae79399..9fab2f395 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -75,6 +75,9 @@ ## various attrs g_attrs = None +## sign convention +g_sign = None + ## debug g_debug_max_l = 2 @@ -443,7 +446,7 @@ def radial_derivative_at_r_chebu(field: np.array, _, len_t, len_n, len_lm = field.shape assert r_1 != r_2 - dx_dr = 2 / (r_2 - r_1) + dx_dr = g_sign * 2 / (r_2 - r_1) if args["debug"] == "y": # populate collocation points, roots of U_i @@ -478,7 +481,7 @@ def radial_derivative_at_r_chebu(field: np.array, dfield = np.zeros(shape=(len([g_re, g_im]), len_t, len_lm)) r = args["radius"] - x = (2 * r - r_1 - r_2) / (r_2 - r_1) + x = g_sign * (2 * r - r_1 - r_2) / (r_2 - r_1) for k in range(len_n): dfield[:, :, :] += field[:, :, k, :] * dUk_dx(k, x) @@ -536,7 +539,7 @@ def Coefficients(self, field: np.array) -> np.array: x_i = self.x_i Uji = self.Uj_xi - coeffs = Uji @ (w_i * field) + coeffs = Uji @ (w_i * field[::-g_sign]) coeffs *= 2.0 / math.pi return coeffs @@ -546,7 +549,7 @@ def __Debug(self): test if the expansion works for different known functions """ - N = self.N+1 # +1 to see the roots of U_N + N = self.N + 1 # +1 to see the roots of U_N x_i = self.x_i # populate funcs using bases themselves @@ -631,7 +634,7 @@ def __init__(self, attrs: dict, args: dict): r_1 = attrs["r_in"][0] r_2 = attrs["r_out"][0] self.r = r = args["radius"] - self.x = (2 * r - r_1 - r_2) / (r_2 - r_1) + self.x = g_sign * (2 * r - r_1 - r_2) / (r_2 - r_1) if args["interpolation"] == "ChebU": self.Uk = np.empty(shape=self.len_n) @@ -801,8 +804,23 @@ def main(args): global g_attrs global g_args + global g_sign g_args = args + """ + note: + in Pitnull: + r = 0.5[(r_2 - r_1)*x + (r_2 + r_1)], r_2 > r_1 & x = [-1,1] + => x = -1 -> r = r_1 + => x = +1 -> r = r_2 + + in AthenaK: + r = 0.5[(r_1 - r_2)*x + (r_2 + r_1)], r_2 > r_1 & x = [1,-1] + => x = -1 -> r = r_2 + => x = +1 -> r = r_1 + + """ + g_sign = int(-1) if args["ftype"] == "bin" else int(1) # check if output dir exist, if not, mkdir if not os.path.exists(g_args["d_out"]): From 3afd56e90ac1256fd5cccf1607d04562f3fe59f2 Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 12:09:36 -0500 Subject: [PATCH 33/72] minors. --- scripts/cce/athk_to_spectre.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 9fab2f395..3d3baa84c 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -429,7 +429,7 @@ def radial_derivative_at_r_chebu(field: np.array, f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) - x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} + x = g_sign*(2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} field(rel/img,t,n,lm) """ @@ -524,14 +524,14 @@ def __init__(self, attrs: dict, args: dict): self.Uj_xi = Uj_xi if args["debug"] == "y": - self.__Debug() + self.__debug() - def Coefficients(self, field: np.array) -> np.array: + def coefficients(self, field: np.array) -> np.array: """ use Gauss Quadrature to expand the given field in ChebU f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) - x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} + x = g_sign*(2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} """ @@ -544,7 +544,7 @@ def Coefficients(self, field: np.array) -> np.array: return coeffs - def __Debug(self): + def __debug(self): """ test if the expansion works for different known functions """ @@ -561,7 +561,7 @@ def __Debug(self): # find coeffs cs = [] for n in range(N): - c = self.Coefficients(fs[n][0]) + c = self.coefficients(fs[n][0]) cs.append((c, fs[n][1])) # expect to see only the n-th entry is 1 for the chebyu(n) of order n @@ -578,7 +578,7 @@ def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, f(x) = sum_{i=0}^{N-1} C_i U_i(x), U_i(x) Chebyshev of 2nd kind collocation points (roots of U_i): x_i = cos(pi*(i+1)/(N+1)) - x = (2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} + x = g_sign*(2*r - r_1 - r_2)/(r_2 - r_1), notes: x != {1 or -1} field(rel/img,t,n,lm) @@ -590,8 +590,8 @@ def radial_expansion_chebu(field: np.array, field_name: str, attrs: dict, for t in range(len_t): for lm in range(len_lm): - field[g_re, t, :, lm] = expand.Coefficients(field[g_re, t, :, lm]) - field[g_im, t, :, lm] = expand.Coefficients(field[g_im, t, :, lm]) + field[g_re, t, :, lm] = expand.coefficients(field[g_re, t, :, lm]) + field[g_im, t, :, lm] = expand.coefficients(field[g_im, t, :, lm]) return field From 70343794e2b7946d5379128743aacb7c95e6e93f Mon Sep 17 00:00:00 2001 From: Alireza Date: Thu, 30 Jan 2025 12:56:40 -0500 Subject: [PATCH 34/72] minor. --- scripts/cce/athk_to_spectre.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/cce/athk_to_spectre.py b/scripts/cce/athk_to_spectre.py index 3d3baa84c..fe409114a 100755 --- a/scripts/cce/athk_to_spectre.py +++ b/scripts/cce/athk_to_spectre.py @@ -636,10 +636,13 @@ def __init__(self, attrs: dict, args: dict): self.r = r = args["radius"] self.x = g_sign * (2 * r - r_1 - r_2) / (r_2 - r_1) + assert -1 < self.x < 1 + if args["interpolation"] == "ChebU": self.Uk = np.empty(shape=self.len_n) for k in range(self.len_n): self.Uk[k] = special.chebyu(k)(self.x) + # print("x,Uk",self.x,self.Uk) self.interp = self.interpolate_at_r_chebu else: raise ValueError("no such option") From 63c47d2ba0ce20a53237bde45bc97dbb92b954f2 Mon Sep 17 00:00:00 2001 From: Hengrui Zhu Date: Sat, 1 Feb 2025 16:00:53 -0500 Subject: [PATCH 35/72] update legendre roots calculation and include unit test for spherical harmonics and integration --- .../unit_tests/gauss_legendre_test.athinput | 48 +++++ src/pgen/unit_tests/gauss_legendre_test.cpp | 118 +++++++++++ src/utils/legendre_roots.hpp | 196 ++++++++++-------- 3 files changed, 280 insertions(+), 82 deletions(-) create mode 100644 inputs/unit_tests/gauss_legendre_test.athinput create mode 100644 src/pgen/unit_tests/gauss_legendre_test.cpp diff --git a/inputs/unit_tests/gauss_legendre_test.athinput b/inputs/unit_tests/gauss_legendre_test.athinput new file mode 100644 index 000000000..9f279a001 --- /dev/null +++ b/inputs/unit_tests/gauss_legendre_test.athinput @@ -0,0 +1,48 @@ +# Athena++ (Kokkos version) input file for spherical blast problem + + +problem = gauss legendre test + + +basename = gauss_legendre_test # problem ID: basename of output filenames + + +nghost = 2 # Number of ghost cells +nx1 = 20 # Number of zones in X1-direction +x1min = -20 # minimum value of X1 +x1max = 20 # maximum value of X1 +ix1_bc = outflow # inner-X1 boundary flag +ox1_bc = outflow # outer-X1 boundary flag + +nx2 = 20 # Number of zones in X2-direction +x2min = -20 # minimum value of X2 +x2max = 20 # maximum value of X2 +ix2_bc = outflow # inner-X2 boundary flag +ox2_bc = outflow # outer-X2 boundary flag + +nx3 = 20 # Number of zones in X3-direction +x3min = -20 # minimum value of X3 +x3max = 20 # maximum value of X3 +ix3_bc = outflow # inner-X3 boundary flag +ox3_bc = outflow # outer-X3 boundary flag + + +nx1 = 20 # Number of cells in each MeshBlock, X1-dir +nx2 = 20 # Number of cells in each MeshBlock, X2-dir +nx3 = 20 # Number of cells in each MeshBlock, X3-dir + +