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

WIP: Uniform Cylindrical and Spherical Coordinates #914

Open
wants to merge 10 commits into
base: develop
Choose a base branch
from
2 changes: 1 addition & 1 deletion src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ 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
set(COORDINATE_TYPE "UniformCartesian" CACHE STRING "UniformCartesian/UniformSpherical/UniformCylindrical")

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

Expand Down
2 changes: 2 additions & 0 deletions src/coordinates/coordinates.hpp
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
47 changes: 47 additions & 0 deletions src/coordinates/uniform_cartesian.hpp
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 @@ -277,13 +294,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