"""Summarize physical time-domain modes, excluding the solver's tail-fit columns. The evolved residual modes already remove A/B. Additional F2/F4/F6 subtraction improves the angular sum without changing the fixed second-order puncture. Reference comparisons always use the same angular cutoff. Resolution, startup relaxation and angular truncation remain separate empirical error diagnostics. """ import argparse import csv import json from pathlib import Path import re import numpy as np def decimal(value): """Read Mathematica's exported arbitrary-precision decimal notation.""" return float(re.sub(r'`[0-9.]*', '', str(value)).replace('*^', 'e')) def read_parameters(path): params = {} for line in path.read_text().splitlines(): if '=' not in line: continue key, value = (part.strip() for part in line.split('=', 1)) if value in ('.true.', '.false.'): value = value == '.true.' elif value.startswith(("'", '"')): value = value[1:-1] else: try: value = int(value) except ValueError: value = float(value.replace('D', 'E').replace('d', 'e')) params[key] = value return params def supported_radius(params): if params['p_orb'] not in (6., 10.): raise ValueError('This benchmark has reference data only at r0/M=6 or 10') return int(params['p_orb']) def regularization(path): result = {} for row in json.loads(path.read_text()): result[int(row['r_over_M'])] = { mode['l']: {order: decimal(value) for order, value in zip(row['orders'], mode['terms'])} for mode in row['modes']} return result def read_reference(path, radius): data = json.loads(path.read_text()) if 'radii' in data: data = data['radii'][str(radius)] return {int(row.get('ell', row.get('l'))): decimal(row['AB']) for row in data['multipoles']} def check_probes(path, coefficients): rows = list(csv.reader(path.open())) b_errors, jump_errors, op_errors, op_relative = [], [], [], [] regularity = {} for row in rows: kind, radius, ell = row[:3] radius, ell = int(radius), int(ell) if kind == 'AB': b_errors.append(abs(float(row[3])-coefficients[radius][ell][0])) jump_errors.append(abs(float(row[4])-2*coefficients[radius][ell][-1])) elif kind == 'operator': error, scale = float(row[4]), float(row[5]) op_errors.append(error) op_relative.append(error/scale) elif kind == 'source': h, left, right = map(float, row[3:6]) regularity.setdefault((radius, ell), []).append((h, left, right)) limits = [] for (radius, ell), samples in sorted(regularity.items()): s0 = next(left for h, left, _ in samples if h == 0.) small = sorted((s for s in samples if s[0] > 0.))[:2] estimates = [] for h, left, right in small: estimates.append(dict(h_over_M=h, left_source=left, right_source=right, max_continuity_difference=max(abs(left-s0), abs(right-s0)), left_derivative=(s0-left)/h, right_derivative=(right-s0)/h)) limits.append(dict(radius=radius, ell=ell, m=ell, source_at_particle=s0, smallest_step_estimates=estimates)) return dict(normalization='G=c=1; M=q=1; covariant radial derivative', max_B_absolute_error=max(b_errors), max_A_jump_absolute_error=max(jump_errors), max_wave_operator_absolute_error=max(op_errors), max_wave_operator_relative_error=max(op_relative), wave_operator_derivative='Fourth-order centered finite difference of the analytic first derivative, h/M=1e-4', regularity_scope='Individual harmonic effective sources, not the unsummed three-dimensional source', regularity=limits) def convergence_summary(results): """Small comparisons reconstructed from the archived raw mode values.""" by_name={row['name']:row for row in results} summary={} spatial=[] for degree in (8,10,14,18): name=f'diag_l16_p{degree}_small_dt_r10' if name in by_name: row=by_name[name] item=dict(degree=degree,dt_over_M=row['dt_over_M'], AB=row['final_AB'][0],last_50M_span=row['last_50M_AB_span'][0]) if 'frequency_domain_comparison' in row: item['absolute_reference_error']=abs(row['frequency_domain_comparison']['AB_error'][0]) spatial.append(item) if spatial: summary['fixed_small_step_spatial_refinement']=spatial name='diag_l16_p10_small_dt_r10' if name in by_name: small_step=by_name[name]['final_AB'][0] temporal=[] for name in ('pilot10','dt025_r10','diag_l16_dt0125_r10'): if name in by_name: row=by_name[name]; index=row['l_values'].index(16) temporal.append(dict(name=name,dt_over_M=row['dt_over_M'], difference_from_small_step=row['final_AB'][index]-small_step, last_50M_span=row['last_50M_AB_span'][index])) summary['fixed_degree10_time_refinement']=temporal differences=[] for left,right,purpose in ( ('wt2_p14_r10','p14_r10','tube width, coarser time step'), ('fine_wt2_r10','fine_r10','tube width, finer time step'), ('fine_boundary_r10','fine_r10','hyperboloidal transition')): if left in by_name and right in by_name: a,b=by_name[left],by_name[right] delta=np.array(a['final_AB'])-np.array(b['final_AB']) differences.append(dict(first=left,second=right,purpose=purpose, force_difference=float(np.sum(delta)), maximum_mode_difference=float(np.max(abs(delta))))) summary['full_force_setting_variations']=differences return summary def analyze(directory, coefficients, reference=None): params = read_parameters(directory/'input.par') if params['mass'] != 1. or params['q_charge'] != 1. or params['ecc'] != 0.: raise ValueError('This analysis requires circular orbits with M=q=1') timing_file = directory/'timing.json' log = (directory/'run.log').read_text() if timing_file.exists(): timing = json.loads(timing_file.read_text()) if timing.get('returncode', 0) != 0: raise RuntimeError('Run did not exit successfully: ' + directory.name) else: # The initial pilot was launched under POSIX time before the common # runner existed. Its completed log records elapsed time directly. elapsed = re.search(r'^real\s+(\d+(?:\.\d+)?)\s*$', log, re.M) if not elapsed: raise RuntimeError('Run has no completion record: ' + directory.name) timing = dict(seconds=float(elapsed.group(1)), threads=None) data = np.loadtxt(directory/'Schwarzschild_wave.frl.asc', ndmin=2) ell = np.arange(params['lmin'], params['lmax']+1) expected = 1+len(ell)+(3 if params['fit_high_l'] else 0) if data.shape[1] != expected or not np.all(np.isfinite(data)): raise RuntimeError('Malformed radial-force output: ' + directory.name) modes = data[:, 1:1+len(ell)] radius = supported_radius(params) correction = np.array([sum(coefficients[radius][int(l)][k] for k in (2, 4, 6)) for l in ell]) reduced = modes-correction drho = float(re.search(r'delta_rho_min\s*=\s*(\S+)', log).group(1)) dt = params['coufac']*drho if params['t_final']-data[-1, 0] > (params['out0d_every']+2)*dt: raise RuntimeError('Incomplete time evolution: ' + directory.name) late = data[:, 0] >= max(0., data[-1, 0]-50.*params['mass']) final = modes[-1] result = dict(name=directory.name, parameters=params, dt_over_M=dt/params['mass'], seconds=timing['seconds'], threads=timing['threads'], last_time_over_M=float(data[-1, 0]/params['mass']), l_values=ell.tolist(), final_AB=final.tolist(), final_AB246=reduced[-1].tolist(), last_50M_AB_span=np.ptp(modes[late], axis=0).tolist(), last_50M_AB246_sum_span=float(np.ptp(np.sum(reduced[late],axis=1))), last_50M_AB246_sum_mean=float(np.mean(np.sum(reduced[late],axis=1))), selected_samples=[]) sample_indices = {len(data)-1} for target in (50., 100., 150., 200., 300., 400., 500., 600., 800., 1000.): if target*params['mass'] <= data[-1, 0]: sample_indices.add(int(np.argmin(abs(data[:, 0]-target*params['mass'])))) for index in sorted(sample_indices): result['selected_samples'].append(dict( time_over_M=float(data[index, 0]/params['mass']), AB=modes[index].tolist(), AB246_sum=float(np.sum(reduced[index])))) if params['lmin'] == 0: result['angular_partial_sums'] = [ dict(lmax=int(l), AB246=float(np.sum(reduced[-1, :i+1]))) for i, l in enumerate(ell) if l in (8, 12, 16)] if reference is not None: reference_modes = np.array([reference[int(l)] for l in ell]) result['frequency_domain_comparison'] = dict( reference_AB=reference_modes.tolist(), AB_error=(final-reference_modes).tolist(), max_AB_error=float(np.max(abs(final-reference_modes))), sum_AB_error=float(np.sum(final-reference_modes)), same_cutoff_reference_AB246=float(np.sum(reference_modes-correction))) return result if __name__ == '__main__': parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--runs', type=Path, required=True) parser.add_argument('--parameters', type=Path, required=True) parser.add_argument('--reference-r6', type=Path) parser.add_argument('--reference-r10', type=Path) parser.add_argument('--output', type=Path, required=True) parser.add_argument('--probes', type=Path) parser.add_argument('--probe-output', type=Path) parser.add_argument('names', nargs='+') args = parser.parse_args() coefficients = regularization(args.parameters) references = {} for radius, path in ((6, args.reference_r6), (10, args.reference_r10)): if path: references[radius] = read_reference(path, radius) results = [] for name in args.names: directory = args.runs/name radius = supported_radius(read_parameters(directory/'input.par')) results.append(analyze(directory, coefficients, references.get(radius))) output = dict(normalization='G=c=1; M=q=1; values denote M^2 Fr/q^2 for the Schwarzschild covariant radial self-force', method='Time-domain discontinuous Galerkin; fixed second-order puncture', angular_postprocessing='Subtract analytic F2, F4, F6 from raw residual modes', shared_data_caveat='The wave solves are independent; their analytic regularization data are shared.', convergence=convergence_summary(results),runs=results) args.output.write_text(json.dumps(output, indent=2)+'\n') if args.probes: if not args.probe_output: parser.error('--probes requires --probe-output') args.probe_output.write_text(json.dumps(check_probes(args.probes, coefficients), indent=2)+'\n') for result in results: comparison = result.get('frequency_domain_comparison', {}) print(result['name'], result['last_time_over_M'], result.get('angular_partial_sums', []), comparison.get('max_AB_error'), comparison.get('sum_AB_error'))