"""Run isolated physical SelfForce-1D parameter variants; no global installs. Upstream numerical code is unmodified. Every run retains its Fortran namelist, log and per-l force output. Puncture order is fixed in the upstream code; the runner varies numerical resolution, worldtube width and the coordinate mapping. """ import argparse import json import os from pathlib import Path import subprocess import sys import shutil import time import re sys.dont_write_bytecode = True ROOT = Path(__file__).resolve().parent PACKAGE = ROOT EXE = ROOT / 'SelfForce-1D/Build/Exe/test.x' BASE = { 'equation_name': 'scalar_schwarzschild', 'n_elems': 32, 'order': 10, 'Sminus': -20.0, 't_size': 5, 'p_orb': 10.0, 'ecc': 0.0, 'mass': 1.0, 'use_particle': True, 'use_constant_acceleration': True, 'omega_ratio': 1.0, 'use_generic_orbit': False, 'evolve_orbit': False, 'lmin': 0, 'lmax': 16, 't_initial': 0.0, 't_final': 1000.0, 'use_world_tube': True, 'world_tube_width': 1, 'use_exact_initial_data': False, 'turn_on_source_smoothly': True, 'tsigma': 10.0, 'torder': 4, 'coufac': 0.5, 'mol_integrator': 'rk4', 'out0d_every': 40, 'out1d_every': -1, 'q_charge': 1.0, 'q_mass': 1.0, 'fit_high_l': True, 'first_l': 10, 'last_l': 16, 'fit_order': 2, 'nfit': 2, } def fortran(value): if isinstance(value, bool): return '.true.' if value else '.false.' if isinstance(value, str): return repr(value) return str(value) def namelist(params): return '¶ms\n' + '\n'.join(key + '=' + fortran(value) for key, value in params.items()) + '\n/\n' def seed_monopole(directory, params, threads): """Exact static l=m=0 retarded field, using upstream's input-data interface. (r^2 f Phi00')'=-4 pi q Y00 delta(r-r0)/ut. The input file is for unit q; the upstream reader applies q itself. Rows on both sides of the particle are separate, so their one-sided derivatives are kept distinct. """ import numpy as np if params['lmin'] != 0: raise ValueError('Monopole initialization requires lmin=0') grid_dir=directory/'initial-grid' grid_dir.mkdir(exist_ok=True) grid_params=dict(params,lmax=0,fit_high_l=False,t_final=0.0, output_coords_for_exact=True,use_exact_initial_data=False, turn_on_source_smoothly=True) (grid_dir/'input.par').write_text(namelist(grid_params)) with (grid_dir/'grid.log').open('w') as log: subprocess.run([str(EXE),'input.par'],cwd=grid_dir,stdout=log, stderr=subprocess.STDOUT,check=True, env=dict(os.environ,OMP_NUM_THREADS=str(threads))) coords=np.loadtxt(grid_dir/'coords.asc')[:-1] mass=params['mass'];r0=params['p_orb']*mass ut=1/(1-3*mass/r0)**0.5 coefficient=-2*np.pi/(np.sqrt(4*np.pi)*mass*ut) middle=params['n_elems']//2*(params['order']+1) radii=coords[:,2]+2*mass field=np.full(len(coords),coefficient*np.log1p(-2*mass/r0)) radial=np.zeros(len(coords)) field[middle:]=coefficient*np.log1p(-2*mass/radii[middle:]) radial[middle:]=coefficient*2*mass/radii[middle:]**2 zeros=np.zeros(len(coords)) np.savetxt(directory/'monopole_l0m0.dat', np.stack([coords[:,0],field,zeros,radial,zeros,zeros,zeros],axis=1), fmt='%.17e') params.update(use_exact_initial_data=True,exact_initial_data_lmax=0, input_directory='.',input_basename='monopole', turn_on_source_smoothly=False) def run(name, changes, threads, monopole=False): directory = ROOT / 'runs' / name directory.mkdir(parents=True, exist_ok=True) if (directory / 'run.log').exists(): raise RuntimeError('Refusing to overwrite an existing run ' + name) params = dict(BASE, **changes) if monopole: seed_monopole(directory,params,threads) (directory / 'input.par').write_text(namelist(params)) start = time.monotonic() with (directory / 'run.log').open('w') as log: result = subprocess.run([str(EXE), 'input.par'], cwd=directory, stdout=log, stderr=subprocess.STDOUT, env=dict(os.environ, OMP_NUM_THREADS=str(threads))) result.check_returncode() import numpy as np data=np.loadtxt(directory/'Schwarzschild_wave.frl.asc',ndmin=2) ncols=1+params['lmax']-params['lmin']+1+(3 if params['fit_high_l'] else 0) if data.shape[1]!=ncols or not np.all(np.isfinite(data)): raise RuntimeError('Invalid or non-finite radial force output') match=re.search(r'delta_rho_min\s*=\s*(\S+)',(directory/'run.log').read_text()) dt=params['coufac']*float(match.group(1)) if params['t_final']-data[-1,0] > (params['out0d_every']+2)*dt: raise RuntimeError('Solver stopped before the requested final time') info = {'seconds': time.monotonic() - start, 'returncode': result.returncode, 'threads': threads, 'parameters': params,'dt':dt, 'last_output_time':float(data[-1,0]),'seeded_monopole':monopole} (directory / 'timing.json').write_text(json.dumps(info, indent=2) + '\n') print(name, json.dumps(info), flush=True) if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('name', nargs='?') parser.add_argument('--changes', default='{}') parser.add_argument('--threads', type=int, default=2) parser.add_argument('--seed-monopole',action='store_true') parser.add_argument('--case', action='store_true', help='Read the named archived case from cases.json') parser.add_argument('--work',type=Path,default=ROOT) parser.add_argument('--executable',type=Path) parser.add_argument('--matrix', action='store_true', help='Run every archived case and regenerate selected results') parser.add_argument('--resume', action='store_true', help='Reuse completed matrix runs with identical parameters') parser.add_argument('--output-directory',type=Path) parser.add_argument('--parameters',type=Path, default=PACKAGE.parent/'scalar-mode-sum/scalar-parameters-60digits.json') parser.add_argument('--reference',type=Path, default=PACKAGE.parent/'scalar-mode-sum/checks.json') args = parser.parse_args() ROOT=args.work.resolve() EXE=(args.executable.resolve() if args.executable else ROOT/'SelfForce-1D/Build/Exe/test.x') if args.matrix: if args.name or args.case or args.changes!='{}' or args.seed_monopole: parser.error('--matrix cannot be combined with a named case or overrides') for path in (args.parameters,args.reference,ROOT/'environment.json', ROOT/'source-probes.csv',EXE): if not path.is_file(): parser.error('Matrix prerequisite missing: '+str(path)) cases=json.loads((PACKAGE/'cases.json').read_text())['cases'] from analyze_runs import read_parameters for case in cases: directory=ROOT/'runs'/case['name'] timing_path=directory/'timing.json' log_path=directory/'run.log' old_timed_log=(log_path.exists() and re.search(r'^real\s+\d+(?:\.\d+)?\s*$',log_path.read_text(),re.M)) if args.resume and (timing_path.exists() or old_timed_log): expected=dict(BASE,**case['changes']) if case.get('seed_monopole'): expected.update(use_exact_initial_data=True,exact_initial_data_lmax=0, input_directory='.',input_basename='monopole', turn_on_source_smoothly=False) if read_parameters(directory/'input.par')!=expected: raise RuntimeError('Existing case parameters differ: '+case['name']) if timing_path.exists() and json.loads(timing_path.read_text())['returncode']!=0: raise RuntimeError('Existing case did not finish: '+case['name']) continue run(case['name'],case['changes'],args.threads,case.get('seed_monopole',False)) destination=(args.output_directory or ROOT).resolve() destination.mkdir(parents=True,exist_ok=True) subprocess.run([sys.executable,str(PACKAGE/'analyze_runs.py'), '--runs',str(ROOT/'runs'),'--parameters',str(args.parameters), '--reference-r6',str(args.reference),'--reference-r10',str(args.reference), '--output',str(destination/'results.json'), '--probes',str(ROOT/'source-probes.csv'), '--probe-output',str(destination/'source-checks.json'), *[case['name'] for case in cases]],check=True) for name in ('environment.json','source-probes.csv'): if (ROOT/name).resolve()!=(destination/name).resolve(): shutil.copy2(ROOT/name,destination/name) sys.exit(0) if not args.name: parser.error('Specify a run name, or --matrix') changes = json.loads(args.changes) if args.case: cases = json.loads((PACKAGE/'cases.json').read_text())['cases'] case = next((row for row in cases if row['name'] == args.name), None) if case is None: parser.error('Unknown archived case: ' + args.name) changes = dict(case['changes'], **changes) args.seed_monopole = args.seed_monopole or case.get('seed_monopole', False) run(args.name, changes, args.threads,args.seed_monopole)