From e721ac56b2b493b0ee9f64a8e50cbab86ca599f4 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Tue, 9 May 2023 20:18:56 +0000 Subject: [PATCH 1/4] feat: drafted GetKSource --- src/radiation/geodesics.hpp | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/src/radiation/geodesics.hpp b/src/radiation/geodesics.hpp index fcfda22f2..c27a9ae49 100644 --- a/src/radiation/geodesics.hpp +++ b/src/radiation/geodesics.hpp @@ -39,7 +39,35 @@ KOKKOS_INLINE_FUNCTION void GetKSource(Real &X0, Real &X1, Real &X2, Real &X3, Real &Kcov0, Real &Kcov1, Real &Kcov2, Real &Kcov3, Real &Kcon0, const Geometry::CoordSysMeshBlock &geom, Real source[4]) { - SPACETIMELOOP(mu) { source[mu] = 0.; } + Real Kcov[NDFULL] = {Kcov0, Kcov1, Kcov2, Kcov3}; + Real Xm[NDFULL], Xp[NDFULL]; + Real Gconm[NDFULL][NDFULL], Gconp[NDFULL][NDFULL]; + Real dG; + + constexpr Real DELTA = 1.0e-6; + + SPACETIMELOOP(mu) { + source[mu] = 0.; + Xp[0] = X0; + Xp[1] = X1; + Xp[2] = X2; + Xp[3] = X3; + SPACETIMELOOP(nu) { Xm[nu] = Xp[nu]; } + + Xm[mu] -= DELTA; + Xp[mu] += DELTA; + + geom.SpacetimeMetricInverse(Xm[0], Xm[1], Xm[2], Xm[3], Gconm); + geom.SpacetimeMetricInverse(Xp[0], Xp[1], Xp[2], Xp[3], Gconp); + + SPACETIMELOOP2(nu, kap) { + dG = (Gconp[nu][kap] - Gconm[nu][kap]) / (Xp[mu] - Xm[mu]); + source[mu] += Kcov[nu] * Kcov[kap] * dG; + } + + source[mu] *= -1.0 / (2.0 * Kcon0); + + } } KOKKOS_INLINE_FUNCTION @@ -80,6 +108,11 @@ void PushParticle(Real &X0, Real &X1, Real &X2, Real &X3, Real &Kcov0, Real &Kco X1 += 0.5 * dt * (c1[1] + c2[1]); X2 += 0.5 * dt * (c1[2] + c2[2]); X3 += 0.5 * dt * (c1[3] + c2[3]); + if ( X1 < 0.0 ) { + std::printf("X1 < 0 in Push %f %f %f %f\n ", X1, c1[1], c2[1], X1-0.5*dt*(c1[1]+c2[1])); + std::printf("X2 in Push %f %f %f %f\n ", X2, c1[2], c2[2], X2-0.5*dt*(c1[2]+c2[2])); + std::printf("X3 in Push %f %f %f %f\n ", X3, c1[3], c2[3], X3-0.5*dt*(c1[3]+c2[3])); + } Kcov0 += 0.5 * dt * (d1[0] + d2[0]); Kcov1 += 0.5 * dt * (d1[1] + d2[1]); From fdf02a2c18a23510b251ceb34847e57b90bf252d Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Tue, 9 May 2023 21:26:07 +0000 Subject: [PATCH 2/4] format and remove unnecessary prints --- scripts/python/plot_snap2d.py | 3 +-- src/radiation/geodesics.hpp | 10 ++-------- 2 files changed, 3 insertions(+), 10 deletions(-) diff --git a/scripts/python/plot_snap2d.py b/scripts/python/plot_snap2d.py index 9602563f5..321467be8 100755 --- a/scripts/python/plot_snap2d.py +++ b/scripts/python/plot_snap2d.py @@ -70,7 +70,6 @@ def plot_dump( coordname = "g.n.coord" coord = data.Get(coordname, False) z = coord[:, 3, :, nG - 1 : -1 - nG, nG - 1 : -1 - nG] - print(np.shape(z)) y = coord[:, 2, :, nG - 1 : -1 - nG, nG - 1 : -1 - nG] x = coord[:, 1, :, nG - 1 : -1 - nG, nG - 1 : -1 - nG] @@ -100,7 +99,7 @@ def plot_dump( fig = plt.figure() p = fig.add_subplot(111, aspect=1) for i in range(NB): - val = q[i, 0, :, :] + val = 100.0 * (q[i, 0, :, :] - 0.225) if len(val.shape) > 2: print("WARNING plotting the 0th index of multidimensional variable!") val = val[:, :, 0] diff --git a/src/radiation/geodesics.hpp b/src/radiation/geodesics.hpp index c27a9ae49..48ecb5898 100644 --- a/src/radiation/geodesics.hpp +++ b/src/radiation/geodesics.hpp @@ -46,7 +46,7 @@ void GetKSource(Real &X0, Real &X1, Real &X2, Real &X3, Real &Kcov0, Real &Kcov1 constexpr Real DELTA = 1.0e-6; - SPACETIMELOOP(mu) { + SPACETIMELOOP(mu) { source[mu] = 0.; Xp[0] = X0; Xp[1] = X1; @@ -66,8 +66,7 @@ void GetKSource(Real &X0, Real &X1, Real &X2, Real &X3, Real &Kcov0, Real &Kcov1 } source[mu] *= -1.0 / (2.0 * Kcon0); - - } + } } KOKKOS_INLINE_FUNCTION @@ -108,11 +107,6 @@ void PushParticle(Real &X0, Real &X1, Real &X2, Real &X3, Real &Kcov0, Real &Kco X1 += 0.5 * dt * (c1[1] + c2[1]); X2 += 0.5 * dt * (c1[2] + c2[2]); X3 += 0.5 * dt * (c1[3] + c2[3]); - if ( X1 < 0.0 ) { - std::printf("X1 < 0 in Push %f %f %f %f\n ", X1, c1[1], c2[1], X1-0.5*dt*(c1[1]+c2[1])); - std::printf("X2 in Push %f %f %f %f\n ", X2, c1[2], c2[2], X2-0.5*dt*(c1[2]+c2[2])); - std::printf("X3 in Push %f %f %f %f\n ", X3, c1[3], c2[3], X3-0.5*dt*(c1[3]+c2[3])); - } Kcov0 += 0.5 * dt * (d1[0] + d2[0]); Kcov1 += 0.5 * dt * (d1[1] + d2[1]); From a9410f64a2a3c4d83a03cedf898718194ec73f43 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Tue, 9 May 2023 21:44:03 +0000 Subject: [PATCH 3/4] remove ye plot hack --- scripts/python/plot_snap2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/python/plot_snap2d.py b/scripts/python/plot_snap2d.py index 321467be8..e6a891620 100755 --- a/scripts/python/plot_snap2d.py +++ b/scripts/python/plot_snap2d.py @@ -99,7 +99,7 @@ def plot_dump( fig = plt.figure() p = fig.add_subplot(111, aspect=1) for i in range(NB): - val = 100.0 * (q[i, 0, :, :] - 0.225) + val = q[i, 0, :, :] if len(val.shape) > 2: print("WARNING plotting the 0th index of multidimensional variable!") val = val[:, :, 0] From 4e8760568b8b3292e948db5299ee121899241735 Mon Sep 17 00:00:00 2001 From: Brandon Barker Date: Tue, 9 May 2023 21:45:06 +0000 Subject: [PATCH 4/4] add comment --- src/radiation/geodesics.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/radiation/geodesics.hpp b/src/radiation/geodesics.hpp index 48ecb5898..30e68f554 100644 --- a/src/radiation/geodesics.hpp +++ b/src/radiation/geodesics.hpp @@ -46,6 +46,7 @@ void GetKSource(Real &X0, Real &X1, Real &X2, Real &X3, Real &Kcov0, Real &Kcov1 constexpr Real DELTA = 1.0e-6; + // TODO(BLB) extend later to exploit spacetime symmetries. Template on geometry? SPACETIMELOOP(mu) { source[mu] = 0.; Xp[0] = X0;