Add 0th order puncture field

src/Evolution/Systems/CurvedScalarWave/Worldtube/CMakeLists.txt

@@ -9,13 +9,15 @@ spectre_target_sources(
   ${LIBRARY}
   PRIVATE
   Tags.cpp
+  PunctureField.cpp
   )
 
 spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   Tags.hpp
+  PunctureField.hpp
   Worldtube.hpp
   )
 

src/Evolution/Systems/CurvedScalarWave/Worldtube/PunctureField.cpp

@@ -0,0 +1,322 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/CurvedScalarWave/Worldtube/PunctureField.hpp"
+
+#include "DataStructures/DataBox/Prefixes.hpp"
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Variables.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+
+void compute_puncture_field(
+    gsl::not_null<Variables<tmpl::list<
+        Tags::Psi, ::Tags::dt<Tags::Psi>,
+        ::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>>*>
+        evolved_vars,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
+    const double orbital_radius, const double bh_mass, const size_t order) {
+  if (order == 0) {
+    compute_puncture_field_0(evolved_vars, coords, time, orbital_radius,
+                             bh_mass);
+  } else {
+    ERROR(
+        "Puncture order is only implemented up to order 0 but requested order "
+        << order);
+  }
+}
+
+void compute_puncture_field_0(
+    gsl::not_null<Variables<tmpl::list<
+        Tags::Psi, ::Tags::dt<Tags::Psi>,
+        ::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>>*>
+        evolved_vars,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
+    const double orbital_radius, const double BH_mass) {
+  evolved_vars->initialize(get<0>(coords).size());
+  const double r0 = orbital_radius;
+  const double M = BH_mass;
+  const double w = pow(r0, -1.5);
+  const double t = time;
+
+  const double charge_pos_x = r0 * cos(w * time);
+  const double charge_pos_y = r0 * sin(w * time);
+  const double charge_pos_z = 0.;
+
+  const DataVector& x = get<0>(coords);
+  const DataVector& y = get<1>(coords);
+  const DataVector& z = get<2>(coords);
+  for (size_t i = 0; i < get<0>(coords).size(); ++i) {
+    const double Dx = x[i] - charge_pos_x;
+    const double Dy = y[i] - charge_pos_y;
+    const double Dz = z[i] - charge_pos_z;
+
+    const double x0 = 1.0 / r0;
+    const double x1 = 3 * M;
+    const double x2 = -r0 + x1;
+    const double x3 = sqrt(-x0 * x2);
+    const double x4 = pow(M, 3.0 / 2.0);
+    const double x5 = 4 * x4;
+    const double x6 = x3 * x5;
+    const double x7 = pow(Dx, 2);
+    const double x8 = pow(Dy, 2);
+    const double x9 = -x8;
+    const double x10 = x7 + x9;
+    const double x11 = t * w;
+    const double x12 = cos(x11);
+    const double x13 = sin(x11);
+    const double x14 = x12 * x13;
+    const double x15 = pow(x13, 2);
+    const double x16 = Dy * x15;
+    const double x17 = Dx * x16;
+    const double x18 = pow(x12, 2);
+    const double x19 = -Dx * Dy * x18 + x17;
+    const double x20 = pow(r0, 2);
+    const double x21 = sqrt(-x2 * x20);
+    const double x22 = x20 * x21;
+    const double x23 = x22 * (-x10 * x14 - x19);
+    const double x24 = pow(x20, 5.0 / 2.0);
+    const double x25 = pow(Dz, 2);
+    const double x26 = x25 + x8;
+    const double x27 = x26 + x7;
+    const double x28 = pow(r0, 3);
+    const double x29 = pow(M, 2);
+    const double x30 = Dx * x12;
+    const double x31 = Dy * x13;
+    const double x32 = x30 * x31;
+    const double x33 = 10 * x32;
+    const double x34 = 6 * x8;
+    const double x35 = 9 * x25;
+    const double x36 = 6 * x7;
+    const double x37 =
+        x29 * (x15 * (x35 + x36 + x8) + x18 * (x34 + x35 + x7) - x33);
+    const double x38 = pow(r0, 4);
+    const double x39 = 2 * x32;
+    const double x40 = 4 * x7;
+    const double x41 = 5 * x8;
+    const double x42 = 6 * x25;
+    const double x43 = 5 * x7;
+    const double x44 = 4 * x8;
+    const double x45 =
+        M * (x15 * (x42 + x43 + x44) + x18 * (x40 + x41 + x42) - x39);
+    const double x46 = x30 + x31;
+    const double x47 = pow(M, 3);
+    const double x48 = 6 * x47;
+    const double x49 = pow(x46, 2) * x48;
+    const double x50 = x20 * x49 + x24 * x27 + x28 * x37 - x38 * x45;
+    const double x51 = 1.0 / x20;
+    const double x52 = pow(x20, 3.0 / 2.0);
+    const double x53 = 6 * M;
+    const double x54 = 9 * x29;
+    const double x55 = r0 * x54 - x20 * x53 + x52;
+    const double x56 = 1.0 / x55;
+    const double x57 = x51 * x56;
+    const double x58 = Dx * x18;
+    const double x59 = Dy * x58;
+    const double x60 = -x10;
+    const double x61 = -x2;
+    const double x62 = sqrt(x20 * x61);
+    const double x63 = x62 * sqrt(x0 * x61);
+    const double x64 = r0 * x63;
+    const double x65 = x5 * x64;
+    const double x66 = pow(x2, -2);
+    const double x67 = r0 * x66 / x52;
+    const double x68 = 2 * x25;
+    const double x69 = -x7;
+    const double x70 = 2 * x8;
+    const double x71 =
+        x15 * (-x68 - 2 * x7 - x9) + x18 * (-x68 - x69 - x70) + 6 * x32;
+    const double x72 = 1.0 / x24;
+    const double x73 = x28 * x72;
+    const double x74 = x46 * x73;
+    const double x75 = 1.0 / x62;
+    const double x76 = Dx * x13;
+    const double x77 = Dy * x12;
+    const double x78 = -x77;
+    const double x79 = x76 + x78;
+    const double x80 = sqrt(M) * r0;
+    const double x81 = x79 * x80;
+    const double x82 = sqrt(x20);
+    const double x83 = 1.0 / x82;
+    const double x84 = sqrt(-x1 * x83 + 1);
+    const double x85 = 1.0 / x84;
+    const double x86 = 2 * M;
+    const double x87 = x0 * x85 * x86;
+    const double x88 = -x46 * x87;
+    const double x89 = x75 * x81 + x88;
+    const double x90 = 4 * x32;
+    const double x91 = x15 * (x10 + x25);
+    const double x92 = x26 + x69;
+    const double x93 = x51 * x85;
+    const double x94 = 2 * x93;
+    const double x95 = x71 * x74 - x89 * x94 * (-x18 * x92 + x90 - x91);
+    const double x96 = M * x20;
+    const double x97 = x95 * x96;
+    const double x98 = 1.0 / x21;
+    const double x99 = x81 * x98 + x88;
+    const double x100 = x83 * x86;
+    const double x101 = x100 * x18 + 1;
+    const double x102 = M * x83;
+    const double x103 =
+        x101 * x7 + x102 * x15 * x70 + x102 * x90 + x26 + pow(x99, 2);
+    const double x104 = 1.0 / x103;
+    const double x105 = x104 * x55;
+    const double x106 = x105 * x97;
+    const double x107 = x23 * x5 * x84 + x50;
+    const double x108 = 1.0 / x107;
+    const double x109 = (1.0 / 2.0) * x106 * x108;
+    const double x110 = x46 * x48;
+    const double x111 = x12 * x13 * x8 - x14 * x7 - x19;
+    const double x112 = r0 * x29;
+    const double x113 = -x90;
+    const double x114 = 2 * x4;
+    const double x115 = x114 * x63;
+    const double x116 = r0 * x53;
+    const double x117 = r0 * x82;
+    const double x118 = 1.0 / (-x116 + x117 + x54);
+    const double x119 = x18 * x36 * x47;
+    const double x120 = x15 * x34 * x47;
+    const double x121 = x18 * x7;
+    const double x122 = x15 * x8;
+    const double x123 = x15 * x96;
+    const double x124 = x18 * x96;
+    const double x125 = x112 * x15;
+    const double x126 = x112 * x18;
+    const double x127 = x42 * x96;
+    const double x128 = 12 * x32 * x47;
+    const double x129 = x21 * x6;
+    const double x130 = x21 * x3;
+    const double x131 = x14 * x4;
+    const double x132 = x130 * x131;
+    const double x133 = pow(
+        x0 * x118 *
+            (x112 * x121 + x112 * x122 - x112 * x33 + x119 + x120 - x123 * x43 -
+             x123 * x44 - x124 * x40 - x124 * x41 + x125 * x35 + x125 * x36 +
+             x126 * x34 + x126 * x35 - x127 * x15 - x127 * x18 + x128 -
+             x129 * x17 + x129 * x59 - x132 * x40 + x132 * x44 + x25 * x52 +
+             x39 * x96 + x52 * x7 + x52 * x8),
+        -3.0 / 2.0);
+    const double x134 = x133 * x56;
+    const double x135 = 5 * x77;
+    const double x136 = 5 * x76;
+    const double x137 = x18 * x76;
+    const double x138 = x15 * x77;
+    const double x139 = x137 - x138;
+    const double x140 = x139 + x15 * (x136 - 4 * x77) + x18 * (-x135 + 4 * x76);
+    const double x141 = M * x38;
+    const double x142 = -x135 * x15 + x136 * x18 + x15 * (6 * x76 + x78) +
+                        x18 * (x76 - 6 * x77);
+    const double x143 = x142 * x29;
+    const double x144 = pow(x12, 3);
+    const double x145 = pow(x13, 3);
+    const double x146 = x76 + x77;
+    const double x147 = -Dx * x12 * x15 + Dx * x144 - Dy * x13 * x18 +
+                        Dy * x145 + 2 * x14 * x146;
+    const double x148 = x115 * x20;
+    const double x149 = x15 + x18;
+    const double x150 = -x113 - x18 * x92 - x91;
+    const double x151 = x20 * x29;
+    const double x152 = M * x52;
+    const double x153 = x15 * x152;
+    const double x154 = x152 * x18;
+    const double x155 = x15 * x151;
+    const double x156 = x151 * x18;
+    const double x157 = x152 * x42;
+    const double x158 = x131 * x64;
+    const double x159 = sqrt(
+        x67 * (r0 * x119 + r0 * x120 + r0 * x128 + x121 * x151 + x122 * x151 -
+               x15 * x157 - x151 * x33 + x152 * x39 - x153 * x43 - x153 * x44 -
+               x154 * x40 - x154 * x41 + x155 * x35 + x155 * x36 + x156 * x34 +
+               x156 * x35 - x157 * x18 - x158 * x40 + x158 * x44 - x17 * x65 +
+               x25 * x38 + x38 * x7 + x38 * x8 + x59 * x65));
+    const double x160 = x108 * x159 * x55;
+    const double x161 = x104 * x160;
+    const double x162 = pow(x103, -2);
+    const double x163 = x75 * x80;
+    const double x164 = x114 * x147;
+    const double x165 = x22 * x84;
+    const double x166 = x159 / pow(x107, 2);
+    const double x167 = 1.0 / x159;
+    const double x168 = Dx * x15;
+    const double x169 = -Dy * x12 * x13 + 5 * x168 + 4 * x58;
+    const double x170 = x13 * x30;
+    const double x171 = -Dy * x18 + x16 + 2 * x170;
+    const double x172 = -x171;
+    const double x173 = x12 * x31;
+    const double x174 = 6 * x168 - 5 * x173 + x58;
+    const double x175 = x28 * x29;
+    const double x176 = x110 * x12;
+    const double x177 = Dx * x24 + x174 * x175 + x176 * x20;
+    const double x178 = x133 * x57;
+    const double x179 = 2 * x74;
+    const double x180 = x71 * x73;
+    const double x181 = -x168 + 2 * x173 + x58;
+    const double x182 = -x12 * x87 + x13 * x80 * x98;
+    const double x183 = x161 * x96;
+    const double x184 = (1.0 / 2.0) * x183;
+    const double x185 = x160 * x162 * x97;
+    const double x186 = x114 * x165;
+    const double x187 = x106 * x166;
+    const double x188 = r0 * x115;
+    const double x189 = Dy * x18;
+    const double x190 = -Dx * x12 * x13 + 4 * x16 + 5 * x189;
+    const double x191 = x16 - 5 * x170 + 6 * x189;
+    const double x192 = x110 * x13;
+    const double x193 = Dy * x24 + x175 * x191 + x192 * x20;
+    const double x194 = x12 * x163 + x13 * x87;
+    const double x195 = x15 * x54;
+    const double x196 = x18 * x54;
+    const double x197 = x38 * x53;
+
+    get(get<Tags::Psi>(*evolved_vars))[i] =
+        x109 * sqrt(x67 * (r0 * x49 + x20 * x37 + x27 * x38 - x45 * x52 +
+                           x65 * (x14 * x60 - x17 + x59))) +
+        pow(x57 * (x23 * x6 + x50), -1.0 / 2.0);
+    get(get<::Tags::dt<Tags::Psi>>(*evolved_vars))[i] =
+        w * (-M * x105 * x166 * x28 * x95 *
+                 (-x140 * x141 + x142 * x28 * x29 - x164 * x165 + x24 * x79) +
+             M * x108 * x159 * x162 * x28 * x55 * x95 *
+                 (x100 * x137 - x101 * x76 - x149 * x163 * x89 + x77) -
+             x0 * x134 * (-x140 * x141 + x143 * x28 - x147 * x148 + x24 * x79) -
+             x109 * x167 * x66 * x83 *
+                 (r0 * x130 * x164 + x140 * x152 - x143 * x20 - x38 * x79) -
+             x134 * (-x110 * x79 - 5 * x111 * x112 + x111 * x96 +
+                     x115 * (x10 * x15 + x113 + x18 * x60)) -
+             x161 * x20 *
+                 (2 * x141 * x46 * x72 * (Dx * x145 - Dy * x144 + x139) +
+                  x149 * x150 * x4 * x85 * x98 -
+                  x87 * x99 * (2 * x137 - 2 * x138 + x146 * x15 - x146 * x18)));
+    get<0>(get<::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>(
+        *evolved_vars))[i] =
+        (1.0 / 2.0) * M * r0 * x104 * x108 * x167 * x55 * x66 * x83 * x95 *
+            (Dx * x38 + r0 * x176 + x151 * x174 - x152 * x169 + x172 * x188) -
+        x178 * (-x141 * x169 + x148 * x172 + x177) -
+        x184 *
+            (-x12 * x180 + 2 * x150 * x182 * x51 * x85 -
+             x179 * (-2 * x168 + 3 * x173 + x58) + 4 * x181 * x51 * x85 * x99) -
+        x185 * (Dx * x101 + x100 * x173 + x182 * x99) -
+        x187 * (-x141 * x169 - x171 * x186 + x177);
+    get<1>(get<::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>(
+        *evolved_vars))[i] =
+        (1.0 / 2.0) * M * r0 * x104 * x108 * x167 * x55 * x66 * x83 * x95 *
+            (Dy * x38 + r0 * x192 + x151 * x191 - x152 * x190 + x181 * x188) -
+        x178 * (-x141 * x190 + x148 * x181 + x193) -
+        x184 * (-x13 * x180 - x150 * x194 * x94 + 4 * x171 * x51 * x85 * x99 -
+                x179 * (x16 + 3 * x170 - 2 * x189)) -
+        x185 * (Dy + x100 * x16 + x100 * x170 - x194 * x99) -
+        x187 * (-x141 * x190 + x181 * x186 + x193);
+    get<2>(get<::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>(
+        *evolved_vars))[i] =
+        Dz *
+        ((1.0 / 2.0) * M * r0 * x104 * x108 * x167 * x55 * x66 * x83 * x95 *
+             (x149 * x20 * x54 - x149 * x52 * x53 + x38) -
+         x118 * x133 * (-x116 * x15 - x116 * x18 + x117 + x195 + x196) -
+         2 * x149 * x183 * (x74 - x93 * x99) - x185 -
+         x187 * (-x15 * x197 - x18 * x197 + x195 * x28 + x196 * x28 + x24));
+  }
+}
+}  // namespace CurvedScalarWave::Worldtube

src/Evolution/Systems/CurvedScalarWave/Worldtube/PunctureField.hpp

@@ -0,0 +1,59 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+
+#include "DataStructures/DataBox/Prefixes.hpp"
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "DataStructures/Variables.hpp"
+#include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
+#include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace CurvedScalarWave::Worldtube {
+/*!
+ * \brief Computes the puncture/singular field \f$\Psi^\mathcal{P}\f$ of a
+ * scalar charge in circular orbit around a Schwarzschild black hole as
+ * described in \ref Detweiler2003.
+ *
+ * \details The field is computed using a Detweiler-Whiting singular
+ * Green's function and perturbatively expanded in the geodesic distance from
+ * the particle. It solves the inhomogeneous wave equation
+ * \f{align*}{
+ * \box \Psi^\mathcal{P} = -4 \pi q \int \sqrt{-g} \delta^4(x^i, z(\tau)) d \tau
+ * \f}
+ * where \f$q\f$ is the scalar charge and \f$z(\tau)\f$ is the worldline of the
+ * particle. The expression is expanded up to order n in geodesic distance and
+ * transformed to Kerr-Schild coordinates.
+ *
+ * The function given here assumes that the particle has scalar charge \f$q=1$\f
+ * and is on a fixed geodesic orbit \f$z(t) = (R \cos{\omega t}, R \sin{\omega
+ * t}, 0)\f$ around a Schwarzschild black hole, where R is the orbital radius
+ * and \f$ \omega = R^(-3/2) \f$. It returns the singular field at the requested
+ * coordinates as well as its time and spatial derivative.
+ *
+ * \note The expressions were computed with Mathematica and optimized by
+ * applying manual common subexpression elimination with sympy.
+ * The calculation is currently not vectorized to avoid a large amount of
+ * temporary allocations but this may be added as a future optimization. In its
+ * current form it does not noticeably impact the simulation time.
+ */
+void compute_puncture_field(
+    gsl::not_null<Variables<tmpl::list<
+        Tags::Psi, ::Tags::dt<Tags::Psi>,
+        ::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>>*>
+        evolved_vars,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
+    const double orbital_radius, const double bh_mass, size_t order);
+
+void compute_puncture_field_0(
+    gsl::not_null<Variables<tmpl::list<
+        Tags::Psi, ::Tags::dt<Tags::Psi>,
+        ::Tags::deriv<Tags::Psi, tmpl::size_t<3>, Frame::Inertial>>>*>
+        evolved_vars,
+    const tnsr::I<DataVector, 3, Frame::Inertial>& coords, const double time,
+    const double orbital_radius, const double bh_mass);
+}  // namespace CurvedScalarWave::Worldtube