Skip to content

Conversation

@JCAurre
Copy link
Member

@JCAurre JCAurre commented Feb 24, 2025

This PR updates the BinaryBH example, whose accuracy has always been highly sensitive to the adaptive mesh refinement. The primary cause of the inaccurate phase evolution of the binary was an artificial growth of the black hole masses, which only converges to zero at very high resolutions. Here we have identified that this mass growth is very sensitive to the Kreiss-Oliger (KO). For example, the growth rate is worse for larger sigma coefficients and the waveforms are more accurate for small values of sigma. However, choosing small sigma does not remove high frequency errors near the extraction spheres, which result in noisy waveforms. Given the dependence of the KO on the timestep Δt, we now introduce a space-dependent coefficient that turns off the KO near the punctures. For now, we have found that lapse^6 * sigma yields optimal results when starting with a pre-collapsed lapse. However, one might want to update this in the future to lapse^4 * chi * sigma if the lapse starts from 1. This PR also implements a tagging criterion closer to the standard 'Moving Boxes' approach, and introduces default params.txt files at two resolutions: N=48, which runs at 100 M/hr on 800 CPUs, and N=64, which runs at 50 M/hr on 1500 CPUs. Hence, these setups generate relatively accurate waveforms in 1 and 2 days, with errors of approximately 1% and 0.5%, respectively. These updates seem to allow us to achieve similar accuracies to the waveforms in arXiv:2112.10567, but with lower resolution (and thus faster) simulations.

image

@JCAurre JCAurre requested a review from KAClough February 24, 2025 20:58
@KAClough KAClough self-assigned this Feb 24, 2025
Copy link
Member

@KAClough KAClough left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This all looks great, just some minor stylistic changes.

Whilst we are here and making this example efficient, can we please define constraint_calculation_level=0 and wrap the constraint calculation into a at_level_timestep_multiple(constraint_calculation_level) as we do for the Weyl data since this is otherwise calculating at every post timestep, and it only outputs every coarse timestep so I can't see why that is needed.

Also it would be good to update the usual params.txt file for the better params - especially to rescale the sigma value.

Another unrelated point that we could fix is that the KerrBH example has an error in the params.txt - the lower boundary should be 2 2 2 not 0 0 0, and the max level should probably be 6 not 3.

We will also need to update the wiki for these changes, and having the speeds there for each set of params will be very useful.

#include "SixthOrderDerivatives.hpp"
#include "SmallDataIO.hpp"
#include "TraceARemoval.hpp"
#include "TwoPuncturesBoxExtractionTaggingCriterion.hpp"
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we rename this MovingBoxesAndExtractionTaggingCriterion.hpp as it is not a TwoPunctures specific method.

double a_dx, //!< The grid spacing
double a_sigma, //!< Kreiss-Oliger dissipation coefficient
int a_formulation = USE_CCZ4, //!< Switches between CCZ4, BSSN,...
int a_rescale_sigma = 0, //!< Allows a space dependent KO coefficient
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be a bool, rather than an int, unless we plan to provide multiple settings?

if (m_rescale_sigma == 1)
{
// rescale KO coefficient with lapse so that it is zero near puncture
sigma = this->m_sigma * pow(vars.lapse, 6.0);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you need the "this->"?

public:
double sigma; // Kreiss-Oliger dissipation parameter

int rescale_sigma;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

bool here too

// of it, which means inner \pm L/4
// we want each level to be double the innermost one in size
const double factor = pow(2.0, m_max_level - m_level - 1);
// const int m_horizon_max_levels = 9;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove old commented code here

#include "SphericalExtraction.hpp"
#include "Tensor.hpp"

class TwoPuncturesBoxExtractionTaggingCriterion
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rename as marked above!

auto regrid = simd_compare_lt(
max_abs_xyz, 2.5 * factor * m_puncture_masses[ipuncture]);
criterion = simd_conditional(regrid, 100.0, criterion);
// if (m_level < m_max_level){
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove commented code

}
else
{
// where am i?
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove irrelevant comment

// }
}
}
else
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we add a comment here about why this is needed - ie for when they are merging/merged we need a larger region.

m_puncture_coords[0][i] - m_puncture_coords[1][i];
puncture_separation_squared += displacement * displacement;
}
// make sure the inner part is regridded around the horizon
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update or remove comment

#TP_swap_xz = false
TP_momentum_plus = -0.000510846 0.0841746 0.0
TP_momentum_minus = 0.000510846 -0.0841746 0.0
TP_spin_plus = 0.0 0.0 0.0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JCAurre as long as you leave this PR open I will keep adding things! Ha ha ha.

It would be a great service to future people if you added a comment here as follows:
# Note that these are the spin values J = chi * M^2 where chi is the dimensionless spin # So e.g. for a dimensionless spin of 0.6, and a target mass of 0.5, the value should be 0.15

@KAClough KAClough added the enhancement Modification of existing feature/general improvement label Nov 19, 2025
@AreefW AreefW force-pushed the bbh_example branch 2 times, most recently from da7b9e5 to 85b25ad Compare November 20, 2025 15:25
@KAClough
Copy link
Member

We should also remove the TwoPunctures README.md file in the BBH example, as this is now duplicated in the wiki here: https://github.com/GRTLCollaboration/GRChombo/wiki/TwoPunctures-Initial-Data.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement Modification of existing feature/general improvement

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants