From 640a7b472e0aa11bd1ca9899a6d2ce2131f85d19 Mon Sep 17 00:00:00 2001
From: Kengo Tomida <tomida@astr.tohoku.ac.jp>
Date: Sat, 2 Mar 2024 21:57:10 +0900
Subject: [PATCH] omega for self-gravity can be set in the input file

---
 src/crdiffusion/crdiffusion.cpp    | 1 -
 src/crdiffusion/mg_crdiffusion.cpp | 7 -------
 src/gravity/mg_gravity.cpp         | 3 ++-
 src/gravity/mg_gravity.hpp         | 6 ++----
 4 files changed, 4 insertions(+), 13 deletions(-)

diff --git a/src/crdiffusion/crdiffusion.cpp b/src/crdiffusion/crdiffusion.cpp
index 1de96672c..b06590758 100644
--- a/src/crdiffusion/crdiffusion.cpp
+++ b/src/crdiffusion/crdiffusion.cpp
@@ -45,7 +45,6 @@ CRDiffusion::CRDiffusion(MeshBlock *pmb, ParameterInput *pin) :
   Dpara_ = pin->GetReal("crdiffusion", "Dpara");
   Dperp_ = pin->GetReal("crdiffusion", "Dperp");
   Lambda_ = pin->GetReal("crdiffusion", "Lambda");
-  zeta_factor_ = pin->GetReal("crdiffusion", "zeta_factor");
 
   output_defect = pin->GetOrAddBoolean("crdiffusion", "output_defect", false);
   if (output_defect)
diff --git a/src/crdiffusion/mg_crdiffusion.cpp b/src/crdiffusion/mg_crdiffusion.cpp
index b0b488b46..5565a9f3d 100644
--- a/src/crdiffusion/mg_crdiffusion.cpp
+++ b/src/crdiffusion/mg_crdiffusion.cpp
@@ -265,13 +265,6 @@ void MGCRDiffusionDriver::Solve(int stage, Real dt) {
 
   crtlist_->DoTaskListOneStage(pmy_mesh_, stage);
 
-  for (auto itr = vmg_.begin(); itr < vmg_.end(); itr++) {
-    Multigrid *pmg = *itr;
-    CRDiffusion *pcrdiff = pmg->pmy_block_->pcrdiff;
-    Hydro *phydro = pmg->pmy_block_->phydro;
-    pcrdiff->CalculateIonizationRate(phydro->w);
-  }
-
   return;
 }
 
diff --git a/src/gravity/mg_gravity.cpp b/src/gravity/mg_gravity.cpp
index cc3a1c4be..062e227fd 100644
--- a/src/gravity/mg_gravity.cpp
+++ b/src/gravity/mg_gravity.cpp
@@ -43,6 +43,7 @@ MGGravityDriver::MGGravityDriver(Mesh *pm, ParameterInput *pin)
                       pm->MGGravitySourceMaskFunction_, pm->MGGravitySourceMaskFunction_,
                       1, 0, 0) {
   four_pi_G_ = pmy_mesh_->four_pi_G_;
+  omega_ = pin->GetOrAddReal("gravity", "omega", 1.15);
   eps_ = pin->GetOrAddReal("gravity", "threshold", -1.0);
   niter_ = pin->GetOrAddInteger("gravity", "niteration", -1);
   ffas_ = pin->GetOrAddBoolean("gravity", "fas", ffas_);
@@ -226,7 +227,7 @@ void MGGravity::Smooth(AthenaArray<Real> &u, const AthenaArray<Real> &src,
   if (rlev <= 0) dx = rdx_*static_cast<Real>(1<<(-rlev));
   else           dx = rdx_/static_cast<Real>(1<<rlev);
   Real dx2 = SQR(dx);
-  Real isix = omega_/6.0;
+  Real isix = static_cast<MGGravityDriver*>(pmy_driver_)->omega_/6.0;
   color ^= pmy_driver_->coffset_;
 
 #pragma omp parallel for num_threads(pmy_driver_->nthreads_) if (th && (ku-kl) >= minth_)
diff --git a/src/gravity/mg_gravity.hpp b/src/gravity/mg_gravity.hpp
index 67dbf1074..e04e3d384 100644
--- a/src/gravity/mg_gravity.hpp
+++ b/src/gravity/mg_gravity.hpp
@@ -41,9 +41,6 @@ class MGGravity : public Multigrid {
   void CalculateFASRHS(AthenaArray<Real> &def, const AthenaArray<Real> &src,
                 const AthenaArray<Real> &coeff, const AthenaArray<Real> &matrix,
                 int rlev, int il, int iu, int jl, int ju, int kl, int ku, bool th) final;
-
- private:
-  static constexpr Real omega_ = 1.15;
 };
 
 
@@ -57,8 +54,9 @@ class MGGravityDriver : public MultigridDriver {
   void Solve(int stage, Real dt = 0.0) final;
   void ProlongateOctetBoundariesFluxCons(AthenaArray<Real> &dst,
                  AthenaArray<Real> &cbuf, const AthenaArray<bool> &ncoarse) final;
+  friend class MGGravity;
  private:
-  Real four_pi_G_;
+  Real four_pi_G_, omega_;
   GravityBoundaryTaskList *gtlist_;
 };