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

WIP: Uniform Cylindrical and Spherical Coordinates #914

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

WIP: Uniform Cylindrical and Spherical Coordinates #914

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
b3d547c
WIP: Uniform cylindrical and spherical added
Jul 28, 2023
c32f395
Add Xs<dir,side>, compilation fixes
Jul 31, 2023
0172270
Turned cached arrays into member functions
Aug 7, 2023
5b6bb61
Always use chuncking when compression is supported even not used
pgrete Aug 23, 2023
0cc1c6c
Small cmake change
Aug 23, 2023
4a01ed5
Added src term helpers for cylindrical
Oct 30, 2023
2cf79a7
WIP:Spherical coordinates
Nov 5, 2023
5bee38e
Merge commit '5b6bb619' into forrestglines/curvilinear
Nov 5, 2023
59d4054
Fixed spherical coords
Nov 7, 2023
f88adb9
GetXmin, GetStartIndex Inlined
Feb 11, 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
9 changes: 8 additions & 1 deletion example/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,19 @@
# the public, perform publicly and display publicly, and to permit others to do so.
#=========================================================================================


add_subdirectory(stochastic_subgrid)
add_subdirectory(advection)
add_subdirectory(calculate_pi)
add_subdirectory(kokkos_pi)
add_subdirectory(particles)
add_subdirectory(particle_leapfrog)
add_subdirectory(particle_tracers)
add_subdirectory(poisson)
add_subdirectory(sparse_advection)


# Not all tests will work with all coordinate systems. Add Cartesian-only tests
# below
if (${COORDINATE_TYPE} STREQUAL "UniformCartesian")
add_subdirectory(poisson)
endif()
5 changes: 3 additions & 2 deletions src/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -85,13 +85,14 @@ else()
set(PAR_LOOP_INNER_LAYOUT_TAG loop_pattern_undefined_tag)
endif()

set(COORDINATE_TYPE "UniformCartesian" CACHE STRING "UniformCartesian/UniformSpherical/UniformCylindrical")
message(STATUS "COORDINATE_TYPE='${COORDINATE_TYPE}' Parthenon Coordinates_t class")

set(EXCEPTION_HANDLING_OPTION ENABLE_EXCEPTIONS) # TODO: Add option to disable exceptions
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

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

add_library(parthenon
Expand Down
2 changes: 2 additions & 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,8 @@
#include "config.hpp"

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

namespace parthenon {

Expand Down
80 changes: 71 additions & 9 deletions src/coordinates/uniform_cartesian.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include "basic_types.hpp"
#include "defs.hpp"
#include "globals.hpp"
#include "parameter_input.hpp"

#include <Kokkos_Macros.hpp>

Expand All @@ -41,6 +42,21 @@ class UniformCartesian {
xmin_[0] = rs.x1min - istart_[0] * dx_[0];
xmin_[1] = rs.x2min - istart_[1] * dx_[1];
xmin_[2] = rs.x3min - istart_[2] * dx_[2];
ndim_ = (rs.nx1 != 1) + (rs.nx2 != 1) + (rs.nx3 != 1);

std::string coord_type_str = pin->GetOrAddString("parthenon/mesh", "coord", "cartesian");
//REMOVE ME
//if (coord_type_str == "cartesian") {
// coord_type = parthenon::uniformOrthMesh::cartesian;
//} else
if (coord_type_str == "cylindrical") {
PARTHENON_THROW(" Please rebuild with -DCOORDINATE_TYPE=UniformCylindrical");
} else if (coord_type_str == "spherical_polar") {
PARTHENON_THROW(" Please rebuild with -DCOORDINATE_TYPE=UniformSpherical");
} else {
PARTHENON_THROW("Invalid coord input in <parthenon/mesh>.");
}
Comment on lines +47 to +58
Copy link
Collaborator

Choose a reason for hiding this comment

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

if coordinate systems are going to compile-time only, maybe we drop coord from the input deck entirely?


}
UniformCartesian(const UniformCartesian &src, int coarsen)
: istart_(src.GetStartIndex()) {
Expand All @@ -56,6 +72,7 @@ class UniformCartesian {
area_[1] = dx_[0] * dx_[2];
area_[2] = dx_[0] * dx_[1];
cell_volume_ = dx_[0] * dx_[1] * dx_[2];
ndim_ = src.ndim_;
}

//----------------------------------------
Expand Down Expand Up @@ -103,11 +120,11 @@ class UniformCartesian {
KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
case X1DIR:
return Xc<dir>(i);
case 2:
case X2DIR:
return Xc<dir>(j);
case 3:
case X3DIR:
return Xc<dir>(k);
default:
PARTHENON_FAIL("Unknown dir");
Expand All @@ -132,11 +149,11 @@ class UniformCartesian {
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
case X1DIR:
return Xf<dir, face>(i);
case 2:
case X2DIR:
return Xf<dir, face>(j);
case 3:
case X3DIR:
return Xf<dir, face>(k);
default:
PARTHENON_FAIL("Unknown dir");
Expand All @@ -154,17 +171,32 @@ class UniformCartesian {
KOKKOS_FORCEINLINE_FUNCTION Real Xf(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4);
switch (dir) {
case 1:
case X1DIR:
return Xf<dir>(i);
case 2:
case X2DIR:
return Xf<dir>(j);
case 3:
case X3DIR:
return Xf<dir>(k);
default:
PARTHENON_FAIL("Unknown dir");
return 0; // To appease compiler
}
}

//----------------------------------------
// Xs: Area averaged positions
//----------------------------------------
template <int dir, int side>
KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int idx) const {
assert(dir > 0 && dir < 4 && side > 0 && side < 4);
return Xc<dir>(idx);
}
template <int dir, int side>
KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int k, const int j, const int i) const {
assert(dir > 0 && dir < 4 && side > 0 && side < 4);
return Xc<dir>(k,j,i);
}
Comment on lines +186 to +198
Copy link
Collaborator

Choose a reason for hiding this comment

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

not sure I like the name here.



template <int dir, TopologicalElement el>
KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const {
Expand Down Expand Up @@ -277,13 +309,43 @@ class UniformCartesian {
return 0.0;
}

//----------------------------------------
// Cylindrical+Spherical Conversions
//----------------------------------------
// TODO
Comment on lines +312 to +315
Copy link
Collaborator

Choose a reason for hiding this comment

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

Hmm... is this actually something we need/want the Parthenon coords tooling to provide? My recommendation is one of the following:

  1. All coords provide functionality to convert between themselves and Cartesian. For UniformCartesian, this is an identity transformation
  2. Just make this not parthenon's problem.

Probably but not necessarily we want to go with option 1, as that's what's needed for vis tooling and initial conditions in some contexts.

Copy link
Collaborator

Choose a reason for hiding this comment

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

longer discussion below. Resolving this comment

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

From what I've done and what Shengtai included in his demo, there's a bigger use case to provide spherical and cylindrical conversions in Parthenon since we're likely to deal with physical systems where those coordinates make sense (i.e. galaxy clusters and stars for spherical setups, jets and disks for cylindrical setups). We could leave it up to each of these systems to implement these conversions on their own but I see them repeated so many times elsewhere it might be better to centralize that code in one place.

I'm not set on including these conversions in Parthenon at the moment but still thinking about it. We might also include conversions for vectors.

Copy link
Collaborator

Choose a reason for hiding this comment

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

I see. No I think that makes sense. As I comment below, we may also want it for vis tools. It would be better if, for example, Paraview could automatically visualize simulations done in spherical coordinates. So we might want to standardize some functions that always map from whatever coordinate system we're in to Cartesian.


//----------------------------------------
// h*v, h*f, dh*vd*, etc.
// These might only be needed for viscocity and conduction and so might exist downstream
//----------------------------------------
// TODO

//----------------------------------------
// Position to Cell index
//----------------------------------------
KOKKOS_INLINE_FUNCTION
void Xtoijk(const Real &x, const Real &y, const Real &z, int &i, int &j, int &k) const {
i = static_cast<int>(
std::floor((x - xmin_[0]) / dx_[0])) +
istart_[0];
j = (ndim_ > 1) ? static_cast<int>(std::floor(
(y - xmin_[1]) / dx_[1])) +
istart_[1]
: istart_[1];
k = (ndim_ > 2) ? static_cast<int>(std::floor(
(z - xmin_[2]) / dx_[2])) +
istart_[2]
: istart_[2];
}
Comment on lines +326 to +339
Copy link
Collaborator

Choose a reason for hiding this comment

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

I have not had good success making the ternary operator vectorize. Maybe we either switch this to using masks or we split it into 3 functions?

Copy link
Collaborator

Choose a reason for hiding this comment

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

Also this appears to be for cell centers only?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This implementation was lifted from an inline function for swarms, I copied it nearly verbatim. I should change the function description for "cell centered index". I think this is sufficient for

My very outdated experience on Kepler GPUs was that ternary operators had less branch splitting, so I also default to using ternary operators.

Copy link
Collaborator

Choose a reason for hiding this comment

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

This implementation was lifted from an inline function for swarms, I copied it nearly verbatim. I should change the function description for "cell centered index". I think this is sufficient for

👍

My very outdated experience on Kepler GPUs was that ternary operators had less branch splitting, so I also default to using ternary operators.

It's not a huge deal either way, but I was thinking of CPUs when I made the comment, not GPUs. My experience is that GPUs are much better at handling ternary/branching than CPUs.


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_;
int ndim_;
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
Real cell_volume_;
constexpr static const char *name_ = "UniformCartesian";

Expand Down
Loading
Loading