From b3d547c790682149607e6e6483f9d54038b7d5bf Mon Sep 17 00:00:00 2001 From: Forrest Glines Date: Fri, 28 Jul 2023 17:44:21 -0600 Subject: [PATCH 1/9] WIP: Uniform cylindrical and spherical added --- src/CMakeLists.txt | 2 +- src/coordinates/coordinates.hpp | 2 + src/coordinates/uniform_cartesian.hpp | 47 ++ src/coordinates/uniform_cylindrical.hpp | 471 ++++++++++++++++++++ src/coordinates/uniform_spherical.hpp | 561 ++++++++++++++++++++++++ src/interface/swarm_device_context.hpp | 12 +- 6 files changed, 1083 insertions(+), 12 deletions(-) create mode 100644 src/coordinates/uniform_cylindrical.hpp create mode 100644 src/coordinates/uniform_spherical.hpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1a3dcb50be83..1258809e7873 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -90,7 +90,7 @@ set(COMPILED_WITH ${CMAKE_CXX_COMPILER}) set(COMPILER_COMMAND "") # TODO: Put something more descriptive here set(COMPILER_FLAGS "") # 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) diff --git a/src/coordinates/coordinates.hpp b/src/coordinates/coordinates.hpp index d1290deea425..83af1cd5de64 100644 --- a/src/coordinates/coordinates.hpp +++ b/src/coordinates/coordinates.hpp @@ -16,6 +16,8 @@ #include "config.hpp" #include "uniform_cartesian.hpp" +#include "uniform_cylindrical.hpp" +#include "uniform_spherical.hpp" namespace parthenon { diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp index 05f8dd9a6f55..af3c14e76446 100644 --- a/src/coordinates/uniform_cartesian.hpp +++ b/src/coordinates/uniform_cartesian.hpp @@ -19,6 +19,7 @@ #include "basic_types.hpp" #include "defs.hpp" #include "globals.hpp" +#include "parameter_input.hpp" #include @@ -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 ."); + } + } UniformCartesian(const UniformCartesian &src, int coarsen) : istart_(src.GetStartIndex()) { @@ -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_; } //---------------------------------------- @@ -277,6 +294,35 @@ class UniformCartesian { return 0.0; } + //---------------------------------------- + // Cylindrical+Spherical Conversions + //---------------------------------------- + // TODO + + //---------------------------------------- + // 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( + std::floor((x - xmin_[0]) / dx_[0])) + + istart_[0]; + j = (ndim_ > 1) ? static_cast(std::floor( + (y - xmin_[1]) / dx_[1])) + + istart_[1] + : istart_[1]; + k = (ndim_ > 2) ? static_cast(std::floor( + (z - xmin_[2]) / dx_[2])) + + istart_[2] + : istart_[2]; + } + const std::array &GetXmin() const { return xmin_; } const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return name_; } @@ -284,6 +330,7 @@ class UniformCartesian { private: std::array istart_; std::array xmin_, dx_, area_; + int ndim_; Real cell_volume_; constexpr static const char *name_ = "UniformCartesian"; diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp new file mode 100644 index 000000000000..e69d129a6de7 --- /dev/null +++ b/src/coordinates/uniform_cylindrical.hpp @@ -0,0 +1,471 @@ +//======================================================================================== +// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== +#ifndef COORDINATES_UNIFORM_CYLINDRICAL_HPP_ +#define COORDINATES_UNIFORM_CYLINDRICAL_HPP_ + +#include +#include + +#include "basic_types.hpp" +#include "defs.hpp" +#include "globals.hpp" +#include "parameter_input.hpp" + +#include + +namespace parthenon { + +class UniformCylindrical { + public: + UniformCylindrical() = default; + UniformCylindrical(const RegionSize &rs, ParameterInput *pin) { + dx_[0] = (rs.x1max - rs.x1min) / rs.nx1; + dx_[1] = (rs.x2max - rs.x2min) / rs.nx2; + dx_[2] = (rs.x3max - rs.x3min) / rs.nx3; + area_[0] = dx_[1] * dx_[2]; + area_[1] = dx_[0] * dx_[2]; + area_[2] = dx_[0] * dx_[1]; + istart_[0] = Globals::nghost; + istart_[1] = (rs.nx2 > 1 ? Globals::nghost : 0); + istart_[2] = (rs.nx3 > 1 ? Globals::nghost : 0); + xmin_[0] = rs.x1min - istart_[0] * dx_[0]; + xmin_[1] = rs.x2min - istart_[1] * dx_[1]; + xmin_[2] = rs.x3min - istart_[2] * dx_[2]; + nx_[0] = rs.nx1; + nx_[1] = rs.nx2; + nx_[2] = rs.nx3; + + std::string coord_type_str = pin->GetOrAddString("parthenon/mesh", "coord", "cylindrical"); + if( coord_type_str != "cylindrical" ){ + if (coord_type_str == "cartesian") { + PARTHENON_THROW(" Please rebuild with -DCOORDINATE_TYPE=UniformCartesian"); + } else if (coord_type_str == "spherical") { + PARTHENON_THROW(" Please rebuild with -DCOORDINATE_TYPE=UniformSpherical"); + } else { + PARTHENON_THROW("Invalid coord input in ."); + } + } + + //Initialized cached coordinates + InitCachedArrays(); + } + UniformCylindrical(const UniformCylindrical &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]; + nx_[0] = src.nx_[0] / ( coarsen ? 2 : 1 ); + nx_[1] = src.nx_[1] / ( coarsen ? 2 : 1 ); + nx_[2] = src.nx_[2] / ( coarsen ? 2 : 1 ); + + //Initialized cached coordinates + InitCachedArrays(); + } + + //---------------------------------------- + // Dxc: Distance between cell centers + // Dxc<1> returns distance between cell centroids in r + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + if( dir == X1DIR ){ + return r_c_(i+1) - r_c_(i); + } else { + return dx_[dir-1]; + } + } + KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch (dir) { + case 1: + return Dxc<1>(k,j,i); + case 2: + return Dxc<2>(k,j,j); + case 3: + return Dxc<3>(k,j,k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // Dxf: Distance between cell faces + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const { + assert(dir > 0 && dir < 4 && face > 0 && face < 4); + return dx_[face - 1]; + } + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const { + assert(dir > 0 && dir < 4); + return dx_[dir - 1]; + } + + //---------------------------------------- + // Xc: Positions at cell centers + // Xc<1> returns r of cell centroid + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { + assert(dir > 0 && dir < 4); + if( dir == X1DIR ){ + return r_c_(idx); + } else { + return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; + } + } + template + 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(i); + case 2: + return Xc(j); + case 3: + return Xc(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // Xf: Positions on Faces + //---------------------------------------- + template + 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(idx); + } + } + template + 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(i); + case 2: + return Xf(j); + case 3: + return Xf(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + template + 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 + 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(i); + case 2: + return Xf(j); + case 3: + return Xf(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { + using TE = TopologicalElement; + bool constexpr X1EDGE = el == TE::F1 || el == TE::E2 || el == TE::E3 || el == TE::NN; + bool constexpr X2EDGE = el == TE::F2 || el == TE::E3 || el == TE::E1 || el == TE::NN; + bool constexpr X3EDGE = el == TE::F3 || el == TE::E1 || el == TE::E2 || el == TE::NN; + if constexpr (dir == X1DIR && X1EDGE) { + return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 + } else if constexpr (dir == X2DIR && X2EDGE) { + return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 + } else if constexpr (dir == X3DIR && X3EDGE) { + 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 + 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(i); + case 2: + return X(j); + case 3: + return X(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // CellWidth: Physical width of cells at cell centers + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + if( dir == X2DIR ){ + return Xf<1>(k, j, i)*dx_[1]; //r*dphi + } else { + return dx_[dir-1]; + } + } + KOKKOS_FORCEINLINE_FUNCTION Real CellWidthFA(const int dir,const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + if( dir == X2DIR ){ + return Xf<1>(k, j, i)*dx_[1]; //r*dphi + } else { + return dx_[dir-1]; + } + } + + //---------------------------------------- + // EdgeLength: Physical length of cell edges + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + if( dir == X2DIR ){ + return Xf<1>(k, j, i)*dx_[1]; //r*dphi + } else { + return dx_[dir-1]; + } + } + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLengthFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + if( dir == X2DIR ){ + return Xf<1>(k, j, i)*dx_[1]; //r*dphi + } else { + return dx_[dir-1]; + } + } + + //---------------------------------------- + // FaceArea: Physical area of cell areas + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir) { + case X1DIR: + return X(i)*area_[0]; //r*dphi*dz + case X2DIR: + return area_[1]; //dr*dz + case X3DIR: + return coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi + } + + } + template + KOKKOS_FORCEINLINE_FUNCTION Real FaceAreaFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir) { + case X1DIR: + return X<1, TE::F1>(i)*area_[0]; //r*dphi*dz + case X2DIR: + return area_[1]; //dr*dz + case X3DIR: + return coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi + } + } + + //---------------------------------------- + // CellVolume + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { + return coord_vol_i_(i)*area_[0]; + } + + //---------------------------------------- + // Generalized physical volume + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Volume(Args... args) const { + using TE = TopologicalElement; + if constexpr (el == TE::CC) { + return CellVolume(args...); + } else if constexpr (el == TE::F1) { + return FaceArea(args...); + } else if constexpr (el == TE::F2) { + return FaceArea(args...); + } else if constexpr (el == TE::F3) { + return FaceArea(args...); + } else if constexpr (el == TE::E1) { + return EdgeLength(args...); + } else if constexpr (el == TE::E2) { + return EdgeLength(args...); + } else if constexpr (el == TE::E3) { + return EdgeLength(args...); + } 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; + } + + //---------------------------------------- + // Cylindrical+Spherical Conversions + //---------------------------------------- + // TODO + + //---------------------------------------- + // Geometric Terms (Find better names for these!) + //---------------------------------------- + KOKKOS_FORCEINLINE_FUNCTION + Real h2v(const int i) const { return r_c_(i); }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h2f(const int i) const { return Xf<1>(i); }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h31v(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h31f(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh2vd1(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh2fd1(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh31vd1(const int i) const { return 0.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh31fd1(const int i) const { return 0.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h32v(const int j) const { return 1.0; } + + KOKKOS_FORCEINLINE_FUNCTION + Real h32f(const int j) const { return 1.0; } + + KOKKOS_FORCEINLINE_FUNCTION + Real dh32vd2(const int j) const { return 0.0; } + + KOKKOS_FORCEINLINE_FUNCTION + Real dh32fd2(const int j) const { return 0.0; } + + //What is this? + KOKKOS_FORCEINLINE_FUNCTION + void GetAcc(const int k, const int j, const int i, const Real &rp, + const Real &cosphip, const Real &sinphip, const Real &zp, + Real &acc1, Real &acc2, Real &acc3) const + { + const Real cosdphi = cosphi_c_(j)*cosphip - sinphi_c_(j)*sinphip; //cos(x3v(k)-phip) + const Real sindphi = sinphi_c_(j)*cosphip - cosphi_c_(j)*sinphip; //sin(x3v(k)-phip) + //Psi = - 1.0/sqrt(r0^2 + rp^2 - 2r0*rp*cos(phi-phip) + (z-zp)^2) + acc1 = (r_c_(i) - rp*cosdphi); //acc1 = dPsi/dr, Psi=-1/dist2 + acc2 = rp*sindphi; //acc2 = dPsi/(r*dphi) + acc3 = Xc<3>(k)-zp; + } + + //---------------------------------------- + // Position to Cell index + //---------------------------------------- + KOKKOS_INLINE_FUNCTION + void Xtoijk(const Real &r, const Real &theta, const Real &z, int &i, int &j, int &k) const { + i = (nx_[0] != 1) ? static_cast( + std::floor((r - xmin_[0]) / dx_[0])) + + istart_[0] : istart_[0] ; + j = (nx_[1] != 1) ? static_cast(std::floor( + (theta - xmin_[1]) / dx_[1])) + + istart_[1] + : istart_[1]; + k = (nx_[2] != 1)? static_cast(std::floor( + (z - xmin_[2]) / dx_[2])) + + istart_[2] + : istart_[2]; + } + + const std::array &GetXmin() const { return xmin_; } + const std::array &GetStartIndex() const { return istart_; } + const char *Name() const { return name_; } + + //private: + std::array istart_; + std::array xmin_, dx_, area_, nx_; + constexpr static const char *name_ = "UniformCylindrical"; + + const std::array &Dx_() const { return dx_; } + + ParArrayND r_c_, coord_vol_i_, cosphi_c_, sinphi_c_; + + void InitCachedArrays() { + int nx1_tot = nx_[0]+2*Globals::nghost+1; + r_c_ = ParArrayND("r_c_", nx1_tot-1); + coord_vol_i_ = ParArrayND("", nx1_tot-1); + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "UniformCylindrical::InitCachedArrays(r)", + parthenon::DevExecSpace(), 0, nx1_tot-2, + KOKKOS_LAMBDA(const int &i) { + Real rm = xmin_[0] + i * dx_[0]; + Real rp = xmin_[0] + (i+1) * dx_[0]; + r_c_(i) = TWO_3RD*(rp*rp*rp - rm*rm*rm)/(SQR(rp) - SQR(rm)); + coord_vol_i_(i) = 0.5*(rp*rp - rm*rm); + }); + int nx2_tot = nx_[1]+1; + + if (istart_[1] > 0) { + nx2_tot += 2*Globals::nghost; + } + sinphi_c_ = ParArrayND("sin(phi_c)", nx2_tot-1); + cosphi_c_ = ParArrayND("cos(phi_c)", nx2_tot-1); + + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "UniformCylindrical::InitCachedArrays(phi)", + parthenon::DevExecSpace(), 0, nx2_tot-2, + KOKKOS_LAMBDA(const int &j) { + Real Phi = xmin_[1] + (j + 0.5) * dx_[1]; + cosphi_c_(j) = std::cos(Phi); + sinphi_c_(j) = std::sin(Phi); + }); + } + +}; + + +} // namespace parthenon + +#endif // COORDINATES_UNIFORM_CYLINDRICAL_HPP_ diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp new file mode 100644 index 000000000000..65c012e11b88 --- /dev/null +++ b/src/coordinates/uniform_spherical.hpp @@ -0,0 +1,561 @@ +//======================================================================================== +// (C) (or copyright) 2023. Triad National Security, LLC. All rights reserved. +// +// This program was produced under U.S. Government contract 89233218CNA000001 for Los +// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC +// for the U.S. Department of Energy/National Nuclear Security Administration. All rights +// in the program are reserved by Triad National Security, LLC, and the U.S. Department +// of Energy/National Nuclear Security Administration. The Government is granted for +// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide +// license in this material to reproduce, prepare derivative works, distribute copies to +// the public, perform publicly and display publicly, and to permit others to do so. +//======================================================================================== +#ifndef COORDINATES_UNIFORM_SPHERICAL_HPP_ +#define COORDINATES_UNIFORM_SPHERICAL_HPP_ + +#include +#include + +#include "basic_types.hpp" +#include "defs.hpp" +#include "globals.hpp" +#include "parameter_input.hpp" + +#include + +namespace parthenon { + +class UniformSpherical { + public: + UniformSpherical() = default; + UniformSpherical(const RegionSize &rs, ParameterInput *pin) { + dx_[0] = (rs.x1max - rs.x1min) / rs.nx1; + dx_[1] = (rs.x2max - rs.x2min) / rs.nx2; + dx_[2] = (rs.x3max - rs.x3min) / rs.nx3; + area_[0] = dx_[1] * dx_[2]; + area_[1] = dx_[0] * dx_[2]; + area_[2] = dx_[0] * dx_[1]; + istart_[0] = Globals::nghost; + istart_[1] = (rs.nx2 > 1 ? Globals::nghost : 0); + istart_[2] = (rs.nx3 > 1 ? Globals::nghost : 0); + xmin_[0] = rs.x1min - istart_[0] * dx_[0]; + xmin_[1] = rs.x2min - istart_[1] * dx_[1]; + xmin_[2] = rs.x3min - istart_[2] * dx_[2]; + + 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 != "spherical" ){ + if (coord_type_str == "cartesian") { + PARTHENON_THROW(" Please rebuild with -DCOORDINATE_TYPE=UniformCartesian"); + } else if (coord_type_str == "cylindrical") { + PARTHENON_THROW(" Please rebuild with -DCOORDINATE_TYPE=UniformCylindrical"); + } else { + PARTHENON_THROW("Invalid coord input in ."); + } + } + + //Initialized cached coordinates + InitCachedArrays(); + } + UniformSpherical(const UniformSpherical &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]; + + //Initialized cached coordinates + InitCachedArrays(); + } + + //---------------------------------------- + // Dxc: Distance between cell centers + // Dxc<1> and Dxc<2> returns distance between cell centroids + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir){ + case X1DIR: + return r_c_(i+1) - r_c_(i); + case X2DIR: + return theta_c_(j+1) - theta_c_(j); + default: + return dx_[2]; + } + } + template + KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir){ + case X1DIR: + return r_c_(i+1) - r_c_(i); + case X2DIR: + return theta_c_(j+1) - theta_c_(j); + default: + return dx_[2]; + } + } + + //---------------------------------------- + // Dxf: Distance between cell faces + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const { + assert(dir > 0 && dir < 4 && face > 0 && face < 4); + //(forrestglines): Double check how dir and face are used + return dx_[face - 1]; + } + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxf(Args... args) const { + assert(dir > 0 && dir < 4); + return dx_[dir - 1]; + } + + //---------------------------------------- + // Xc: Positions at cell centers + // Xc<1> returns r of cell centroid + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { + assert(dir > 0 && dir < 4); + switch (dir) { + case 1: + return r_c_(idx); + case 2: + return theta_c_(idx); + case 3: + return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + template + 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(i); + case 2: + return Xc(j); + case 3: + return Xc(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // Xf: Positions on Faces + //---------------------------------------- + template + 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(idx); + } + } + template + 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(i); + case 2: + return Xf(j); + case 3: + return Xf(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + template + 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 + 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(i); + case 2: + return Xf(j); + case 3: + return Xf(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + template + KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { + using TE = TopologicalElement; + bool constexpr X1EDGE = el == TE::F1 || el == TE::E2 || el == TE::E3 || el == TE::NN; + bool constexpr X2EDGE = el == TE::F2 || el == TE::E3 || el == TE::E1 || el == TE::NN; + bool constexpr X3EDGE = el == TE::F3 || el == TE::E1 || el == TE::E2 || el == TE::NN; + if constexpr (dir == X1DIR && X1EDGE) { + return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 + } else if constexpr (dir == X2DIR && X2EDGE) { + return xmin_[dir - 1] + idx * dx_[dir - 1]; // idx - 1/2 + } else if constexpr (dir == X3DIR && X3EDGE) { + 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 + 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(i); + case 2: + return X(j); + case 3: + return X(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // CellWidth: Physical width of cells at cell centers + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real CellWidth(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir){ + case X1DIR: + return Dxf<1>(k,j,i); + case X2DIR: + return Xc<1>(k,j,i)*Dxf<2>(k,j,i);//r*dphi; + case X3DIR: + return Xc<1>(k,j,i)*sintht_c_(j)*Dxf<3>(k,j,i); //r*sin(th)*dphi + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + KOKKOS_FORCEINLINE_FUNCTION Real CellWidthFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + if( dir == X2DIR ){ + return Xf<1>( k, j, i)*dx_[1]; //r*dphi + } else { + return dx_[dir-1]; + } + } + + //---------------------------------------- + // EdgeLength: Physical length of cell edges + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLength(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir){ + case X1DIR: + return Dxf<1>(k,j,i); + case X2DIR: + return Xf<1>(k,j,i)*Dxf<2>(k,j,i);//r*dphi; + case X3DIR: + return Xf<1>(k,j,i)*sintht_f_(j)*Dxf<3>(k,j,i); //r*sin(th)*dphi + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + KOKKOS_FORCEINLINE_FUNCTION Real EdgeLengthFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir){ + case X1DIR: + return EdgeLength(k, j, i); + case X2DIR: + return EdgeLength(k, j, i); + case X3DIR: + return EdgeLength(k, j, i); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // FaceArea: Physical area of cell areas + //---------------------------------------- + //(How is this different from Area(dir, k, j, i)? + template + KOKKOS_FORCEINLINE_FUNCTION Real FaceArea(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir) { + case X1DIR: + return (Xf<1>(i)*Xf<1>(i)*(costht_f_(j)-costht_f_(j+1))*dx_[2]); //r^2*d[-cos(tht)]*dph + case X2DIR: + return (coord_area2_i_(i)*sintht_f_(j)*dx_[2]);//rdr*sin(th)*dph + case X3DIR: + return (coord_area2_i_(i)*dx_[1]); //d(r^2/2)*dtheta; + } + + } + KOKKOS_FORCEINLINE_FUNCTION Real FaceAreaFA(const int dir, const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch(dir){ + case X1DIR: + return FaceArea(k, j, i); + case X2DIR: + return FaceArea(k, j, i); + case X3DIR: + return FaceArea(k, j, i); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + //---------------------------------------- + // CellVolume + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { + return (coord_vol_i_(i)*coord_area1_j_(j)*dx_[2]); + } + + //---------------------------------------- + // Generalized physical volume + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Volume(Args... args) const { + using TE = TopologicalElement; + if constexpr (el == TE::CC) { + return CellVolume(args...); + } else if constexpr (el == TE::F1) { + return FaceArea(args...); + } else if constexpr (el == TE::F2) { + return FaceArea(args...); + } else if constexpr (el == TE::F3) { + return FaceArea(args...); + } else if constexpr (el == TE::E1) { + return EdgeLength(args...); + } else if constexpr (el == TE::E2) { + return EdgeLength(args...); + } else if constexpr (el == TE::E3) { + return EdgeLength(args...); + } 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; + } + + //---------------------------------------- + // Cylindrical+Spherical Conversions + //---------------------------------------- + // TODO + + //---------------------------------------- + // Area Averaged Positions (for CT MHD?) AMR + //---------------------------------------- + // TODO + // x1s2 + // x1s3 + // x2s1 + // etc. + + //---------------------------------------- + // Coordinate source terms? + //---------------------------------------- + // TODO + // coord_src1_i + // coord_src2_i + // coord_area1_j + // coord_area2_j + + //---------------------------------------- + // Geometric Terms (Find better names for these!) + //---------------------------------------- + KOKKOS_FORCEINLINE_FUNCTION + Real h2v(const int i) const { return r_c_(i); }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h2f(const int i) const { return Xf<1>(i); }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h31v(const int i) const { return r_c_(i); }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h31f(const int i) const { return Xf<1>(i); }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh2vd1(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh2fd1(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh31vd1(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real dh31fd1(const int i) const { return 1.0; }; + + KOKKOS_FORCEINLINE_FUNCTION + Real h32v(const int j) const { return sintht_c_(j); } + + KOKKOS_FORCEINLINE_FUNCTION + Real h32f(const int j) const { return sintht_f_(j); } + + KOKKOS_FORCEINLINE_FUNCTION + Real dh32vd2(const int j) const { return costht_c_(j); } + + KOKKOS_FORCEINLINE_FUNCTION + Real dh32fd2(const int j) const { return costht_f_(j); } + + //What is this? + KOKKOS_FORCEINLINE_FUNCTION + void GetAcc(const int k, const int j, const int i, const Real &rp, + const Real &cosphip, const Real &sinphip, const Real &zp, + Real &acc1, Real &acc2, Real &acc3) const + { + const Real cosdphi = cosphi_c_(j)*cosphip - sinphi_c_(j)*sinphip; //cos(x3v(k)-phip) + const Real sindphi = sinphi_c_(j)*cosphip - cosphi_c_(j)*sinphip; //sin(x3v(k)-phip) + //Psi = - 1.0/sqrt(r0^2 + rp^2 - 2r0*rp*sin(tht)cos(phi-phip) + zp^2- 2*r0*cos(tht)*zp + acc1 = (r_c_(i) - rp*sintht_c_(j)*cosdphi - costht_c_(j)*zp); //dPsi/dr0 + acc2 = (-rp*costht_c_(j)*cosdphi + sintht_c_(j)*zp); //dPsi/(r0*dtht) + acc3 = (rp*sindphi); //dPsi/(r0*sintht*dphi) + } + + //---------------------------------------- + // Position to Cell index + //---------------------------------------- + KOKKOS_INLINE_FUNCTION + void Xtoijk(const Real &r, const Real &theta, const Real &phi, int &i, int &j, int &k) const { + i = (nx_[0] != 1) ? static_cast( + std::floor((r - xmin_[0]) / dx_[0])) + + istart_[0] : istart_[0]; + j = (nx_[1] != 1) ? static_cast(std::floor( + (theta - xmin_[1]) / dx_[1])) + + istart_[1] + : istart_[1]; + k = (nx_[2] != 1)? static_cast(std::floor( + (phi - xmin_[2]) / dx_[2])) + + istart_[2] + : istart_[2]; + } + + const std::array &GetXmin() const { return xmin_; } + const std::array &GetStartIndex() const { return istart_; } + const char *Name() const { return name_; } + + //private: + std::array istart_, nx_; + std::array xmin_, dx_, area_; + constexpr static const char *name_ = "UniformCylindrical"; + + const std::array &Dx_() const { return dx_; } + + ParArrayND r_c_, theta_c_, x1s2_, coord_vol_i_, coord_area2_i_, costht_f_, sintht_f_, sintht_c_, costht_c_, coord_area1_j_, cosphi_c_, sinphi_c_; + void InitCachedArrays() { + int nx1_tot = nx_[0]+2*Globals::nghost+1; + r_c_ = ParArrayND("centroid(r)", nx1_tot-1); + x1s2_ = ParArrayND("area1(r)", nx1_tot-1); + coord_vol_i_ = ParArrayND("volume(r)", nx1_tot-1); + coord_area2_i_ = ParArrayND("area2(r)", nx1_tot-1); + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "initializeArraySph_r", + parthenon::DevExecSpace(), 0, nx1_tot-2, + KOKKOS_LAMBDA(const int &i) { + //for (int i = 0; i < nx1_tot-1; i++) { + Real rm = xmin_[0] + i * dx_[0]; + Real rp = rm + dx_[0]; + r_c_(i) = 0.75*(std::pow(rp, 4) - std::pow(rm, 4)) / + (std::pow(rp, 3) - std::pow(rm, 3)); + x1s2_(i) = TWO_3RD*(std::pow(rp,3) - std::pow(rm,3))/(SQR(rp) - SQR(rm)); + coord_vol_i_(i) = (ONE_3RD)*(rp*rp*rp - rm*rm*rm); + coord_area2_i_(i) = 0.5*(rp*rp - rm*rm); + }); + + int nx2_tot = nx_[1]+1; + if (istart_[1] > 0) { + nx2_tot += 2*Globals::nghost; + } + costht_f_ = ParArrayND("cos(theta_f)", nx2_tot); + sintht_f_ = ParArrayND("sin(theta_f)", nx2_tot); + sintht_c_ = ParArrayND("sin(theta_c)", nx2_tot-1); + costht_c_ = ParArrayND("cos(theta_c)", nx2_tot-1); + coord_area1_j_ = ParArrayND("dcos(theta)", nx2_tot-1); + theta_c_ = ParArrayND("centroid-theta", std::max(2,nx2_tot-1)); + + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "initializeArraySph_th1", + parthenon::DevExecSpace(), 0, nx2_tot-1, + KOKKOS_LAMBDA(const int &j) { + //for (int j = 0; j < nx2_tot; j++) { + Real theta = xmin_[1] + j * dx_[1]; + costht_f_(j) = std::cos(theta); + sintht_f_(j) = std::sin(theta); + }); + + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "initializeArraySph_th2", + parthenon::DevExecSpace(), 0, nx2_tot-2, + KOKKOS_LAMBDA(const int &j) { + //for (int j = 0; j < nx2_tot-1; j++) { + coord_area1_j_(j) = std::abs(costht_f_(j ) - costht_f_(j+1)); + Real tm = xmin_[1] + j * dx_[1]; + Real tp = tm + dx_[1]; + theta_c_(j) = (((sintht_f_(j+1) - tp*costht_f_(j+1)) - + (sintht_f_(j ) - tm*costht_f_(j )))/ + (costht_f_(j ) - costht_f_(j+1))); + sintht_c_(j) = std::sin(theta_c_(j)); + costht_c_(j) = std::cos(theta_c_(j)); + }); + + if (nx2_tot==2) { + theta_c_(1) = theta_c_(0) + dx_[1]; + } + + //phi-direction + int nx3_tot = nx_[2]+1; + if (istart_[2] > 0) { + nx3_tot += 2*Globals::nghost; + } + sinphi_c_ = ParArrayND("sin(phi_c)", nx3_tot-1); + cosphi_c_ = ParArrayND("cos(phi_c)", nx3_tot-1); + parthenon::par_for( + parthenon::loop_pattern_flatrange_tag, "initializeArraySph_ph", + parthenon::DevExecSpace(), 0, nx3_tot-2, + KOKKOS_LAMBDA(const int &k) { + //for (int k = 0; k < nx3_tot-1; k++) { + Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; + cosphi_c_(k) = std::cos(Phi); + sinphi_c_(k) = std::sin(Phi); + }); + } + +}; + +} // namespace parthenon + +#endif // COORDINATES_UNIFORM_SPHERICAL_HPP_ diff --git a/src/interface/swarm_device_context.hpp b/src/interface/swarm_device_context.hpp index e01b542b4925..4f89e0afdc8d 100644 --- a/src/interface/swarm_device_context.hpp +++ b/src/interface/swarm_device_context.hpp @@ -86,17 +86,7 @@ class SwarmDeviceContext { // TODO(BRR) This logic will change for non-uniform cartesian meshes KOKKOS_INLINE_FUNCTION void Xtoijk(const Real &x, const Real &y, const Real &z, int &i, int &j, int &k) const { - i = static_cast( - std::floor((x - x_min_) / coords_.Dxc())) + - ib_s_; - j = (ndim_ > 1) ? static_cast(std::floor( - (y - y_min_) / coords_.Dxc())) + - jb_s_ - : jb_s_; - k = (ndim_ > 2) ? static_cast(std::floor( - (z - z_min_) / coords_.Dxc())) + - kb_s_ - : kb_s_; + coords_.Xtoijk(x,y,z,i,j,k); } KOKKOS_INLINE_FUNCTION From c32f395c50191c88d3e493bc592c01a7597b72c7 Mon Sep 17 00:00:00 2001 From: Forrest Glines Date: Mon, 31 Jul 2023 16:15:38 -0600 Subject: [PATCH 2/9] Add Xs, compilation fixes --- example/CMakeLists.txt | 9 +- src/coordinates/uniform_cartesian.hpp | 33 ++++++-- src/coordinates/uniform_cylindrical.hpp | 80 ++++++++++++++---- src/coordinates/uniform_spherical.hpp | 108 ++++++++++++++++++------ 4 files changed, 179 insertions(+), 51 deletions(-) diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index eb01d46e51bb..d016154e2d9f 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -11,6 +11,7 @@ # 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) @@ -18,5 +19,11 @@ 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() diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp index af3c14e76446..4c655115a9a6 100644 --- a/src/coordinates/uniform_cartesian.hpp +++ b/src/coordinates/uniform_cartesian.hpp @@ -120,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(i); - case 2: + case X2DIR: return Xc(j); - case 3: + case X3DIR: return Xc(k); default: PARTHENON_FAIL("Unknown dir"); @@ -149,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(i); - case 2: + case X2DIR: return Xf(j); - case 3: + case X3DIR: return Xf(k); default: PARTHENON_FAIL("Unknown dir"); @@ -171,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(i); - case 2: + case X2DIR: return Xf(j); - case 3: + case X3DIR: return Xf(k); default: PARTHENON_FAIL("Unknown dir"); return 0; // To appease compiler } } + + //---------------------------------------- + // Xs: Area averaged positions + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int idx) const { + assert(dir > 0 && dir < 4 && side > 0 && side < 4); + return Xc(idx); + } + template + 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(k,j,i); + } + template KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp index e69d129a6de7..5d5b6c188582 100644 --- a/src/coordinates/uniform_cylindrical.hpp +++ b/src/coordinates/uniform_cylindrical.hpp @@ -85,23 +85,54 @@ class UniformCylindrical { // Dxc<1> returns distance between cell centroids in r //---------------------------------------- template - KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const { assert(dir > 0 && dir < 4); if( dir == X1DIR ){ - return r_c_(i+1) - r_c_(i); + return r_c_(idx+1) - r_c_(idx); } else { return dx_[dir-1]; } } + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch (dir) { + case X1DIR: + return Dxc(i); + case X2DIR: + return Dxc(j); + case X3DIR: + return Dxc(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + + KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, const int idx) const { + assert(dir > 0 && dir < 4); + switch (dir) { + case X1DIR: + return Dxc(idx); + case X2DIR: + return Dxc(idx); + case X3DIR: + return Dxc(idx); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, const int k, const int j, const int i) const { assert(dir > 0 && dir < 4); switch (dir) { - case 1: - return Dxc<1>(k,j,i); - case 2: - return Dxc<2>(k,j,j); - case 3: - return Dxc<3>(k,j,k); + case X1DIR: + return Dxc(i); + case X2DIR: + return Dxc(j); + case X3DIR: + return Dxc(k); default: PARTHENON_FAIL("Unknown dir"); return 0; // To appease compiler @@ -139,11 +170,11 @@ class UniformCylindrical { 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(i); - case 2: + case X2DIR: return Xc(j); - case 3: + case X3DIR: return Xc(k); default: PARTHENON_FAIL("Unknown dir"); @@ -151,6 +182,7 @@ class UniformCylindrical { } } + //---------------------------------------- // Xf: Positions on Faces //---------------------------------------- @@ -168,11 +200,11 @@ class UniformCylindrical { 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(i); - case 2: + case X2DIR: return Xf(j); - case 3: + case X3DIR: return Xf(k); default: PARTHENON_FAIL("Unknown dir"); @@ -190,11 +222,11 @@ class UniformCylindrical { 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(i); - case 2: + case X2DIR: return Xf(j); - case 3: + case X3DIR: return Xf(k); default: PARTHENON_FAIL("Unknown dir"); @@ -202,6 +234,20 @@ class UniformCylindrical { } } + //---------------------------------------- + // Xs: Area averaged positions + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int idx) const { + assert(dir > 0 && dir < 4 && side > 0 && side < 4); + return Xc(idx); + } + template + 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(k,j,i); + } + template KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { using TE = TopologicalElement; diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp index 65c012e11b88..28bb2c95f0a7 100644 --- a/src/coordinates/uniform_spherical.hpp +++ b/src/coordinates/uniform_spherical.hpp @@ -82,28 +82,59 @@ class UniformSpherical { // Dxc: Distance between cell centers // Dxc<1> and Dxc<2> returns distance between cell centroids //---------------------------------------- - template - KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const { assert(dir > 0 && dir < 4); switch(dir){ case X1DIR: - return r_c_(i+1) - r_c_(i); + return r_c_(idx+1) - r_c_(idx); case X2DIR: - return theta_c_(j+1) - theta_c_(j); + return theta_c_(idx+1) - theta_c_(idx); default: return dx_[2]; } } - template + template + KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4); + switch (dir) { + case X1DIR: + return Dxc(i); + case X2DIR: + return Dxc(j); + case X3DIR: + return Dxc(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } + KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, const int idx) const { + assert(dir > 0 && dir < 4); + switch (dir) { + case X1DIR: + return Dxc(idx); + case X2DIR: + return Dxc(idx); + case X3DIR: + return Dxc(idx); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } KOKKOS_FORCEINLINE_FUNCTION Real DxcFA(const int dir, const int k, const int j, const int i) const { assert(dir > 0 && dir < 4); - switch(dir){ - case X1DIR: - return r_c_(i+1) - r_c_(i); - case X2DIR: - return theta_c_(j+1) - theta_c_(j); - default: - return dx_[2]; + switch (dir) { + case X1DIR: + return Dxc(i); + case X2DIR: + return Dxc(j); + case X3DIR: + return Dxc(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler } } @@ -130,11 +161,11 @@ class UniformSpherical { KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { assert(dir > 0 && dir < 4); switch (dir) { - case 1: + case X1DIR: return r_c_(idx); - case 2: + case X2DIR: return theta_c_(idx); - case 3: + case X3DIR: return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; default: PARTHENON_FAIL("Unknown dir"); @@ -145,11 +176,11 @@ class UniformSpherical { 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(i); - case 2: + case X2DIR: return Xc(j); - case 3: + case X3DIR: return Xc(k); default: PARTHENON_FAIL("Unknown dir"); @@ -174,11 +205,11 @@ class UniformSpherical { 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(i); - case 2: + case X2DIR: return Xf(j); - case 3: + case X3DIR: return Xf(k); default: PARTHENON_FAIL("Unknown dir"); @@ -196,17 +227,46 @@ class UniformSpherical { 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(i); - case 2: + case X2DIR: return Xf(j); - case 3: + case X3DIR: return Xf(k); default: PARTHENON_FAIL("Unknown dir"); return 0; // To appease compiler } } + + //---------------------------------------- + // Xs: Area averaged positions + //---------------------------------------- + template + KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int idx) const { + assert(dir > 0 && dir < 4 && side > 0 && side < 4); + if( dir == X1DIR){ + (2.0/3.0)*(std::pow(Xf<1>(idx+1),3) - std::pow(Xf<1>(idx),3)) + /(SQR(Xf<1>(idx+1)) - SQR(Xf<1>(idx))); + } else { + return Xc(idx); + } + } + template + KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int k, const int j, const int i) const { + assert(dir > 0 && dir < 4 && side > 0 && side < 4); + switch (dir) { + case X1DIR: + return Xs(i); + case X2DIR: + return Xs(j); + case X3DIR: + return Xs(k); + default: + PARTHENON_FAIL("Unknown dir"); + return 0; // To appease compiler + } + } template KOKKOS_FORCEINLINE_FUNCTION Real X(const int idx) const { From 0172270198d3ecd196d24cfbaa340bc02a47007b Mon Sep 17 00:00:00 2001 From: Forrest Glines Date: Mon, 7 Aug 2023 10:58:28 -0600 Subject: [PATCH 3/9] Turned cached arrays into member functions --- src/coordinates/uniform_cylindrical.hpp | 76 +++++------ src/coordinates/uniform_spherical.hpp | 165 ++++++++++++------------ 2 files changed, 109 insertions(+), 132 deletions(-) diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp index 5d5b6c188582..64ba9b1e5e97 100644 --- a/src/coordinates/uniform_cylindrical.hpp +++ b/src/coordinates/uniform_cylindrical.hpp @@ -55,9 +55,6 @@ class UniformCylindrical { PARTHENON_THROW("Invalid coord input in ."); } } - - //Initialized cached coordinates - InitCachedArrays(); } UniformCylindrical(const UniformCylindrical &src, int coarsen) : istart_(src.GetStartIndex()) { @@ -75,9 +72,6 @@ class UniformCylindrical { nx_[0] = src.nx_[0] / ( coarsen ? 2 : 1 ); nx_[1] = src.nx_[1] / ( coarsen ? 2 : 1 ); nx_[2] = src.nx_[2] / ( coarsen ? 2 : 1 ); - - //Initialized cached coordinates - InitCachedArrays(); } //---------------------------------------- @@ -88,7 +82,7 @@ class UniformCylindrical { KOKKOS_FORCEINLINE_FUNCTION Real Dxc(const int idx) const { assert(dir > 0 && dir < 4); if( dir == X1DIR ){ - return r_c_(idx+1) - r_c_(idx); + return R_c_(idx+1) - R_c_(idx); } else { return dx_[dir-1]; } @@ -161,7 +155,7 @@ class UniformCylindrical { KOKKOS_FORCEINLINE_FUNCTION Real Xc(const int idx) const { assert(dir > 0 && dir < 4); if( dir == X1DIR ){ - return r_c_(idx); + return R_c_(idx); } else { return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; } @@ -336,7 +330,7 @@ class UniformCylindrical { case X2DIR: return area_[1]; //dr*dz case X3DIR: - return coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi + return Coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi } } @@ -349,7 +343,7 @@ class UniformCylindrical { case X2DIR: return area_[1]; //dr*dz case X3DIR: - return coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi + return Coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi } } @@ -358,7 +352,7 @@ class UniformCylindrical { //---------------------------------------- template KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { - return coord_vol_i_(i)*area_[0]; + return Coord_vol_i_(i)*area_[0]; } //---------------------------------------- @@ -398,7 +392,7 @@ class UniformCylindrical { // Geometric Terms (Find better names for these!) //---------------------------------------- KOKKOS_FORCEINLINE_FUNCTION - Real h2v(const int i) const { return r_c_(i); }; + Real h2v(const int i) const { return R_c_(i); }; KOKKOS_FORCEINLINE_FUNCTION Real h2f(const int i) const { return Xf<1>(i); }; @@ -439,10 +433,10 @@ class UniformCylindrical { const Real &cosphip, const Real &sinphip, const Real &zp, Real &acc1, Real &acc2, Real &acc3) const { - const Real cosdphi = cosphi_c_(j)*cosphip - sinphi_c_(j)*sinphip; //cos(x3v(k)-phip) - const Real sindphi = sinphi_c_(j)*cosphip - cosphi_c_(j)*sinphip; //sin(x3v(k)-phip) + const Real cosdphi = Cos_phi_c_(j)*cosphip - Sin_phi_c_(j)*sinphip; //cos(x3v(k)-phip) + const Real sindphi = Sin_phi_c_(j)*cosphip - Cos_phi_c_(j)*sinphip; //sin(x3v(k)-phip) //Psi = - 1.0/sqrt(r0^2 + rp^2 - 2r0*rp*cos(phi-phip) + (z-zp)^2) - acc1 = (r_c_(i) - rp*cosdphi); //acc1 = dPsi/dr, Psi=-1/dist2 + acc1 = (R_c_(i) - rp*cosdphi); //acc1 = dPsi/dr, Psi=-1/dist2 acc2 = rp*sindphi; //acc2 = dPsi/(r*dphi) acc3 = Xc<3>(k)-zp; } @@ -469,44 +463,32 @@ class UniformCylindrical { const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return name_; } - //private: + private: std::array istart_; std::array xmin_, dx_, area_, nx_; constexpr static const char *name_ = "UniformCylindrical"; const std::array &Dx_() const { return dx_; } - ParArrayND r_c_, coord_vol_i_, cosphi_c_, sinphi_c_; - - void InitCachedArrays() { - int nx1_tot = nx_[0]+2*Globals::nghost+1; - r_c_ = ParArrayND("r_c_", nx1_tot-1); - coord_vol_i_ = ParArrayND("", nx1_tot-1); - parthenon::par_for( - parthenon::loop_pattern_flatrange_tag, "UniformCylindrical::InitCachedArrays(r)", - parthenon::DevExecSpace(), 0, nx1_tot-2, - KOKKOS_LAMBDA(const int &i) { - Real rm = xmin_[0] + i * dx_[0]; - Real rp = xmin_[0] + (i+1) * dx_[0]; - r_c_(i) = TWO_3RD*(rp*rp*rp - rm*rm*rm)/(SQR(rp) - SQR(rm)); - coord_vol_i_(i) = 0.5*(rp*rp - rm*rm); - }); - int nx2_tot = nx_[1]+1; - - if (istart_[1] > 0) { - nx2_tot += 2*Globals::nghost; - } - sinphi_c_ = ParArrayND("sin(phi_c)", nx2_tot-1); - cosphi_c_ = ParArrayND("cos(phi_c)", nx2_tot-1); - - parthenon::par_for( - parthenon::loop_pattern_flatrange_tag, "UniformCylindrical::InitCachedArrays(phi)", - parthenon::DevExecSpace(), 0, nx2_tot-2, - KOKKOS_LAMBDA(const int &j) { - Real Phi = xmin_[1] + (j + 0.5) * dx_[1]; - cosphi_c_(j) = std::cos(Phi); - sinphi_c_(j) = std::sin(Phi); - }); + //ParArrayND r_c_, coord_vol_i_, cos_phi_c_, sin_phi_c_; + + KOKKOS_INLINE_FUNCTION Real R_c_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = xmin_[0] + (i+1) * dx_[0]; + return TWO_3RD*(rp*rp*rp - rm*rm*rm)/(SQR(rp) - SQR(rm)); + } + KOKKOS_INLINE_FUNCTION Real Coord_vol_i_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = xmin_[0] + (i+1) * dx_[0]; + return 0.5*(rp*rp - rm*rm); + } + KOKKOS_INLINE_FUNCTION Real Cos_phi_c_(const int j) const { + const Real Phi = xmin_[1] + (j + 0.5) * dx_[1]; + return std::cos(Phi); + } + KOKKOS_INLINE_FUNCTION Real Sin_phi_c_(const int j) const { + const Real Phi = xmin_[1] + (j + 0.5) * dx_[1]; + return std::sin(Phi); } }; diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp index 28bb2c95f0a7..72b16863508d 100644 --- a/src/coordinates/uniform_spherical.hpp +++ b/src/coordinates/uniform_spherical.hpp @@ -56,9 +56,6 @@ class UniformSpherical { PARTHENON_THROW("Invalid coord input in ."); } } - - //Initialized cached coordinates - InitCachedArrays(); } UniformSpherical(const UniformSpherical &src, int coarsen) : istart_(src.GetStartIndex()) { @@ -73,9 +70,6 @@ class UniformSpherical { area_[0] = dx_[1] * dx_[2]; area_[1] = dx_[0] * dx_[2]; area_[2] = dx_[0] * dx_[1]; - - //Initialized cached coordinates - InitCachedArrays(); } //---------------------------------------- @@ -87,9 +81,9 @@ class UniformSpherical { assert(dir > 0 && dir < 4); switch(dir){ case X1DIR: - return r_c_(idx+1) - r_c_(idx); + return R_c_(idx+1) - R_c_(idx); case X2DIR: - return theta_c_(idx+1) - theta_c_(idx); + return Theta_c_(idx+1) - Theta_c_(idx); default: return dx_[2]; } @@ -162,9 +156,9 @@ class UniformSpherical { assert(dir > 0 && dir < 4); switch (dir) { case X1DIR: - return r_c_(idx); + return R_c_(idx); case X2DIR: - return theta_c_(idx); + return Theta_c_(idx); case X3DIR: return xmin_[dir - 1] + (idx + 0.5) * dx_[dir - 1]; default: @@ -246,8 +240,9 @@ class UniformSpherical { KOKKOS_FORCEINLINE_FUNCTION Real Xs(const int idx) const { assert(dir > 0 && dir < 4 && side > 0 && side < 4); if( dir == X1DIR){ - (2.0/3.0)*(std::pow(Xf<1>(idx+1),3) - std::pow(Xf<1>(idx),3)) + return (2.0/3.0)*(std::pow(Xf<1>(idx+1),3) - std::pow(Xf<1>(idx),3)) /(SQR(Xf<1>(idx+1)) - SQR(Xf<1>(idx))); + //FIXME: Shouldn't we use X1s2_? } else { return Xc(idx); } @@ -457,13 +452,13 @@ class UniformSpherical { // Geometric Terms (Find better names for these!) //---------------------------------------- KOKKOS_FORCEINLINE_FUNCTION - Real h2v(const int i) const { return r_c_(i); }; + Real h2v(const int i) const { return R_c_(i); }; KOKKOS_FORCEINLINE_FUNCTION Real h2f(const int i) const { return Xf<1>(i); }; KOKKOS_FORCEINLINE_FUNCTION - Real h31v(const int i) const { return r_c_(i); }; + Real h31v(const int i) const { return R_c_(i); }; KOKKOS_FORCEINLINE_FUNCTION Real h31f(const int i) const { return Xf<1>(i); }; @@ -501,7 +496,7 @@ class UniformSpherical { const Real cosdphi = cosphi_c_(j)*cosphip - sinphi_c_(j)*sinphip; //cos(x3v(k)-phip) const Real sindphi = sinphi_c_(j)*cosphip - cosphi_c_(j)*sinphip; //sin(x3v(k)-phip) //Psi = - 1.0/sqrt(r0^2 + rp^2 - 2r0*rp*sin(tht)cos(phi-phip) + zp^2- 2*r0*cos(tht)*zp - acc1 = (r_c_(i) - rp*sintht_c_(j)*cosdphi - costht_c_(j)*zp); //dPsi/dr0 + acc1 = (R_c_(i) - rp*sintht_c_(j)*cosdphi - costht_c_(j)*zp); //dPsi/dr0 acc2 = (-rp*costht_c_(j)*cosdphi + sintht_c_(j)*zp); //dPsi/(r0*dtht) acc3 = (rp*sindphi); //dPsi/(r0*sintht*dphi) } @@ -536,82 +531,82 @@ class UniformSpherical { const std::array &Dx_() const { return dx_; } ParArrayND r_c_, theta_c_, x1s2_, coord_vol_i_, coord_area2_i_, costht_f_, sintht_f_, sintht_c_, costht_c_, coord_area1_j_, cosphi_c_, sinphi_c_; - void InitCachedArrays() { - int nx1_tot = nx_[0]+2*Globals::nghost+1; - r_c_ = ParArrayND("centroid(r)", nx1_tot-1); - x1s2_ = ParArrayND("area1(r)", nx1_tot-1); - coord_vol_i_ = ParArrayND("volume(r)", nx1_tot-1); - coord_area2_i_ = ParArrayND("area2(r)", nx1_tot-1); - parthenon::par_for( - parthenon::loop_pattern_flatrange_tag, "initializeArraySph_r", - parthenon::DevExecSpace(), 0, nx1_tot-2, - KOKKOS_LAMBDA(const int &i) { - //for (int i = 0; i < nx1_tot-1; i++) { - Real rm = xmin_[0] + i * dx_[0]; - Real rp = rm + dx_[0]; - r_c_(i) = 0.75*(std::pow(rp, 4) - std::pow(rm, 4)) / + + KOKKOS_INLINE_FUNCTION + Real R_c_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = rm + dx_[0]; + return 0.75*(std::pow(rp, 4) - std::pow(rm, 4)) / (std::pow(rp, 3) - std::pow(rm, 3)); - x1s2_(i) = TWO_3RD*(std::pow(rp,3) - std::pow(rm,3))/(SQR(rp) - SQR(rm)); - coord_vol_i_(i) = (ONE_3RD)*(rp*rp*rp - rm*rm*rm); - coord_area2_i_(i) = 0.5*(rp*rp - rm*rm); - }); - - int nx2_tot = nx_[1]+1; - if (istart_[1] > 0) { - nx2_tot += 2*Globals::nghost; - } - costht_f_ = ParArrayND("cos(theta_f)", nx2_tot); - sintht_f_ = ParArrayND("sin(theta_f)", nx2_tot); - sintht_c_ = ParArrayND("sin(theta_c)", nx2_tot-1); - costht_c_ = ParArrayND("cos(theta_c)", nx2_tot-1); - coord_area1_j_ = ParArrayND("dcos(theta)", nx2_tot-1); - theta_c_ = ParArrayND("centroid-theta", std::max(2,nx2_tot-1)); - - parthenon::par_for( - parthenon::loop_pattern_flatrange_tag, "initializeArraySph_th1", - parthenon::DevExecSpace(), 0, nx2_tot-1, - KOKKOS_LAMBDA(const int &j) { - //for (int j = 0; j < nx2_tot; j++) { - Real theta = xmin_[1] + j * dx_[1]; - costht_f_(j) = std::cos(theta); - sintht_f_(j) = std::sin(theta); - }); - - parthenon::par_for( - parthenon::loop_pattern_flatrange_tag, "initializeArraySph_th2", - parthenon::DevExecSpace(), 0, nx2_tot-2, - KOKKOS_LAMBDA(const int &j) { - //for (int j = 0; j < nx2_tot-1; j++) { - coord_area1_j_(j) = std::abs(costht_f_(j ) - costht_f_(j+1)); + } + + KOKKOS_INLINE_FUNCTION + Real X1s2_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = rm + dx_[0]; + return TWO_3RD*(std::pow(rp,3) - std::pow(rm,3))/(SQR(rp) - SQR(rm)); + } + + KOKKOS_INLINE_FUNCTION + Real Coord_vol_i_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = rm + dx_[0]; + return (ONE_3RD)*(rp*rp*rp - rm*rm*rm); + } + + KOKKOS_INLINE_FUNCTION + Real Coord_area2_i_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = rm + dx_[0]; + return (ONE_3RD)*(rp*rp*rp - rm*rm*rm); + } + + KOKKOS_INLINE_FUNCTION + Real Coord_area1_j_(const int j) const { + return std::abs( Cos_theta_f_(j) - Cos_theta_f_(j+1)); + } + + KOKKOS_INLINE_FUNCTION + Real Cos_theta_f_(const int j) const { + const Real theta = xmin_[1] + j * dx_[1]; + return std::cos(theta); + } + + KOKKOS_INLINE_FUNCTION + Real Sin_theta_f_(const int j) const { + const Real theta = xmin_[1] + j * dx_[1]; + return std::sin(theta); + } + + KOKKOS_INLINE_FUNCTION + Real Theta_c_(const int j) const { Real tm = xmin_[1] + j * dx_[1]; Real tp = tm + dx_[1]; - theta_c_(j) = (((sintht_f_(j+1) - tp*costht_f_(j+1)) - - (sintht_f_(j ) - tm*costht_f_(j )))/ - (costht_f_(j ) - costht_f_(j+1))); - sintht_c_(j) = std::sin(theta_c_(j)); - costht_c_(j) = std::cos(theta_c_(j)); - }); - - if (nx2_tot==2) { - theta_c_(1) = theta_c_(0) + dx_[1]; - } + return (((Sin_theta_f_(j+1) - tp*Cos_theta_f_(j+1)) - + (Sin_theta_f_(j ) - tm*Cos_theta_f_(j )))/ + (Cos_theta_f_(j ) - Cos_theta_f_(j+1))); + } - //phi-direction - int nx3_tot = nx_[2]+1; - if (istart_[2] > 0) { - nx3_tot += 2*Globals::nghost; - } - sinphi_c_ = ParArrayND("sin(phi_c)", nx3_tot-1); - cosphi_c_ = ParArrayND("cos(phi_c)", nx3_tot-1); - parthenon::par_for( - parthenon::loop_pattern_flatrange_tag, "initializeArraySph_ph", - parthenon::DevExecSpace(), 0, nx3_tot-2, - KOKKOS_LAMBDA(const int &k) { - //for (int k = 0; k < nx3_tot-1; k++) { + KOKKOS_INLINE_FUNCTION + Real Cos_theta_c_(const int j) const { + return std::cos(Theta_c_(j)); + } + + KOKKOS_INLINE_FUNCTION + Real Sin_theta_c_(const int j) const { + return std::sin(Theta_c_(j)); + } + + KOKKOS_INLINE_FUNCTION + Real Cos_phi_c_(const int k) const { + Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; + return std::cos(Phi); + } + + KOKKOS_INLINE_FUNCTION + Real Sin_phi_c_(const int k) const { Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; - cosphi_c_(k) = std::cos(Phi); - sinphi_c_(k) = std::sin(Phi); - }); + return std::sin(Phi); } }; From 5b6bb61906f7c278f9724ee9f38e79dee8707098 Mon Sep 17 00:00:00 2001 From: Philipp Grete Date: Wed, 23 Aug 2023 15:53:41 -0400 Subject: [PATCH 4/9] Always use chuncking when compression is supported even not used --- src/outputs/parthenon_hdf5.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp index 66c1a5405490..d01cd8fea13c 100644 --- a/src/outputs/parthenon_hdf5.cpp +++ b/src/outputs/parthenon_hdf5.cpp @@ -444,11 +444,11 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm local_count[vinfo.tensor_rank + 3] = global_count[vinfo.tensor_rank + 3] = nx1; #ifndef PARTHENON_DISABLE_HDF5_COMPRESSION - if (output_params.hdf5_compression_level > 0) { + // if (output_params.hdf5_compression_level > 0) { for (int i = ndim - 3; i < ndim; i++) { chunk_size[i] = local_count[i]; } - } + // } #endif } else if (vinfo.where == MetadataFlag(Metadata::None)) { ndim = vinfo.tensor_rank + 1; @@ -457,12 +457,12 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm } #ifndef PARTHENON_DISABLE_HDF5_COMPRESSION - if (output_params.hdf5_compression_level > 0) { + // if (output_params.hdf5_compression_level > 0) { int nchunk_indices = std::min(vinfo.tensor_rank, 3); for (int i = ndim - nchunk_indices; i < ndim; i++) { chunk_size[i] = alldims[6 - nchunk_indices + i]; } - } + // } #endif } else { PARTHENON_THROW("Only Cell and None locations supported!"); From 0cc1c6c018176c92fa19c013cc75d63b4be1c130 Mon Sep 17 00:00:00 2001 From: Forrest Glines Date: Wed, 23 Aug 2023 15:13:42 -0600 Subject: [PATCH 5/9] Small cmake change --- src/CMakeLists.txt | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1258809e7873..70b84f027bf8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -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 "") # TODO: Put something more descriptive here set(COMPILER_FLAGS "") # TODO: Put something more descriptive here -set(COORDINATE_TYPE "UniformCartesian" CACHE STRING "UniformCartesian/UniformSpherical/UniformCylindrical") - configure_file(config.hpp.in generated/config.hpp @ONLY) add_library(parthenon From 4a01ed5b330b3e0b432c0ea5917efe46c0b2c106 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Mon, 30 Oct 2023 12:38:48 -0600 Subject: [PATCH 6/9] Added src term helpers for cylindrical --- src/coordinates/uniform_cylindrical.hpp | 37 +++++++++++++++++++++---- 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp index 64ba9b1e5e97..ecfb547e56ab 100644 --- a/src/coordinates/uniform_cylindrical.hpp +++ b/src/coordinates/uniform_cylindrical.hpp @@ -332,7 +332,7 @@ class UniformCylindrical { case X3DIR: return Coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi } - + return NAN; //To appease compiler } template KOKKOS_FORCEINLINE_FUNCTION Real FaceAreaFA(const int dir, const int k, const int j, const int i) const { @@ -345,6 +345,7 @@ class UniformCylindrical { case X3DIR: return Coord_vol_i_(i)*dx_[1]; //d(r^2/2)*dphi } + return NAN; //To appease compiler } //---------------------------------------- @@ -459,6 +460,35 @@ class UniformCylindrical { : istart_[2]; } + //---------------------------------------- + // Terms for Source Terms + //---------------------------------------- + + KOKKOS_INLINE_FUNCTION Real Coord_vol_i_(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = xmin_[0] + (i+1) * dx_[0]; + return 0.5*(rp*rp - rm*rm); + } + + KOKKOS_INLINE_FUNCTION Real CoordSrc1i(const int i) const { + return Dxf<1,1>(i)/Coord_vol_i_(i); + } + KOKKOS_INLINE_FUNCTION Real CoordSrc2i(const int i) const { + const Real rm = xmin_[0] + i * dx_[0]; + const Real rp = xmin_[0] + (i+1) * dx_[0]; + return Dxf<1,1>(i)/( (rm + rp)*Coord_vol_i_(i)); + } + + KOKKOS_INLINE_FUNCTION + Real PhySrc1i(const int i) const { + return 1./( Xc<1>(i)*Xf<1>(i) ); + } + + KOKKOS_INLINE_FUNCTION + Real PhySrc2i(const int i) const { + return 1./( Xc<1>(i)*Xf<1>(i+1) ); + } + const std::array &GetXmin() const { return xmin_; } const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return name_; } @@ -477,11 +507,6 @@ class UniformCylindrical { const Real rp = xmin_[0] + (i+1) * dx_[0]; return TWO_3RD*(rp*rp*rp - rm*rm*rm)/(SQR(rp) - SQR(rm)); } - KOKKOS_INLINE_FUNCTION Real Coord_vol_i_(const int i) const { - const Real rm = xmin_[0] + i * dx_[0]; - const Real rp = xmin_[0] + (i+1) * dx_[0]; - return 0.5*(rp*rp - rm*rm); - } KOKKOS_INLINE_FUNCTION Real Cos_phi_c_(const int j) const { const Real Phi = xmin_[1] + (j + 0.5) * dx_[1]; return std::cos(Phi); From 2cf79a798552d3e4b15dc39fcfb1ab1e32be85e9 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Sun, 5 Nov 2023 14:07:33 -0700 Subject: [PATCH 7/9] WIP:Spherical coordinates --- src/coordinates/uniform_spherical.hpp | 40 ++++++++++++++------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp index 72b16863508d..4d34a31d8d03 100644 --- a/src/coordinates/uniform_spherical.hpp +++ b/src/coordinates/uniform_spherical.hpp @@ -519,6 +519,27 @@ class UniformSpherical { : istart_[2]; } + //---------------------------------------- + // Terms for Source Terms + //---------------------------------------- + + KOKKOS_INLINE_FUNCTION Real Coord_vol_i_(const int i) const { + const Real rm = Xf(i); + const Real rp = Xf(i+1); + return ONE_3RD*( std::pow(rp,3) - std::pow(rm,3)); + } + + KOKKOS_INLINE_FUNCTION Real CoordSrc1i(const int i) const { + const Real rm = Xf(i); + const Real rp = Xf(i+1); + return 0.5*(SQR(rp) - SQR(rm))/Coord_vol_i_(i); + } + KOKKOS_INLINE_FUNCTION Real CoordSrc2i(const int i) const { + const Real rm = Xf(i); + const Real rp = Xf(i+1); + return Dxf(i)/( (rm + rp)*Coord_vol_i_(i)); + } + const std::array &GetXmin() const { return xmin_; } const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return name_; } @@ -547,25 +568,6 @@ class UniformSpherical { return TWO_3RD*(std::pow(rp,3) - std::pow(rm,3))/(SQR(rp) - SQR(rm)); } - KOKKOS_INLINE_FUNCTION - Real Coord_vol_i_(const int i) const { - const Real rm = xmin_[0] + i * dx_[0]; - const Real rp = rm + dx_[0]; - return (ONE_3RD)*(rp*rp*rp - rm*rm*rm); - } - - KOKKOS_INLINE_FUNCTION - Real Coord_area2_i_(const int i) const { - const Real rm = xmin_[0] + i * dx_[0]; - const Real rp = rm + dx_[0]; - return (ONE_3RD)*(rp*rp*rp - rm*rm*rm); - } - - KOKKOS_INLINE_FUNCTION - Real Coord_area1_j_(const int j) const { - return std::abs( Cos_theta_f_(j) - Cos_theta_f_(j+1)); - } - KOKKOS_INLINE_FUNCTION Real Cos_theta_f_(const int j) const { const Real theta = xmin_[1] + j * dx_[1]; From 59d4054012e989c5a03df4b5e9bbbbc62e386d00 Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Tue, 7 Nov 2023 15:11:56 -0700 Subject: [PATCH 8/9] Fixed spherical coords --- src/coordinates/uniform_spherical.hpp | 85 +++++++++++++++------------ 1 file changed, 48 insertions(+), 37 deletions(-) diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp index 4d34a31d8d03..4b3127195a31 100644 --- a/src/coordinates/uniform_spherical.hpp +++ b/src/coordinates/uniform_spherical.hpp @@ -309,7 +309,7 @@ class UniformSpherical { case X2DIR: return Xc<1>(k,j,i)*Dxf<2>(k,j,i);//r*dphi; case X3DIR: - return Xc<1>(k,j,i)*sintht_c_(j)*Dxf<3>(k,j,i); //r*sin(th)*dphi + return Xc<1>(k,j,i)*Sin_theta_c_(j)*Dxf<3>(k,j,i); //r*sin(th)*dphi default: PARTHENON_FAIL("Unknown dir"); return 0; // To appease compiler @@ -336,7 +336,7 @@ class UniformSpherical { case X2DIR: return Xf<1>(k,j,i)*Dxf<2>(k,j,i);//r*dphi; case X3DIR: - return Xf<1>(k,j,i)*sintht_f_(j)*Dxf<3>(k,j,i); //r*sin(th)*dphi + return Xf<1>(k,j,i)*Sin_theta_f_(j)*Dxf<3>(k,j,i); //r*sin(th)*dphi default: PARTHENON_FAIL("Unknown dir"); return 0; // To appease compiler @@ -366,11 +366,11 @@ class UniformSpherical { assert(dir > 0 && dir < 4); switch(dir) { case X1DIR: - return (Xf<1>(i)*Xf<1>(i)*(costht_f_(j)-costht_f_(j+1))*dx_[2]); //r^2*d[-cos(tht)]*dph + return (Xf<1>(i)*Xf<1>(i)*(Cos_theta_f_(j)-Cos_theta_f_(j+1))*dx_[2]); //r^2*d[-cos(tht)]*dph case X2DIR: - return (coord_area2_i_(i)*sintht_f_(j)*dx_[2]);//rdr*sin(th)*dph + return (Coord_area2_i_(i)*Sin_theta_f_(j)*dx_[2]);//rdr*sin(th)*dph case X3DIR: - return (coord_area2_i_(i)*dx_[1]); //d(r^2/2)*dtheta; + return (Coord_area2_i_(i)*dx_[1]); //d(r^2/2)*dtheta; } } @@ -394,7 +394,7 @@ class UniformSpherical { //---------------------------------------- template KOKKOS_FORCEINLINE_FUNCTION Real CellVolume(const int k, const int j, const int i) const { - return (coord_vol_i_(i)*coord_area1_j_(j)*dx_[2]); + return (Coord_vol_i_(i)*Coord_area1_j_(j)*dx_[2]); } //---------------------------------------- @@ -476,16 +476,16 @@ class UniformSpherical { Real dh31fd1(const int i) const { return 1.0; }; KOKKOS_FORCEINLINE_FUNCTION - Real h32v(const int j) const { return sintht_c_(j); } + Real h32v(const int j) const { return Sin_theta_c_(j); } KOKKOS_FORCEINLINE_FUNCTION - Real h32f(const int j) const { return sintht_f_(j); } + Real h32f(const int j) const { return Sin_theta_f_(j); } KOKKOS_FORCEINLINE_FUNCTION - Real dh32vd2(const int j) const { return costht_c_(j); } + Real dh32vd2(const int j) const { return Cos_theta_c_(j); } KOKKOS_FORCEINLINE_FUNCTION - Real dh32fd2(const int j) const { return costht_f_(j); } + Real dh32fd2(const int j) const { return Cos_theta_f_(j); } //What is this? KOKKOS_FORCEINLINE_FUNCTION @@ -493,11 +493,11 @@ class UniformSpherical { const Real &cosphip, const Real &sinphip, const Real &zp, Real &acc1, Real &acc2, Real &acc3) const { - const Real cosdphi = cosphi_c_(j)*cosphip - sinphi_c_(j)*sinphip; //cos(x3v(k)-phip) - const Real sindphi = sinphi_c_(j)*cosphip - cosphi_c_(j)*sinphip; //sin(x3v(k)-phip) + const Real cosdphi = Cos_phi_c_(j)*cosphip - Sin_phi_c_(j)*sinphip; //cos(x3v(k)-phip) + const Real sindphi = Sin_phi_c_(j)*cosphip - Cos_phi_c_(j)*sinphip; //sin(x3v(k)-phip) //Psi = - 1.0/sqrt(r0^2 + rp^2 - 2r0*rp*sin(tht)cos(phi-phip) + zp^2- 2*r0*cos(tht)*zp - acc1 = (R_c_(i) - rp*sintht_c_(j)*cosdphi - costht_c_(j)*zp); //dPsi/dr0 - acc2 = (-rp*costht_c_(j)*cosdphi + sintht_c_(j)*zp); //dPsi/(r0*dtht) + acc1 = (R_c_(i) - rp*Sin_theta_c_(j)*cosdphi - Cos_theta_c_(j)*zp); //dPsi/dr0 + acc2 = (-rp*Cos_theta_c_(j)*cosdphi + Sin_theta_c_(j)*zp); //dPsi/(r0*dtht) acc3 = (rp*sindphi); //dPsi/(r0*sintht*dphi) } @@ -547,68 +547,79 @@ class UniformSpherical { //private: std::array istart_, nx_; std::array xmin_, dx_, area_; - constexpr static const char *name_ = "UniformCylindrical"; + constexpr static const char *name_ = "UniformSpherical"; const std::array &Dx_() const { return dx_; } - ParArrayND r_c_, theta_c_, x1s2_, coord_vol_i_, coord_area2_i_, costht_f_, sintht_f_, sintht_c_, costht_c_, coord_area1_j_, cosphi_c_, sinphi_c_; + //ParArrayND r_c_, theta_c_, x1s2_, coord_vol_i_, coord_area2_i_, Cos_theta_f_, Sin_theta_f_, Sin_theta_c_, Cos_theta_c_, coord_area1_j_, Cos_phi_c_, Sin_phi_c_; KOKKOS_INLINE_FUNCTION Real R_c_(const int i) const { - const Real rm = xmin_[0] + i * dx_[0]; - const Real rp = rm + dx_[0]; - return 0.75*(std::pow(rp, 4) - std::pow(rm, 4)) / + const Real rm = Xf(i); + const Real rp = Xf(i+1); + return 0.75*(std::pow(rp, 4) - std::pow(rm, 4)) / (std::pow(rp, 3) - std::pow(rm, 3)); } KOKKOS_INLINE_FUNCTION Real X1s2_(const int i) const { - const Real rm = xmin_[0] + i * dx_[0]; - const Real rp = rm + dx_[0]; - return TWO_3RD*(std::pow(rp,3) - std::pow(rm,3))/(SQR(rp) - SQR(rm)); + const Real rm = Xf(i); + const Real rp = Xf(i+1); + return TWO_3RD*(std::pow(rp,3) - std::pow(rm,3))/(SQR(rp) - SQR(rm)); } KOKKOS_INLINE_FUNCTION Real Cos_theta_f_(const int j) const { - const Real theta = xmin_[1] + j * dx_[1]; - return std::cos(theta); + const Real theta = xmin_[1] + j * dx_[1]; + return std::cos(theta); } KOKKOS_INLINE_FUNCTION Real Sin_theta_f_(const int j) const { - const Real theta = xmin_[1] + j * dx_[1]; - return std::sin(theta); + const Real theta = xmin_[1] + j * dx_[1]; + return std::sin(theta); } KOKKOS_INLINE_FUNCTION Real Theta_c_(const int j) const { - Real tm = xmin_[1] + j * dx_[1]; - Real tp = tm + dx_[1]; - return (((Sin_theta_f_(j+1) - tp*Cos_theta_f_(j+1)) - - (Sin_theta_f_(j ) - tm*Cos_theta_f_(j )))/ - (Cos_theta_f_(j ) - Cos_theta_f_(j+1))); + Real tm = xmin_[1] + j * dx_[1]; + Real tp = tm + dx_[1]; + return (((Sin_theta_f_(j+1) - tp*Cos_theta_f_(j+1)) - + (Sin_theta_f_(j ) - tm*Cos_theta_f_(j )))/ + (Cos_theta_f_(j ) - Cos_theta_f_(j+1))); } KOKKOS_INLINE_FUNCTION Real Cos_theta_c_(const int j) const { - return std::cos(Theta_c_(j)); + return std::cos(Theta_c_(j)); } KOKKOS_INLINE_FUNCTION Real Sin_theta_c_(const int j) const { - return std::sin(Theta_c_(j)); + return std::sin(Theta_c_(j)); } KOKKOS_INLINE_FUNCTION Real Cos_phi_c_(const int k) const { - Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; - return std::cos(Phi); + Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; + return std::cos(Phi); } KOKKOS_INLINE_FUNCTION Real Sin_phi_c_(const int k) const { - Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; - return std::sin(Phi); + Real Phi = xmin_[2] + (k + 0.5) * dx_[2]; + return std::sin(Phi); + } + + KOKKOS_INLINE_FUNCTION Real Coord_area2_i_(const int i) const { + const Real rm = Xf(i); + const Real rp = Xf(i+1); + return 0.5*( SQR(rp) - SQR(rm) ); + } + KOKKOS_INLINE_FUNCTION Real Coord_area1_j_(const int j) const { + const Real cm = std::cos(Xf(j)); + const Real cp = std::cos(Xf(j+1)); + return std::abs(cm - cp); } }; From f88adb907afe0451162656bdd4e3a4d4ee8539ff Mon Sep 17 00:00:00 2001 From: Forrest Wolfgang Glines Date: Sun, 11 Feb 2024 15:41:22 -0700 Subject: [PATCH 9/9] GetXmin, GetStartIndex Inlined --- src/coordinates/uniform_cylindrical.hpp | 2 ++ src/coordinates/uniform_spherical.hpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/coordinates/uniform_cylindrical.hpp b/src/coordinates/uniform_cylindrical.hpp index ecfb547e56ab..85dc2f8e9693 100644 --- a/src/coordinates/uniform_cylindrical.hpp +++ b/src/coordinates/uniform_cylindrical.hpp @@ -489,7 +489,9 @@ class UniformCylindrical { return 1./( Xc<1>(i)*Xf<1>(i+1) ); } + KOKKOS_INLINE_FUNCTION const std::array &GetXmin() const { return xmin_; } + KOKKOS_INLINE_FUNCTION const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return name_; } diff --git a/src/coordinates/uniform_spherical.hpp b/src/coordinates/uniform_spherical.hpp index 4b3127195a31..7e9478c12334 100644 --- a/src/coordinates/uniform_spherical.hpp +++ b/src/coordinates/uniform_spherical.hpp @@ -540,7 +540,9 @@ class UniformSpherical { return Dxf(i)/( (rm + rp)*Coord_vol_i_(i)); } + KOKKOS_INLINE_FUNCTION const std::array &GetXmin() const { return xmin_; } + KOKKOS_INLINE_FUNCTION const std::array &GetStartIndex() const { return istart_; } const char *Name() const { return name_; }