diff --git a/src/output/load_diag_output_data.cpp b/src/output/load_diag_output_data.cpp index 320f79d0..ff909e3f 100644 --- a/src/output/load_diag_output_data.cpp +++ b/src/output/load_diag_output_data.cpp @@ -265,12 +265,14 @@ void OutputType::loadDiagOutputData(MeshBlockImpl* pmb, Variables const& vars) { auto area = pcoord->face_area1(); auto ny = peos->nvar() - 5; int il = pcoord->il(); + int iu = pcoord->iu(); if (ny > 0) { pod = new OutputData; pod->type = "VECTORS"; pod->name = get_hydro_names(pmb, "path_"); - auto u_sum = (u * vol).narrow(0, ICY, ny).sum(-1) / area.select(-1, il); + auto u_sum = (u * vol).narrow(0, ICY, ny).slice(-1, il, iu + 1).sum(-1) / + area.select(-1, il); pod->data.CopyFromTensor(u_sum); AppendOutputDataNode(pod); diff --git a/src/recon/cp3.cpp b/src/recon/cp3.cpp index 05c7a579..b1e9894a 100644 --- a/src/recon/cp3.cpp +++ b/src/recon/cp3.cpp @@ -11,7 +11,7 @@ namespace snap { void Center3InterpImpl::reset() { cm = register_buffer( - "cm", torch::tensor({-1. / 3., 5. / 6., -1. / 6.}, torch::kFloat64)); + "cm", torch::tensor({1. / 3., 5. / 6., -1. / 6.}, torch::kFloat64)); cp = register_buffer("cp", cm.flip({0})); } diff --git a/src/recon/interp_impl.cuh b/src/recon/interp_impl.cuh index 3b2f197d..f7d2306f 100644 --- a/src/recon/interp_impl.cuh +++ b/src/recon/interp_impl.cuh @@ -55,8 +55,8 @@ __device__ T interp_shared_weno3_coeff_impl(T const *line, T const *coeff, T beta0 = SQR(_vvdot<3>(phi, c3)); T beta1 = SQR(_vvdot<3>(phi, c4)); - T alpha0 = (1.0 / 3.0) / SQR(beta0 + 1e-6); - T alpha1 = (2.0 / 3.0) / SQR(beta1 + 1e-6); + T alpha0 = (2.0 / 3.0) / SQR(beta0 + 1e-6); + T alpha1 = (1.0 / 3.0) / SQR(beta1 + 1e-6); return (alpha0 * p0 + alpha1 * p1) / (alpha0 + alpha1) * vscale; } @@ -109,7 +109,7 @@ __device__ T interp_shared_fused_impl(T const *line, int v, int start, int axis_size, FusedReconScheme scheme, bool right, bool scale) { if (scheme == FusedReconScheme::CP3) { - constexpr T cm[3] = {-1. / 3., 5. / 6., -1. / 6.}; + constexpr T cm[3] = {1. / 3., 5. / 6., -1. / 6.}; T c[3]; for (int k = 0; k < 3; ++k) c[k] = right ? cm[2 - k] : cm[k]; diff --git a/src/recon/interp_impl.h b/src/recon/interp_impl.h index 68fd0c2a..c920df25 100644 --- a/src/recon/interp_impl.h +++ b/src/recon/interp_impl.h @@ -10,7 +10,7 @@ namespace snap { template -DISPATCH_MACRO inline T _vvdot(T *v1, T *v2) { +DISPATCH_MACRO inline T _vvdot(T* v1, T* v2) { T out = 0.; for (int i = 0; i < N; ++i) { out += v1[i] * v2[i]; @@ -20,7 +20,7 @@ DISPATCH_MACRO inline T _vvdot(T *v1, T *v2) { // polynomial template -DISPATCH_MACRO void interp_poly_impl(T *out, T *inp, T *coeff, int stride1, +DISPATCH_MACRO void interp_poly_impl(T* out, T* inp, T* coeff, int stride1, int stride2, int stride_out, int nvar) { for (int j = 0; j < nvar; ++j) { OUT(j) = 0.; @@ -32,13 +32,13 @@ DISPATCH_MACRO void interp_poly_impl(T *out, T *inp, T *coeff, int stride1, // WENO 3 interpolation template -DISPATCH_MACRO void interp_weno3_impl(T *out, T *inp, T *coeff, int stride1, +DISPATCH_MACRO void interp_weno3_impl(T* out, T* inp, T* coeff, int stride1, int stride2, int stride_out, int nvar, double scale) { - T *c1 = coeff; - T *c2 = c1 + 3; - T *c3 = c2 + 3; - T *c4 = c3 + 3; + T* c1 = coeff; + T* c2 = c1 + 3; + T* c3 = c2 + 3; + T* c4 = c3 + 3; T phi[3]; @@ -56,20 +56,14 @@ DISPATCH_MACRO void interp_weno3_impl(T *out, T *inp, T *coeff, int stride1, continue; } - if (vscale != 0.0) { - phi[0] /= vscale; - phi[1] /= vscale; - phi[2] /= vscale; - } - T p0 = _vvdot<3>(phi, c1); T p1 = _vvdot<3>(phi, c2); T beta0 = SQR(_vvdot<3>(phi, c3)); T beta1 = SQR(_vvdot<3>(phi, c4)); - T alpha0 = (1.0 / 3.0) / SQR(beta0 + 1e-6); - T alpha1 = (2.0 / 3.0) / SQR(beta1 + 1e-6); + T alpha0 = (2.0 / 3.0) / SQR(beta0 + 1e-6); + T alpha1 = (1.0 / 3.0) / SQR(beta1 + 1e-6); OUT(j) = (alpha0 * p0 + alpha1 * p1) / (alpha0 + alpha1) * vscale; } @@ -77,18 +71,18 @@ DISPATCH_MACRO void interp_weno3_impl(T *out, T *inp, T *coeff, int stride1, // WENO 5 interpolation template -DISPATCH_MACRO void interp_weno5_impl(T *out, T *inp, T *coeff, int stride1, +DISPATCH_MACRO void interp_weno5_impl(T* out, T* inp, T* coeff, int stride1, int stride2, int stride_out, int nvar, double scale) { - T *c1 = coeff; - T *c2 = c1 + 5; - T *c3 = c2 + 5; - T *c4 = c3 + 5; - T *c5 = c4 + 5; - T *c6 = c5 + 5; - T *c7 = c6 + 5; - T *c8 = c7 + 5; - T *c9 = c8 + 5; + T* c1 = coeff; + T* c2 = c1 + 5; + T* c3 = c2 + 5; + T* c4 = c3 + 5; + T* c5 = c4 + 5; + T* c6 = c5 + 5; + T* c7 = c6 + 5; + T* c8 = c7 + 5; + T* c9 = c8 + 5; T phi[5]; diff --git a/src/recon/interp_simple.hpp b/src/recon/interp_simple.hpp index 0ccf1a4d..4a68d0a2 100644 --- a/src/recon/interp_simple.hpp +++ b/src/recon/interp_simple.hpp @@ -99,8 +99,8 @@ inline DISPATCH_MACRO T interp_weno3(T const& phim1, T const& phi, T beta0 = (phim1 - phi) * (phim1 - phi); T beta1 = (phi - phip1) * (phi - phip1); - T alpha0 = (1.0 / 3.0) / sqr(beta0 + 1e-6); - T alpha1 = (2.0 / 3.0) / sqr(beta1 + 1e-6); + T alpha0 = (2.0 / 3.0) / sqr(beta0 + 1e-6); + T alpha1 = (1.0 / 3.0) / sqr(beta1 + 1e-6); return (alpha0 * p0 + alpha1 * p1) / (alpha0 + alpha1); }; diff --git a/src/recon/recon_dispatch.cpp b/src/recon/recon_dispatch.cpp index e4624811..0afdf56c 100644 --- a/src/recon/recon_dispatch.cpp +++ b/src/recon/recon_dispatch.cpp @@ -87,8 +87,8 @@ void call_weno3_mps(at::TensorIterator& iter, torch::Tensor coeff, int dim, wu /= wscale.unsqueeze(-1); } - auto alpha1 = 1. / 3. / (wu.matmul(c3).square() + 1e-6).square(); - auto alpha2 = 2. / 3. / (wu.matmul(c4).square() + 1e-6).square(); + auto alpha1 = 2. / 3. / (wu.matmul(c3).square() + 1e-6).square(); + auto alpha2 = 1. / 3. / (wu.matmul(c4).square() + 1e-6).square(); torch::add_out(result, alpha1 * wu.matmul(c1), alpha2 * wu.matmul(c2)); result /= alpha1 + alpha2; diff --git a/tests/run_shallow_splash.cmake b/tests/run_shallow_splash.cmake index f810c7ca..d7147bf1 100644 --- a/tests/run_shallow_splash.cmake +++ b/tests/run_shallow_splash.cmake @@ -17,7 +17,13 @@ if(NOT _status EQUAL 0) message(FATAL_ERROR "Failed to download reference file with exit code ${_status}") endif() -execute_process(COMMAND ln -sf ../bin/shallow_splash.yaml shallow_splash.yaml) +file(READ "../bin/shallow_splash.yaml" config) +file(READ "../configure.h" configure_h) +if(configure_h MATCHES "NO_PNETCDFOUTPUT") + string(REPLACE "type: pnetcdf" "type: netcdf" config "${config}") +endif() +file(REMOVE "shallow_splash.yaml") +file(WRITE "shallow_splash.yaml" "${config}") execute_process( COMMAND torchrun --no-python --nproc-per-node=6 ../bin/shallow_splash.${buildl} diff --git a/tests/run_straka.cmake b/tests/run_straka.cmake index c832222f..362c0ee6 100644 --- a/tests/run_straka.cmake +++ b/tests/run_straka.cmake @@ -17,7 +17,13 @@ if(NOT _status EQUAL 0) message(FATAL_ERROR "Failed to download reference file with exit code ${_status}") endif() -execute_process(COMMAND ln -sf ../bin/straka.yaml straka.yaml) +file(READ "../bin/straka.yaml" config) +file(READ "../configure.h" configure_h) +if(configure_h MATCHES "NO_PNETCDFOUTPUT") + string(REPLACE "type: pnetcdf" "type: netcdf" config "${config}") +endif() +file(REMOVE "straka.yaml") +file(WRITE "straka.yaml" "${config}") execute_process( COMMAND torchrun --no-python --nproc-per-node=2 ../bin/straka.${buildl} diff --git a/tests/test_fused_cs_sync_smoke.cu b/tests/test_fused_cs_sync_smoke.cu index aaf45bad..7e11135f 100644 --- a/tests/test_fused_cs_sync_smoke.cu +++ b/tests/test_fused_cs_sync_smoke.cu @@ -47,20 +47,6 @@ torch::Device select_cuda_device() { return torch::Device(torch::kCUDA, local_rank); } -void require_cuda_6rank_or_skip() { - if (!torch::cuda::is_available()) { - GTEST_SKIP() << "CUDA runtime is unavailable"; - } - int world_size = env_int("WORLD_SIZE", 1); - if (world_size != 6) { - GTEST_SKIP() << "test_fused_cs_sync_smoke requires torchrun with 6 ranks"; - } - int device_count = c10::cuda::device_count(); - if (device_count < world_size) { - GTEST_SKIP() << "test_fused_cs_sync_smoke requires 6 CUDA devices"; - } -} - LayoutOptions make_layout_options() { auto opts = LayoutOptionsImpl::create(); opts->type("cubed-sphere"); @@ -134,6 +120,21 @@ bool fused_exchange_uses_right_interp_for_side(int side) { return !fused_exchange_uses_right_state_start_for_side(side); } +bool cuda_6rank_available_or_skip() { + if (!torch::cuda::is_available()) { + return false; + } + int world_size = env_int("WORLD_SIZE", 1); + if (world_size != 6) { + return false; + } + int device_count = c10::cuda::device_count(); + if (device_count < world_size) { + return false; + } + return true; +} + void initialize_symmetric_memory_group(LayoutImpl const &layout, std::string const &group_name) { static std::set initialized; @@ -318,7 +319,8 @@ __global__ void verify_hydro_remote_constant_kernel(void **buffer_ptrs, } // namespace TEST(FusedCubedSphereSymmetricMemory, RendezvousCompletes) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); auto ctx = make_smoke_context(); std::string group_name = "snapy:test:fused-cs-rendezvous"; initialize_symmetric_memory_group(*ctx.layout, group_name); @@ -331,7 +333,8 @@ TEST(FusedCubedSphereSymmetricMemory, RendezvousCompletes) { } TEST(FusedCubedSphereSymmetricMemory, PreviousKernelSyncCompletes) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); auto ctx = make_smoke_context(); std::string group_name = "snapy:test:fused-cs-sync"; initialize_symmetric_memory_group(*ctx.layout, group_name); @@ -357,7 +360,8 @@ TEST(FusedCubedSphereSymmetricMemory, PreviousKernelSyncCompletes) { } TEST(FusedCubedSphereSymmetricMemory, OrientationMetadataMatchesStaged) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); auto ctx = make_smoke_context(); auto meta = ctx.edge_meta.cpu(); @@ -378,7 +382,8 @@ TEST(FusedCubedSphereSymmetricMemory, OrientationMetadataMatchesStaged) { } TEST(FusedCubedSphereSymmetricMemory, BoundaryStateConventionMatchesStaged) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); // Staged hydro sends the local right state to lower-side neighbors and the // local left state to upper-side neighbors. EXPECT_TRUE(fused_exchange_uses_right_state_start_for_side(SIDE_L)); @@ -392,7 +397,8 @@ TEST(FusedCubedSphereSymmetricMemory, BoundaryStateConventionMatchesStaged) { } TEST(FusedCubedSphereSymmetricMemory, RemoteEdgePayloadMatches) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); auto ctx = make_smoke_context(); std::string group_name = "snapy:test:fused-cs-remote-read"; initialize_symmetric_memory_group(*ctx.layout, group_name); @@ -427,7 +433,8 @@ TEST(FusedCubedSphereSymmetricMemory, RemoteEdgePayloadMatches) { } TEST(FusedCubedSphereSymmetricMemory, RemoteStateSelectionMatchesStaged) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); auto ctx = make_smoke_context(); std::string group_name = "snapy:test:fused-cs-remote-state"; initialize_symmetric_memory_group(*ctx.layout, group_name); @@ -463,7 +470,8 @@ TEST(FusedCubedSphereSymmetricMemory, RemoteStateSelectionMatchesStaged) { } TEST(FusedCubedSphereSymmetricMemory, ConstantStateExchangeHasZeroMassFlux) { - require_cuda_6rank_or_skip(); + if (!cuda_6rank_available_or_skip()) + GTEST_SKIP(); auto ctx = make_smoke_context(); std::string group_name = "snapy:test:fused-cs-constant-state"; initialize_symmetric_memory_group(*ctx.layout, group_name); diff --git a/tests/test_riemann.cpp b/tests/test_riemann.cpp index a57f91ce..bb54c4e8 100644 --- a/tests/test_riemann.cpp +++ b/tests/test_riemann.cpp @@ -27,7 +27,7 @@ enum { DIM3 = 1, }; -const char *block_config = R"( +const char* block_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -66,7 +66,7 @@ boundary-condition: x3-outer: reflecting )"; -const char *small_ideal_gas_config = R"( +const char* small_ideal_gas_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -106,7 +106,7 @@ boundary-condition: x3-outer: reflecting )"; -const char *small_ideal_gas_gravity_config = R"( +const char* small_ideal_gas_gravity_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -151,7 +151,7 @@ boundary-condition: x3-outer: reflecting )"; -const char *small_ideal_gas_cubed_sphere_hllc_config = R"( +const char* small_ideal_gas_cubed_sphere_hllc_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -202,7 +202,7 @@ boundary-condition: x3-outer: custom )"; -const char *small_ideal_gas_implicit_config = R"( +const char* small_ideal_gas_implicit_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -250,7 +250,7 @@ boundary-condition: x3-outer: reflecting )"; -const char *small_ideal_moist_sedimentation_config = R"( +const char* small_ideal_moist_sedimentation_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -319,7 +319,7 @@ boundary-condition: x3-outer: reflecting )"; -const char *small_ideal_moist_lmars_cloud_config = R"( +const char* small_ideal_moist_lmars_cloud_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -377,7 +377,7 @@ boundary-condition: x3-outer: reflecting )"; -const char *small_shallow_water_config = R"( +const char* small_shallow_water_config = R"( dynamics: equation-of-state: type: shallow-water @@ -411,9 +411,15 @@ using namespace snap; namespace { +int env_int(char const* name, int fallback) { + char const* value = std::getenv(name); + if (value == nullptr || value[0] == '\0') return fallback; + return std::stoi(value); +} + struct ScopedEnv { - explicit ScopedEnv(char const *name_) : name(name_) { - auto *value = std::getenv(name); + explicit ScopedEnv(char const* name_) : name(name_) { + auto* value = std::getenv(name); if (value) { had_value = true; old_value = value; @@ -427,7 +433,7 @@ struct ScopedEnv { } } - char const *name; + char const* name; bool had_value = false; std::string old_value; }; @@ -984,7 +990,7 @@ TEST_P(DeviceTest, fused_recon_riemann_honors_weno_scale_and_shock) { GTEST_SKIP() << "fused reconstruction/Riemann path is CUDA-only"; } - for (auto const &options : {std::pair{true, false}, std::pair{false, true}, + for (auto const& options : {std::pair{true, false}, std::pair{false, true}, std::pair{true, true}}) { SCOPED_TRACE(testing::Message() << "scale=" << options.first << " shock=" << options.second); @@ -1051,6 +1057,10 @@ TEST_P(DeviceTest, fused_recon_riemann_matches_staged_cubed_sphere_hllc) { if (device.type() != torch::kCUDA) { GTEST_SKIP() << "fused reconstruction/Riemann path is CUDA-only"; } + if (env_int("WORLD_SIZE", 1) != 6) { + GTEST_SKIP() + << "cubed-sphere fused reconstruction/Riemann test requires 6 ranks"; + } char staged_name[80] = "/tmp/tempfile.XXXXXX"; mkstemp(staged_name); @@ -1476,7 +1486,7 @@ TEST_P(DeviceTest, test_hllc) { std::remove(fname); } -int main(int argc, char **argv) { +int main(int argc, char** argv) { testing::InitGoogleTest(&argc, argv); return RUN_ALL_TESTS(); diff --git a/tests/test_weno.cpp b/tests/test_weno.cpp index a53fd5a2..507aa555 100644 --- a/tests/test_weno.cpp +++ b/tests/test_weno.cpp @@ -151,6 +151,27 @@ TEST_P(DeviceTest, interp_weno3a) { EXPECT_NEAR(result, expected_result, 1.E-10); } +TEST_P(DeviceTest, interp_weno3_smooth_limit) { + Weno3Interp interp; + interp->to(device, torch::kFloat64); + + constexpr double eps = 1.e-7; + auto phi = torch::tensor({1. + eps, 1. + 2. * eps, 1. + 4. * eps}, + torch::device(device).dtype(torch::kFloat64)); + + auto left_expected = interp_cp3(phi[0].item(), phi[1].item(), + phi[2].item()); + auto left_result = torch::zeros({1}, phi.options()); + interp->left(phi, 0, left_result); + EXPECT_NEAR(left_expected, left_result.item(), 1.E-6); + + auto right_expected = interp_cp3(phi[2].item(), phi[1].item(), + phi[0].item()); + auto right_result = torch::zeros({1}, phi.options()); + interp->right(phi, 0, right_result); + EXPECT_NEAR(right_expected, right_result.item(), 1.E-6); +} + TEST_P(DeviceTest, interp_weno5b) { double phim2 = 1.0; double phim1 = 2.0; @@ -633,7 +654,7 @@ TEST(weno5, interp_weno5b) { std::cout << "1) " << result << std::endl; } -int main(int argc, char **argv) { +int main(int argc, char** argv) { testing::InitGoogleTest(&argc, argv); return RUN_ALL_TESTS();