Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update Parthenon and move to type-based variables #191

Merged
merged 15 commits into from
Dec 11, 2023
Merged

Conversation

Yurlungur
Copy link
Collaborator

@Yurlungur Yurlungur commented Dec 7, 2023

PR Summary

This PR was originally going to be #180. However, after making the changes described below, I decided this was a big enough changeset that it should go in first, before continuing to full sparse packs.

In this PR I move all our variables to the type-based variables. What this means is:

  1. I updated to current parthenon develop, which necessitated changing some calls to GetParentPointer (as this API was cleaned up a little bit).
  2. I added macros that declare variable types per variable and use them to declare variables in phoebus_utils/variables.hpp
  3. Because variable names are now types, not strings, the strings must be evoked via the variable::name() method.

After this goes through, I will start porting each variable pack to a sparse pack, which should eliminate much of the index gymnastics going on throughout the code.

@curiousmiah @mari2895 this is the first step in the refactor we discussed forever ago, and it should make said refactor much simpler to pursue. I also realize it'll be necessary to add face-centered constrained transport, as the topological elements are only exposed through sparse packs.

PR Checklist

  • Adds a test for any bugs fixed. Adds tests for new features.
  • Format your changes by calling scripts/bash/format.sh.
  • Explain what you did.

@Yurlungur Yurlungur added the enhancement New feature or request label Dec 7, 2023
@Yurlungur Yurlungur self-assigned this Dec 7, 2023
Copy link
Collaborator

@mari2895 mari2895 left a comment

Choose a reason for hiding this comment

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

@Yurlungur this is a lot of work, I appreciate it a lot! Looks good to me!

@Yurlungur
Copy link
Collaborator Author

Yurlungur commented Dec 7, 2023

This is an example of how the new typed variables and sparse packs can simplify code:

                                      const IndexRange &jb, const IndexRange &kb) {
  namespace p = fluid_prim;
  namespace c = fluid_cons;
  namespace impl = internal_variables;
  using parthenon::MakePackDescriptor;

  auto *pmb = rc->GetParentPointer();
  auto &resolved_pkgs = pmb->resolved_packages;

  static auto desc =
      MakePackDescriptor<p::density, c::density, p::velocity, c::momentum, p::energy,
                         c::energy, p::bfield, c::bfield, p::ye, c::ye, p::pressure,
                         p::gamma1, impl::cell_signal_speed>(resolved_pkgs.get());
  auto v = desc.GetPack(rc);

  // We need these to check whether or not these variables are present
  // in the pack. They are -1 if not present.
  const int pb_hi = v.GetUpperBoundHost(0, p::bfield());
  const int pye = v.GetUpperBoundHost(0, p::ye()); 

  auto geom = Geometry::GetCoordinateSystem(rc);
  const int nblocks = v.GetNBlocks();

  parthenon::par_for(
      DEFAULT_LOOP_PATTERN, "PrimToCons", DevExecSpace(), 0, nblocks - 1, kb.s, kb.e,
      jb.s, jb.e, ib.s, ib.e,
      KOKKOS_LAMBDA(const int b, const int k, const int j, const int i) {
        auto &coords = v.GetCoordinates(b);
        Real gcov4[4][4];
        geom.SpacetimeMetric(CellLocation::Cent, k, j, i, gcov4);
        Real gcon[3][3];
        geom.MetricInverse(CellLocation::Cent, k, j, i, gcon);
        Real gdet = geom.DetGamma(CellLocation::Cent, k, j, i);
        Real lapse = geom.Lapse(CellLocation::Cent, k, j, i);
        Real shift[3];
        geom.ContravariantShift(CellLocation::Cent, k, j, i, shift);

        Real S[3];
        const Real vel[] = {v(b, p::velocity(0), k, j, i), v(b, p::velocity(1), k, j, i),
                            v(b, p::velocity(2), k, j, i)};
        Real bcons[3];
        Real bp[3] = {0.0, 0.0, 0.0};
        if (pb_hi > 0) {
          for (int d = 0; d < 3; ++d) {
            bp[d] = v(b, p::bfield(d), k, j, i);
          }
        }
        Real ye_cons;
        Real ye_prim = 0.0;
        if (pye > 0) {
          ye_prim = v(b, p::ye(), k, j, i);
        }

        Real sig[3];
        prim2con::p2c(v(b, p::density(), k, j, i), vel, bp, v(b, p::energy(), k, j, i),
                      ye_prim, v(b, p::pressure(), k, j, i), v(b, p::gamma1(), k, j, i),
                      gcov4, gcon, shift, lapse, gdet, v(b, c::density(), k, j, i), S,
                      bcons, v(b, c::energy(), k, j, i), ye_cons, sig);

        for (int d = 0; d < 3; ++d) {
          v(b, c::momentum(d), k, j, i) = S[d];
        }

        if (pb_hi > 0) {
          for(int d = 0; d < 3; ++d) {
            v(b, c::bfield(d), k, j, i) = bcons[d];
          }
        }

        if (pye > 0) {
          v(b, c::ye(), k, j, i) = ye_cons;
        }

        for (int d = 0; d < 3; d++) {
          v(b, impl::cell_signal_speed(d), k, j, i) = sig[d];
        }
      });

  return TaskStatus::complete;
}

Copy link
Collaborator

@brryan brryan left a comment

Choose a reason for hiding this comment

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

A huge amount of changes, a great improvement to the code though, thanks for doing this, I don't see any problems, approved.

@Yurlungur Yurlungur merged commit 0aaac0c into main Dec 11, 2023
4 checks passed
@Yurlungur Yurlungur deleted the jmm/go-to-new-packs branch December 11, 2023 18:40
@Yurlungur Yurlungur mentioned this pull request Dec 12, 2023
3 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants