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.

Sign up for GitHub

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

Jump to bottom

WIP: Another stab at more general coordinate support #1137

Open
wants to merge 9 commits into
base: develop
Choose a base branch
from
Open

WIP: Another stab at more general coordinate support #1137

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
f966020
Initial refactor and addition of spherical. Compiles at least.
jdolence Jul 10, 2024
ad63ab7
small update to make COORDINATE_TYPE settable
jdolence Jul 10, 2024
d3045bf
fix bug and add using statements to pick up overloads in base
jdolence Jul 10, 2024
eeaf8ee
CellWidth
jdolence Jul 10, 2024
0175b54
Scale
jdolence Jul 10, 2024
2250dc5
add accessors for i and j in index split
jdolence Jul 30, 2024
f231157
missing const
jdolence Jul 30, 2024
17fdc73
another missing const
jdolence Jul 30, 2024
5199495
typos
jdolence Jul 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions example/poisson/poisson_package.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,9 +160,9 @@ TaskStatus SumMass(T *u, Real *reduce_sum) {

const parthenon::Coordinates_t &coords = GetCoords(pm);
const int ndim = v.GetNdim();
const Real dx = coords.Dxc<X1DIR>();
const Real dx = coords.Dxc<X1DIR>(0, 0, 0);
for (int i = X2DIR; i <= ndim; i++) {
const Real dy = coords.DxcFA(i);
const Real dy = coords.Dxc(i, 0, 0, 0);
PARTHENON_REQUIRE_THROWS(dx == dy,
"SumMass requires that DX be equal in all directions.");
}
Expand Down Expand Up @@ -231,9 +231,9 @@ TaskStatus UpdatePhi(T *u, T *du) {

const parthenon::Coordinates_t &coords = GetCoords(pm);
const int ndim = v.GetNdim();
const Real dx = coords.Dxc<X1DIR>();
const Real dx = coords.Dxc<X1DIR>(0, 0, 0);
for (int i = X2DIR; i <= ndim; i++) {
const Real dy = coords.DxcFA(i);
const Real dy = coords.Dxc(i, 0, 0, 0);
PARTHENON_REQUIRE_THROWS(dx == dy,
"UpdatePhi requires that DX be equal in all directions.");
}
Expand Down
7 changes: 6 additions & 1 deletion src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,10 @@ set(COMPILED_WITH ${CMAKE_CXX_COMPILER})
set(COMPILER_COMMAND "<not-implemented>") # TODO: Put something more descriptive here
set(COMPILER_FLAGS "<not-implemented>") # TODO: Put something more descriptive here

set(COORDINATE_TYPE UniformCartesian) # TODO: Make this an option when more are available
if (NOT COORDINATE_TYPE)
set(COORDINATE_TYPE UniformCartesian)
endif()
message(STATUS "COORDINATE_TYPE = " ${COORDINATE_TYPE})

configure_file(config.hpp.in generated/config.hpp @ONLY)

Expand All @@ -115,6 +118,8 @@ add_library(parthenon

coordinates/coordinates.hpp
coordinates/uniform_cartesian.hpp
coordinates/uniform_coordinates.hpp
coordinates/uniform_spherical.hpp

driver/driver.cpp
driver/driver.hpp
Expand Down
1 change: 1 addition & 0 deletions src/coordinates/coordinates.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "config.hpp"

#include "uniform_cartesian.hpp"
#include "uniform_spherical.hpp"

namespace parthenon {

Expand Down
311 changes: 6 additions & 305 deletions src/coordinates/uniform_cartesian.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,319 +13,20 @@
#ifndef COORDINATES_UNIFORM_CARTESIAN_HPP_
#define COORDINATES_UNIFORM_CARTESIAN_HPP_

#include <array>
#include <cassert>

#include "basic_types.hpp"
#include "defs.hpp"
#include "globals.hpp"

#include <Kokkos_Macros.hpp>
#include "uniform_coordinates.hpp"

namespace parthenon {

class UniformCartesian {
class UniformCartesian : public UniformCoordinates<UniformCartesian> {
public:
UniformCartesian() = default;
UniformCartesian(const RegionSize &rs, ParameterInput *pin) {
for (auto &dir : {X1DIR, X2DIR, X3DIR}) {
dx_[dir - 1] = (rs.xmax(dir) - rs.xmin(dir)) / rs.nx(dir);
istart_[dir - 1] = (!rs.symmetry(dir) ? Globals::nghost : 0);
xmin_[dir - 1] = rs.xmin(dir) - istart_[dir - 1] * dx_[dir - 1];
}
area_[0] = dx_[1] * dx_[2];
area_[1] = dx_[0] * dx_[2];
area_[2] = dx_[0] * dx_[1];
cell_volume_ = dx_[0] * dx_[1] * dx_[2];
}
UniformCartesian(const RegionSize &rs, ParameterInput *pin)
: UniformCoordinates<UniformCartesian>(rs, pin) {}
UniformCartesian(const UniformCartesian &src, int coarsen)
: istart_(src.GetStartIndex()) {
dx_ = src.Dx_();
xmin_ = src.GetXmin();
xmin_[0] += istart_[0] * dx_[0] * (1 - coarsen);
xmin_[1] += istart_[1] * dx_[1] * (1 - coarsen);
xmin_[2] += istart_[2] * dx_[2] * (1 - coarsen);
dx_[0] *= coarsen;
dx_[1] *= (istart_[1] > 0 ? coarsen : 1);
dx_[2] *= (istart_[2] > 0 ? coarsen : 1);
area_[0] = dx_[1] * dx_[2];
area_[1] = dx_[0] * dx_[2];
area_[2] = dx_[0] * dx_[1];
cell_volume_ = dx_[0] * dx_[1] * dx_[2];
}

//----------------------------------------
// Dxc: Distance between cell centers
//----------------------------------------
template <int dir>
KOKKOS_FORCEINLINE_FUNCTION Real Dxc() const {
assert(dir > 0 && dir < 4);
return dx_[dir - 1];
}
template <int dir, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real Dxc(Args... args) const {
assert(dir > 0 && dir < 4);
return dx_[dir - 1];
}
template <class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, Args... args) const {
assert(dir > 0 && dir < 4);
return dx_[dir - 1];
}

//----------------------------------------
// Dxf: Distance between cell faces
//----------------------------------------
template <int dir, int face, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const {
assert(dir > 0 && dir < 4 && face > 0 && face < 4);
return dx_[face - 1];
}
template <int dir, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const {
assert(dir > 0 && dir < 4);
return dx_[dir - 1];
}

//----------------------------------------
// Xf: Positions at cell centers
//----------------------------------------
template <int dir>
KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const {
assert(dir > 0 && dir < 4);
return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1];
}
template <int dir>
KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
return Xc<dir>(i);
case 2:
return Xc<dir>(j);
case 3:
return Xc<dir>(k);
default:
PARTHENON_FAIL("Unknown dir");
return 0; // To appease compiler
}
}

//----------------------------------------
// Xf: Positions on Faces
//----------------------------------------
template <int dir, int face>
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const {
assert(dir > 0 && dir < 4 && face > 0 && face < 4);
// Return position in direction "dir" along index "idx" on face "face"
if constexpr (dir == face) {
return xmin_[dir - 1] + idx * dx_[dir - 1];
} else {
return Xc<dir>(idx);
}
}
template <int dir, int face>
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
return Xf<dir, face>(i);
case 2:
return Xf<dir, face>(j);
case 3:
return Xf<dir, face>(k);
default:
PARTHENON_FAIL("Unknown dir");
return 0; // To appease compiler
}
}

template <int dir>
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int idx) const {
assert(dir > 0 && dir < 4);
// Return position in direction "dir" along index "idx" on face "dir"
return xmin_[dir - 1] + idx * dx_[dir - 1];
}
template <int dir>
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
return Xf<dir>(i);
case 2:
return Xf<dir>(j);
case 3:
return Xf<dir>(k);
default:
PARTHENON_FAIL("Unknown dir");
return 0; // To appease compiler
}
}

template <int dir, TopologicalElement el>
KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const {
if constexpr (dir == X1DIR && TopologicalOffsetI(el)) {
return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
} else if constexpr (dir == X2DIR && TopologicalOffsetJ(el)) {
return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
} else if constexpr (dir == X3DIR && TopologicalOffsetK(el)) {
return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2
} else {
return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; // idx
}
return 0; // This should never be reached, but w/o it some compilers generate warnings
}

template <int dir, TopologicalElement el>
KOKKOS_FORCEINLINE_FUNCTION Real X(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
return X<dir, el>(i);
case 2:
return X<dir, el>(j);
case 3:
return X<dir, el>(k);
default:
PARTHENON_FAIL("Unknown dir");
return 0; // To appease compiler
}
}

//----------------------------------------
// CellWidth: Width of cells at cell centers
//----------------------------------------
template <int dir, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(Args... args) const {
assert(dir > 0 && dir < 4);
return dx_[dir - 1];
}
template <class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real CellWidthFA(const int dir, Args... args) const {
assert(dir > 0 && dir < 4);
return dx_[dir - 1];
}
: UniformCoordinates<UniformCartesian>(src, coarsen) {}

//----------------------------------------
// EdgeLength: Length of cell edges
//----------------------------------------
template <int dir, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(Args... args) const {
assert(dir > 0 && dir < 4);
return CellWidth<dir>();
}
template <class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real EdgeLengthFA(const int dir, Args... args) const {
return CellWidthFA(dir);
}

//----------------------------------------
// FaceArea: Area of cell areas
//----------------------------------------
template <int dir, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(Args... args) const {
assert(dir > 0 && dir < 4);
return area_[dir - 1];
}
template <class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real FaceAreaFA(const int dir, Args... args) const {
assert(dir > 0 && dir < 4);
return area_[dir - 1];
}

//----------------------------------------
// CellVolume
//----------------------------------------
template <class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(Args... args) const {
return cell_volume_;
}

//----------------------------------------
// Generalized volume
//----------------------------------------
template <TopologicalElement el, class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real Volume(Args... args) const {
using TE = TopologicalElement;
if constexpr (el == TE::CC) {
return cell_volume_;
} else if constexpr (el == TE::F1) {
return area_[X1DIR - 1];
} else if constexpr (el == TE::F2) {
return area_[X2DIR - 1];
} else if constexpr (el == TE::F3) {
return area_[X3DIR - 1];
} else if constexpr (el == TE::E1) {
return dx_[X1DIR - 1];
} else if constexpr (el == TE::E2) {
return dx_[X2DIR - 1];
} else if constexpr (el == TE::E3) {
return dx_[X3DIR - 1];
} else if constexpr (el == TE::NN) {
return 1.0;
}
PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
"TopologicalElement enum.");
return 0.0;
}

template <class... Args>
KOKKOS_FORCEINLINE_FUNCTION Real Volume(CellLevel cl, TopologicalElement el,
Args... args) const {
using TE = TopologicalElement;
if (cl == CellLevel::same) {
if (el == TE::CC) {
return cell_volume_;
} else if (el == TE::F1) {
return area_[X1DIR - 1];
} else if (el == TE::F2) {
return area_[X2DIR - 1];
} else if (el == TE::F3) {
return area_[X3DIR - 1];
} else if (el == TE::E1) {
return dx_[X1DIR - 1];
} else if (el == TE::E2) {
return dx_[X2DIR - 1];
} else if (el == TE::E3) {
return dx_[X3DIR - 1];
} else if (el == TE::NN) {
return 1.0;
}
} else if (cl == CellLevel::fine) {
if (el == TE::CC) {
return cell_volume_ / 8.0;
} else if (el == TE::F1) {
return area_[X1DIR - 1] / 4.0;
} else if (el == TE::F2) {
return area_[X2DIR - 1] / 4.0;
} else if (el == TE::F3) {
return area_[X3DIR - 1] / 4.0;
} else if (el == TE::E1) {
return dx_[X1DIR - 1] / 2.0;
} else if (el == TE::E2) {
return dx_[X2DIR - 1] / 2.0;
} else if (el == TE::E3) {
return dx_[X3DIR - 1] / 2.0;
} else if (el == TE::NN) {
return 1.0;
}
}
PARTHENON_FAIL("If you reach this point, someone has added a new value to the the "
"TopologicalElement enum.");
return 0.0;
}

const std::array<Real, 3> &GetXmin() const { return xmin_; }
const std::array<int, 3> &GetStartIndex() const { return istart_; }
const char *Name() const { return name_; }

private:
std::array<int, 3> istart_;
std::array<Real, 3> xmin_, dx_, area_;
Real cell_volume_;
constexpr static const char *name_ = "UniformCartesian";

const std::array<Real, 3> &Dx_() const { return dx_; }
private:
};

} // namespace parthenon
Expand Down
Loading
Loading