From 83b6898036608da56fe70d8654019a38f9dfca81 Mon Sep 17 00:00:00 2001 From: schnetter Date: Fri, 21 Sep 2007 15:47:13 +0000 Subject: Rewrite puncture tracker git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/PunctureTracker/trunk@6 a2659f00-0f4f-0410-9214-a4596bbb8a4f --- README | 9 +- interface.ccl | 19 +- par/qc0-punctrac.par | 638 +++++++++++++++++++++++++++++++++++++++++++++++++ param.ccl | 55 +---- schedule.ccl | 38 ++- src/adjust.c | 90 ------- src/make.code.defn | 5 +- src/punc_track.c | 170 ------------- src/puncture_tracker.c | 285 ++++++++++++++++++++++ 9 files changed, 962 insertions(+), 347 deletions(-) create mode 100644 par/qc0-punctrac.par delete mode 100644 src/adjust.c delete mode 100644 src/punc_track.c create mode 100644 src/puncture_tracker.c diff --git a/README b/README index 14af522..dfaeac4 100644 --- a/README +++ b/README @@ -1,11 +1,10 @@ -CVS info : $Header$ - Cactus Code Thorn PunctureTracker Thorn Author(s) : Michael Koppitz -Thorn Maintainer(s) : Michael Koppitz + Erik Schnetter +Thorn Maintainer(s) : Erik Schnetter -------------------------------------------------------------------------- Purpose of the thorn: -time integrate shift at location of puncture to follow location of puncture - +Time integrate the shift at the location of the punctures to follow the +location of the punctures. diff --git a/interface.ccl b/interface.ccl index 8c79d3f..beac81a 100644 --- a/interface.ccl +++ b/interface.ccl @@ -1,20 +1,17 @@ # Interface definition for thorn PunctureTracker -# $Header$ IMPLEMENTS: PunctureTracker -INHERITS: ADMBase, TwoPunctures,SphericalSurface +INHERITS: ADMBase CarpetRegrid2 -CCTK_REAL pt_loc[pt_num_tracked] TYPE=scalar timelevels=2 -{ - pt_loc_x, pt_loc_y, pt_loc_z - -} "location of punctures" -CCTK_REAL pt_shift[pt_num_tracked] TYPE=scalar +CCTK_REAL pt_loc[10] TYPE=scalar { - pt_shiftx, pt_shifty, pt_shiftz - -} "shift at location of punctures" + pt_loc_t pt_loc_x pt_loc_y pt_loc_z +} "Location of punctures" +CCTK_REAL pt_loc_p[10] TYPE=scalar +{ + pt_loc_t_p pt_loc_x_p pt_loc_y_p pt_loc_z_p +} "Previous location of punctures" diff --git a/par/qc0-punctrac.par b/par/qc0-punctrac.par new file mode 100644 index 0000000..19ba635 --- /dev/null +++ b/par/qc0-punctrac.par @@ -0,0 +1,638 @@ +Cactus::cctk_run_title = "QC-0" + +Cactus::cctk_full_warnings = yes +Cactus::highlight_warning_messages = no +Cactus::cctk_timer_output = "full" + +Cactus::terminate = "time" +#Cactus::cctk_final_time = 1.0 +#Cactus::cctk_final_time = 10.0 +Cactus::cctk_final_time = 100.0 + + + +ActiveThorns = "IOUtil" + +IO::out_dir = $parfile + + + +ActiveThorns = "ManualTermination" + +ManualTermination::max_walltime = @WALLTIME_HOURS@ # hours +ManualTermination::on_remaining_walltime = 60 # minutes + +ManualTermination::termination_from_file = yes +ManualTermination::create_termination_file = yes +ManualTermination::termination_file = "../TERMINATE" + + + +ActiveThorns = "AEILocalInterp" + +ActiveThorns = "Fortran" + +ActiveThorns = "LocalInterp" + +ActiveThorns = "MPIClock" + +#ActiveThorns = "NaNCatcher" + +ActiveThorns = "Slab" + +ActiveThorns = "TGRtensor" + + + +ActiveThorns = "SphericalSurface" + +SphericalSurface::nsurfaces = 6 +SphericalSurface::maxntheta = 39 +SphericalSurface::maxnphi = 76 + +SphericalSurface::ntheta [0] = 39 +SphericalSurface::nphi [0] = 76 +SphericalSurface::nghoststheta[0] = 2 +SphericalSurface::nghostsphi [0] = 2 + +SphericalSurface::ntheta [1] = 39 +SphericalSurface::nphi [1] = 76 +SphericalSurface::nghoststheta[1] = 2 +SphericalSurface::nghostsphi [1] = 2 + +SphericalSurface::ntheta [2] = 39 +SphericalSurface::nphi [2] = 76 +SphericalSurface::nghoststheta[2] = 2 +SphericalSurface::nghostsphi [2] = 2 + +SphericalSurface::ntheta [3] = 39 +SphericalSurface::nphi [3] = 76 +SphericalSurface::nghoststheta[3] = 2 +SphericalSurface::nghostsphi [3] = 2 + +SphericalSurface::ntheta [4] = 39 +SphericalSurface::nphi [4] = 76 +SphericalSurface::nghoststheta[4] = 2 +SphericalSurface::nghostsphi [4] = 2 + +SphericalSurface::ntheta [5] = 39 +SphericalSurface::nphi [5] = 76 +SphericalSurface::nghoststheta[5] = 2 +SphericalSurface::nghostsphi [5] = 2 + + + +ActiveThorns = "SummationByParts" + +SummationByParts::order = 4 + + + +ActiveThorns = "InitBase" + + + +ActiveThorns = "Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab" + +Carpet::verbose = no +Carpet::schedule_barriers = no +Carpet::veryverbose = no +CarpetLib::output_bboxes = no + +Carpet::domain_from_coordbase = yes +Carpet::max_refinement_levels = 10 + +driver::ghost_size = 3 +Carpet::use_buffer_zones = yes + +Carpet::prolongation_order_space = 5 +Carpet::prolongation_order_time = 2 + +Carpet::convergence_level = 0 + +Carpet::init_3_timelevels = yes + +Carpet::poison_new_timelevels = yes +CarpetLib::poison_new_memory = yes + +Carpet::output_timers_every = 512 +CarpetLib::print_timestats_every = 512 +CarpetLib::timestat_timer = "rdtsc" +CarpetLib::print_memstats_every = 512 + + + +ActiveThorns = "NaNChecker" + +NaNChecker::check_every = 512 +NaNChecker::action_if_found = "just warn" +NaNChecker::check_vars = " + ADM_BSSN::ADM_BSSN_phi + ADM_BSSN::ADM_BSSN_metric + ADM_BSSN::ADM_BSSN_curv + ADM_BSSN::ADM_BSSN_K + ADM_BSSN::ADM_BSSN_gamma + ADMBase::lapse + ADMBase::shift + ADM_BSSN::ADM_BSSN_B +" + + + +ActiveThorns = "Boundary CartGrid3D CoordBase ReflectionSymmetry RotatingSymmetry180 SymBase" + +CoordBase::domainsize = "minmax" + +CoordBase::xmin = 0.00 +CoordBase::ymin = -128.00 +CoordBase::zmin = 0.00 +CoordBase::xmax = +128.00 +CoordBase::ymax = +128.00 +CoordBase::zmax = +128.00 +CoordBase::dx = 2.56 +CoordBase::dy = 2.56 +CoordBase::dz = 2.56 + +CoordBase::boundary_size_x_lower = 3 +CoordBase::boundary_size_z_lower = 3 +CoordBase::boundary_shiftout_x_lower = 1 +CoordBase::boundary_shiftout_z_lower = 1 + +CartGrid3D::type = "coordbase" + +ReflectionSymmetry::reflection_z = yes +ReflectionSymmetry::avoid_origin_z = no + + + +ActiveThorns = "CarpetMask" + +CarpetMask::verbose = yes + +CarpetMask::excluded_surface [0] = 0 +CarpetMask::excluded_surface_factor[0] = 1.0 + +CarpetMask::excluded_surface [1] = 1 +CarpetMask::excluded_surface_factor[1] = 1.0 + +CarpetMask::excluded_surface [2] = 2 +CarpetMask::excluded_surface_factor[2] = 1.0 + + + +ActiveThorns = "CarpetRegrid2" + +CarpetRegrid2::regrid_every = 128 + +CarpetRegrid2::symmetry_rotating180 = yes + +CarpetRegrid2::num_centres = 2 + +CarpetRegrid2::num_levels_1 = 7 +CarpetRegrid2::position_x_1 = +1.168642873 +CarpetRegrid2::radius_1[ 1] = 64.0 # 1.28 +CarpetRegrid2::radius_1[ 2] = 16.0 # 0.64 +CarpetRegrid2::radius_1[ 3] = 8.0 # 0.32 +CarpetRegrid2::radius_1[ 4] = 4.0 # 0.16 +CarpetRegrid2::radius_1[ 5] = 2.0 # 0.08 +CarpetRegrid2::radius_1[ 6] = 1.0 # 0.04 +CarpetRegrid2::movement_threshold_1 = 0.16 + +CarpetRegrid2::num_levels_2 = 7 +CarpetRegrid2::position_x_2 = -1.168642873 +CarpetRegrid2::radius_2[ 1] = 64.0 # 1.28 +CarpetRegrid2::radius_2[ 2] = 16.0 # 0.64 +CarpetRegrid2::radius_2[ 3] = 8.0 # 0.32 +CarpetRegrid2::radius_2[ 4] = 4.0 # 0.16 +CarpetRegrid2::radius_2[ 5] = 2.0 # 0.08 +CarpetRegrid2::radius_2[ 6] = 1.0 # 0.04 +CarpetRegrid2::movement_threshold_2 = 0.16 + + + +ActiveThorns = "MoL Time" + +MoL::ODE_Method = "RK4" +MoL::MoL_Intermediate_Steps = 4 +MoL::MoL_Num_Scratch_Levels = 1 + +Time::dtfac = 0.40 + + + +ActiveThorns = "ADMBase ADMCoupling ADMMacros CoordGauge SpaceMask StaticConformal TmunuBase" + +ADMMacros::spatial_order = 4 + + + +ActiveThorns = "TwoPunctures" + +ADMBase::metric_type = "physical" + +ADMBase::initial_data = "twopunctures" +ADMBase::initial_lapse = "twopunctures-averaged" +ADMBase::initial_shift = "zero" + +TwoPunctures::par_b = 1.168642873 +TwoPunctures::par_m_plus = 0.453 +TwoPunctures::par_m_minus = 0.453 +TwoPunctures::par_P_plus [1] = +0.3331917498 +TwoPunctures::par_P_minus[1] = -0.3331917498 + +#TODO# TwoPunctures::grid_setup_method = "evaluation" + +#TwoPunctures::TP_epsilon = 1.0e-4 +TwoPunctures::TP_epsilon = 1.0e-2 +TwoPunctures::TP_Tiny = 1.0e-2 + +TwoPunctures::verbose = yes + + + +ActiveThorns = "BSSN_MoL" + +ADMBase::evolution_method = "ADM_BSSN" + +ADM_BSSN::timelevels = 3 +ADM_BSSN::stencil_size = 3 +ADM_BSSN::advection = "upwind4" +ADM_BSSN::bound = "newrad" + +ADMBase::lapse_evolution_method = "1+log" +ADM_BSSN::lapsesource = "straight" +ADM_BSSN::harmonic_f = 2.0 +ADM_BSSN::lapse_advection_coeff = 1.0 +ADM_BSSN::force_lapse_positive = yes + +ADMBase::shift_evolution_method = "gamma0" +ADM_BSSN::ShiftGammaCoeff = 0.75 +ADM_BSSN::BetaDriver = 0.5 +ADM_BSSN::gamma_driver_advection_coeff = 1.0 +ADM_BSSN::ApplyShiftBoundary = yes + + + +ActiveThorns = "Dissipation" + +Dissipation::order = 5 +Dissipation::vars = " + ADM_BSSN::ADM_BSSN_phi + ADM_BSSN::ADM_BSSN_metric + ADM_BSSN::ADM_BSSN_curv + ADM_BSSN::ADM_BSSN_K + ADM_BSSN::ADM_BSSN_gamma + ADMBase::lapse + ADMBase::shift + ADM_BSSN::ADM_BSSN_B +" + + + +ActiveThorns = "ADMConstraints" + +ADMConstraints::constraints_persist = yes +ADMConstraints::constraints_timelevels = 3 + + + +ActiveThorns = "AHFinderDirect" + +AHFinderDirect::find_every = 128 + +AHFinderDirect::run_at_CCTK_POST_RECOVER_VARIABLES = no + +AHFinderDirect::move_origins = yes +AHFinderDirect::reshape_while_moving = yes +AHFinderDirect::predict_origin_movement = yes + +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::output_h_every = 0 + +AHFinderDirect::N_horizons = 6 + +AHFinderDirect::origin_x [1] = +1.168642873 +AHFinderDirect::initial_guess__coord_sphere__x_center[1] = +1.168642873 +AHFinderDirect::initial_guess__coord_sphere__radius [1] = 0.25 +AHFinderDirect::which_surface_to_store_info [1] = 0 +AHFinderDirect::reset_horizon_after_not_finding [1] = no +AHFinderDirect::dont_find_after_individual [1] = 24000 # t=30.0 + +AHFinderDirect::origin_x [2] = -1.168642873 +AHFinderDirect::initial_guess__coord_sphere__x_center[2] = -1.168642873 +AHFinderDirect::initial_guess__coord_sphere__radius [2] = 0.25 +AHFinderDirect::which_surface_to_store_info [2] = 1 +AHFinderDirect::reset_horizon_after_not_finding [2] = no +AHFinderDirect::dont_find_after_individual [2] = 24000 + +AHFinderDirect::initial_guess__coord_sphere__radius [3] = 1.0 +AHFinderDirect::which_surface_to_store_info [3] = 2 +AHFinderDirect::reset_horizon_after_not_finding [3] = no +AHFinderDirect::find_after_individual [3] = 12000 # t=15.0 + +AHFinderDirect::surface_definition [4] = "expansion product" +AHFinderDirect::surface_selection [4] = "areal radius" +AHFinderDirect::desired_value [4] = 50.0 +AHFinderDirect::initial_guess__coord_sphere__radius [4] = 50.0 +AHFinderDirect::which_surface_to_store_info [4] = 3 +AHFinderDirect::reset_horizon_after_not_finding [4] = no + +AHFinderDirect::depends_on [5] = 1 +AHFinderDirect::desired_value_offset [5] = 0.001 +AHFinderDirect::which_surface_to_store_info [5] = 4 +AHFinderDirect::reset_horizon_after_not_finding [5] = no +AHFinderDirect::dont_find_after_individual [5] = 24000 # t=30.0 + +AHFinderDirect::depends_on [6] = 3 +AHFinderDirect::desired_value_offset [6] = 0.001 +AHFinderDirect::which_surface_to_store_info [6] = 5 +AHFinderDirect::reset_horizon_after_not_finding [6] = no +AHFinderDirect::find_after_individual [6] = 12000 # t=15.0 + + + +ActiveThorns = "IsolatedHorizon" + +IsolatedHorizon::verbose = yes +IsolatedHorizon::interpolator = "Lagrange polynomial interpolation" +IsolatedHorizon::interpolator_options = "order=4" +IsolatedHorizon::spatial_order = 6 + +IsolatedHorizon::num_horizons = 6 + +IsolatedHorizon::surface_index [0] = 0 +IsolatedHorizon::companion_index[0] = 4 + +IsolatedHorizon::surface_index [1] = 1 + +IsolatedHorizon::surface_index [2] = 2 + +IsolatedHorizon::surface_index [3] = 3 +IsolatedHorizon::companion_index[3] = 5 + +IsolatedHorizon::surface_index [4] = 4 + +IsolatedHorizon::surface_index [5] = 5 + + + +ActiveThorns = "ProperDistance ProperTime" + +ProperDistance::number_geodesics = 512 +ProperDistance::direction = "x" +ProperDistance::plane = "xy" +ProperDistance::opening_angle = 180.0 +ProperDistance::step_size = 0.04 +ProperDistance::horizon_number = 1 +ProperDistance::interpolation_order = 3 +ProperDistance::integration_method = "rk4" +ProperDistance::eps = 0.00001 +ProperDistance::max_proper_distance = 15 +ProperDistance::calc_every = 128 +ProperDistance::use_second_horizon = yes +ProperDistance::second_horizon_number = 2 + + + +ActiveThorns = "PsiKadelia" + +PsiKadelia::psikadelia_persists = yes +PsiKadelia::calc_every = 256 # only on levels 0 and 1 +PsiKadelia::ricci_prolongation_type = "none" +PsiKadelia::weyl_timelevels = 3 +PsiKadelia::psif_vec = "standard-radial" + + + +ActiveThorns = "PunctureTracker" + +PunctureTracker::verbose = yes + +PunctureTracker::track [0] = yes +PunctureTracker::initial_x[0] = +1.168642873 + +PunctureTracker::track [1] = yes +PunctureTracker::initial_x[1] = -1.168642873 + + + +ActiveThorns = "SphericalHarmonics" + +SphericalHarmonics::grid_type = "cart3d" +SphericalHarmonics::interp_integration_order = 4 +SphericalHarmonics::InterpPointsTheta = 41 +SphericalHarmonics::InterpPointsPhi = 76 + +SphericalHarmonics::number_of_radii = 5 +SphericalHarmonics::ex_radii[0] = 30 +SphericalHarmonics::ex_radii[1] = 40 +SphericalHarmonics::ex_radii[2] = 50 +SphericalHarmonics::ex_radii[3] = 60 +SphericalHarmonics::ex_radii[4] = 70 + +SphericalHarmonics::number_of_vars = 10 +SphericalHarmonics::vars[0] = "PsiKadelia::psi0re" +SphericalHarmonics::vars[1] = "PsiKadelia::psi0im" +SphericalHarmonics::vars[2] = "PsiKadelia::psi1re" +SphericalHarmonics::vars[3] = "PsiKadelia::psi1im" +SphericalHarmonics::vars[4] = "PsiKadelia::psi2re" +SphericalHarmonics::vars[5] = "PsiKadelia::psi2im" +SphericalHarmonics::vars[6] = "PsiKadelia::psi3re" +SphericalHarmonics::vars[7] = "PsiKadelia::psi3im" +SphericalHarmonics::vars[8] = "PsiKadelia::psi4re" +SphericalHarmonics::vars[9] = "PsiKadelia::psi4im" +SphericalHarmonics::SH_spin_weight[0] = +2 +SphericalHarmonics::SH_spin_weight[1] = +2 +SphericalHarmonics::SH_spin_weight[2] = +1 +SphericalHarmonics::SH_spin_weight[3] = +1 +SphericalHarmonics::SH_spin_weight[4] = 0 +SphericalHarmonics::SH_spin_weight[5] = 0 +SphericalHarmonics::SH_spin_weight[6] = -1 +SphericalHarmonics::SH_spin_weight[7] = -1 +SphericalHarmonics::SH_spin_weight[8] = -2 +SphericalHarmonics::SH_spin_weight[9] = -2 + + + +ActiveThorns = "WaveExtract" + +WaveExtract::out_every = 128 +WaveExtract::maximum_detector_number = 5 +WaveExtract::switch_output_format = 100 +WaveExtract::rsch2_computation = "average Schwarzschild metric" +WaveExtract::l_mode = 8 +WaveExtract::m_mode = 8 +WaveExtract::detector_radius[0] = 30 +WaveExtract::detector_radius[1] = 40 +WaveExtract::detector_radius[2] = 50 +WaveExtract::detector_radius[3] = 60 +WaveExtract::detector_radius[4] = 70 + + + +ActiveThorns = "CarpetIOBasic" + +IOBasic::outInfo_every = 128 +IOBasic::outInfo_reductions = "norm2" +IOBasic::outInfo_vars = " + Carpet::timing + ADMConstraints::ham + IsolatedHorizon::ih_spin[0] + IsolatedHorizon::ih_radius[0] +" + + + +ActiveThorns = "CarpetIOScalar" + +IOScalar::one_file_per_group = yes + +IOScalar::outScalar_every = 128 +IOScalar::outScalar_vars = " + CarpetReduce::weight + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMConstraints::hamiltonian + ADMConstraints::momentum + SphericalSurface::sf_radius + IsolatedHorizon::ih_shapes + IsolatedHorizon::ih_coordinates + IsolatedHorizon::ih_tetrad_l + IsolatedHorizon::ih_tetrad_n + IsolatedHorizon::ih_tetrad_m + IsolatedHorizon::ih_newman_penrose + IsolatedHorizon::ih_weyl_scalars + IsolatedHorizon::ih_twometric + IsolatedHorizon::ih_killing_vector + IsolatedHorizon::ih_killed_twometric + IsolatedHorizon::ih_invariant_coordinates + IsolatedHorizon::ih_multipole_moments + IsolatedHorizon::ih_3determinant + PsiKadelia::WeylComponents + PunctureTracker::pt_loc + SphericalHarmonics::decomposed_vars +" + + + +ActiveThorns = "CarpetIOASCII" + +IOASCII::one_file_per_group = yes + +IOASCII::output_symmetry_points = no +IOASCII::out3D_ghosts = no + +IOASCII::out0D_every = 128 +IOASCII::out0D_vars = " + Carpet::timing + CarpetReduce::weight + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMConstraints::hamiltonian + ADMConstraints::momentum + SphericalSurface::sf_active + SphericalSurface::sf_valid + SphericalSurface::sf_info + SphericalSurface::sf_radius + SphericalSurface::sf_origin + SphericalSurface::sf_coordinate_descriptors + IsolatedHorizon::ih_shapes + IsolatedHorizon::ih_state + IsolatedHorizon::ih_grid_int + IsolatedHorizon::ih_grid_real + IsolatedHorizon::ih_shapes + IsolatedHorizon::ih_coordinates + IsolatedHorizon::ih_tetrad_l + IsolatedHorizon::ih_tetrad_n + IsolatedHorizon::ih_tetrad_m + IsolatedHorizon::ih_newman_penrose + IsolatedHorizon::ih_weyl_scalars + IsolatedHorizon::ih_twometric + IsolatedHorizon::ih_killing_vector + IsolatedHorizon::ih_killed_twometric + IsolatedHorizon::ih_scalars + IsolatedHorizon::ih_invariant_coordinates + IsolatedHorizon::ih_multipole_moments + IsolatedHorizon::ih_3determinant + PsiKadelia::WeylComponents + SphericalHarmonics::decomposed_vars + ProperDistance::pdistance + ProperTime::ptime +" + +IOASCII::out1D_every = 128 +IOASCII::out1D_vars = " + CarpetReduce::weight + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMConstraints::hamiltonian + ADMConstraints::momentum + SphericalSurface::sf_radius + IsolatedHorizon::ih_shapes + IsolatedHorizon::ih_weyl_scalars + IsolatedHorizon::ih_killing_vector + IsolatedHorizon::ih_killed_twometric + IsolatedHorizon::ih_3determinant + PsiKadelia::WeylComponents + SphericalHarmonics::decomposed_vars +" + +IOASCII::out2D_every = 128 +IOASCII::out2D_vars = " + SphericalSurface::sf_radius + SphericalHarmonics::decomposed_vars +" + + + +Activethorns = "CarpetIOHDF5" + +IOHDF5::out_every = 512 +IOHDF5::one_file_per_group = yes +IOHDF5::compression_level = 1 +IOHDF5::out_vars = " + CarpetReduce::weight + ADMBase::metric + ADMBase::curv + ADMBase::lapse + ADMBase::shift + ADMConstraints::hamiltonian + ADMConstraints::momentum + PsiKadelia::WeylComponents +" + +IOHDF5::checkpoint = yes +IO::checkpoint_dir = $parfile +IO::checkpoint_ID = yes +IO::checkpoint_every = 512 +IO::checkpoint_on_terminate = yes + +IO::recover = "autoprobe" +IO::recover_dir = $parfile + + + +ActiveThorns = "Formaline" + +Formaline::collect_metadata = no + +#Formaline::send_as_rdf = yes +#Formaline::rdf_hostname[0] = "numrel07.cct.lsu.edu" +#Formaline::rdf_port [0] = 24000 + + + +ActiveThorns = "TimerReport" + +TimerReport::out_every = 512 +TimerReport::out_filename = "TimerReport" diff --git a/param.ccl b/param.ccl index cfa152d..1475ebb 100644 --- a/param.ccl +++ b/param.ccl @@ -1,61 +1,24 @@ # Parameter definitions for thorn PunctureTracker -# $Header$ -INT pt_verbose "speak up?" +BOOLEAN verbose "speak up?" STEERABLE=always { -0:10 :: "eloquence" -} 1 - -INT pt_num_tracked "Number of punctures that should be tracked" -{ - 0:10 :: "" -} 0 - -# not to be used -# because the puncture has to be tracked every timestep -#INT pt_every "how often to calculate new position" -#{ -# 0:*::"" -#} 0 +} "no" -BOOLEAN pt_track_punctures "do it?" +BOOLEAN track[10] "Track this puncture" { } "no" -REAL pt_initial_x[10] "Initial x coordinate positions of BHs" -#This overwrites parameter par_b from Twopunctures if nonzero +REAL initial_x[10] "Initial x coordinate positions of punctures" { - : :: "Anything goes" + *:* :: "" } 0.0 -REAL pt_initial_y[10] "Initial y coordinate positions of BHs" -#This overwrites parameter par_b from Twopunctures if nonzero +REAL initial_y[10] "Initial y coordinate positions of punctures" { - : :: "Anything goes" + *:* :: "" } 0.0 -REAL pt_initial_z[10] "Initial z coordinate positions of BHs" -#This overwrites parameter par_b from Twopunctures if nonzero +REAL initial_z[10] "Initial z coordinate positions of punctures" { - : :: "Anything goes" + *:* :: "" } 0.0 - -BOOLEAN pt_regrid "should we use PT to Regrid?" -{ -} "no" - -REAL pt_regrid_buffer "regrid if new location-old_location > buffer" -{ - : ::"anything greater zero" -} 0.5 - -INT pt_which_surface_to_take[10] "which surface" -{ -0:42 :: "" -} 0 - - -SHARES: TwoPunctures - -USES REAL par_b - diff --git a/schedule.ccl b/schedule.ccl index ab3b151..cce26af 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -1,27 +1,21 @@ # Schedule definitions for thorn PunctureTracker -# $Header$ +STORAGE: pt_loc pt_loc_p -storage: pt_loc[2] -storage: pt_shift +SCHEDULE PunctureTracker_Init AT initial +{ + LANG: C + OPTIONS: GLOBAL +} "Calculate initial location of punctures" -if(pt_track_punctures){ +SCHEDULE PunctureTracker_Track AT evol AFTER MoL_Evolution +{ + LANG: C + OPTIONS: GLOBAL +} "Calculate new location of punctures" - SCHEDULE PunctureTracker_init at CCTK_INITIAL - { - LANG: C - OPTIONS: GLOBAL - } "calculate new location of punctures" - - SCHEDULE PunctureTracker at CCTK_POSTSTEP before SphericalSurface_HasBeenSet - { - LANG: C - OPTIONS: GLOBAL - } "calculate new location of punctures" - - SCHEDULE update_punc_loc at CCTK_POSTSTEP after PunctureTracker before SphericalSurface_HasBeenSet - { - LANG: C - OPTIONS: GLOBAL - } "update new location of punctures" -} +SCHEDULe PunctureTracker_SetPositions AT preregrid +{ + LANG: C + OPTIONS: global +} "Set positions of refined regions" diff --git a/src/adjust.c b/src/adjust.c deleted file mode 100644 index 43a2df0..0000000 --- a/src/adjust.c +++ /dev/null @@ -1,90 +0,0 @@ -/* $Header$ */ - -#include -#include - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" - - - -void -BHTracker_Adjust (CCTK_ARGUMENTS); - -static void -setpar (char const * name, int index, char const * thorn, CCTK_REAL value); - - - -void -BHTracker_Adjust (CCTK_ARGUMENTS) -{ - DECLARE_CCTK_PARAMETERS; - DECLARE_CCTK_ARGUMENTS; - - int n; - int si, ri; - char name[100], val[100]; - int icnt, ierr; - - CCTK_REAL xpos,ypos,zpos; - - if ((cctk_iteration-1) < track_first) return; - if ((cctk_iteration-1) % track_every != 0) return; - - - - CCTK_INFO ("Tracking..."); - - for (n=0; n=0 && si=0 && ri 0) { - - CCTK_VInfo (CCTK_THORNSTRING, - " Object #%d has moved to [%g,%g,%g]", n, - (double)sf_centroid_x[si], (double)sf_centroid_y[si], - (double)sf_centroid_z[si]); - - /* calculate new positions */ - xpos = (sf_centroid_x[si] - initial_x[si]); - ypos = (sf_centroid_y[si] - initial_y[si]); - zpos = (sf_centroid_z[si] - initial_z[si]); - - setpar ("offsetx", ri, "CarpetRegrid", xpos); - setpar ("offsety", ri, "CarpetRegrid", ypos); - setpar ("offsetz", ri, "CarpetRegrid", zpos); - - } else { - - CCTK_VInfo (CCTK_THORNSTRING, - " Object #%d cannot be detected at this time", n); - - } - } -} - - - -static void -setpar (char const * const name, int const index, char const * const thorn, - CCTK_REAL const value) -{ - char nambuf[100], valbuf[100]; - int icnt, ierr; - - icnt = snprintf (nambuf, sizeof nambuf, "%s[%d]", name, index); - assert (icnt>=0 && icnt=0 && icnt -#include - -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" -#include "util_Table.h" - -#define REFLEVEL round(log10((CCTK_REAL)(cctkGH->cctk_levfac[0]))/log10(2.0)) - -void PunctureTracker (CCTK_ARGUMENTS); -void PunctureTracker_init (CCTK_ARGUMENTS); - - -void PunctureTracker_init (CCTK_ARGUMENTS) -{ - DECLARE_CCTK_PARAMETERS; - DECLARE_CCTK_ARGUMENTS; - int i; - - if(pt_verbose>0) - CCTK_INFO("Initializing PunctureTracker"); - - if((pt_initial_x[0]!=0.0)||(pt_initial_y[0]!=0.0)||(pt_initial_z[0]!=0.0)) - for (i=0;i0) - CCTK_INFO("done initializing PunctureTracker"); -} - -void PunctureTracker (CCTK_ARGUMENTS) -{ - DECLARE_CCTK_PARAMETERS; - DECLARE_CCTK_ARGUMENTS; - - int n,ierror,i,numpoints; - // portland doesnt like that - //CCTK_REAL pt_betax[pt_num_tracked], pt_betay[pt_num_tracked], pt_betaz[pt_num_tracked]; - //so I do that - - CCTK_REAL * restrict const pt_betax = malloc(pt_num_tracked*sizeof(CCTK_REAL)); - CCTK_REAL * restrict const pt_betay = malloc(pt_num_tracked*sizeof(CCTK_REAL)); - CCTK_REAL * restrict const pt_betaz = malloc(pt_num_tracked*sizeof(CCTK_REAL)); - - int operator_handle, param_table_handle, coordsys_handle; - int input_array_indices[3]; - const void *interp_coords[3]; - void *output_arrays[3]; - - const CCTK_INT input_array_type_codes[3] = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL}; - const CCTK_INT output_array_type_codes[3]= { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL}; - - if(pt_verbose >0) - CCTK_INFO("tracking Punctures"); - - interp_coords[0] = (const void *) pt_loc_x_p; - interp_coords[1] = (const void *) pt_loc_y_p; - interp_coords[2] = (const void *) pt_loc_z_p; - - operator_handle = -1; - coordsys_handle = -1; - param_table_handle = -1; - - // set interpolation handles - operator_handle = CCTK_InterpHandle("Lagrange polynomial interpolation"); - if (operator_handle < 0) CCTK_WARN(0, "can’t get interpolation handle!"); - - param_table_handle = Util_TableCreateFromString("order=4"); - if (param_table_handle < 0) CCTK_WARN(0, "can't ge parameter table!"); - - coordsys_handle = CCTK_CoordSystemHandle("cart3d"); - if (coordsys_handle < 0) CCTK_WARN(0,"can't get coordsys handle!"); - if(pt_verbose>0) - for (i=0; i1) - printf("at the location of puncture %i:\n is pt_betax=%f, pt_betay=%f,pt_betaz=%f \n",n,pt_betax[n],pt_betay[n],pt_betaz[n]); - pt_shiftx[n] = pt_betax[n]; - pt_shifty[n] = pt_betay[n]; - pt_shiftz[n] = pt_betaz[n]; - } - - free(pt_betax,pt_betay,pt_betaz); - if(pt_verbose>0) - CCTK_INFO ("done"); -} - -void update_punc_loc(CCTK_ARGUMENTS){ - - DECLARE_CCTK_PARAMETERS; - DECLARE_CCTK_ARGUMENTS; - - CCTK_REAL pt_dt; - int n; - - pt_dt = CCTK_DELTA_TIME; - - for (n=0; n< pt_num_tracked; n++) { - //correct small errors - // if (pt_shiftx[n] < 10e-14) pt_shiftx[n] = 0.0; - // if (pt_shifty[n] < 10e-14) pt_shifty[n] = 0.0; - // if (pt_shiftz[n] < 10e-14) pt_shiftz[n] = 0.0; - - if(pt_verbose > 2){ - printf("updating puncture %i\n",n); - printf(" pt_dt is %f\n",pt_dt); - printf(" pt_loc_x is %f\n",pt_loc_x[n]); - printf(" pt_loc_x_p is %f\n",pt_loc_x_p[n]); - printf(" pt_shiftx is %f\n",pt_shiftx[n]); - printf(" pt_dt is %f\n",pt_dt); - fflush(stdout); - } - - printf("WRITING TO SURFACE %d\n",pt_which_surface_to_take[n]); - pt_loc_x[n] = pt_loc_x_p[n] - pt_shiftx[n] * pt_dt; - pt_loc_y[n] = pt_loc_y_p[n] - pt_shifty[n] * pt_dt; - pt_loc_z[n] = pt_loc_z_p[n] - pt_shiftz[n] * pt_dt; - - sf_centroid_x[pt_which_surface_to_take[n]]=pt_loc_x_p[n]; - sf_centroid_y[pt_which_surface_to_take[n]]=pt_loc_y_p[n]; - sf_centroid_z[pt_which_surface_to_take[n]]=pt_loc_z_p[n]; - - sf_valid[pt_which_surface_to_take[n]]=1.0; - sf_valid[pt_which_surface_to_take[n]]=1.0; - sf_valid[pt_which_surface_to_take[n]]=1.0; - } -} diff --git a/src/puncture_tracker.c b/src/puncture_tracker.c new file mode 100644 index 0000000..d3da795 --- /dev/null +++ b/src/puncture_tracker.c @@ -0,0 +1,285 @@ +#include +#include + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_Table.h" + + + +static int const max_num_tracked = 10; + + + +void +PunctureTracker_Init (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (verbose) { + CCTK_INFO ("Initializing PunctureTracker"); + } + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + pt_loc_t[n] = cctk_time; + pt_loc_x[n] = initial_x[n]; + pt_loc_y[n] = initial_y[n]; + pt_loc_z[n] = initial_z[n]; + } else { + // Initialise to some sensible but unimportant values + pt_loc_t[n] = 0.0; + pt_loc_x[n] = 0.0; + pt_loc_y[n] = 0.0; + pt_loc_z[n] = 0.0; + } + pt_loc_t_p[n] = 0.0; + pt_loc_x_p[n] = 0.0; + pt_loc_y_p[n] = 0.0; + pt_loc_z_p[n] = 0.0; + } +} + + + +void +PunctureTracker_Track (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + // Do not track while setting up initial data; + // time interpolation may fail + + if (cctk_iteration == 0) { + return; + } + + // Some output + + if (verbose) { + CCTK_INFO ("Tracking punctures..."); + } + + if (verbose) { + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + CCTK_VInfo (CCTK_THORNSTRING, + "Puncture #%d is at (%g,%g,%g)", + n, + (double)pt_loc_x[n], + (double)pt_loc_y[n], + (double)pt_loc_z[n]); + } + } + } + + // Manual time level cycling + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + pt_loc_t_p[n] = pt_loc_t[n]; + pt_loc_x_p[n] = pt_loc_x[n]; + pt_loc_y_p[n] = pt_loc_y[n]; + pt_loc_z_p[n] = pt_loc_z[n]; + + pt_loc_t[n] = cctk_time; + } + } + + // Interpolate + + // Dimensions + int const dim = 3; + + // Interpolation operator + int const operator_handle = + CCTK_InterpHandle ("Lagrange polynomial interpolation"); + if (operator_handle < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can’t get interpolation handle"); + goto label_return; + } + + // Interpolation parameter table + int const param_table_handle = Util_TableCreateFromString ("order=4"); + if (param_table_handle < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can't create parameter table"); + goto label_return; + } + + { + + // Interpolation coordinate system + int const coordsys_handle = CCTK_CoordSystemHandle ("cart3d"); + if (coordsys_handle < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can't get coordinate system handle"); + goto label_free_param_table; + } + + // Only processor 0 interpolates + int const num_points = CCTK_MyProc(cctkGH) == 0 ? max_num_tracked : 0; + + // Interpolation coordinates + assert (dim == 3); + CCTK_POINTER_TO_CONST interp_coords[3]; + interp_coords[0] = pt_loc_x_p; + interp_coords[1] = pt_loc_y_p; + interp_coords[2] = pt_loc_z_p; + + // Number of interpolation variables + int const num_vars = 3; + + // Interpolated variables + assert (num_vars == 3); + int input_array_indices[3]; + input_array_indices[0] = CCTK_VarIndex ("ADMBase::betax"); + input_array_indices[1] = CCTK_VarIndex ("ADMBase::betay"); + input_array_indices[2] = CCTK_VarIndex ("ADMBase::betaz"); + + // Interpolation result types + assert (num_vars == 3); + CCTK_INT output_array_type_codes[3]; + output_array_type_codes[0] = CCTK_VARIABLE_REAL; + output_array_type_codes[1] = CCTK_VARIABLE_REAL; + output_array_type_codes[2] = CCTK_VARIABLE_REAL; + + // Interpolation result + CCTK_REAL pt_betax[max_num_tracked]; + CCTK_REAL pt_betay[max_num_tracked]; + CCTK_REAL pt_betaz[max_num_tracked]; + + assert (num_vars == 3); + CCTK_POINTER output_arrays[3]; + output_arrays[0] = pt_betax; + output_arrays[1] = pt_betay; + output_arrays[2] = pt_betaz; + + // Interpolate + int const ierr = CCTK_InterpGridArrays + (cctkGH, dim, + operator_handle, param_table_handle, coordsys_handle, + num_points, + CCTK_VARIABLE_REAL, + interp_coords, + num_vars, input_array_indices, + num_vars, output_array_type_codes, output_arrays); + if (ierr < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Interpolation error"); + goto label_free_param_table; + } + + if (CCTK_MyProc(cctkGH) == 0) { + + // Some more output + + if (verbose && CCTK_MyProc(cctkGH) == 0) { + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + CCTK_VInfo (CCTK_THORNSTRING, + "Shift at puncture #%d is at (%g,%g,%g)", + n, + (double)pt_betax[n], + (double)pt_betay[n], + (double)pt_betaz[n]); + } + } + } + + // Time evolution + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + CCTK_REAL const dt = pt_loc_t[n] - pt_loc_t_p[n]; + // First order time integrator + // Michael Koppitz says this works... + // if it doesn't, we can make it second order accurate + pt_loc_x[n] = pt_loc_x_p[n] + dt * (- pt_betax[n]); + pt_loc_y[n] = pt_loc_y_p[n] + dt * (- pt_betay[n]); + pt_loc_z[n] = pt_loc_z_p[n] + dt * (- pt_betaz[n]); + } + } + + } + + // Broadcast result + + CCTK_REAL loc_local[3*max_num_tracked]; + if (CCTK_MyProc(cctkGH) == 0) { + for (int n = 0; n < max_num_tracked; ++ n) { + loc_local[ n] = pt_loc_x[n]; + loc_local[ max_num_tracked+n] = pt_loc_y[n]; + loc_local[2*max_num_tracked+n] = pt_loc_z[n]; + } + } else { + for (int n = 0; n < max_num_tracked; ++ n) { + loc_local[ n] = 0.0; + loc_local[ max_num_tracked+n] = 0.0; + loc_local[2*max_num_tracked+n] = 0.0; + } + } + + CCTK_REAL loc_global[3*max_num_tracked]; + + int const handle_sum = CCTK_ReductionHandle ("sum"); + if (handle_sum < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Can't get redunction handle"); + goto label_free_param_table; + } + + int const ierr2 = CCTK_ReduceLocArrayToArray1D + (cctkGH, -1, handle_sum, + loc_local, loc_global, 3*max_num_tracked, CCTK_VARIABLE_REAL); + if (ierr2 < 0) { + CCTK_WARN (CCTK_WARN_ALERT, "Reduction error"); + goto label_free_param_table; + } + + for (int n = 0; n < max_num_tracked; ++ n) { + pt_loc_x[n] = loc_global[ n]; + pt_loc_y[n] = loc_global[ max_num_tracked+n]; + pt_loc_z[n] = loc_global[2*max_num_tracked+n]; + } + + } + + // Done + + // Poor man's exception handling + label_free_param_table: + Util_TableDestroy (param_table_handle); + + label_return: + return; +} + + + +void +PunctureTracker_SetPositions (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + for (int n = 0; n < max_num_tracked; ++ n) { + if (track[n]) { + + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Setting position of refined region from puncture #%d to (%g,%g,%g)", + n, + (double)pt_loc_x[n], + (double)pt_loc_y[n], + (double)pt_loc_z[n]); + } + + // Set position in CarpetRegrid2 + position_x[n] = pt_loc_x[n]; + position_y[n] = pt_loc_y[n]; + position_z[n] = pt_loc_z[n]; + + } + } +} -- cgit v1.2.3