diff --git a/src/coord/gnomonic_equiangle.cpp b/src/coord/gnomonic_equiangle.cpp index cecf38c4..72394703 100644 --- a/src/coord/gnomonic_equiangle.cpp +++ b/src/coord/gnomonic_equiangle.cpp @@ -384,11 +384,15 @@ torch::Tensor GnomonicEquiangleImpl::forward(torch::Tensor prim, if (face_pressure1.defined()) { int si = il(); int ei = iu() + 1; + auto pressure_gradient = (face_pressure1.slice(-1, si + 1, ei + 1) - + face_pressure1.slice(-1, si, ei)) / + dx1f.slice(0, si, ei); src1.slice(-1, si, ei) += (CoordinateImpl::face_area1(si + 1, ei + 1) * face_pressure1.slice(-1, si + 1, ei + 1) - CoordinateImpl::face_area1(si, ei) * face_pressure1.slice(-1, si, ei)) / - cell_volume().slice(-1, si, ei); + cell_volume().slice(-1, si, ei) - + pressure_gradient; } else { src1 += 2.0 * pr / radius; } diff --git a/src/coord/spherical_polar.cpp b/src/coord/spherical_polar.cpp index caabb360..6e8a3ece 100644 --- a/src/coord/spherical_polar.cpp +++ b/src/coord/spherical_polar.cpp @@ -233,11 +233,15 @@ torch::Tensor SphericalPolarImpl::forward(torch::Tensor prim, coord_src1_i * prim[IDN] * (prim[IVY].square() + prim[IVZ].square()); if (eos_type != "shallow-water") { if (face_pressure1.defined()) { + auto pressure_gradient = (face_pressure1.slice(-1, si + 1, ei + 1) - + face_pressure1.slice(-1, si, ei)) / + dx1f.slice(0, si, ei); src1.slice(-1, si, ei) += (CoordinateImpl::face_area1(si + 1, ei + 1) * face_pressure1.slice(-1, si + 1, ei + 1) - CoordinateImpl::face_area1(si, ei) * face_pressure1.slice(-1, si, ei)) / - vol.slice(-1, si, ei); + vol.slice(-1, si, ei) - + pressure_gradient; } else { src1 += 2.0 * coord_src1_i * prim[IPR]; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 21055086..9ff6d6af 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -19,6 +19,7 @@ setup_test(test_slab) setup_test(test_cubed) setup_test(test_cubed_sphere) setup_test(test_coordinate) +setup_test(test_hydrostatic) setup_test(test_output_dir) setup_test(test_netcdf_utils) setup_test(test_user_output) diff --git a/tests/test_coordinate.cpp b/tests/test_coordinate.cpp index d6297963..4343db64 100644 --- a/tests/test_coordinate.cpp +++ b/tests/test_coordinate.cpp @@ -27,7 +27,7 @@ using namespace snap; namespace { -const char *gnomonic_radial_config = R"( +const char* gnomonic_radial_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -63,7 +63,7 @@ boundary-condition: x3-outer: custom )"; -const char *spherical_polar_config = R"( +const char* spherical_polar_config = R"( reference-state: Tref: 300. Pref: 1.e5 @@ -105,7 +105,7 @@ boundary-condition: x3-outer: periodic )"; -std::string write_temp_config(char const *config) { +std::string write_temp_config(char const* config) { char fname[] = "/tmp/test-coordinate-XXXXXX"; int fd = mkstemp(fname); EXPECT_NE(fd, -1); @@ -553,6 +553,78 @@ TEST_P(DeviceTest, radial_source_uses_face_pressure_in_x1_momentum) { EXPECT_TRUE(torch::allclose(div_lo[IVX], div_hi[IVX], 1.e-8, 1.e-8)); } +TEST_P(DeviceTest, radial_source_preserves_face_pressure_gradient) { + auto fname = write_temp_config(gnomonic_radial_config); + auto op = MeshBlockOptionsImpl::from_yaml(fname); + auto block = MeshBlock(op); + block->to(device, dtype); + std::remove(fname.c_str()); + + auto pcoord = block->pcoord; + int nc1 = pcoord->options->nc1(); + int nc2 = pcoord->options->nc2(); + int nc3 = pcoord->options->nc3(); + int si = pcoord->il(); + int ei = pcoord->iu() + 1; + auto opts = torch::TensorOptions().dtype(dtype).device(device); + + auto prim = torch::zeros({5, nc3, nc2, nc1}, opts); + prim[IDN].fill_(1.0); + auto face_pressure = pcoord->x1f.slice(0, 0, nc1) + .to(opts) + .view({1, 1, nc1}) + .expand({nc3, nc2, nc1}); + auto flux1 = torch::zeros_like(prim); + flux1[IVX].copy_(face_pressure); + + auto div = pcoord->forward(prim, flux1, torch::Tensor(), torch::Tensor(), + face_pressure); + auto radial_div = div[IVX].slice(-1, si, ei); + auto expected = (face_pressure.slice(-1, si + 1, ei + 1) - + face_pressure.slice(-1, si, ei)) / + pcoord->dx1f.slice(0, si, ei); + + EXPECT_TRUE(torch::allclose(radial_div, expected, 1.e-6, 1.e-6)) + << "radial_div=" << radial_div << "\nexpected=" << expected; +} + +TEST_P(DeviceTest, + spherical_polar_radial_source_preserves_face_pressure_gradient) { + auto fname = write_temp_config(spherical_polar_config); + auto op = MeshBlockOptionsImpl::from_yaml(fname); + auto block = MeshBlock(op); + block->to(device, dtype); + std::remove(fname.c_str()); + + auto pcoord = block->pcoord; + int nc1 = pcoord->options->nc1(); + int nc2 = pcoord->options->nc2(); + int nc3 = pcoord->options->nc3(); + int si = pcoord->il(); + int ei = pcoord->iu() + 1; + auto opts = torch::TensorOptions().dtype(dtype).device(device); + + auto prim = torch::zeros({5, nc3, nc2, nc1}, opts); + prim[IDN].fill_(1.0); + auto face_pressure = pcoord->x1f.slice(0, 0, nc1) + .to(opts) + .view({1, 1, nc1}) + .expand({nc3, nc2, nc1}); + auto flux1 = torch::zeros_like(prim); + flux1[IVX].copy_(face_pressure); + auto flux2 = torch::zeros_like(prim); + auto flux3 = torch::zeros_like(prim); + + auto div = pcoord->forward(prim, flux1, flux2, flux3, face_pressure); + auto radial_div = div[IVX].slice(-1, si, ei); + auto expected = (face_pressure.slice(-1, si + 1, ei + 1) - + face_pressure.slice(-1, si, ei)) / + pcoord->dx1f.slice(0, si, ei); + + EXPECT_TRUE(torch::allclose(radial_div, expected, 1.e-6, 1.e-6)) + << "radial_div=" << radial_div << "\nexpected=" << expected; +} + TEST_P(DeviceTest, spherical_polar_radial_source_uses_face_pressure_in_x1_momentum) { auto fname = write_temp_config(spherical_polar_config); @@ -592,7 +664,7 @@ TEST_P(DeviceTest, torch::allclose(div_face_lo[IVX], div_face_hi[IVX], 1.e-6, 1.e-6)); } -int main(int argc, char **argv) { +int main(int argc, char** argv) { testing::InitGoogleTest(&argc, argv); int result = RUN_ALL_TESTS(); diff --git a/tests/test_hydrostatic.cpp b/tests/test_hydrostatic.cpp new file mode 100644 index 00000000..ef449aec --- /dev/null +++ b/tests/test_hydrostatic.cpp @@ -0,0 +1,168 @@ +// C/C++ +#include +#include +#include +#include +#include +#include +#include + +// POSIX +#include + +// gtest +#include + +// snap +#include + +#include + +using namespace snap; + +namespace { + +const char* cubed_sphere_hydrostatic_config = R"( +reference-state: + Tref: 300. + Pref: 100. + +species: + - name: dry + composition: {O: 0.42, N: 1.56, Ar: 0.01} + cv_R: 2.5 + +dynamics: + equation-of-state: + type: ideal-gas + gammad: 1.4 + density-floor: 1.e-12 + pressure-floor: 1.e-12 + limiter: false + + reconstruct: + vertical: {type: weno5, scale: false, shock: false} + horizontal: {type: weno5, scale: false, shock: false} + + riemann-solver: + type: lmars + +integration: + type: rk3 + cfl: 0.4 + implicit-scheme: 0 + +forcing: + const-gravity: + grav1: -1. + non-hydrostatic: 0. + +distribute: + layout: cubed-sphere + nb2: 1 + nb3: 1 + blocks_per_process: 6 + verbose: false + +geometry: + type: gnomonic-equiangle + cells: {nx1: 24, nx2: 8, nx3: 8, nghost: 3} + bounds: + x1min: 10. + x1max: 11. + x2min_pi: -0.25 + x2max_pi: 0.25 + x3min_pi: -0.25 + x3max_pi: 0.25 + +boundary-condition: + external: + x1-inner: reflecting + x1-outer: reflecting + x2-inner: custom + x2-outer: custom + x3-inner: custom + x3-outer: custom +)"; + +std::string write_temp_config() { + char fname[] = "/tmp/test-hydrostatic-XXXXXX"; + int fd = mkstemp(fname); + EXPECT_NE(fd, -1); + if (fd != -1) close(fd); + + std::ofstream outfile(fname); + outfile << cubed_sphere_hydrostatic_config; + outfile.close(); + return fname; +} + +void initialize_hydrostatic_column(MeshBlock const& block, Variables& vars) { + constexpr double gamma = 1.4; + constexpr double pressure0 = 100.; + constexpr double density0 = 1.; + constexpr double gravity = 1.; + constexpr double radius0 = 10.; + + auto pcoord = block->pcoord; + auto entropy = pressure0 / std::pow(density0, gamma); + auto exponent = (gamma - 1.) / gamma; + auto pressure_power = std::pow(pressure0, exponent) - + exponent * gravity * (pcoord->x1v - radius0) / + std::pow(entropy, 1. / gamma); + auto pressure = pressure_power.pow(1. / exponent); + auto density = (pressure / entropy).pow(1. / gamma); + + auto nc1 = pcoord->options->nc1(); + auto nc2 = pcoord->options->nc2(); + auto nc3 = pcoord->options->nc3(); + auto w = torch::zeros({block->phydro->peos->nvar(), nc3, nc2, nc1}, + torch::kFloat64); + w[IDN].copy_(density.view({1, 1, nc1}).expand({nc3, nc2, nc1})); + w[IPR].copy_(pressure.view({1, 1, nc1}).expand({nc3, nc2, nc1})); + vars["hydro_w"] = w; +} + +} // namespace + +TEST(HydrostaticAtmosphere, cubed_sphere_remains_at_rest) { + setenv("FUSED", "OFF", 1); + torch::set_num_threads(1); + torch::set_num_interop_threads(1); + + auto fname = write_temp_config(); + auto mesh = Mesh(MeshOptionsImpl::from_yaml(fname)); + std::remove(fname.c_str()); + mesh->to(torch::kCPU, torch::kFloat64); + ASSERT_EQ(mesh->blocks.size(), 6); + + MeshVariables vars(mesh->blocks.size()); + for (size_t i = 0; i < mesh->blocks.size(); ++i) { + initialize_hydrostatic_column(mesh->blocks[i], vars[i]); + } + mesh->initialize(vars); + + constexpr int nsteps = 100; + double max_vertical_velocity = 0.; + for (int step = 0; step < nsteps; ++step) { + auto dt = mesh->max_time_step(vars); + for (int stage = 0; stage < mesh->blocks.front()->pintg->stages.size(); + ++stage) { + mesh->forward(vars, dt, stage); + } + + for (size_t i = 0; i < mesh->blocks.size(); ++i) { + auto interior = mesh->blocks[i]->part( + {0, 0, 0}, PartOptions().exterior(false).ndim(3)); + auto vertical_velocity = + vars[i].at("hydro_w")[IVX].index(interior).abs().max().item(); + max_vertical_velocity = + std::max(max_vertical_velocity, vertical_velocity); + } + } + + std::cout << "maximum spurious vertical velocity after " << nsteps + << " steps: " << max_vertical_velocity << std::endl; + EXPECT_TRUE(std::isfinite(max_vertical_velocity)); + EXPECT_LT(max_vertical_velocity, 1.e-8); +}