From caee775c79e09fd8f6ef7e1b696e8fe9730b1748 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 12 Aug 2022 16:17:35 +0200
Subject: [PATCH] Advance entropic time step implementation

- Remains lots of work in 2d and 3d (great progresses in 2d)
- Remains to improve some constants evaluations (for Mie-Grunisen)
- The case Gamma_infty<1 is not started
---
 src/analysis/Polynomial1D.hpp             |  24 +-
 src/analysis/Polynomial1DRootComputer.cpp |   5 +-
 src/language/modules/SchemeModule.cpp     |  43 +-
 src/scheme/AcousticSolver.cpp             | 497 ++++++++++++++++++++--
 src/scheme/AcousticSolver.hpp             |  52 ++-
 5 files changed, 544 insertions(+), 77 deletions(-)

diff --git a/src/analysis/Polynomial1D.hpp b/src/analysis/Polynomial1D.hpp
index db033a4e8..1e0ba6ead 100644
--- a/src/analysis/Polynomial1D.hpp
+++ b/src/analysis/Polynomial1D.hpp
@@ -24,18 +24,18 @@ class [[nodiscard]] Polynomial1D
   PUGS_INLINE
   friend Polynomial1D _simplify(Polynomial1D && p)
   {
-    return std::move(p);
-    // size_t real_degree = p._getRealDegree();
-
-    // if (real_degree != p.degree()) {
-    //   Polynomial1D q(real_degree);
-    //   for (size_t i = 0; i <= real_degree; ++i) {
-    //     q.coefficient(i) = p.coefficient(i);
-    //   }
-    //   return q;
-    // } else {
-    //   return std::move(p);
-    // }
+    // return std::move(p);
+    size_t real_degree = p._getRealDegree();
+
+    if (real_degree != p.degree()) {
+      Polynomial1D q(real_degree);
+      for (size_t i = 0; i <= real_degree; ++i) {
+        q.coefficient(i) = p.coefficient(i);
+      }
+      return q;
+    } else {
+      return std::move(p);
+    }
   }
 
  public:
diff --git a/src/analysis/Polynomial1DRootComputer.cpp b/src/analysis/Polynomial1DRootComputer.cpp
index a2c54f88c..59b7fdc6c 100644
--- a/src/analysis/Polynomial1DRootComputer.cpp
+++ b/src/analysis/Polynomial1DRootComputer.cpp
@@ -7,7 +7,7 @@ Polynomial1DRootComputer::_computeRootBetween(const PositionSignChanges& left_bo
   double left  = left_bound.position;
   double right = right_bound.position;
 
-  const double interval_size = right - left;
+  const double interval_size = m_right - m_left;
 
   const Polynomial1D& p0 = m_sturm_sequence.p(0);
   double left_value      = p0(left);
@@ -49,6 +49,9 @@ Polynomial1DRootComputer::getFirstRoot() const
         } else if (l.nb_sign_changes - rk.nb_sign_changes == 0) {
           a = rk;
         }
+        if (b.position - a.position < 1E-14 * (m_right - m_left)) {
+          return a.position;
+        }
       } while (a.nb_sign_changes - b.nb_sign_changes != 1);
       return _computeRootBetween(a, b);
     } else {
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 186304326..6e234910d 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -297,16 +297,16 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("acoustic_epsilon_dt",
+  this->_addBuiltinFunction("acoustic_pg_epsilon_dt",
                             std::function(
 
                               [](const std::shared_ptr<const IDiscreteFunction>& rho,
                                  const std::shared_ptr<const IDiscreteFunction>& u,
-                                 const std::shared_ptr<const IDiscreteFunction>& E,
+                                 const std::shared_ptr<const IDiscreteFunction>& epsilon,
                                  const std::shared_ptr<const IDiscreteFunction>& p,
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr,
                                  const std::shared_ptr<const ItemValueVariant>& ur, const double& dt_max) -> double {
-                                return acoustic_epsilon_dt(rho, u, E, p, Ajr, ur, dt_max);
+                                return acoustic_pg_epsilon_dt(rho, u, epsilon, p, Ajr, ur, dt_max);
                               }
 
                               ));
@@ -326,6 +326,43 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("acoustic_mg_epsilon_dt",
+                            std::function(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                 const std::shared_ptr<const IDiscreteFunction>& u,
+                                 const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                 const std::shared_ptr<const IDiscreteFunction>& p_k,
+                                 const std::shared_ptr<const IDiscreteFunction>& rho0k0_expN0p1,
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                                 const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                                 const double& dt_max) -> double {
+                                return acoustic_mg_epsilon_dt(rho, u, epsilon, p, p_k, rho0k0_expN0p1, Ajr_variant,
+                                                              ur_variant, dt_max);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("acoustic_mg_entropy_dt",
+                            std::function(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                 const std::shared_ptr<const IDiscreteFunction>& u,
+                                 const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                 const std::shared_ptr<const IDiscreteFunction>& p_k,
+                                 const std::shared_ptr<const IDiscreteFunction>& Gamma_tau,
+                                 const std::shared_ptr<const IDiscreteFunction>& Gamma_inf,
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                                 const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                                 const double& dt_max) -> double {
+                                return acoustic_mg_entropy_dt(rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf,
+                                                              Ajr_variant, ur_variant, dt_max);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("glace_fluxes", std::function(
 
                                               [](const std::shared_ptr<const IDiscreteFunction>& rho,
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index f63eef8ae..afeabedb9 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -113,7 +113,7 @@ acoustic_rho_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       for (size_t R = 0; R < cell_nodes.size(); ++R) {
         const NodeId r   = cell_nodes[R];
         const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
-        G -= dot(ur[rp1], Rd{ur[r][1], -ur[r][0]});
+        G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
       }
       Gj[j] = 0.5 * G;
     }
@@ -143,13 +143,13 @@ acoustic_rho_dt(const DiscreteFunctionP0<Dimension, double>& rho,
 
 template <size_t Dimension>
 double
-acoustic_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
-                    const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
-                    const DiscreteFunctionP0<Dimension, double>& E,
-                    const DiscreteFunctionP0<Dimension, double>& p,
-                    const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
-                    const std::shared_ptr<const ItemValueVariant>& ur_variant,
-                    const double& dt_max)
+acoustic_pg_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
+                       const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
+                       const DiscreteFunctionP0<Dimension, double>& epsilon,
+                       const DiscreteFunctionP0<Dimension, double>& p,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                       const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                       const double& dt_max)
 {
   using Rd           = TinyVector<Dimension>;
   using Rdxd         = TinyMatrix<Dimension>;
@@ -223,8 +223,7 @@ acoustic_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
   for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
     const double inv_mj = 1 / (rho[j] * Vj[j]);
 
-    Polynomial1D p_e(
-      {E[j] - 0.5 * dot(u[j], u[j]), inv_mj * (dSj[j] - p[j] * dVj[j]), -0.5 * inv_mj * inv_mj * dot(dPj[j], dPj[j])});
+    Polynomial1D p_e({epsilon[j], inv_mj * (dSj[j] - p[j] * dVj[j]), -0.5 * inv_mj * inv_mj * dot(dPj[j], dPj[j])});
     Polynomial1DRootComputer computer{p_e, 0, dt};
     std::optional dt_e = computer.getFirstRoot();
     if (dt_e.has_value()) {
@@ -296,7 +295,7 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       for (size_t R = 0; R < cell_nodes.size(); ++R) {
         const NodeId r   = cell_nodes[R];
         const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
-        G -= dot(ur[rp1], Rd{ur[r][1], -ur[r][0]});
+        G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
       }
       Gj[j] = 0.5 * G;
     }
@@ -341,8 +340,6 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
   DiscreteFunctionP0<Dimension, const double> tauj_gamma_1 = pow(rho, 1 - gamma);
   DiscreteFunctionP0<Dimension, const double> tauj_gamma_2 = tauj_gamma_1 * rho;
 
-  std::vector<Polynomial1D> q_Sj(mesh.numberOfCells(), Polynomial1D(std::vector<double>{0}));
-
   Polynomial1D DT({0, 1});
   for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
     if (dSj[j] == 0)
@@ -357,7 +354,8 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       } else if constexpr (Dimension == 2) {
         return mj * p[j] * Gj[j];
       } else {
-        throw NotImplementedError("3d");
+#warning 3d
+        // throw NotImplementedError("3d");
       }
     }();
 
@@ -370,20 +368,25 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       q_tau = q_tau + Polynomial1D({0, 0, inv_mj * Gj[j]});
     }
 
+    Polynomial1D dV(std::vector<double>{dVj[j]});
+    if constexpr (Dimension == 2) {
+      dV = dV + Gj[j] * DT;
+    }
+
     Polynomial1D S1 = Polynomial1D(std::vector<double>{(mj_pi_Gj - 0.5 * dPj2) * tauj_gamma_1[j]}) +
-                      ((gamma[j] - 1) * (dVj[j] /* +DT*G[j] */) * tauj_gamma_2[j]) * dS1;
+                      ((gamma[j] - 1) * tauj_gamma_2[j]) * dV * dS1;
 
     if (gamma[j] < 2) {
+#warning Not valid for d>1
       const double tauj_min = std::min(q_tau(0), q_tau(dt_max));
-      S1 = S1 + Polynomial1D(std::vector<double>{0.5 * (gamma[j] - 1) * (gamma[j] - 2) * (dVj[j] /* +DT*G[j]  */) *
-                                                 (dVj[j] /* +DT*G[j] */) * std::pow(tauj_min, gamma[j] - 3)}) *
-                  (Polynomial1D(std::vector<double>{epsilon[j]}) + (inv_mj * DT) * dS1);
+
+      S1 = S1 +
+           Polynomial1D(std::vector<double>{0.5 * (gamma[j] - 1) * (gamma[j] - 2) * std::pow(tauj_min, gamma[j] - 3)}) *
+             dV * dV * (Polynomial1D(std::vector<double>{epsilon[j]}) + (inv_mj * DT) * dS1);
     }
 
     Polynomial1D p_S = Polynomial1D(std::vector<double>{delta_S}) + (inv_mj * DT) * S1;
 
-    q_Sj[j] = p_S;
-
     Polynomial1DRootComputer computer{p_S, 0, dt};
     std::optional dt_S = computer.getFirstRoot();
     if (dt_S.has_value()) {
@@ -395,6 +398,284 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
   return parallel::allReduceMin(dt);
 }
 
+template <size_t Dimension>
+double
+acoustic_mg_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
+                       const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
+                       const DiscreteFunctionP0<Dimension, double>& epsilon,
+                       const DiscreteFunctionP0<Dimension, double>& p,
+                       const DiscreteFunctionP0<Dimension, double>& p_k,
+                       const DiscreteFunctionP0<Dimension, double>& rho0k0_expN0p1,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                       const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                       const double& dt_max)
+{
+  using Rd           = TinyVector<Dimension>;
+  using Rdxd         = TinyMatrix<Dimension>;
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  double dt = dt_max;
+
+  const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
+
+  const NodeValue<const Rd> ur           = ur_variant->get<NodeValue<const Rd>>();
+  const NodeValuePerCell<const Rdxd> Ajr = Ajr_variant->get<NodeValuePerCell<const Rdxd>>();
+  MeshDataType& mesh_data                = MeshDataManager::instance().getMeshData(mesh);
+
+  const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+  const CellValue<const double> Vj     = mesh_data.Vj();
+
+  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+  CellValue<double> dVj{mesh.connectivity()};
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dV = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+        dV += dot(Cjr(j, R), ur[r]);
+      }
+
+      dVj[j] = dV;
+    });
+
+  CellValue<double> dSj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dS = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+
+        Rd du = u[j] - ur[r];
+        dS += dot(Ajr(j, R) * du, du);
+      }
+
+      dSj[j] = dS;
+    });
+
+  CellValue<Rd> dPj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      Rd dP = zero;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+
+        const Rd du = u[j] - ur[r];
+        dP += Ajr(j, R) * du;
+      }
+
+      dPj[j] = dP;
+    });
+
+  Polynomial1D DT({0, 1});
+
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    const double inv_mj = 1 / (rho[j] * Vj[j]);
+
+    Polynomial1D dV(std::vector<double>{dVj[j]});
+    if constexpr (Dimension == 2) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double G = 0;
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r   = cell_nodes[R];
+        const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+        G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
+      }
+
+      dV = dV + 0.5 * G * DT;
+      // throw NotImplementedError("Gj for d>1");
+    }
+
+    Polynomial1D p_e({epsilon[j], inv_mj * (dSj[j] - p[j] * dVj[j]), -0.5 * inv_mj * inv_mj * dot(dPj[j], dPj[j])});
+    p_e = p_e + inv_mj * DT * (p_k[j] * dV - (0.5 * inv_mj * rho0k0_expN0p1[j]) * dV * dV * DT);
+
+    Polynomial1DRootComputer computer{p_e, 0, dt};
+    std::optional dt_e = computer.getFirstRoot();
+    if (dt_e.has_value()) {
+      Assert(dt_e.value() <= dt);
+      dt = dt_e.value();
+    }
+  }
+
+  if (dt < dt_max) {
+    std::cout << "energy variation imposes time step\n";
+    dt *= 0.95;
+  }
+
+  return parallel::allReduceMin(dt);
+}
+
+template <size_t Dimension>
+double
+acoustic_mg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
+                       const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
+                       const DiscreteFunctionP0<Dimension, double>& epsilon,
+                       const DiscreteFunctionP0<Dimension, double>& p,
+                       const DiscreteFunctionP0<Dimension, double>& p_k,
+                       const DiscreteFunctionP0<Dimension, double>& Gamma_tau,
+                       const DiscreteFunctionP0<Dimension, double>& Gamma_inf,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                       const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                       const double& dt_max)
+{
+  double dt = dt_max;
+
+  using Rd           = TinyVector<Dimension>;
+  using Rdxd         = TinyMatrix<Dimension>;
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
+
+  const NodeValue<const Rd> ur           = ur_variant->get<NodeValue<const Rd>>();
+  const NodeValuePerCell<const Rdxd> Ajr = Ajr_variant->get<NodeValuePerCell<const Rdxd>>();
+  MeshDataType& mesh_data                = MeshDataManager::instance().getMeshData(mesh);
+
+  const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+  const CellValue<const double> Vj     = mesh_data.Vj();
+
+  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+  CellValue<double> dVj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dV = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+        dV += dot(Cjr(j, R), ur[r]);
+      }
+
+      dVj[j] = dV;
+    });
+
+  CellValue<double> Gj{mesh.connectivity()};
+  if constexpr (Dimension == 2) {
+    for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double G = 0;
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r   = cell_nodes[R];
+        const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+        G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
+      }
+      Gj[j] = 0.5 * G;
+    }
+  }
+
+  CellValue<double> dSj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dS = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+
+        const Rd du = u[j] - ur[r];
+        dS += dot(Ajr(j, R) * du, du);
+      }
+
+      dSj[j] = dS;
+    });
+
+  CellValue<Rd> dPj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      Rd dP = zero;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+
+        Rd du = u[j] - ur[r];
+        dP += Ajr(j, R) * du;
+      }
+
+      dPj[j] = dP;
+    });
+
+  Polynomial1D DT({0, 1});
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    if (dSj[j] == 0)
+      continue;
+    const double mj     = (rho[j] * Vj[j]);
+    const double inv_mj = 1 / mj;
+
+    const double mj_pi_Gj = [&]() -> double {
+      if constexpr (Dimension == 1) {
+        return 0;
+      } else if constexpr (Dimension == 2) {
+        return mj * p[j] * Gj[j];
+      } else {
+#warning 3D
+        //        throw NotImplementedError("3d");
+      }
+    }();
+
+    Polynomial1D dV(std::vector<double>{dVj[j]});
+    if constexpr (Dimension == 2) {
+      dV = dV + Gj[j] * DT;
+    }
+
+    const double dPj2 = dot(dPj[j], dPj[j]);
+
+    Polynomial1D dS1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dPj2});
+    dS1 = dS1 + p_k[j] * dV;
+
+    Polynomial1D q_tau({1. / rho[j], inv_mj * dVj[j]});
+    if constexpr (Dimension == 2) {
+      q_tau = q_tau + Polynomial1D({0, 0, inv_mj * Gj[j]});
+    }
+
+    if (Gamma_inf[j] >= 1) {
+      Polynomial1D S1 = Polynomial1D(std::vector<double>{(mj_pi_Gj - 0.5 * dPj2)});
+
+      // if (gamma[j] < 2) {
+      //   const double tauj_min = std::min(q_tau(0), q_tau(dt_max));
+
+      S1 = S1 + Gamma_tau[j] * rho[j] * dV * dS1;
+
+#warning APPROX
+      S1 = S1 - 0.5 * 5e+17 * inv_mj * dV * dV * DT;
+
+      Polynomial1D p_S = Polynomial1D(std::vector<double>{dSj[j]}) + (inv_mj * DT) * S1;
+
+      Polynomial1DRootComputer computer{p_S, 0, dt};
+      std::optional dt_S = computer.getFirstRoot();
+      if (dt_S.has_value()) {
+        Assert(dt_S.value() <= dt);
+        dt = dt_S.value();
+      }
+    } else {
+      throw NotImplementedError("Gamma_inf<1");
+    }
+  }
+
+  return parallel::allReduceMin(dt);
+}
+
 double
 acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                 const std::shared_ptr<const ItemValueVariant>& ur,
@@ -423,44 +704,44 @@ acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
 }
 
 double
-acoustic_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
-                    const std::shared_ptr<const IDiscreteFunction>& u,
-                    const std::shared_ptr<const IDiscreteFunction>& E,
-                    const std::shared_ptr<const IDiscreteFunction>& p,
-                    const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
-                    const std::shared_ptr<const ItemValueVariant>& ur_variant,
-                    const double& dt_max)
+acoustic_pg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                       const std::shared_ptr<const IDiscreteFunction>& u,
+                       const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                       const std::shared_ptr<const IDiscreteFunction>& p,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                       const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                       const double& dt_max)
 {
-  if (not checkDiscretizationType({rho, u, E, p}, DiscreteFunctionType::P0)) {
+  if (not checkDiscretizationType({rho, u, epsilon, p}, DiscreteFunctionType::P0)) {
     throw NormalError("acoustic solver expects P0 functions");
   }
 
-  std::shared_ptr mesh = getCommonMesh({rho, u, E, p});
+  std::shared_ptr mesh = getCommonMesh({rho, u, epsilon, p});
 
   switch (mesh->dimension()) {
   case 1: {
     constexpr size_t Dimension = 1;
-    return acoustic_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
-                               dt_max);
+    return acoustic_pg_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant,
+                                  ur_variant, dt_max);
   }
   case 2: {
     constexpr size_t Dimension = 2;
-    return acoustic_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
-                               dt_max);
+    return acoustic_pg_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant,
+                                  ur_variant, dt_max);
   }
   case 3: {
     constexpr size_t Dimension = 3;
-    return acoustic_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
-                               dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
-                               dt_max);
+    return acoustic_pg_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant,
+                                  ur_variant, dt_max);
   }
   default: {
     throw UnexpectedError("invalid mesh dimension");
@@ -518,6 +799,118 @@ acoustic_pg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
   }
 }
 
+double
+acoustic_mg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                       const std::shared_ptr<const IDiscreteFunction>& u,
+                       const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                       const std::shared_ptr<const IDiscreteFunction>& p,
+                       const std::shared_ptr<const IDiscreteFunction>& p_k,
+                       const std::shared_ptr<const IDiscreteFunction>& rho0k0_expN0p1,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                       const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                       const double& dt_max)
+{
+  if (not checkDiscretizationType({rho, u, epsilon, p, p_k, rho0k0_expN0p1}, DiscreteFunctionType::P0)) {
+    throw NormalError("acoustic solver expects P0 functions");
+  }
+
+  std::shared_ptr mesh = getCommonMesh({rho, u, epsilon, p, p_k, rho0k0_expN0p1});
+
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    return acoustic_mg_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho0k0_expN0p1),
+                                  Ajr_variant, ur_variant, dt_max);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    return acoustic_mg_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho0k0_expN0p1),
+                                  Ajr_variant, ur_variant, dt_max);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    return acoustic_mg_epsilon_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho0k0_expN0p1),
+                                  Ajr_variant, ur_variant, dt_max);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+double
+acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                       const std::shared_ptr<const IDiscreteFunction>& u,
+                       const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                       const std::shared_ptr<const IDiscreteFunction>& p,
+                       const std::shared_ptr<const IDiscreteFunction>& p_k,
+                       const std::shared_ptr<const IDiscreteFunction>& Gamma_tau,
+                       const std::shared_ptr<const IDiscreteFunction>& Gamma_inf,
+                       const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                       const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                       const double& dt_max)
+{
+  if (not checkDiscretizationType({rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf}, DiscreteFunctionType::P0)) {
+    throw NormalError("acoustic solver expects P0 functions");
+  }
+
+  std::shared_ptr mesh = getCommonMesh({rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf});
+
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    return acoustic_mg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_tau),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf), Ajr_variant,
+                                  ur_variant, dt_max);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    return acoustic_mg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_tau),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf), Ajr_variant,
+                                  ur_variant, dt_max);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    return acoustic_mg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_tau),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf), Ajr_variant,
+                                  ur_variant, dt_max);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
 template <size_t Dimension>
 class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
 {
@@ -942,8 +1335,22 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
           const NodeId r = cell_nodes[R];
           dV += dot(Cjr(j, R), ur[r]);
         }
+        if constexpr (Dimension == 1) {
+          new_rho[j] = 1. / (1. / rho[j] + dt / mj * dV);
+        } else if constexpr (Dimension == 2) {
+          double G = 0;
+          for (size_t R = 0; R < cell_nodes.size(); ++R) {
+            const NodeId r   = cell_nodes[R];
+            const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+            G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
+          }
+          G *= 0.5;
 
-        new_rho[j] = 1. / (1. / rho[j] + dt / mj * dV);
+          new_rho[j] = 1. / (1. / rho[j] + dt / mj * (dV + dt * G));
+        } else {
+#warning 3d
+          // throw NotImplementedError("3D with GCL error");
+        }
       });
 
     return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 1fc8e0355..e7bd26fa8 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -17,22 +17,21 @@ double acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                        const std::shared_ptr<const ItemValueVariant>& ur,
                        const double& dt_max);
 
-double acoustic_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
-                           const std::shared_ptr<const IDiscreteFunction>& u,
-                           const std::shared_ptr<const IDiscreteFunction>& E,
-                           const std::shared_ptr<const IDiscreteFunction>& p,
-                           const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
-                           const std::shared_ptr<const ItemValueVariant>& ur_variant,
-                           const double& dt_max);
-
-double acoustic_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
-                           const std::shared_ptr<const IDiscreteFunction>& u,
-                           const std::shared_ptr<const IDiscreteFunction>& epsilon,
-                           const std::shared_ptr<const IDiscreteFunction>& p,
-                           const std::shared_ptr<const IDiscreteFunction>& gamma,
-                           const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
-                           const std::shared_ptr<const ItemValueVariant>& ur_variant,
-                           const double& dt_max);
+double acoustic_pg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                              const std::shared_ptr<const IDiscreteFunction>& u,
+                              const std::shared_ptr<const IDiscreteFunction>& E,
+                              const std::shared_ptr<const IDiscreteFunction>& p,
+                              const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                              const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                              const double& dt_max);
+
+double acoustic_pg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                              const std::shared_ptr<const IDiscreteFunction>& u,
+                              const std::shared_ptr<const IDiscreteFunction>& E,
+                              const std::shared_ptr<const IDiscreteFunction>& p,
+                              const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                              const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                              const double& dt_max);
 
 double acoustic_pg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                               const std::shared_ptr<const IDiscreteFunction>& u,
@@ -43,6 +42,27 @@ double acoustic_pg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rh
                               const std::shared_ptr<const ItemValueVariant>& ur_variant,
                               const double& dt_max);
 
+double acoustic_mg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                              const std::shared_ptr<const IDiscreteFunction>& u,
+                              const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                              const std::shared_ptr<const IDiscreteFunction>& p,
+                              const std::shared_ptr<const IDiscreteFunction>& p_k,
+                              const std::shared_ptr<const IDiscreteFunction>& rho0k0_expN0p1,
+                              const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                              const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                              const double& dt_max);
+
+double acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                              const std::shared_ptr<const IDiscreteFunction>& u,
+                              const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                              const std::shared_ptr<const IDiscreteFunction>& p,
+                              const std::shared_ptr<const IDiscreteFunction>& p_k,
+                              const std::shared_ptr<const IDiscreteFunction>& Gamma_tau,
+                              const std::shared_ptr<const IDiscreteFunction>& Gamma_inf,
+                              const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                              const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                              const double& dt_max);
+
 class AcousticSolverHandler
 {
  public:
-- 
GitLab