Updating your branches

We make regular updates to the code that add new features/enhancements or fix bugs. It is good practice to keep your fork of this repository up to date in order to take advantage of these new features/enhancements and bugfixes. We recommend updating your code at the start of new research projects unless you are sure of what you are doing.

Most of the changes we make are backwards compatible meaning that if you merge the main branch of this repository into the branches of your fork, there shouldn't be a need to fix merge conflicts so long as you haven't modified any files in the Source subdirectory and there shouldn't be any significant changes to function prototypes that break compilation of your own examples.

Below we will try to list new features/enhancements and bugfixes, how to use them/take advantage of them and link to the relevant pull requests which will have more information in order from newest to oldest. We will also try to describe any more significant changes that require modifying existing examples in order to work.

Update GitHub actions #221

• The main user-facing change in this PR is that the version of clang-format used to enforce the style of the code has been updated from v10 to v15 which formats very slightly differently.
• You should upgrade your local version of clang-format to the new version. Instructions on how to install clang-format can be found in the wiki.

Diagnostics variable tagging #219

• This just adds the ability to use diagnostics variables in the tagging criterion.
• There is a new virtual function with the signature

  virtual void
  GRAMRLevel::computeTaggingCriterion(FArrayBox &tagging_criterion,
                                      const FArrayBox &current_state,
                                      const FArrayBox &current_state_diagnostics);

  which can you overload in your examples' PhysicsLevel class.
• You can continue using the old overload so there is no need to change your examples if you don't want to use diagnostic variable tagging.

Clang Format pre-commit hook #198

• There are now instructions on the wiki on how to set-up a pre-commit hook that will check your code is formatted correctly according to our Clang Format configuration before allowing you to commit changes.
• In order to use it, you will need to be up-to-date with this PR.
• There is no action required but we recommend setting this up to keep your code tidy.

Add Matter and Z4 terms to Weyl4 #174

Action strongly recommended for users of the Weyl4 class. See below.

• The main change in this PR is the addition of matter and Z4 terms to the Weyl scalar. The calculation of the electric and magnetic parts of the Weyl scalar has also been modified to ensure they are explicitly trace-free and symmetric respectively rather than relying on satisfaction of the constraints in order to do so.
• If you evolving with CCZ4 and wish to include the Z4 terms to the Weyl scalar, simply pass the formulation parameter as the final argument to the constructor of Weyl4 e.g. instead of doing (in SpecificExampleLevel.cpp)

  BoxLoops::loop(Weyl4(m_p.extraction_params.center, m_dx),
                 m_state_new, m_state_diagnostics,
                 EXCLUDE_GHOST_CELLS);

  you would do

  BoxLoops::loop(
      Weyl4(m_p.extraction_params.center, m_dx, m_p.formulation),
      m_state_new, m_state_diagnostics, EXCLUDE_GHOST_CELLS);

  Note that this formulation argument defaults to CCZ4RHS<>::USE_BSSN which does not include the Z4 terms which is similar to existing behaviour.
• If you are evolving with matter and wish to include the matter terms when calculating the Weyl scalar, you should switch to using the new MatterWeyl4 class defined in the MatterWeyl4.hpp header with constructor

  MatterWeyl4(matter_t a_matter,
              const std::array<double, CH_SPACEDIM> a_center,
              const double a_dx,
              const int a_formulation = CCZ4RHS<>::USE_CCZ4,
              double a_G_Newton = 1.0);

  where matter_t is a template parameter as for the MatterCCZ4RHS class.
• There are also some minor changes to the calculation of the Ricci tensor though the changes should be analytically identical to the old calculation. For details, please read the comments in the PR.
• A new EMTensor class has been added to calculate and store components of the EM tensor. This is defined in the EMTensor.hpp header file

Improve folder structure of hdf5 / pouts / .dat #166

Action required for PunctureTracker. See below.

• This improves organization of directories for all outputs in a given run.
• You can specify a global directory where all files are located using output_path. This will be prepended to all internal paths. Defaults to "".
• hdf5 files, pout files and .dat files (puncture files, Weyl4 integrated modes, etc.) can now be put in their own specific directory using hdf5_subpath, pout_subpath and data_subpath. These are created at runtime if the path doesn't exist. They default to "". Additionally extraction_subpath can be used to specify a directory for extraction dumped data when write_extraction = 1.
• Change the name prefix of pout files using pout_prefix.
• pout verbosity can be reduced with print_progress_only_to_rank_0, which will print the number of boxes per rank on ranks other than rank 0 only at regrid time steps. For large number of nodes this improves performance.
• Example: Adding to the BinaryBH Example params_very_cheap.txt the following:

  output_path = "results"
  hdf5_subpath = "hdf5"
  pout_subpath = "pout"
  data_subpath = "data"
  pout_prefix = "poutChanged"
  

  would put all the files in a folder results/ with files in subdirectories: hdf5's in results/hdf5, pout's in results/pout and .dats files in results/data. It would also change the name of the pout files to poutChanged.
• Action required: The signature of PunctureTracker::initial_setup was changed to force replacing checkpoint_prefix by data_path. In your Example's Main:

  bh_amr.m_puncture_tracker.initial_setup(
              {sim_params.bh1_params.center, sim_params.bh2_params.center},
              sim_params.checkpoint_prefix, puncture_tracker_min_level);

  Should be changed to

  bh_amr.m_puncture_tracker.initial_setup(
              {sim_params.bh1_params.center, sim_params.bh2_params.center},
              "punctures", sim_params.data_path, puncture_tracker_min_level);

  (notice how the file name is customizable, above set to punctures.dat)

Add TwoPunctures integration #160

• This allows the use of the TwoPunctures initial data code in the BinaryBH example.
• The TwoPunctures code is kept in a separate repository for licensing reasons. Instructions are provided in the README explaining how to build the BinaryBH example with TwoPunctures.
• The TwoPunctures code creates high quality Bowen-York puncture data for two black holes using a spectral method. The constraint violations are much smaller than the existing approximate initial data that is normally used in this example. Additionally spinning BHs are now possible.
• All the new code is conditionally compiled on the USE_TWOPUNCTURES macro so there should be no changes if you do not wish to use it.

Add sixth order spatial derivatives #162

• This adds a new SixthOrderDerivatives class which is similar to the existing FourthOrderDerivatives class except that it uses sixth order stencils (but still sixth order Kreiss-Oliger dissipation).
• Note that in order to evaluate sixth order stencils, you need 4 ghost cells (as opposed to 3 for fourth order stencils).
• There is a new max_spatial_derivative_order parameter in ChomboParameters which defaults to 4 but can be set to 6. If set to 6, it will change the default num_ghosts parameter to 4. This parameter can be used in examples to decide which stencils to use. For example:

  if (m_p.max_spatial_derivative_order == 4)
  {
      BoxLoops::loop(CCZ4RHS<MovingPunctureGauge, FourthOrderDerivatives>(
                         m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
                     a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
  }
  else if (m_p.max_spatial_derivative_order == 6)
  {
      BoxLoops::loop(CCZ4RHS<MovingPunctureGauge, SixthOrderDerivatives>(
                         m_p.ccz4_params, m_dx, m_p.sigma, m_p.formulation),
                     a_soln, a_rhs, EXCLUDE_GHOST_CELLS);
  }

• Note that this parameter does not automatically change all existing derivatives from fourth to sixth order and it is up to individual users to add this to their examples as above.

Separate Out Gauge #154

• This refactors the CCZ4 and MatterCCZ4 classes by templating them over a gauge and derivative class.
• The new template classes are called CCZ4RHS and MatterCCZ4RHS respectively.
• The computation of the gauge RHS is now done in the gauge class. The old gauge which is still the default can be found in MovingPunctureGauge.hpp.
• A new gauge has also been added which is similar to the existing moving puncture gauge but has a first-order Gamma-driver shift condition. It can be found in IntegratedMovingPunctureGauge.hpp. It uses the same parameters.
• An alias in the form

  using CCZ4 = CCZ4RHS<>;

  (and similarly for MatterCCZ4) is provided as well as a symbolic link such that the old CCZ4.hpp points to CCZ4RHS.hpp so existing examples will continue to compile without changes. However this is marked deprecated so the compiler will issue a warning. You should update your examples to use the new CCZ4RHS or MatterCCZ4RHS template classes as has been done for all of the provided examples.

Add more flexibility when filling ghosts #152

• This allows the user to be more specific when filling/exchanging ghosts by choosing which variables to apply it to.
• The reason why you might want to do this is because filling ghosts is necessary but expensive so it should be minimised.
• There should be no need to modify existing example code.
• If you wish to take advantage of the speedup, you can now choose to manually fill ghosts rather than it being done automatically in AMRInterpolator::refresh(). For example, previously in BinaryBHLevel::specifiPostTimeStep, we called

  m_gr_amr.m_interpolator->refresh();

  before doing WeylExtraction. This exchanges all ghosts for all variables (both evolution and diagnostic) on all levels. This was changed to

  bool fill_ghosts = false;
  m_gr_amr.m_interpolator->refresh(fill_ghosts);
  m_gr_amr.fill_multilevel_ghosts(
      VariableType::diagnostic, Interval(c_Weyl4_Re, c_Weyl4_Im),
      min_level);

  Now, the ghosts are only exchanged for Weyl4_Re and Weyl4_Im, from the minimun extraction level to max_level, since this is all that is necessary for WeylExtraction
• Similarly, you can manually fill ghosts before tagging cells by overriding the GRAMRLevel::preTagCells function (see any of the provided examples e.g. BinaryBHLevel.cpp. By default, this function will fill all evolution variable ghosts.

Add class to call specificPostTimeStep at t=0 #142

• This allows the user to call specificPostTimeStep at t=0. In general, it allows the user to call any GRAMRLevel function at regular intervals.
• To call specificPostTimeStep at t=0, you should add to your Main[Example].cpp:

    auto task = [](GRAMRLevel *level) { if (level->time() == 0.) level->specificPostTimeStep(); };
    MultiLevelTaskPtr<> call_task(task);
    call_task.execute(bh_amr);

  Note this will also be called at restarts, so the if (level->time() == 0.) is important.
• If this is done, inside specificPostTimeStep, recall to set the argument bool first_step = (m_time == 0.); that is typically passed to Extraction classes or the SmallDataIO class (opposed to the usual bool first_step = (m_time == m_dt);).
• Alternatively, to call a personalised function GRAMRLevel::my_function every 2 timesteps:

    auto task = [](GRAMRLevel *level) { level->my_function(); };
    int some_interval = 2;
    bool reverse_levels = true;
    MultiLevelTaskPtr<> call_task(task, reverse_levels, some_interval);
    bh_amr.schedule(call_task);

  Note that this function is called post-PlotLevel.
• The default behaviour is the use of reverse_levels = true, as specificPostTimeStep is called after each level concludes its timestep in reverse order (e.g. for GW extraction, the Weyl4 scalar is filled starting from max_level and extraction is performed when min_extraction_level is reached).

Add grid_buffer_size parameter #153

• This adds a new grid_buffer_size parameter which specifies the number of cells on level l between the boundary of level l+1 and the boundary of level l. Before, this was hardcoded to ceil(num_ghosts/ref_ratio) + 3 = 2 + 3 = 5. Now, the default values is 8.
• Increasing this parameter will help keep refinement boundaries further apart which will reduce reflections of noise.
• There should be no need to modify existing example code.

Speed up SurfaceExtraction #147

• This speeds up the SurfaceExtraction class (and derived classes such as SphericalExtraction and WeylExtraction) by making only rank 0 request points for interpolation. This should be particularly helpful in the case of a large number of MPI ranks/processes.
• In most cases there should be no need to modify existing example code. However in the unlikely circumstance you used the values of calculated surface integrals other than writing them to a file e.g. by doing something like

  SphericalExtraction spherical_extraction(...);
  spherical_extraction.extract(m_gr_amr.interpolator);
  std::vector<double> integral;
  spherical_extraction.add_integrand(integrand, integral);
  spherical_extraction.integrate();
  <code that uses the value of integral>

  they will no longer be calculated on ranks > 0 by default. However, the previous behaviour can be implemented by passing a_broadcast_integral = true to the very final argument of add_integrand (and similarly to add_var_integrand, etc.). See SphericalExtractionTest.cpp for an example.

Allow building GRChombo with NAMESPACE=TRUE #151

• This adds #include "UsingNamespace.H" to source files that directly include Chombo header files (after all other header files).
• There should be no need to modify existing example code.
• Any PRs with new source files to this repository should include UsingNamespace.H if other Chombo header files are directly included.
• Building Chombo and GRChombo with NAMESPACE=TRUE places most Chombo functions/classes in the Chombo namespace and the file UsingNamespace.H adds the line

  using namespace Chombo;

  in this case. This is currently unnecessary but may be used in future to link with libraries that have definitions of functions/classes that conflict with the Chombo ones (e.g. In #91 as both the SLURM (a commonly used job scheduler) library and Chombo have structs/classes called List)

Parameter Checking #149

• This adds the checking of parameters after they are loaded at the start of simulations in order to prevent a crash later on due to invalid parameters. It also provides a readable message that can make it easier to figure out which parameter is invalid and why.
• Additionally, the provided examples can now be run to just check parameters without starting the whole simulation (e.g. on the login nodes of a cluster before job submission) by passing the -check_params flag on the command line. e.g.

  mpiexec /path/to/exec.ex /path/to/params.txt -check_params
  

  To add this feature to your own examples, simply add the following two lines after the line SimulationParameters sim_params(pp); in Main_example.cpp:runGRChombo(int argc, char *argv[]):

      if (sim_params.just_check_params)
        return 0;

  e.g. in the BinaryBH example: Main_BinaryBH.cpp.
• You can also add the checking of your own examples' parameters. See the BinaryBH/SimulationParameters.hpp file for examples. The check_parameter() and warn_parameter() functions are defined in ChomboParameters.hpp.

Update for extrapolating BC and update SF example to use them #140

• This updated added a new class of boundary condition that extrapolated the values of variables radially from within the grid to the outer ghost cells. This is particularly useful where one wants a "floating" boundary condition, e.g. in the case where there is an oscillating (in time, but spatially homogeneous) field.

• No updates to code should be required if one does not want to use these updates.

• One can activate the condition by choosing option 3 in the params.txt file for the boundary in each direction, and specifying the extrapolation order. So far only order 0 or 1 extrapolation is supported. e.g.

  # if not periodic, then specify the boundary type
  # 0 = static, 1 = sommerfeld, 2 = reflective
  # 3 = extrapolating, 4 = mixed
  # (see BoundaryConditions.hpp for details)
  hi_boundary = 3 3 3
  lo_boundary = 2 2 2
  extrapolation_order = 1
  

• there is also the option to use extrapolating boundaries for a subset of vars and Sommerfeld for the rest using option 4. See the ScalarField Example params.txt file for a model of this.

• The way in which Sommerfeld boundary params is read in was streamlined so one can input only non zero values, e.g.

  num_nonzero_asymptotic_vars = 5
  nonzero_asymptotic_vars = chi h11 h22 h33 lapse
  nonzero_asymptotic_values = 1.0 1.0 1.0 1.0 1.0
  

  This should be preferred over the previous method of inputting the full UserVariables vector, although that remains supported for now.

Update for Diagnostic Variables #119

Note that this update on 2 July 2020 exceptionally broke old code - please see this page for details.