wip

examples/2024-YZhou-jupiter/jupiter_dry.inp

@@ -6,7 +6,7 @@ problem_id  = jupiter2d     # problem ID: basename of output filenames
 
 <output0>
 file_type   = rst
-dt          = 100.E5
+dt          = 10.E5
 
 <output1>
 file_type   = hst           # History data dump
@@ -35,33 +35,33 @@ dt          = 4.E4          # time increment
 <time>
 cfl_number  = 0.9
 nlim        = -1            # cycle limit
-tlim        = 200.E5
+tlim        = 1000.E5
 xorder      = 5             # horizontal reconstruction order
 integrator  = rk3           # integration method
 
 <mesh>
-nx1         = 80          # Number of zones in X1-direction
+nx1         = 40          # Number of zones in X1-direction
 x1min       = -320.E3     # minimum value of X1
 x1max       = 40.E3       # maximum value of X1
 ix1_bc      = reflecting  # Inner-X1 boundary condition flag
 ox1_bc      = reflecting  # Outer-X1 boundary condition flag
 
-nx2         = 320         # Number of zones in X2-direction
+nx2         = 80         # Number of zones in X2-direction
 x2min       = 0.          # minimum value of X2
-x2max       = 3200.E3     # maximum value of X2
+x2max       = 1600.E3     # maximum value of X2
 ix2_bc      = periodic    # Inner-X2 boundary condition flag
 ox2_bc      = periodic    # Outer-X2 boundary condition flag
 
-nx3         = 1           # Number of zones in X3-direction
-x3min       = -0.5        # minimum value of X3
-x3max       = 0.5         # maximum value of X3
+nx3         = 80         # Number of zones in X3-direction
+x3min       = 0.        	# minimum value of X3
+x3max       = 1600.E3    	# maximum value of X3
 ix3_bc      = periodic    # Inner-X3 boundary condition flag
 ox3_bc      = periodic    # Outer-X3 boundary condition flag
 
 <meshblock>
-nx1         = 80
-nx2         = 40
-nx3         = 1
+nx1         = 40
+nx2         = 20
+nx3         = 20
 
 <hydro>
 gamma       = 1.42
@@ -81,7 +81,7 @@ Rd          = 3777.   # mu = 2.3175 g/mol
 P0          = 1.E5
 T0          = 169.
 Tmin        = 110.
-hrate       = -100.
+hrate       = -20.
 prad        = 2.e5
 init_Ttol   = 2.
-diagnostics = div, curl, b, mean, div_h, tempa, presa, radflux, hydroflux
+diagnostics = div, curl, b, mean, div_h, anomaly, radflux, hydroflux

src/diagnostics/anomaly.cpp

@@ -0,0 +1,91 @@
+// athena
+#include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
+
+// snap
+#include <snap/thermodynamics/thermodynamics.hpp>
+
+// exchanger
+#include <exchanger/exchanger.hpp>
+
+// diagnostics
+#include "diagnostics.hpp"
+
+Anomaly::Anomaly(MeshBlock *pmb) : Diagnostics(pmb, "rhoa,tempa,v1a,presa") {
+  type = "VECTORS";
+  data.NewAthenaArray(4, ncells3_, ncells2_, ncells1_);
+  mean_.resize(4 * ncells1_);
+
+  pexh = std::make_shared<PlanarExchanger<Real, 2>>("diag/anomaly");
+
+  // volumn
+  pexh->send_buffer[0].resize(total_vol_.size());
+  pexh->recv_buffer[0].resize(total_vol_.size());
+
+  // value
+  pexh->send_buffer[1].resize(4 * ncells1_);
+  pexh->recv_buffer[1].resize(4 * ncells1_);
+}
+
+void Anomaly::Finalize(MeshBlock *pmb) {
+  Coordinates *pcoord = pmb->pcoord;
+  auto const &w = pmb->phydro->w;
+
+  auto pthermo = Thermodynamics::GetInstance();
+
+  int is = pmb->is, js = pmb->js, ks = pmb->ks;
+  int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
+
+  std::fill(total_vol_.begin(), total_vol_.end(), 0.);
+  std::fill(mean_.begin(), mean_.end(), 0.);
+
+  // calculate horizontal mean
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j) {
+      pcoord->CellVolume(k, j, is, ie, vol_);
+
+      for (int i = is; i <= ie; ++i) {
+        total_vol_[i] += vol_(i);
+        // density anomaly
+        mean_[i] += vol_(i) * w(IDN, k, j, i);
+        // temperature anomaly
+        mean_[ncells1_ + i] += vol_(i) * pthermo->GetTemp(pmb, k, j, i);
+        // vertical velocity anomaly
+        mean_[2 * ncells1_ + i] += vol_(i) * w(IVX, k, j, i);
+        // pressure anomaly
+        mean_[3 * ncells1_ + i] += vol_(i) * w(IPR, k, j, i);
+      }
+    }
+
+  pexh->Regroup(pmb, X1DIR);
+  packData(pmb);
+  pexh->ReduceAll(pmb, MPI_SUM);
+  unpackData(pmb);
+
+  // pressure anomaly
+  for (int k = ks; k <= ke; ++k)
+    for (int j = js; j <= je; ++j)
+      for (int i = is; i <= ie; ++i) {
+        // density anomaly
+        data(0, k, j, i) = w(IDN, k, j, i) - mean_[i] / total_vol_[i];
+        // temperature anomaly
+        data(1, k, j, i) = pthermo->GetTemp(pmb, k, j, i) -
+                           mean_[ncells1_ + i] / total_vol_[i];
+        // vertical velocity anomaly
+        data(2, k, j, i) =
+            w(IVX, k, j, i) - mean_[2 * ncells1_ + i] / total_vol_[i];
+        // pressure anomaly
+        data(3, k, j, i) =
+            w(IPR, k, j, i) - mean_[3 * ncells1_ + i] / total_vol_[i];
+      }
+}
+
+void Anomaly::packData(MeshBlock const *pmb) {
+  pexh->send_buffer[0].swap(total_vol_);
+  pexh->send_buffer[1].swap(mean_);
+}
+
+void Anomaly::unpackData(MeshBlock const *pmb) {
+  total_vol_.swap(pexh->recv_buffer[0]);
+  mean_.swap(pexh->recv_buffer[1]);
+}

src/diagnostics/diagnostics.hpp

@@ -12,12 +12,12 @@
 #include <configure.hpp>
 #include <virtual_groups.hpp>
 
-// exchanger
-#include <exchanger/exchanger.hpp>
-
 class ParameterInput;
 class Meshlock;
 
+template <typename T, int N>
+class PlanarExchanger;
+
 class Diagnostics : public NamedGroup {
  public:
   // data
@@ -128,37 +128,17 @@ class HorizontalDivergence : public Diagnostics {
   AthenaArray<Real> v2f2_, v3f3_;
 };
 
-// 6. temperature anomaly
-class TemperatureAnomaly : public Diagnostics {
+// 6. anomaly
+class Anomaly : public Diagnostics {
  public:
   std::shared_ptr<PlanarExchanger<Real, 2>> pexh;
 
  public:
-  TemperatureAnomaly(MeshBlock *pmb);
-  virtual ~TemperatureAnomaly() {}
+  Anomaly(MeshBlock *pmb);
+  virtual ~Anomaly() {}
 
   void Finalize(MeshBlock *pmb) override;
-  int GetNumVars() const override { return 1; }
-
- protected:
-  void packData(MeshBlock const *pmb);
-  void unpackData(MeshBlock const *pmb);
-
-  //! mean component
-  std::vector<Real> mean_;
-};
-
-// 7. pressure anomaly
-class PressureAnomaly : public Diagnostics {
- public:
-  std::shared_ptr<PlanarExchanger<Real, 2>> pexh;
-
- public:
-  PressureAnomaly(MeshBlock *pmb);
-  virtual ~PressureAnomaly() {}
-
-  void Finalize(MeshBlock *pmb) override;
-  int GetNumVars() const override { return 1; }
+  int GetNumVars() const override { return 4; }
 
  protected:
   void packData(MeshBlock const *pmb);
@@ -208,6 +188,23 @@ class HydroFlux : public Diagnostics {
   int ncycle_;
 };
 
+// 10. w_avg
+class V1Moments : public Diagnostics {
+ public:
+  std::shared_ptr<PlanarExchanger<Real, 2>> pexh;
+
+ public:
+  V1Moments(MeshBlock *pmb);
+  virtual ~V1Moments() {}
+
+  void Finalize(MeshBlock *pmb) override;
+  int GetNumVars() const override { return 3; }
+
+ protected:
+  void packData(MeshBlock const *pmb);
+  void unpackData(MeshBlock const *pmb);
+};
+
 /* 6. eddy flux
 class EddyFlux : public Diagnostics {
  public:

src/diagnostics/diagnostics_factory.cpp

@@ -28,14 +28,14 @@ DiagnosticsContainer DiagnosticsFactory::CreateFrom(MeshBlock *pmb,
       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 == "tempa") {  // 6.
-      diag.push_back(std::make_shared<TemperatureAnomaly>(pmb));
-    } else if (name == "presa") {  // 7.
-      diag.push_back(std::make_shared<PressureAnomaly>(pmb));
+    } else if (name == "anomaly") {  // 6.
+      diag.push_back(std::make_shared<Anomaly>(pmb));
     } else if (name == "radflux") {  // 8.
       diag.push_back(std::make_shared<RadiativeFlux>(pmb));
     } else if (name == "hydroflux") {  // 9.
       diag.push_back(std::make_shared<HydroFlux>(pmb));
+    } else if (name == "w_avg") {  // 10.
+      diag.push_back(std::make_shared<V1Moments>(pmb));
       /*} else if (name == "eddyflux") { // 6.
         diag.push_back(std::make_shared<EddyFlux>(pmb));
       } else if (name == "am") { // 11.

src/diagnostics/hydro_flux.cpp

@@ -2,6 +2,9 @@
 #include <athena/coordinates/coordinates.hpp>
 #include <athena/hydro/hydro.hpp>
 
+// exchanger
+#include <exchanger/exchanger.hpp>
+
 // snap
 #include <snap/thermodynamics/thermodynamics.hpp>
 

src/diagnostics/radiative_flux.cpp

@@ -1,6 +1,9 @@
 // athena
 #include <athena/coordinates/coordinates.hpp>
 
+// exchanger
+#include <exchanger/exchanger.hpp>
+
 // harp
 #include <harp/radiation.hpp>