Skip to content

Commit

Permalink
Output metric in embedding/KS coordinates as well as native
Browse files Browse the repository at this point in the history
  • Loading branch information
Ben Prather committed Dec 7, 2023
1 parent ce23cbb commit 946562d
Show file tree
Hide file tree
Showing 4 changed files with 55 additions and 28 deletions.
37 changes: 26 additions & 11 deletions kharma/coord_output/coord_output.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,23 +40,18 @@ std::shared_ptr<KHARMAPackage> CoordinateOutput::Initialize(ParameterInput *pin,
auto pkg = std::make_shared<KHARMAPackage>("CoordinateOutput");
Params &params = pkg->AllParams();

// Options
//Real n = pin->GetOrAddReal("wind", "ne", 2.e-4);
//params.Add("ne", n);
// Any options? Which fields to output is determined in outputs

// Fields: cell-center values for geometry only
// TODO test faces when available, optional lists of locations?
Metadata::AddUserFlag("Geometry");
// TODO add values elsewhere e.g. faces, edges?
// Note these are all marked with the "CoordinateOutput" flag
std::vector<int> s_4vector({GR_DIM});
std::vector<int> s_4tensor({GR_DIM, GR_DIM});
std::vector<int> s_4conn({GR_DIM, GR_DIM, GR_DIM});
std::vector<MetadataFlag> flags_geom = {Metadata::Real, Metadata::Cell, Metadata::Derived,
Metadata::OneCopy, Metadata::GetUserFlag("Geometry")};
std::vector<MetadataFlag> flags_geom_face = {Metadata::Real, Metadata::Face, Metadata::Derived,
Metadata::OneCopy, Metadata::GetUserFlag("Geometry")};
Metadata::OneCopy};
auto m0 = Metadata(flags_geom);
auto m1 = Metadata(flags_geom, s_4vector);
auto m1f = Metadata(flags_geom_face, s_4vector);
auto m2 = Metadata(flags_geom, s_4tensor);
auto m3 = Metadata(flags_geom, s_4conn);

Expand All @@ -83,6 +78,14 @@ std::shared_ptr<KHARMAPackage> CoordinateOutput::Initialize(ParameterInput *pin,
pkg->AddField("coords.lapse", m0);
pkg->AddField("coords.conn", m3);

// Metric, embedding (i.e., KS) coordinates
pkg->AddField("coords.gcon_embed", m2);
pkg->AddField("coords.gcov_embed", m2);
pkg->AddField("coords.gdet_embed", m0);
// TODO?
//pkg->AddField("coords.lapse_embed", m0);
//pkg->AddField("coords.conn_embed", m3);

// Register our output. This will be called before *any* output,
// but we will only fill the fields before the first.
// This is all that's needed unless:
Expand All @@ -100,7 +103,7 @@ TaskStatus CoordinateOutput::BlockUserWorkBeforeOutput(MeshBlock *pmb, Parameter
auto rc = pmb->meshblock_data.Get();

PackIndexMap geom_map;
auto Geom = rc->PackVariables({Metadata::GetUserFlag("Geometry")}, geom_map);
auto Geom = rc->PackVariables({Metadata::GetUserFlag("CoordinateOutput")}, geom_map);

const auto& G = pmb->coords;

Expand All @@ -125,6 +128,10 @@ TaskStatus CoordinateOutput::BlockUserWorkBeforeOutput(MeshBlock *pmb, Parameter
const int mlapse = geom_map["coords.lapse"].first;
const int mconn = geom_map["coords.conn"].first;

const int mgcov_embed = geom_map["coords.gcov_embed"].first;
const int mgcon_embed = geom_map["coords.gcon_embed"].first;
const int mgdet_embed = geom_map["coords.gdet_embed"].first;

IndexRange3 b = KDomain::GetRange(rc, IndexDomain::entire);
pmb->par_for("set_geometry", b.ks, b.ke, b.js, b.je, b.is, b.ie,
KOKKOS_LAMBDA (const int &k, const int &j, const int &i) {
Expand All @@ -143,7 +150,7 @@ TaskStatus CoordinateOutput::BlockUserWorkBeforeOutput(MeshBlock *pmb, Parameter
Geom(mXsph+2, k, j, i) = Geom(mth, k, j, i) = G.th(k, j, i);
Geom(mXsph+3, k, j, i) = Geom(mphi, k, j, i) = G.phi(k, j, i);

// Metric
// Metric, native
DLOOP2 Geom(mgcov+GR_DIM*mu+nu, k, j, i) = G.gcov(Loci::center, j, i, mu, nu);
DLOOP2 Geom(mgcon+GR_DIM*mu+nu, k, j, i) = G.gcon(Loci::center, j, i, mu, nu);
Geom(mgdet, k, j, i) = G.gdet(Loci::center, j, i);
Expand All @@ -152,6 +159,14 @@ TaskStatus CoordinateOutput::BlockUserWorkBeforeOutput(MeshBlock *pmb, Parameter
// Connection
DLOOP3 Geom(mconn+GR_DIM*GR_DIM*mu+GR_DIM*nu+lam, k, j, i) = G.conn(j, i, mu, nu, lam);

// Metric, embedding
GReal Xembed[GR_DIM], gcov_embed[GR_DIM][GR_DIM], gcon_embed[GR_DIM][GR_DIM];
G.coord_embed(k, j, i, Loci::center, Xembed);
G.coords.gcov_embed(Xembed, gcov_embed);
GReal gdet = G.coords.gcon_embed(Xembed, gcon_embed); // Save a tiny bit of time
DLOOP2 Geom(mgcov_embed+GR_DIM*mu+nu, k, j, i) = gcov_embed[mu][nu];
DLOOP2 Geom(mgcon_embed+GR_DIM*mu+nu, k, j, i) = gcon_embed[mu][nu];
Geom(mgdet_embed, k, j, i) = gdet;
}
);

Expand Down
35 changes: 23 additions & 12 deletions kharma/coordinates/coordinate_embedding.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -271,15 +271,32 @@ class CoordinateEmbedding {
return mpark::holds_alternative<CartMinkowskiCoords>(base) && mpark::holds_alternative<NullTransform>(transform);
}

// Spell out the interface we take from BaseCoords
// TODO add a gcon_embed, gdet_embed
// Note this is the one thing we need from BaseCoords
KOKKOS_INLINE_FUNCTION void gcov_embed(const GReal Xembed[GR_DIM], Real gcov[GR_DIM][GR_DIM]) const
{
mpark::visit( [&Xembed, &gcov](const auto& self) {
self.gcov_embed(Xembed, gcov);
}, base);
}
// and from the Transform
// All the quantities we can derive from that
KOKKOS_INLINE_FUNCTION Real gcon_from_gcov(const Real gcov[GR_DIM][GR_DIM], Real gcon[GR_DIM][GR_DIM]) const
{
Real gdet = invert(&gcov[0][0], &gcon[0][0]);
return m::sqrt(m::abs(gdet));
}
KOKKOS_INLINE_FUNCTION Real gcon_embed(const GReal Xembed[GR_DIM], Real gcon[GR_DIM][GR_DIM]) const
{
GReal gcov[GR_DIM][GR_DIM];
gcov_embed(Xembed, gcov);
return gcon_from_gcov(gcov, gcon);
}
KOKKOS_INLINE_FUNCTION Real gdet_embed(const GReal Xembed[GR_DIM]) const
{
GReal gcon[GR_DIM][GR_DIM];
return gcon_embed(Xembed, gcon);
}

// Now, everything we take from CoordinateTransform
KOKKOS_INLINE_FUNCTION void coord_to_embed(const GReal Xnative[GR_DIM], GReal Xembed[GR_DIM]) const
{
mpark::visit( [&Xnative, &Xembed](const auto& self) {
Expand Down Expand Up @@ -484,18 +501,12 @@ class CoordinateEmbedding {
{
Real gcov[GR_DIM][GR_DIM];
gcov_native(X, gcov);
return gcon_native(gcov, gcon);
}
KOKKOS_INLINE_FUNCTION Real gcon_native(const Real gcov[GR_DIM][GR_DIM], Real gcon[GR_DIM][GR_DIM]) const
{
Real gdet = invert(&gcov[0][0], &gcon[0][0]);
return m::sqrt(m::abs(gdet));
return gcon_from_gcov(gcov, gcon);
}
KOKKOS_INLINE_FUNCTION Real gdet_native(const GReal X[GR_DIM]) const
{
Real gcov[GR_DIM][GR_DIM], gcon[GR_DIM][GR_DIM];
gcov_native(X, gcov);
return gcon_native(gcov, gcon);
Real gcon[GR_DIM][GR_DIM];
return gcon_native(X, gcon);
}

KOKKOS_INLINE_FUNCTION void conn_native(const GReal X[GR_DIM], const GReal delta, Real conn[GR_DIM][GR_DIM][GR_DIM]) const
Expand Down
6 changes: 3 additions & 3 deletions kharma/coordinates/gr_coordinates.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ void init_GRCoordinates(GRCoordinates& G) {
// Get geometry at points
GReal gcov_loc[GR_DIM][GR_DIM], gcon_loc[GR_DIM][GR_DIM];
G.coords.gcov_native(X, gcov_loc);
const GReal gdet = G.coords.gcon_native(gcov_loc, gcon_loc);
const GReal gdet = G.coords.gcon_from_gcov(gcov_loc, gcon_loc);
// Add to running averages
gdet_local(loc, j, i) += gdet / square;
DLOOP2 {
Expand Down Expand Up @@ -177,7 +177,7 @@ void init_GRCoordinates(GRCoordinates& G) {
// Get geometry at the point
GReal gcov_loc[GR_DIM][GR_DIM], gcon_loc[GR_DIM][GR_DIM];
G.coords.gcov_native(X, gcov_loc);
const GReal gdet = G.coords.gcon_native(gcov_loc, gcon_loc);
const GReal gdet = G.coords.gcon_from_gcov(gcov_loc, gcon_loc);
// Add to running averages
gdet_local(loc, j, i) += gdet / diameter;
DLOOP2 {
Expand All @@ -192,7 +192,7 @@ void init_GRCoordinates(GRCoordinates& G) {
// Get geometry
GReal gcov_loc[GR_DIM][GR_DIM], gcon_loc[GR_DIM][GR_DIM];
G.coords.gcov_native(X, gcov_loc);
const GReal gdet = G.coords.gcon_native(gcov_loc, gcon_loc);
const GReal gdet = G.coords.gcon_from_gcov(gcov_loc, gcon_loc);
// Set geometry
gdet_local(loc, j, i) = gdet;
DLOOP2 {
Expand Down
5 changes: 3 additions & 2 deletions pars/tori_2d/sane2d.par
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,10 +98,11 @@ dt = 0.1
# a bunch of other contexts.
<parthenon/output3>
file_type = hdf5
analysis_output = true
dt = 1e20
single_precision_output = false
# If you want the convenience and hate disk space, you can add these
# variables to "normal" dump files like output0 too.
variables = coords.Xnative, coords.Xsph, &
coords.gcon, coords.gcov, coords.gdet, coords.lapse, &
coords.conn
coords.gcon, coords.gcov, coords.gdet, coords.lapse, coords.conn, &
coords.gcon_embed, coords.gcov_embed, coords.gdet_embed

0 comments on commit 946562d

Please sign in to comment.