lanl / phoebus

Phifty One Ergs Blows Up A Star
BSD 3-Clause "New" or "Revised" License
32 stars 0 forks source link

Update Parthenon and move to type-based variables #191

Closed Yurlungur closed 9 months ago

Yurlungur commented 9 months ago

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

Yurlungur commented 9 months ago

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;
}