Add BnsInitialData System #6054

isaaclegred · 2024-06-03T18:51:08Z

Proposed changes

This PR adds the system for generating BNS Irrotational Initial data. To do this, an elliptic system is introduced to solve for the velocity potential using the framework of Baumgarte + Shapiro Ch. 15. The background consists of a matter profile which is not yet implemented (as it is contained in the initial data).

Upgrade instructions

Code review checklist

The code is documented and the documentation renders correctly. Run
make doc to generate the documentation locally into BUILD_DIR/docs/html.
Then open index.html.
The code follows the stylistic and code quality guidelines listed in the
code review guide.
The PR lists upgrade instructions and is labeled bugfix or
new feature if appropriate.

Further comments

isaaclegred · 2024-06-03T18:51:28Z

@nilsvu whenever you get a chance!

src/Elliptic/Systems/IrrotationalBns/BoundaryConditions/Neumann.cpp

src/Elliptic/Systems/IrrotationalBns/BoundaryConditions/Neumann.hpp

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

nilsvu · 2024-06-06T00:10:10Z

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

+ *  formulated as a set of coupled first-order PDEs.
+ *
+ * \details This system formulates the Irrotational Bns HSE equations for the
+ *  velocity potential \f$\Phi\f$. For a background matter distribution (given


(optional) note that you can write standard latex $\Phi$ and \begin{align}...\end{align}; you don't need those \f.

Oh huh, that makes things easier. I'll start doing that from now on

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

src/Elliptic/Systems/IrrotationalBns/Equations.cpp

nilsvu · 2024-06-06T20:30:03Z

src/Elliptic/Systems/IrrotationalBns/Equations.cpp

+  ::tenex::update(source_for_potential,
+                  (*source_for_potential)() +
+                      flux_for_potential(ti::I) *
+                          (log_deriv_of_lapse_over_specific_enthalpy(ti::i) -


Does this include the factor of \rho as well in D_i \ln(\alpha\rho / h)?

Also I think the Christoffels should be added, not subtracted? The equation is -\partial_i F^i + S = f please check.

Yes, it does include the factor of \rho, I'll change the variable name to make that obvious but maybe we can condense it in the future.

Hm, I would think the Christoffel should be be subtracted, we have the equation which looks like -D_i F^i + S_flat = f and then we expand the term D_i F^i = \partial_i F^i + Gamma^i_{ij}F^j so -Gamma^i_{ij} F^j gets added to the sources

Yes you're right it should be subtracted, but the other term as well if I did the math correctly

src/Elliptic/Systems/IrrotationalBns/IrrotationalBns.hpp

src/Elliptic/Systems/IrrotationalBns/Tags.hpp

nilsvu · 2024-06-06T20:40:14Z

src/Elliptic/Systems/IrrotationalBns/Tags.hpp

+ * \f]
+ */
+template <typename DataType>
+struct AuxiliaryVelocity : db::SimpleTag {


You can delete this tag, and also some or all of the tags below this (at least for now).

So basically If I deleted these now they would get added back in when the executable is added?

You can keep the ones you surely will need, like the Euler constant. But e.g. the FixedSources tag you don't need because there's already one in DataStructures/DataBox/Prefixes.hpp.

nilsvu · 2024-06-13T19:26:08Z

src/Elliptic/Systems/IrrotationalBns/BoundaryConditions/Neumann.cpp

+// The constant term on the RHS is nonlinear with respect to the
+// operator, so the linearized form of the RHS is zero.


"The term on the RHS is constant w.r.t. the variables, so the linearized form of the RHS is zero."

nilsvu · 2024-06-13T19:31:22Z

src/Elliptic/Systems/IrrotationalBns/BoundaryConditions/Neumann.hpp

-class Neumann : public elliptic::BoundaryConditions::BoundaryCondition<3> {
+/// Impose StarSurface boundary conditions
+/// \f[
+///    n_i F^i = \frac{1}{\sqrt{\partial_i \rho \partial_j \rho \gamma^{ij}}


This should be just n_i F^i = C/square(alpha) B^i n_i, right?

nilsvu · 2024-06-13T19:32:34Z

src/Elliptic/Systems/IrrotationalBns/BoundaryConditions/Neumann.cpp

  set_number_of_grid_points(n_dot_auxiliary_velocity, *velocity_potential);
-  get(*n_dot_auxiliary_velocity) = 0.0;
+  tenex::evaluate<>(n_dot_auxiliary_velocity,
+                    -euler_enthalpy_constant / square(lapse()) *


I think the minus sign is incorrect?

nilsvu · 2024-06-13T19:35:47Z

src/Elliptic/Systems/IrrotationalBns/Equations.cpp

@@ -47,6 +47,8 @@ void fluxes_on_face(
    const tnsr::I<DataVector, 3>& face_normal_vector,
    const tnsr::II<DataVector, 3>& rotational_shift_stress,
    const Scalar<DataVector>& velocity_potential) {
+  // potential optimization : precompute the produce of the `face_normal`


produce -> product

🥦 -> ✖️ 👍

nilsvu · 2024-06-13T19:42:40Z

src/Elliptic/Systems/IrrotationalBns/Equations.cpp

+      source_for_potential,
+      (*source_for_potential)() +
+          flux_for_potential(ti::I) *
+              (log_deriv_lapse_times_density_over_specific_enthalpy(ti::i) -


I think both of these terms should have a minus sign. Please check.

still a minus sign to fix (unless I got the math wrong)

nilsvu · 2024-06-13T22:18:24Z

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

+ * velocity potential \f$\Phi\f$. For a background matter distribution (given
+ * by the specific enthalpy h) and a background metric \f$\gamma_{ij}\f$.
+ * The velocity potential is defined by \f$D_i \Phi = h u_i\f$ with \f$u_i\f$
+ (the spatial part of) the four velocity


fix this line (and more below that don't start with *)

nilsvu · 2024-06-13T22:20:37Z

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

 *
 * The system can be formulated in terms of these fluxes and sources (see
 * `elliptic::protocols::FirstOrderSystem`):
 *
+ *
+ * \f{align*}
+ * D_i F^i + S = f


-\partial_i F^i + S = f (note the minus sign and \partial_i). This means some minus signs in the equations below need to be fixed too.

nilsvu · 2024-06-13T22:21:13Z

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

- * F^i_{v_j} &= u \delta^i_j \\
- * S_{v_j} &= v_j \\
- * f_{v_j} &= 0 \text{.}
+ * F^i &=  h u^i - \frac{B^j h u_j}{\alpha^2}B^i  \\


h u^i -> \partial_i \phi to make clear that the flux depends on the gradient of the potential

isaaclegred · 2024-06-14T05:41:51Z

I think everything should be done except removing the "flat" things. I'll do that in a separate commit after squashing all of this because I plan on keeping it around on a branch just in case it becomes useful for checking that the executable is working
.

nilsvu

You can squash everything

nilsvu · 2024-06-14T23:45:57Z

src/Elliptic/Systems/IrrotationalBns/Equations.cpp

+      source_for_potential,
+      (*source_for_potential)() +
+          flux_for_potential(ti::I) *
+              (log_deriv_lapse_times_density_over_specific_enthalpy(ti::i) -


still a minus sign to fix (unless I got the math wrong)

nilsvu · 2024-06-14T23:47:42Z

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

+ *
+ * \f{align*}
+ * F^i &=  D_i \phi - \frac{B^j D_j \phi}{\alpha^2}B^i  \\
+ * S &= -F^iD_i \left( \ln \frac{\alpha \rho}{h}\right) -\Gamma^i_{ij}F^j \\


Can you add spin terms already? This probably also means we should name the executable something other than "irrotational".

I'd like to have the irrotational case merged first, at that point it should be straightforward to extend to the spin case but also at that point things stop being "exact". I'd like to go through the sources and understand what's going on and write up notes before I add that though.

The minus sign is my bad I have confused myself thoroughly with signs.

Ok that's fine. Probably still we should not name things "irrotational", don't you think? That seems confusing once we add the spin terms (and we don't want to rename everything later).

src/Elliptic/Systems/IrrotationalBns/BoundaryConditions/CMakeLists.txt

nilsvu · 2024-06-18T03:55:50Z

src/Elliptic/Systems/IrrotationalBns/CMakeLists.txt

+  HEADERS
+  Equations.hpp
+  FirstOrderSystem.hpp
+  Geometry.hpp


remove Geometry.hpp

nilsvu · 2024-06-18T03:56:30Z

src/Elliptic/Systems/IrrotationalBns/Equations.cpp

@@ -47,6 +47,8 @@ void fluxes_on_face(
    const tnsr::I<DataVector, 3>& face_normal_vector,
    const tnsr::II<DataVector, 3>& rotational_shift_stress,
    const Scalar<DataVector>& velocity_potential) {
+  // potential optimization : precompute the produce of the `face_normal`


src/Elliptic/Systems/IrrotationalBns/Equations.cpp

src/Elliptic/Systems/IrrotationalBns/Equations.hpp

src/Elliptic/Systems/IrrotationalBns/Tags.hpp

tests/Unit/Elliptic/Systems/IrrotationalBns/BoundaryConditions/CMakeLists.txt

nilsvu · 2024-06-18T04:07:47Z

tests/Unit/Elliptic/Systems/IrrotationalBns/BoundaryConditions/Test_Neumann.cpp

+        const double euler_enthalpy_constant = 1.0;
+        apply_neumann_boundary_condition<false>(
+            make_not_null(&n_dot_auxilliary_velocity));
+        CHECK(get(n_dot_auxilliary_velocity) ==


You aren't actually comparing to the Python implementation I think?

tests/Unit/Elliptic/Systems/IrrotationalBns/Test_Tags.cpp

nilsvu · 2024-06-18T04:14:55Z

src/Elliptic/Systems/IrrotationalBns/FirstOrderSystem.hpp

+ *
+ * \f{align*}
+ * F^i &=  D_i \phi - \frac{B^j D_j \phi}{\alpha^2}B^i  \\
+ * S &= -F^iD_i \left( \ln \frac{\alpha \rho}{h}\right) -\Gamma^i_{ij}F^j \\


Ok that's fine. Probably still we should not name things "irrotational", don't you think? That seems confusing once we add the spin terms (and we don't want to rename everything later).

isaaclegred · 2024-06-21T18:19:50Z

@nilsvu I renamed everything to BnsInitialData in order to make it more general. Also fixed up comments.

nilsvu

You can squash everything

src/Elliptic/Systems/BnsInitialData/BoundaryConditions/Factory.hpp

tests/Unit/Elliptic/Systems/BnsInitialData/BoundaryConditions/Test_StarSurface.cpp

nilsvu · 2024-06-21T18:36:32Z

src/Elliptic/Systems/BnsInitialData/BoundaryConditions/StarSurface.hpp

+/*!
+ * Impose StarSurface boundary conditions:
+ * \f[
+ *   n_i F^i = \frac{C}{\alpha^2} B^i n_i.


Question: Since this is a Neumann boundary condition, it leaves the freedom for a constant offset in \Phi. I think the equations also leave this freedom. Have you thought about how to fix this freedom?

I haven't, I would think a reasonable choice would be that Phi should be zero at the center of the star or something like that, we probably can't set Phi=0 at the surface because I'm almost sure it's not "equipotential"

src/Elliptic/Systems/BnsInitialData/BoundaryConditions/StarSurface.hpp

src/Elliptic/Systems/BnsInitialData/Equations.hpp

nilsvu

Just have to appease CI :)

nilsvu self-requested a review June 4, 2024 19:24

nilsvu requested changes Jun 6, 2024

View reviewed changes

nilsvu requested changes Jun 13, 2024

View reviewed changes

nilsvu requested changes Jun 14, 2024

View reviewed changes

isaaclegred force-pushed the BNSInitialDataSystems branch 2 times, most recently from d7efda3 to a3d4259 Compare June 18, 2024 00:56

nilsvu requested changes Jun 18, 2024

View reviewed changes

nilsvu requested changes Jun 21, 2024

View reviewed changes

isaaclegred force-pushed the BNSInitialDataSystems branch 2 times, most recently from 5b4b7e0 to 5a3fd8a Compare June 27, 2024 18:29

nilsvu previously approved these changes Jun 28, 2024

View reviewed changes

isaaclegred dismissed nilsvu’s stale review via 393e308 July 2, 2024 19:57

isaaclegred force-pushed the BNSInitialDataSystems branch 2 times, most recently from 393e308 to 7af735e Compare July 5, 2024 16:41

Add BnsInitialData System

5981b40

isaaclegred force-pushed the BNSInitialDataSystems branch from 7af735e to 5981b40 Compare July 5, 2024 16:45

nilsvu approved these changes Jul 5, 2024

View reviewed changes

nilsvu changed the title ~~Add IrrotationalBns System~~ Add BnsInitialData System Jul 5, 2024

nilsvu merged commit b085f46 into sxs-collaboration:develop Jul 5, 2024
21 of 22 checks passed

		// The constant term on the RHS is nonlinear with respect to the
		// operator, so the linearized form of the RHS is zero.

Add BnsInitialData System #6054

Add BnsInitialData System #6054

Conversation

isaaclegred commented Jun 3, 2024

Proposed changes

Upgrade instructions

Code review checklist

Further comments

isaaclegred commented Jun 3, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nilsvu Jun 13, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

isaaclegred commented Jun 14, 2024

nilsvu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nilsvu Jun 14, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

isaaclegred commented Jun 21, 2024

nilsvu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nilsvu left a comment

Choose a reason for hiding this comment

nilsvu Jun 13, 2024 •

edited

Loading

nilsvu Jun 14, 2024 •

edited

Loading