From 93b5c2a58a10bffdb1d63150b7f477642b96f424 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Mon, 15 Dec 2025 12:44:09 +0100
Subject: [PATCH] Scale energy norm with time and porevolume

---
 opm/simulators/flow/TemperatureModel.hpp | 55 +++++++++++++++++++-----
 1 file changed, 45 insertions(+), 10 deletions(-)

diff --git a/opm/simulators/flow/TemperatureModel.hpp b/opm/simulators/flow/TemperatureModel.hpp
index 501cc525a9f..7fdeeb95ad7 100644
--- a/opm/simulators/flow/TemperatureModel.hpp
+++ b/opm/simulators/flow/TemperatureModel.hpp
@@ -36,6 +36,7 @@
 #include <opm/simulators/flow/BlackoilModelParameters.hpp>
 #include <opm/simulators/flow/GenericTemperatureModel.hpp>
 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
+#include <opm/simulators/aquifers/AquiferGridUtils.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -217,14 +218,21 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
     {
         const int max_iter = 20;
         const int min_iter = 1;
+        bool is_converged = false;
         // solve using Newton
         for (int iter = 0; iter < max_iter; ++iter) {
             assembleEquations();
             if (iter >= min_iter && converged(iter)) {
+                is_converged = true;
                 break;
             }
             solveAndUpdate();
         }
+        if (!is_converged) {
+            const auto msg = fmt::format("Temperature model (TEMP): Newton did not converge after {} iterations. \n"
+                                         "The Simulator will continue to the next step with an unconverged solution.",max_iter);
+            OpmLog::debug(msg);
+        }
     }
 
     void solveAndUpdate()
@@ -248,24 +256,51 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
 
     bool converged(const int iter)
     {
-        const unsigned int numCells = simulator_.model().numTotalDof();
+        Scalar dt = simulator_.timeStepSize();
         Scalar maxNorm = 0.0;
         Scalar sumNorm = 0.0;
+        auto tolerance_cnv_energy = Parameters::Get<Parameters::ToleranceCnvEnergy<Scalar>>();
         const auto& elemMapper = simulator_.model().elementMapper();
+        const IsNumericalAquiferCell isNumericalAquiferCell(simulator_.gridView().grid());
+        Scalar sum_pv = 0.0;
+        Scalar sum_pv_not_converged = 0.0;
         for (const auto& elem : elements(simulator_.gridView(), Dune::Partitions::interior)) {
             unsigned globI = elemMapper.index(elem);
-            maxNorm = max(maxNorm, std::abs(this->energyVector_[globI]));
-            sumNorm += std::abs(this->energyVector_[globI]);
+            const auto pvValue =  simulator_.problem().referencePorosity(globI, /*timeIdx=*/0)
+            *  simulator_.model().dofTotalVolume(globI);
+
+            const Scalar scaled_norm = dt * std::abs(this->energyVector_[globI])/ pvValue;
+            maxNorm = max(maxNorm, scaled_norm);
+            sumNorm += scaled_norm;
+            if (!isNumericalAquiferCell(elem)) {
+                if (scaled_norm > tolerance_cnv_energy) {
+                    sum_pv_not_converged += pvValue;
+                }
+                sum_pv += pvValue;
+            }
         }
         maxNorm = simulator_.gridView().comm().max(maxNorm);
         sumNorm = simulator_.gridView().comm().sum(sumNorm);
-        const int globalNumCells = simulator_.gridView().comm().sum(numCells);
-        sumNorm /= globalNumCells;
-        const auto tolerance_cnv_energy = Parameters::Get<Parameters::ToleranceCnvEnergy<Scalar>>();
-        const auto tolerance_energy_balance = Parameters::Get<Parameters::ToleranceEnergyBalance<Scalar>>();
-        if (maxNorm < tolerance_cnv_energy || sumNorm < tolerance_energy_balance) {
-            const auto msg = fmt::format("Temperature model (TEMP): Newton converged after {} iterations", iter);
-            OpmLog::debug(msg);
+        sum_pv = simulator_.gridView().comm().sum(sum_pv);
+        sumNorm /= sum_pv;
+
+        // Use relaxed tolerance if the fraction of unconverged cells porevolume is less than relaxed_max_pv_fraction
+        sum_pv_not_converged = simulator_.gridView().comm().sum(sum_pv_not_converged);
+        Scalar relaxed_max_pv_fraction = Parameters::Get<Parameters::RelaxedMaxPvFraction<Scalar>>();
+        const bool relax = (sum_pv_not_converged / sum_pv) <  relaxed_max_pv_fraction;
+        const auto tolerance_energy_balance = relax? Parameters::Get<Parameters::ToleranceEnergyBalanceRelaxed<Scalar>>():
+                                            Parameters::Get<Parameters::ToleranceEnergyBalance<Scalar>>();
+        tolerance_cnv_energy = relax? Parameters::Get<Parameters::ToleranceCnvEnergyRelaxed<Scalar>>():
+                                            Parameters::Get<Parameters::ToleranceCnvEnergy<Scalar>>();
+
+        const auto msg = fmt::format("Temperature model (TEMP): Newton iter {}: "
+                                     "CNV(E): {:.1e}, EB: {:.1e}",
+                                     iter, maxNorm, sumNorm);
+        OpmLog::debug(msg);
+        if (maxNorm < tolerance_cnv_energy && sumNorm < tolerance_energy_balance) {
+            const auto msg2 = fmt::format("Temperature model (TEMP): Newton converged after {} iterations"
+                                         , iter);
+            OpmLog::debug(msg2);
             return true;
         }
         return false;