Add tracer particles to phoebus #190

AstroBarker · 2023-12-06T22:45:34Z

I believe that the tracer treatment is now mature enough to have more eyes on it. It's a pretty big PR and has been awhile coming, so interact in whatever level you feel.

Summary

Here I add tracer particle support in phoebus. We create a tracer package that, if enabled, registers the appropriate swarm variables. See src/tracers/tracers.cpp. At present, the tracers are coupled in first order operator split, though in the future I'd like to explore improving this, after some phoebus task list refactoring. Tracers are advected in AdvectTracers(...), interpolating the appropriate quantities to the tracer positions and then pushing (with RK2). The "filling" of the tracers is handled in FillTracers registered with UserWorkBeforeOutput. To enable this, I added phoebus::UserWorkBeforeOutput to the driver.
It's a lot of quantities.. and may desire to shorten later, or otherwise make what is stored runtime modular?

To deal with black hole accretion, when FMKS geometry is enabled, we use PurgeParticles in src/fixup/fixup_particles.cpp If FMKS is not enabled, the constexpr if resolves to a no op. In principle, here, future "purge particles if.." conditions can be supplied.

At present, I have implemented tracers in the torus and advection pgens. The advection test yields relative errors on tracer positions of ~ 1e-14 after one period. Tracer orbits in a non-magnetized torus work as desired.

To run with tracers (when the pgen supports it), simply

<physics>
tracers = true

<tracers>
num_tracers = ...

TODO:

Add tracers into advection, at least. Good to have a simpler-than-torus test problem.
Do 3D, non-magnetized torus orbit test.
Unit & regression test.

cc: @bprather

Yurlungur

Overall this looks great! Thanks for your hard work @AstroBarker this is awesome. A few nitpicks below. And also tests seem to be not passing.

Yurlungur · 2023-12-07T00:55:45Z

inputs/blast_wave_2d.pin

@@ -93,6 +93,7 @@ hydro = true
 he = false
 3t = false
 rad = false
+tracers = false


do you need to add this to all of the pgens if the default is false?

Was using GetBoolean not GetOrAddBoolean(_, _, false). Fixed!

Yurlungur · 2023-12-07T01:03:27Z

src/pgen/torus.cpp


    auto ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
    auto jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
    auto kb = pmb->cellbounds.GetBoundsK(IndexDomain::interior);

    PackIndexMap imap;
-    auto v = rc->PackVariables({fluid_prim::bfield}, imap);
+    std::vector<std::string> vars = {fluid_prim::bfield, fluid_prim::density};


note that in #180 this syntax will change slightly. Should be easy to fix though.

Yurlungur · 2023-12-07T01:05:03Z

src/pgen/torus.cpp

+      // TODO (BLB): better way to construct & access indices array?
+      // Access as indices( flatten( k-kb.s, j-jb.s, i-ib.s ) )


I like @brryan 's idea of using a "volume-weighted" mask, as discussed on MM. We should look at the tracers example parthenon.

Note we don't have to fix this now. PR can go in even if the initialization is rough around the edges.

src/tracers/tracers.cpp

AstroBarker · 2023-12-07T14:15:34Z

re: tests failing, I had some old testing code in radiation_advection and radiation_equilibration problems that didn't work with tracers disabled. I've just reverted those files to pre-tracers, and will put tracers in the fluid advection problem.

AstroBarker · 2023-12-14T21:25:58Z

src/pgen/torus.cpp

+      // TODO (BLB): better way to construct & access indices array?
+      // Access as indices( flatten( k-kb.s, j-jb.s, i-ib.s ) )
+      const int mysize = (kb.e - kb.s + 1) * (jb.e - jb.s + 1) * (ib.e - ib.s + 1);
+      ParArray1D<int> indices("map from k,j,i to 1D particle array", mysize);


Weirdness on GPU accessing the indices gymnastics @Yurlungur . In Phoebus::ProblemGenerator::Torus::SumDiskCells kernel, n_cells_disk is always 0.

AH! I see the problem. indices is not guarded with an atomic. You need a Kokkos::atomic on the indices assignment. Alternatively, you need to make it a raw kokkos view with the atomic trait.

Yurlungur · 2023-12-14T22:44:48Z

src/pgen/torus.cpp

+            Real uphi;
+            if (r > rin) lnh = log_enthalpy(r, th, a, rin, angular_mom, uphi);
+            if (lnh > 0.0 && x1 > xh) {
+              indices(flatten(k - kb.s, j - jb.s, i - ib.s)) = n_cells_disk;


n_cells_disk is not guaranteed to be assigned in any particular order... so making indices have atomic writes may not solve this problem. But let's try that first.

AstroBarker · 2024-01-22T20:13:01Z

Some updates:

I have refactored the torus pgen to use rejection sampling. It is cleaner now, and doesn't fail on GPU. As implemented, it should be uniform in volume.
I have split the tracer rhs calculation into a separate function
Tracer integration is now RK2, as in nubhlight
Added Bernoulli parameter - u^0 h - 1 to tracer output.
Tracers are implemented into the advection pgen
- Implementation here could be improved. As is, we get nblocks * num_tracers particles. Should just move to post init...

In the advection test, relative errors on tracer positions was ~ 1e-14 or lower after one period.
I also ran a 3D, non-magnetized torus to see tracer orbits (notes: below only goes to t=500.. longer is running. And this was before the upgrade to RK2)

Yurlungur · 2024-01-22T20:19:15Z

Looks awesome! When do you think you'll be ready to remove the WIP markings? :)

AstroBarker · 2024-01-22T23:42:28Z

One last change: I've moved the FillTracers call from PostFillDerivedBlock into UserWorkBeforeOutput. To accommodate this, I created phoebus::UserWorkBeforeOutput in the driver and registered it in main. At present it just unpacks MeshBlockData and calls FillTracers, but it could accommodate extension in the future. Should cut down on unnecessary work. Let me know how this bit looks @Yurlungur.

(Also fixed something I forgot to commit in the torus init..).

Yurlungur

This looks great, @AstroBarker ! I'm happy to merge this, but I would like @brryan 's eyes on it if he gets the chance to look.

Also this is a question for Ben: Are particle IDs integers? Or longs? We might want to make them longs if possible (which might be a parthenon change). I recall hitting integer overflow in nubhlight.

src/tracers/tracers.cpp

Yurlungur · 2024-01-24T00:20:22Z

src/phoebus_driver.cpp

+void UserWorkBeforeOutput(MeshBlock *pmb, ParameterInput *pin) {
+  auto tracer_pkg = pmb->packages.Get("tracers");
+  bool do_tracers = tracer_pkg->Param<bool>("active");
+  if (do_tracers) {
+    auto &mbd = pmb->meshblock_data.Get();
+    tracers::FillTracers(mbd.get());
+  }
+}


👍 When particle meshblock packs work, we should refactor to use those :)

Yurlungur · 2024-01-24T00:21:49Z

src/pgen/torus.cpp

+            Real uphi;
+            if (r > rin) lnh = log_enthalpy(r, th, a, rin, angular_mom, uphi);
+            if (lnh > 0.0) {
+              const Real vol = coords.CellVolume(k, j, i);


This should have a factor of the metric in it I think. (It probably doesn't matter for problems we run, but we should include it.)

Yurlungur · 2024-01-24T00:23:00Z

src/pgen/torus.cpp

+            if (r > rin) lnh = log_enthalpy(r, th, a, rin, angular_mom, uphi);
+            if (lnh > 0.0 && x1 > xh) {
+              const Real vol = coords.CellVolume(k, j, i);
+              number_block_reduce += vol;


factor of metric determinant missing (see above)

Yurlungur · 2024-01-24T00:23:41Z

src/pgen/torus.cpp

+            if (swarm_d.IsActive(n)) {
+              int k, j, i;
+              swarm_d.Xtoijk(x(n), y(n), z(n), i, j, k);
+              const Real vol = coords.CellVolume(k, j, i);


missing detgamma

Yurlungur · 2024-01-24T00:24:18Z

src/pgen/advection.cpp

+            // sample in ye ball
+            Real r2 = 1.0 + rin * rin; // init > rin^2
+            while (r2 > rin * rin) {
+              x(n) = x_min + rng_gen.drand() * (x_max - x_min);
+              y(n) = y_min + rng_gen.drand() * (y_max - y_min);
+              z(n) = z_min + rng_gen.drand() * (z_max - z_min);
+              r2 = x(n) * x(n) + y(n) * y(n) + z(n) * z(n);
+            }
+            id(n) = num_tracers_total * gid + n;


Yurlungur · 2024-01-24T00:26:16Z

Actually @AstroBarker before merge, please correct the missing metric factor in the volumes in the torus calculation.

Yurlungur · 2024-02-16T00:52:55Z

Since @brryan never had a chance to look. I'm going to click the button. Ben whenever you get the chance, I'd still like your thoughts.

Yurlungur · 2024-02-16T00:53:25Z

Also thanks for all your hard work on this @AstroBarker !

AstroBarker added 28 commits October 21, 2022 15:38

feat(wip): begin adding tracer particles

7af4981

Merge branch 'main' into blb/tracers

71f2b67

feat: tracers loaded in driver. begin AdvectTracers.

a908fae

feat: draft of AdvectTracers

4f6f46f

feat(wip): working on tracer initialization for torus. bugs.

4dca3f7

re-comment out meshblock section

1d4c85f

fix: a few small fixes in tracers.cpp. was setting num_tracers as a Real

24500c2

feat: simple tracer setup working in rad_eq and rad_advect

74ac2c7

feat: now interp geometry quantities. Next is Lorentz factor?

a76f8fb

Merge branch 'main' into blb/tracers

577f1bc

Merge branch 'main' into blb/tracers

2474f22

conflict resolve

9a1dc95

updates to tracer interpolation, sme pgen edits.

8344a7a

missed files

2459720

merge main into blb/tracers

abec1ed

progress initializing tracers in torus. need to work locally

14d689f

bug typo fix

618389b

feat: drafted tracer init for torus

e031306

feat: update tracers in driver, draft FillTracers, bugged on torus

cc9eb33

feat: draft PurgeParticles to remove particles beneath event horizon

f255df3

merging main

9548463

fix: reformat so that tracers are properly disabled when tracers = false

129b378

fix torus tracer indexing

b2f2c9d

small renamings to driver

1447e1f

tracer fixes and quality of life updates

4b798b8

final B logic, preambles

d3cfd80

Merge branch 'main' into blb/tracers

93048ea

add tracers = false to all input decks

660dafd

AstroBarker added the enhancement New feature or request label Dec 6, 2023

AstroBarker requested a review from Yurlungur December 6, 2023 22:45

Yurlungur approved these changes Dec 7, 2023

View reviewed changes

AstroBarker added 2 commits December 7, 2023 13:36

revert rad advection and equilbration

c24d73d

fix: use GetOrAddBoolean for tracers pkg

7c39420

restore test default values for rad_eq

692c465

AstroBarker commented Dec 14, 2023

View reviewed changes

Yurlungur reviewed Dec 14, 2023

View reviewed changes

AstroBarker added 7 commits January 10, 2024 20:59

feat: refactor torus tracer init. rejection sample. todo: gpu

e6dce7d

fix: merge main, refeactor, fix bad host access on device

0af88eb

fmt

5d8fbf1

feat: add bernoulli parameter to tracers

923a97e

fmt...

e6aa218

feat: move to RK2. add tracers to advection pgen

a0eaca3

remove debug print

209a956

AstroBarker added 3 commits January 22, 2024 22:50

feat: move FillTracers to UserWorkBeforeOutput

a273147

doc: update some tracer comments

989b262

fix: forgot to commit fix for torus rejection sample loop

8d3dd80

AstroBarker changed the title ~~WIP: Add tracer particles to phoebus~~ Add tracer particles to phoebus Jan 22, 2024

AstroBarker added 2 commits January 22, 2024 23:45

fmt

6b79480

fix: do NOT fill tracers when there are no tracers...

a00d395

Yurlungur approved these changes Jan 24, 2024

View reviewed changes

fix: add missing detgammas for volume terms

baca991

Yurlungur merged commit 4e3d2c7 into main Feb 16, 2024
2 checks passed

Yurlungur deleted the blb/tracers branch February 16, 2024 00:53

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add tracer particles to phoebus #190

Add tracer particles to phoebus #190

AstroBarker commented Dec 6, 2023 •

edited

Loading

Yurlungur left a comment

Yurlungur Dec 7, 2023

AstroBarker Dec 7, 2023

Yurlungur Dec 7, 2023

Yurlungur Dec 7, 2023

Yurlungur Dec 7, 2023

AstroBarker commented Dec 7, 2023

AstroBarker Dec 14, 2023

Yurlungur Dec 14, 2023

Yurlungur Dec 14, 2023

AstroBarker commented Jan 22, 2024

Yurlungur commented Jan 22, 2024

AstroBarker commented Jan 22, 2024

Yurlungur left a comment

Yurlungur Jan 24, 2024

Yurlungur Jan 24, 2024

Yurlungur Jan 24, 2024

Yurlungur Jan 24, 2024

Yurlungur Jan 24, 2024

Yurlungur commented Jan 24, 2024

Yurlungur commented Feb 16, 2024

Yurlungur commented Feb 16, 2024

		// TODO (BLB): better way to construct & access indices array?
		// Access as indices( flatten( k-kb.s, j-jb.s, i-ib.s ) )

Add tracer particles to phoebus #190

Add tracer particles to phoebus #190

Conversation

AstroBarker commented Dec 6, 2023 • edited Loading

Summary

TODO:

Yurlungur left a 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

AstroBarker commented Dec 7, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

AstroBarker commented Jan 22, 2024

Yurlungur commented Jan 22, 2024

AstroBarker commented Jan 22, 2024

Yurlungur left a 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

Yurlungur commented Jan 24, 2024

Yurlungur commented Feb 16, 2024

Yurlungur commented Feb 16, 2024

AstroBarker commented Dec 6, 2023 •

edited

Loading