Fix pressure grad

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;
     }

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];
     }

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)

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();

tests/test_hydrostatic.cpp

@@ -0,0 +1,168 @@
+// C/C++
+#include <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <string>
+
+// POSIX
+#include <unistd.h>
+
+// gtest
+#include <gtest/gtest.h>
+
+// snap
+#include <snap/snap.h>
+
+#include <snap/mesh/mesh.hpp>
+
+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<double>();
+      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);
+}