Refactoring Exchanger (#128)

- simplify exchanger
- add back more complex diagnostics

patches/24.time_integrator.patch

@@ -123,7 +123,7 @@ index 2bfc5178..e0fca29d 100644
 +  auto all_diags = pmb->pimpl->all_diags;
 +
 +  for (auto& diag : all_diags) {
-+    diag->Progress(pmb->phydro->w);
++    diag->Progress(pmb);
 +  }
 +
    pmb->UserWorkInLoop();

src/diagnostics/buoyancy.cpp

@@ -14,8 +14,10 @@ Buoyancy::Buoyancy(MeshBlock* pmb) : Diagnostics(pmb, "b") {
   grav_ = pmb->phydro->hsrc.GetG1();
 }
 
-void Buoyancy::Finalize(AthenaArray<Real> const& w) {
-  MeshBlock* pmb = pmy_block_;
+void Buoyancy::Finalize(MeshBlock* pmb) {
+  Coordinates* pcoord = pmb->pcoord;
+  auto const& w = pmb->phydro->w;
+
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
   int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
 

src/diagnostics/curl.cpp

@@ -1,5 +1,6 @@
 // athena
 #include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
 #include <athena/reconstruct/interpolation.hpp>
 
 // canoe
@@ -23,9 +24,10 @@ Curl::Curl(MeshBlock *pmb) : Diagnostics(pmb, "curl") {
   }
 }
 
-void Curl::Finalize(AthenaArray<Real> const &w) {
-  MeshBlock *pmb = pmy_block_;
+void Curl::Finalize(MeshBlock *pmb) {
   Coordinates *pcoord = pmb->pcoord;
+  auto const &w = pmb->phydro->w;
+
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
   int ie = pmb->ie, je = pmb->je, ke = pmb->ke;
 

src/diagnostics/diagnostics.cpp

@@ -12,16 +12,21 @@
 // diagnostics
 #include "diagnostics.hpp"
 
-Diagnostics::Diagnostics(MeshBlock *pmb, std::string name)
-    : NamedGroup(name), pmy_block_(pmb), ncycle_(0) {
+Diagnostics::Diagnostics(MeshBlock *pmb, std::string name) : NamedGroup(name) {
   Application::Logger app("main");
   app->Log("Initialize Diagnostics");
 
   ncells1_ = pmb->block_size.nx1 + 2 * (NGHOST);
   ncells2_ = 1;
   ncells3_ = 1;
-  if (pmb->pmy_mesh->f2) ncells2_ = pmb->block_size.nx2 + 2 * (NGHOST);
-  if (pmb->pmy_mesh->f3) ncells3_ = pmb->block_size.nx3 + 2 * (NGHOST);
+
+  if (pmb->pmy_mesh->f2) {
+    ncells2_ = pmb->block_size.nx2 + 2 * (NGHOST);
+  }
+
+  if (pmb->pmy_mesh->f3) {
+    ncells3_ = pmb->block_size.nx3 + 2 * (NGHOST);
+  }
 
   x1edge_.NewAthenaArray(ncells1_ + 1);
   x1edge_p1_.NewAthenaArray(ncells1_);
@@ -38,6 +43,7 @@ Diagnostics::Diagnostics(MeshBlock *pmb, std::string name)
 
   vol_.NewAthenaArray(ncells1_);
   total_vol_.NewAthenaArray(ncells1_);
+  total_area_.NewAthenaArray(ncells1_);
 }
 
 Diagnostics::~Diagnostics() {

src/diagnostics/diagnostics.hpp

@@ -16,6 +16,7 @@
 #include <exchanger/exchanger.hpp>
 
 class ParameterInput;
+class Meshlock;
 
 class Diagnostics : public NamedGroup {
  public:
@@ -28,19 +29,12 @@ class Diagnostics : public NamedGroup {
   virtual ~Diagnostics();
 
   virtual int GetNumVars() const = 0;
-  virtual void Progress(AthenaArray<Real> const &w) {}
-  virtual void Finalize(AthenaArray<Real> const &w) {}
+  virtual void Progress(MeshBlock *pmb) {}
+  virtual void Finalize(MeshBlock *pmb) {}
 
  protected:
-  MeshBlock *pmy_block_;
-
   int ncells1_, ncells2_, ncells3_;
 
-  int ncycle_;
-
-  //! mean and eddy component
-  AthenaArray<Real> mean_, eddy_;
-
   //! MPI color of each rank
   std::vector<int> color_;
 
@@ -53,10 +47,7 @@ class Diagnostics : public NamedGroup {
 
   AthenaArray<Real> x1area_, x2area_, x2area_p1_, x3area_, x3area_p1_;
 
-  AthenaArray<Real> vol_, total_vol_;
-
-  void GatherVolumnData(AthenaArray<Real> &total_vol,
-                        AthenaArray<Real> &total_data);
+  AthenaArray<Real> vol_, total_vol_, total_area_;
 };
 
 using DiagnosticsPtr = std::shared_ptr<Diagnostics>;
@@ -74,7 +65,7 @@ class Divergence : public Diagnostics {
   Divergence(MeshBlock *pmb);
   virtual ~Divergence() {}
 
-  void Finalize(AthenaArray<Real> const &w) override;
+  void Finalize(MeshBlock *pmb) override;
   int GetNumVars() const override { return 1; }
 
  protected:
@@ -87,7 +78,7 @@ class Curl : public Diagnostics {
   Curl(MeshBlock *pmb);
   virtual ~Curl() {}
 
-  void Finalize(AthenaArray<Real> const &w) override;
+  void Finalize(MeshBlock *pmb) override;
   int GetNumVars() const override { return 1; }
 
  protected:
@@ -100,7 +91,7 @@ class Buoyancy : public Diagnostics {
   Buoyancy(MeshBlock *pmb);
   virtual ~Buoyancy() {}
 
-  void Finalize(AthenaArray<Real> const &w) override;
+  void Finalize(MeshBlock *pmb) override;
   int GetNumVars() const override { return 1; }
 
  protected:
@@ -112,71 +103,99 @@ class Buoyancy : public Diagnostics {
 class HydroMean : public Diagnostics {
  public:
   HydroMean(MeshBlock *pmb);
-  virtual ~HydroMean();
+  virtual ~HydroMean() {}
 
-  void Progress(AthenaArray<Real> const &w) override;
-  void Finalize(AthenaArray<Real> const &w) override;
+  void Progress(MeshBlock *pmb) override;
+  void Finalize(MeshBlock *pmb) override;
   int GetNumVars() const override { return NHYDRO; }
+
+ protected:
+  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_;
 };
 
-// 5. temperature anomaly
-class TemperatureAnomaly : public Diagnostics {
+// 6. temperature anomaly
+class TemperatureAnomaly : public Diagnostics, public PlanarExchanger<Real, 0> {
  public:
   TemperatureAnomaly(MeshBlock *pmb);
   virtual ~TemperatureAnomaly() {}
-  void Finalize(AthenaArray<Real> const &w);
+
+  void Finalize(MeshBlock *pmb) override;
+  int GetNumVars() const override { return 1; }
+
+ protected:
+  //! mean component
+  AthenaArray<Real> mean_;
 };
 
-// 6. pressure anomaly
-class PressureAnomaly : public Diagnostics {
+// 7. pressure anomaly
+class PressureAnomaly : public Diagnostics, public PlanarExchanger<Real, 0> {
  public:
   PressureAnomaly(MeshBlock *pmb);
   virtual ~PressureAnomaly() {}
-  void Finalize(AthenaArray<Real> const &w);
+
+  void Finalize(MeshBlock *pmb) override;
+  int GetNumVars() const override { return 1; }
 
  protected:
-  AthenaArray<Real> pf_;
+  //! mean component
+  AthenaArray<Real> mean_;
 };
 
-/* 6. eddy flux
-class EddyFlux : public Diagnostics {
+// 8. total radiative flux
+class RadiativeFlux : public Diagnostics, public PlanarExchanger<Real, 0> {
  public:
-  EddyFlux(MeshBlock *pmb);
-  virtual ~EddyFlux();
-  void Progress(AthenaArray<Real> const &w);
-  void Finalize(AthenaArray<Real> const &w);
+  RadiativeFlux(MeshBlock *pmb);
+  virtual ~RadiativeFlux() {}
+
+  void Progress(MeshBlock *pmb) override;
+  void Finalize(MeshBlock *pmb) override;
+  int GetNumVars() const override { return 2; }
+
+ protected:
+  int ncycle_;
 };
 
-// 7. hydro flux
-class HydroFlux : public Diagnostics {
+// 9. hydro flux
+class HydroFlux : public Diagnostics, public PlanarExchanger<Real, 0> {
  public:
   HydroFlux(MeshBlock *pmb);
   virtual ~HydroFlux() {}
-  void Progress(AthenaArray<Real> const &w);
-  void Finalize(AthenaArray<Real> const &w);
-};
 
-// 8. horizontal divergence
-class HorizontalDivergence : public Diagnostics {
- public:
-  HorizontalDivergence(MeshBlock *pmb);
-  virtual ~HorizontalDivergence();
-  void Finalize(AthenaArray<Real> const &w);
+  void Progress(MeshBlock *pmb) override;
+  void Finalize(MeshBlock *pmb) override;
+  int GetNumVars() const override { return NHYDRO; }
 
  protected:
-  AthenaArray<Real> v2f2_, v3f3_;
+  int ncycle_;
 };
 
-
-// 10. total radiative flux
-class RadiativeFlux : public Diagnostics {
+/* 6. eddy flux
+class EddyFlux : public Diagnostics {
  public:
-  RadiativeFlux(MeshBlock *pmb);
-  virtual ~RadiativeFlux() {}
+  EddyFlux(MeshBlock *pmb);
+  virtual ~EddyFlux();
   void Progress(AthenaArray<Real> const &w);
   void Finalize(AthenaArray<Real> const &w);
+
+ protected:
+  //! mean and eddy component
+  AthenaArray<Real> mean_, eddy_;
 };
 
+
 // 11. total angular momentum
 class SphericalAngularMomentum : public Diagnostics {
  public:

src/diagnostics/diagnostics_factory.cpp

@@ -26,18 +26,18 @@ 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 == "tempa") { // 4.
-        diag.push_back(std::make_shared<TemperatureAnomaly>(pmb));
-      } else if (name == "presa") { // 5.
-        diag.push_back(std::make_shared<PressureAnomaly>(pmb));
-      } else if (name == "eddyflux") { // 6.
+    } 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 == "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 == "eddyflux") { // 6.
         diag.push_back(std::make_shared<EddyFlux>(pmb));
-      } else if (name == "hydroflux") { // 7.
-        diag.push_back(std::make_shared<HydroFlux>(pmb));
-      } else if (name == "div_h") { // 8.
-        diag.push_back(std::make_shared<HorizontalDivergence>(pmb));
-      } else if (name == "radflux") { // 10.
-        diag.push_back(std::make_shared<RadiativeFlux>(pmb));
       } else if (name == "am") { // 11.
         diag.push_back(std::make_shared<AngularMomentum>(pmb));
       } else if (name == "eke") { // 12.

src/diagnostics/divergence.cpp

@@ -1,5 +1,6 @@
 // athena
 #include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
 #include <athena/reconstruct/interpolation.hpp>
 
 // canoe
@@ -20,9 +21,9 @@ Divergence::Divergence(MeshBlock *pmb) : Diagnostics(pmb, "div") {
   }
 }
 
-void Divergence::Finalize(AthenaArray<Real> const &w) {
-  MeshBlock *pmb = pmy_block_;
+void Divergence::Finalize(MeshBlock *pmb) {
   Coordinates *pcoord = pmb->pcoord;
+  auto const &w = pmb->phydro->w;
 
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
   int ie = pmb->ie, je = pmb->je, ke = pmb->ke;

src/diagnostics/todo/horizontal_divergence.cpp → src/diagnostics/horizontal_divergence.cpp

@@ -1,5 +1,6 @@
 // athena
 #include <athena/coordinates/coordinates.hpp>
+#include <athena/hydro/hydro.hpp>
 #include <athena/reconstruct/interpolation.hpp>
 
 // canoe
@@ -9,21 +10,17 @@ HorizontalDivergence::HorizontalDivergence(MeshBlock *pmb)
     : Diagnostics(pmb, "div_h") {
   type = "SCALARS";
   data.NewAthenaArray(ncells3_, ncells2_, ncells1_);
+
   if (pmb->block_size.nx2 > 1)
     v2f2_.NewAthenaArray(ncells3_, ncells2_ + 1, ncells1_);
   if (pmb->block_size.nx3 > 1)
     v3f3_.NewAthenaArray(ncells3_ + 1, ncells2_, ncells1_);
 }
 
-HorizontalDivergence::~HorizontalDivergence() {
-  data.DeleteAthenaArray();
-  if (pmy_block_->block_size.nx2 > 1) v2f2_.DeleteAthenaArray();
-  if (pmy_block_->block_size.nx3 > 1) v3f3_.DeleteAthenaArray();
-}
-
-void HorizontalDivergence::Finalize(AthenaArray<Real> const &w) {
-  MeshBlock *pmb = pmy_block_;
+void HorizontalDivergence::Finalize(MeshBlock *pmb) {
   Coordinates *pcoord = pmb->pcoord;
+  auto const &w = pmb->phydro->w;
+
   int is = pmb->is, js = pmb->js, ks = pmb->ks;
   int ie = pmb->ie, je = pmb->je, ke = pmb->ke;