#!/usr/bin/env python # Copyright Barry Wardell, Ian Hinder, Eloisa Bentivegna # We ask that if you make use of the parameter file or the example # data, then please cite # Simulation of GW150914 binary black hole merger using the # Einstein Toolkit - https://doi.org/10.5281/zenodo.155394 # as well as the Einstein Toolkit, the Llama multi-block # infrastructure, the Carpet mesh-refinement driver, the apparent # horizon finder AHFinderDirect, the TwoPunctures initial data code, # QuasiLocalMeasures, Cactus, and the McLachlan spacetime evolution # code, the Kranc code generation package, and the Simulation Factory. # An appropriate bibtex file, etgw150914.bib, is provided with this # parameter file. from math import * import sys import re from string import Template ################################################################################ # Binary black hole configuration ################################################################################ # BHs labeled as '+' and '-' ('p' and 'm') for their initial position # on the x axis. The more massive BH is '+'. D = 10.0 # Separation q = 36.0/29.0 # Mass ratio: q = mp/mm >= 1 M = 1.0 # Total mass chip = [0, 0, 0.31] # Dimensionsless spin of + BH (x0 > 0, more massive) chim = [0, 0, -0.46] # Dimensionsless spin of - BH (x0 < 0, less massive) Pr = -0.00084541526517121; # Radial linear momentum Pphi = 0.09530152296974252; # Azimuthal linear momentum ################################################################################ # Initial data ################################################################################ mp = M * q/(1+q) # Heavier, larger BH, AH1, SS 0 mm = M * 1/(1+q) # Lighter, smaller BH, AH2, SS 1 xp = D * mm xm = -D * mp half_D = D/2.0 center_offset = xp - half_D Spx = chip[0] * mp**2; Spy = chip[1] * mp**2; Spz = chip[2] * mp**2; Smx = chim[0] * mm**2; Smy = chim[1] * mm**2; Smz = chim[2] * mm**2; Ppx = Pr Ppy = Pphi Ppz = 0 Pmx = -Ppx Pmy = -Ppy Pmz = -Ppz ################################################################################ # Grid structure ################################################################################ sphere_inner_radius = 45 expected_merger = 1000.0 # Number of cells across finest grid radius n = int("@N@") if "@N@"[0] != "@" else 28 i = n/4 n_angular = 2*i*2 # Estimated eventual AH radii ahrp = mp * 1.0 ahrm = mm * 1.0 # Coordinate size of finest boxes around each BH rp = ahrp * 1.2 rm = ahrm * 1.2 # Minimum acceptable radial grid spacing hr_min = 2*0.96*M maxrls = 9 # Essentially determines iteration counting # Smaller '-' BH # Decisions are made independent of n, to avoid roundoff problems. # This is achieved by using nmin for the calculations and scaling by # n/nmin at the end. n_min = 24 # Cells across AHm radius hfm_min = rm/n_min # Fine grid spacing for '-' BH rlsm = 1 + int(log(hr_min/hfm_min,2)) # Number of refinements to attain hr_min h0_min = hfm_min * 2**(rlsm-1) # RL0 Cartesian spacing # Radii of each level for the centre around the '-' BH levelsm = "[0," + ",".join(["%f" %(rm*2**(rlsm-l-1)) for l in range(1,rlsm)])+"]" # '+' BH # Number of levels for '+' BH estimated to keep error in angular # velocity the same as for the '-' BH, assuming vErr \propto v * # (h0/2**(rls-1)/rAH)**8 rlsp = rlsm + log(ahrm/ahrp * (xp/-xm)**(1./8))/log(2) rlsp = int(round(rlsp)) levelsp = "[0," + ",".join(["%f" %(rp*2**(rlsp-l-1)) for l in range(1,rlsp)])+"]" hr = h0_min * float(n_min)/float(n) # This could be changed (h0_min -> # hr_min) to get the minimum # radial resolution for efficiency h0 = h0_min * float(n_min)/float(n) time_after_merger = 200.0 waveform_length = expected_merger + time_after_merger outermost_detector = 500.0 final_time = waveform_length + outermost_detector sphere_outer_radius = int((outermost_detector + final_time)/(i*hr))*i*hr sphere_outer_radius = int(sphere_outer_radius / hr) * hr + hr # round up to a multiple of hr sphere_inner_radius = int(ceil(sphere_inner_radius/(i*h0))) * h0 * i ################################################################################ # Frequencies ################################################################################ dtdx = 0.45 # Courant factor dtdx0 = dtdx * 0.5 # RL0 is evolved with the same frequency as RL1 dt0 = dtdx0 * h0 # Time step on RL0 rl0_every = 2**(maxrls-2) rl1_every = rl0_every rl_max = maxrls-1 dt_it = dt0/2.0**(rl_max-1) # Time step of one iteration find_cah = max(int((expected_merger - 200)/dt_it), 0) wave_extract_every = rl0_every * 2 # Every other coarse grid step # (TODO: should this be every # coarse grid step?) horizon_every = rl0_every * 2 out_every = rl0_every # out3d_every = rl0_every * 2 out3d_every = 0 out2d_every = rl0_every * 8 ################################################################################ # Record all script variables in generated parameter file ################################################################################ local_vars = locals() col_width = 0 for key,val in sorted(local_vars.items()): if isinstance(val, (int, float, complex, str)): col_width = max(len(str(key)), col_width) var_settings = [] for key,val in sorted(local_vars.items()): if isinstance(val, (int, float, complex, str)): var_settings = var_settings + ["# {0}{1} = {2}".format(key," "*(col_width-len(key)),val)] var_settings_str = "\n".join(var_settings) lines = """ ################################################################################ # Script variables ################################################################################ $var_settings_str ################################################################################ # Active thorns ################################################################################ ActiveThorns = " ADMBase ML_ADMConstraints ADMCoupling ADMMacros AEILocalInterp AHFinderDirect Boundary Carpet CarpetIOASCII CarpetIOBasic CarpetIOHDF5 CarpetIOScalar CarpetInterp CarpetInterp2 CarpetLib CarpetReduce CarpetRegrid2 CarpetTracker CartGrid3D CoordBase CoordGauge Coordinates CoordinatesSymmetry Formaline GlobalDerivative hwloc IOUtil InitBase Interpolate2 QuasiLocalMeasures LocalInterp LoopControl MoL NaNChecker PunctureTracker Slab SpaceMask SphericalSurface StaticConformal SummationByParts SymBase SystemStatistics SystemTopology TerminationTrigger TensorTypes Time TmunuBase TwoPunctures Vectors ML_BSSN ML_BSSN_Helper NewRad GenericFD WeylScal4 Multipole # WaveExtractCPM ADMDerivatives " ################################################################################ # Grid structure ################################################################################ Carpet::domain_from_multipatch = yes CartGrid3D::type = "multipatch" CartGrid3D::set_coordinate_ranges_on = "all maps" Coordinates::coordinate_system = "Thornburg04" Coordinates::h_cartesian = $h0 Coordinates::h_radial = $hr Coordinates::sphere_inner_radius = $sphere_inner_radius Coordinates::sphere_outer_radius = $sphere_outer_radius Coordinates::n_angular = $n_angular Driver::ghost_size = 5 Coordinates::patch_boundary_size = 5 Coordinates::additional_overlap_size = 3 Coordinates::outer_boundary_size = 5 CoordinatesSymmetry::reflection_z = yes CoordinatesSymmetry::stagger = no Coordinates::symmetry = "+z bitant" Coordinates::additional_symmetry_size = 1 Coordinates::verbose = no Time::timestep_method = "given" Time::timestep = $dt0 Carpet::time_refinement_factors = "[1,1,2,4,8,16,32,64,128,256]" ################################################################################ # Mesh refinement ################################################################################ Carpet::max_refinement_levels = $maxrls CarpetRegrid2::num_centres = 2 CarpetRegrid2::num_levels_1 = $rlsp CarpetRegrid2::position_x_1 = $xp CarpetRegrid2::radius_1 = $levelsp CarpetRegrid2::num_levels_2 = $rlsm CarpetRegrid2::position_x_2 = $xm CarpetRegrid2::radius_2 = $levelsm Carpet::use_buffer_zones = yes Carpet::prolongation_order_space = 5 Carpet::prolongation_order_time = 2 CarpetRegrid2::regrid_every = $rl1_every CarpetRegrid2::verbose = no Carpet::grid_coordinates_filename = "carpet-grid.asc" ################################################################################ # Initial data ################################################################################ ADMBase::initial_data = "twopunctures" ADMBase::initial_lapse = "twopunctures-averaged" ADMBase::initial_shift = "zero" ADMBase::initial_dtlapse = "zero" ADMBase::initial_dtshift = "zero" TwoPunctures::target_M_plus = $mp TwoPunctures::target_M_minus = $mm # Initial guesses for iterative solve TwoPunctures::par_m_plus = $mp TwoPunctures::par_m_minus = $mm TwoPunctures::par_b = $half_D TwoPunctures::center_offset[0] = $center_offset TwoPunctures::par_P_plus[0] = $Ppx TwoPunctures::par_P_plus[1] = $Ppy TwoPunctures::par_P_plus[2] = $Ppz TwoPunctures::par_P_minus[0] = $Pmx TwoPunctures::par_P_minus[1] = $Pmy TwoPunctures::par_P_minus[2] = $Pmz TwoPunctures::par_S_plus[0] = $Spx TwoPunctures::par_S_plus[1] = $Spy TwoPunctures::par_S_plus[2] = $Spz TwoPunctures::par_S_minus[0] = $Smx TwoPunctures::par_S_minus[1] = $Smy TwoPunctures::par_S_minus[2] = $Smz TwoPunctures::grid_setup_method = "evaluation" TwoPunctures::give_bare_mass = no TwoPunctures::TP_epsilon = 1e-6 Carpet::init_fill_timelevels = yes InitBase::initial_data_setup_method = "init_single_level" ################################################################################ # Evolution and boundary ################################################################################ GenericFD::jacobian_group = "Coordinates::jacobian" GenericFD::jacobian_derivative_group = "Coordinates::jacobian2" GenericFD::jacobian_identity_map = 0 ADMBase::evolution_method = "ML_BSSN" ADMBase::lapse_evolution_method = "ML_BSSN" ADMBase::shift_evolution_method = "ML_BSSN" ADMBase::dtlapse_evolution_method = "ML_BSSN" ADMBase::dtshift_evolution_method = "ML_BSSN" ML_BSSN::fdOrder = 8 # 1+log slicing requires harmonicN = 1 and harmonicF = 2.0 ML_BSSN::harmonicN = 1 ML_BSSN::harmonicF = 2.0 ML_BSSN::ShiftGammaCoeff = 0.75 ML_BSSN::BetaDriver = 1.0 ML_BSSN::advectLapse = 1 ML_BSSN::advectShift = 1 ML_BSSN::MinimumLapse = 1.0e-8 # conformalaMethod = 1 for W, 0 for phi ML_BSSN::conformalMethod = 1 # We apply dissipation using GlobalDerivaitive so disable it here ML_BSSN::epsDiss = 0.0 ML_BSSN::dt_lapse_shift_method = "noLapseShiftAdvection" ML_BSSN::initial_boundary_condition = "extrapolate-gammas" ML_BSSN::rhs_boundary_condition = "scalar" Boundary::radpower = 2 ################################################################################ # BH tracking ################################################################################ CarpetTracker::surface [0] = 0 CarpetTracker::surface [1] = 1 PunctureTracker::track [0] = yes PunctureTracker::initial_x [0] = $xp PunctureTracker::which_surface_to_store_info[0] = 0 PunctureTracker::track [1] = yes PunctureTracker::initial_x [1] = $xm PunctureTracker::which_surface_to_store_info[1] = 1 ################################################################################ # Spatial finite differencing ################################################################################ SummationByParts::order = 8 # Drop order instead of using upwinded stencils, only for advection derivatives SummationByParts::sbp_upwind_deriv = no SummationByParts::sbp_1st_deriv = yes SummationByParts::sbp_2nd_deriv = no SummationByParts::onesided_interpatch_boundaries = no SummationByParts::onesided_outer_boundaries = yes SummationByParts::use_dissipation = no GlobalDerivative::use_dissipation = yes SummationByParts::scale_with_h = yes SummationByParts::dissipation_type = "Kreiss-Oliger" SummationByParts::epsdis = 0.15 # Because dt/dx is half that of the other levels we also need half the dissipation GlobalDerivative::epsdis_for_level [0] = 0.075 # Variables for dissipation SummationByParts::vars = " ML_BSSN::ML_log_confac ML_BSSN::ML_metric ML_BSSN::ML_trace_curv ML_BSSN::ML_curv ML_BSSN::ML_Gamma ML_BSSN::ML_lapse ML_BSSN::ML_shift ML_BSSN::ML_dtlapse ML_BSSN::ML_dtshift " ################################################################################ # Time integration ################################################################################ MoL::ODE_Method = "rk4" MoL::MoL_Intermediate_Steps = 4 MoL::MoL_Num_Scratch_Levels = 1 ################################################################################ # Interpolation ################################################################################ CarpetInterp::check_tree_search = no CarpetInterp::tree_search = yes # Use 5-th order interpatch interpolation on the Llama grid Interpolate::interpolator_order = 5 ################################################################################ # Psi4 computation ################################################################################ WeylScal4::fdOrder = 8 WeylScal4::calc_scalars = "psis" WeylScal4::calc_invariants = "always" ################################################################################ # Psi4 mode decomposition ################################################################################ # Radii are chosen to be evenly spaced in 1/r as that is the variable # extrapolation is performed in Multipole::nradii = 7 Multipole::radius[0] = 100 Multipole::radius[1] = 115 Multipole::radius[2] = 136 Multipole::radius[3] = 167 Multipole::radius[4] = 214 Multipole::radius[5] = 300 Multipole::radius[6] = 500 Multipole::ntheta = 120 Multipole::nphi = 240 Multipole::variables = "WeylScal4::Psi4r{sw=-2 cmplx='WeylScal4::Psi4i' name='psi4'}" Multipole::out_every = $wave_extract_every Multipole::l_max = 8 Multipole::output_hdf5 = yes # Disable ASCII output to avoid creating a large number of files Multipole::output_ascii = no ################################################################################ # Gauge invariant perturbations of Schwarzschild (ZM-CPM variables) ################################################################################ # WaveExtractCPM::out_every = $wave_extract_every # WaveExtractCPM::use_carpetinterp2 = no # WaveExtractCPM::calc_when_necessary = no # WaveExtractCPM::verbose = 0 # WaveExtractCPM::maximum_detector_number = 7 # WaveExtractCPM::switch_output_format = 100 # WaveExtractCPM::rsch2_computation = "average Schwarzschild metric" # WaveExtractCPM::l_mode = 8 # WaveExtractCPM::m_mode = 8 # WaveExtractCPM::detector_radius [0] = 100 # WaveExtractCPM::detector_radius [1] = 115 # WaveExtractCPM::detector_radius [2] = 136 # WaveExtractCPM::detector_radius [3] = 167 # WaveExtractCPM::detector_radius [4] = 214 # WaveExtractCPM::detector_radius [5] = 300 # WaveExtractCPM::detector_radius [6] = 500 # WaveExtractCPM::maxntheta = 120 # WaveExtractCPM::maxnphi = 240 # WaveExtractCPM::output_hdf5 = yes # WaveExtractCPM::output_ascii = no # WaveExtractCPM::output_h = yes # WaveExtractCPM::output_Psi = yes ################################################################################ # Apparent Horizons ################################################################################ AHFinderDirect::N_horizons = 3 AHFinderDirect::find_every = $horizon_every AHFinderDirect::output_h_every = 0 AHFinderDirect::max_Newton_iterations__initial = 50 AHFinderDirect::max_Newton_iterations__subsequent = 50 AHFinderDirect::max_allowable_Theta_growth_iterations = 10 AHFinderDirect::max_allowable_Theta_nonshrink_iterations = 10 AHFinderDirect::geometry_interpolator_name = "Lagrange polynomial interpolation" AHFinderDirect::geometry_interpolator_pars = "order=4" AHFinderDirect::surface_interpolator_name = "Lagrange polynomial interpolation" AHFinderDirect::surface_interpolator_pars = "order=4" AHFinderDirect::verbose_level = "physics details" AHFinderDirect::move_origins = yes AHFinderDirect::origin_x [1] = $xp AHFinderDirect::initial_guess__coord_sphere__x_center[1] = $xp AHFinderDirect::initial_guess__coord_sphere__radius [1] = $rp AHFinderDirect::which_surface_to_store_info [1] = 2 AHFinderDirect::set_mask_for_individual_horizon [1] = no AHFinderDirect::reset_horizon_after_not_finding [1] = no AHFinderDirect::track_origin_from_grid_scalar [1] = yes AHFinderDirect::track_origin_source_x [1] = "PunctureTracker::pt_loc_x[0]" AHFinderDirect::track_origin_source_y [1] = "PunctureTracker::pt_loc_y[0]" AHFinderDirect::track_origin_source_z [1] = "PunctureTracker::pt_loc_z[0]" AHFinderDirect::max_allowable_horizon_radius [1] = 3 AHFinderDirect::origin_x [2] = $xm AHFinderDirect::initial_guess__coord_sphere__x_center[2] = $xm AHFinderDirect::initial_guess__coord_sphere__radius [2] = $rm AHFinderDirect::which_surface_to_store_info [2] = 3 AHFinderDirect::set_mask_for_individual_horizon [2] = no AHFinderDirect::reset_horizon_after_not_finding [2] = no AHFinderDirect::track_origin_from_grid_scalar [2] = yes AHFinderDirect::track_origin_source_x [2] = "PunctureTracker::pt_loc_x[1]" AHFinderDirect::track_origin_source_y [2] = "PunctureTracker::pt_loc_y[1]" AHFinderDirect::track_origin_source_z [2] = "PunctureTracker::pt_loc_z[1]" AHFinderDirect::max_allowable_horizon_radius [2] = 3 AHFinderDirect::origin_x [3] = 0 AHFinderDirect::find_after_individual [3] = $find_cah AHFinderDirect::initial_guess__coord_sphere__x_center[3] = 0 AHFinderDirect::initial_guess__coord_sphere__radius [3] = 1.0 AHFinderDirect::which_surface_to_store_info [3] = 4 AHFinderDirect::set_mask_for_individual_horizon [3] = no AHFinderDirect::max_allowable_horizon_radius [3] = 6 ################################################################################ # Spherical surfaces ################################################################################ SphericalSurface::nsurfaces = 5 SphericalSurface::maxntheta = 66 SphericalSurface::maxnphi = 124 SphericalSurface::verbose = no # Surfaces 0 and 1 are used by PunctureTracker # Horizon 1 SphericalSurface::ntheta [2] = 41 SphericalSurface::nphi [2] = 80 SphericalSurface::nghoststheta [2] = 2 SphericalSurface::nghostsphi [2] = 2 # Horizon 2 SphericalSurface::ntheta [3] = 41 SphericalSurface::nphi [3] = 80 SphericalSurface::nghoststheta [3] = 2 SphericalSurface::nghostsphi [3] = 2 # Common horizon SphericalSurface::ntheta [4] = 41 SphericalSurface::nphi [4] = 80 SphericalSurface::nghoststheta [4] = 2 SphericalSurface::nghostsphi [4] = 2 ################################################################################ # Isolated Horizons ################################################################################ QuasiLocalMeasures::verbose = no QuasiLocalMeasures::veryverbose = no QuasiLocalMeasures::interpolator = "Lagrange polynomial interpolation" QuasiLocalMeasures::interpolator_options = "order=4" QuasiLocalMeasures::spatial_order = 4 QuasiLocalMeasures::num_surfaces = 3 QuasiLocalMeasures::surface_index [0] = 2 QuasiLocalMeasures::surface_index [1] = 3 QuasiLocalMeasures::surface_index [2] = 4 QuasiLocalMeasures::output_vtk_every = $out3d_every ################################################################################ # Correctness checking ################################################################################ Carpet::poison_new_timelevels = no Carpet::check_for_poison = no NaNChecker::check_every = 256 NanChecker::check_after = 0 NaNChecker::report_max = 10 NaNChecker::verbose = "all" NaNChecker::action_if_found = terminate NaNChecker::out_NaNmask = yes NaNChecker::check_vars = " ML_BSSN::ML_log_confac " ################################################################################ # Timers ################################################################################ Carpet::output_timer_tree_every = 1024 Carpet::output_initialise_timer_tree = yes ################################################################################ # Output ################################################################################ IO::out_dir = "@SIMULATION_NAME@" IOScalar::one_file_per_group = yes IOASCII::one_file_per_group = yes IOBasic::outInfo_every = 1 IOBasic::outInfo_reductions = "minimum maximum" IOBasic::outInfo_vars = " ML_BSSN::ML_log_confac Carpet::physical_time_per_hour SystemStatistics::maxrss_mb SystemStatistics::swap_used_mb " IOScalar::outScalar_every = 256 IOScalar::outScalar_reductions = "minimum maximum average" IOScalar::outScalar_vars = "SystemStatistics::process_memory_mb" IOASCII::out0D_every = 256 IOASCII::out0D_vars = " #Carpet::timing #PunctureTracker::pt_loc #QuasiLocalMeasures::qlm_scalars{out_every = $horizon_every} " IOASCII::out1D_every = 0 IOASCII::out1D_d = no IOASCII::out1D_x = yes IOASCII::out1D_y = no IOASCII::out1D_z = yes IOASCII::out1D_vars = " #ML_BSSN::ML_log_confac #ML_BSSN::ML_trace_curv #WeylScal4::Psi4r " IOASCII::out2D_every = 0 IOASCII::out2D_vars = " " IOHDF5::out_every = 1#$out3d_every IOHDF5::out_vars = " #Grid::Coordinates{out_every=1000000000 refinement_levels={0}} QuasiLocalMeasures::qlm_twometric #ML_BSSN::ML_log_confac #WeylScal4::Psi4r #WeylScal4::Psi4i #WeylScal4::curvIr{refinement_levels={3 5}} #WeylScal4::curvIi{refinement_levels={3 5}} #WeylScal4::curvJr{refinement_levels={3 5}} #WeylScal4::curvJi{refinement_levels={3 5}} " IOHDF5::out2D_every = $out2d_every IOHDF5::out2D_vars = " #Grid::Coordinates{out_every=1000000000 refinement_levels={0}} #ML_BSSN::alpha #ML_BSSN::phi #ML_BSSN::trK #WeylScal4::Psi4r #WeylScal4::Psi4i " ################################################################################ # Checkpointing and recovery ################################################################################ CarpetIOHDF5::checkpoint = yes IO::checkpoint_ID = no IO::recover = "autoprobe" IO::out_proc_every = 2 IO::checkpoint_on_terminate = yes IO::checkpoint_dir = "../checkpoints" IO::recover_dir = "../checkpoints" IO::abort_on_io_errors = yes CarpetIOHDF5::open_one_input_file_at_a_time = yes CarpetIOHDF5::compression_level = 7 ################################################################################ # Run termination ################################################################################ TerminationTrigger::max_walltime = @WALLTIME_HOURS@ # Trigger termination 30 minutes before the walltime is reached TerminationTrigger::on_remaining_walltime = 30 TerminationTrigger::output_remtime_every_minutes = 30 TerminationTrigger::termination_from_file = yes TerminationTrigger::termination_file = "terminate.txt" TerminationTrigger::create_termination_file = yes #Cactus::terminate = time #Cactus::cctk_final_time = $final_time Cactus::cctk_itlast = 4 """ open(re.sub(r'(.*)\.rpar$', r'\1.par', sys.argv[0]), 'w').write(re.sub(r'\n *',r'\n',Template(Template(lines).substitute(locals())).substitute(locals())))