wip

chengcli · Apr 11, 2024 · 8659446 · 8659446
1 parent 6bb06f4
commit 8659446
Show file tree

Hide file tree

Showing 5 changed files with 23 additions and 94 deletions.
diff --git a/src/diagnostics/diagnostics.hpp b/src/diagnostics/diagnostics.hpp
@@ -68,7 +68,7 @@ class Divergence : public Diagnostics {
   virtual ~Divergence() {}
 
   void Finalize(MeshBlock *pmb) override;
-  int GetNumVars() const override { return 1; }
+  int GetNumVars() const override { return 2; }
 
  protected:
   AthenaArray<Real> v1f1_, v2f2_, v3f3_;
@@ -115,19 +115,6 @@ class HydroMean : public Diagnostics {
   int ncycle_;
 };
 
-// 5. horizontal divergence
-class HorizontalDivergence : public Diagnostics {
- public:
-  HorizontalDivergence(MeshBlock *pmb);
-  virtual ~HorizontalDivergence() {}
-
-  void Finalize(MeshBlock *pmb) override;
-  int GetNumVars() const override { return 1; }
-
- protected:
-  AthenaArray<Real> v2f2_, v3f3_;
-};
-
 // 6. anomaly
 class Anomaly : public Diagnostics {
  public:

diff --git a/src/diagnostics/diagnostics_factory.cpp b/src/diagnostics/diagnostics_factory.cpp
@@ -26,15 +26,13 @@ DiagnosticsContainer DiagnosticsFactory::CreateFrom(MeshBlock *pmb,
       diag.push_back(std::make_shared<Buoyancy>(pmb));
     } else if (name == "mean") {  // 4.
       diag.push_back(std::make_shared<HydroMean>(pmb));
-    } else if (name == "div_h") {  // 5.
-      diag.push_back(std::make_shared<HorizontalDivergence>(pmb));
-    } else if (name == "anomaly") {  // 6.
+    } else if (name == "anomaly") {  // 5.
       diag.push_back(std::make_shared<Anomaly>(pmb));
-    } else if (name == "radflux") {  // 8.
+    } else if (name == "radflux") {  // 6.
       diag.push_back(std::make_shared<RadiativeFlux>(pmb));
-    } else if (name == "hydroflux") {  // 9.
+    } else if (name == "hydroflux") {  // 7.
       diag.push_back(std::make_shared<HydroFlux>(pmb));
-    } else if (name == "w_avg") {  // 10.
+    } else if (name == "w_avg") {  // 8.
       diag.push_back(std::make_shared<V1Moments>(pmb));
       /*} else if (name == "eddyflux") { // 6.
         diag.push_back(std::make_shared<EddyFlux>(pmb));

diff --git a/src/diagnostics/divergence.cpp b/src/diagnostics/divergence.cpp
@@ -6,10 +6,10 @@
 // canoe
 #include "diagnostics.hpp"
 
-Divergence::Divergence(MeshBlock *pmb) : Diagnostics(pmb, "div") {
-  type = "SCALARS";
+Divergence::Divergence(MeshBlock *pmb) : Diagnostics(pmb, "div,div_h") {
+  type = "VECTORS";
 
-  data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
+  data.NewAthenaArray(2, ncells3_, ncells2_, ncells1_);
   v1f1_.NewAthenaArray(ncells3_, ncells2_, ncells1_ + 1);
 
   if (pmb->block_size.nx2 > 1) {
@@ -56,24 +56,29 @@ void Divergence::Finalize(MeshBlock *pmb) {
       pcoord->CellVolume(k, j, is, ie, vol_);
       for (int i = is; i <= ie; ++i) {
         pcoord->Face1Area(k, j, is, ie + 1, x1area_);
-        data(k, j, i) =
+        data(0, k, j, i) =
             x1area_(i + 1) * v1f1_(k, j, i + 1) - x1area_(i) * v1f1_(k, j, i);
 
         if (pmb->block_size.nx2 > 1) {
           pcoord->Face2Area(k, j, is, ie, x2area_);
           pcoord->Face2Area(k, j + 1, is, ie, x2area_p1_);
-          data(k, j, i) +=
+          data(0, k, j, i) +=
+              x2area_p1_(i) * v2f2_(k, j + 1, i) - x2area_(i) * v2f2_(k, j, i);
+          data(1, k, j, i) =
               x2area_p1_(i) * v2f2_(k, j + 1, i) - x2area_(i) * v2f2_(k, j, i);
         }
 
         if (pmb->block_size.nx3 > 1) {
           pcoord->Face3Area(k, j, is, ie, x3area_);
           pcoord->Face3Area(k + 1, j, is, ie, x3area_p1_);
-          data(k, j, i) +=
+          data(0, k, j, i) +=
+              x3area_p1_(i) * v3f3_(k + 1, j, i) - x3area_(i) * v3f3_(k, j, i);
+          data(1, k, j, i) +=
               x3area_p1_(i) * v3f3_(k + 1, j, i) - x3area_(i) * v3f3_(k, j, i);
         }
 
-        data(k, j, i) /= vol_(i);
+        data(0, k, j, i) /= vol_(i);
+        data(1, k, j, i) /= vol_(i);
       }
     }
 }
diff --git a/src/diagnostics/horizontal_divergence.cpp b/src/diagnostics/horizontal_divergence.cpp
diff --git a/src/outputs/output_utils.cpp b/src/outputs/output_utils.cpp
@@ -36,7 +36,9 @@ __attribute__((weak)) MetadataTable::MetadataTable() {
       {"x3f", "distance at cell boundary", "m", "F--"},
       {"rho", "density", "kg/m^3", "CCC"},
       {"press", "pressure", "pa", "CCC"},
-      {"vel", "velocity", "m/s", "CCC"},
+      {"vel1", "vertical velocity", "m/s", "CCC"},
+      {"vel2", "horizontal velocity", "m/s", "CCC"},
+      {"vel3", "horizontal velocity", "m/s", "CCC"},
       {"vapor", "mass mixing ratio of vapor", "kg/kg", "CCC"},
       {"temp", "temperature", "K", "CCC"},
       {"theta", "potential temperature", "K", "CCC"},
@@ -51,7 +53,9 @@ __attribute__((weak)) MetadataTable::MetadataTable() {
       // radiation
       {"radiance", "top-of-atmosphere radiance", "K", "RCC"},
       // curl
-      {"curl", "curl", "1/s", "CCC"},
+      {"curl1", "curl in the vertical direction", "1/s", "CCC"},
+      {"curl2", "curl in the horizontal direction", "1/s", "CCC"},
+      {"curl3", "curl in the horizontal direction", "1/s", "CCC"},
       // divergence
       {"div", "divergence", "1/s", "CCC"},
       {"div_h", "horizontal divergence", "1/s", "CCC"},