From 81549a04774ffc67232e658e6d20fc085f5f20a1 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Fri, 20 Oct 2023 13:56:33 +0200
Subject: [PATCH] Add output of the residual at the end of a report step

To activte the output you will need to use
--enable-opm-rst-file=true
and
RPTRST
 RESIDUAL /
---
 ebos/eclgenericoutputblackoilmodule.cc | 32 ++++++++++++++++++++++++++
 ebos/eclgenericoutputblackoilmodule.hh |  3 +++
 ebos/ecloutputblackoilmodule.hh        | 10 ++++++++
 3 files changed, 45 insertions(+)

diff --git a/ebos/eclgenericoutputblackoilmodule.cc b/ebos/eclgenericoutputblackoilmodule.cc
index 99d91c1d1..2fd5b7daf 100644
--- a/ebos/eclgenericoutputblackoilmodule.cc
+++ b/ebos/eclgenericoutputblackoilmodule.cc
@@ -551,6 +551,28 @@ assignToSolution(data::Solution& sol)
                    data::TargetType::RESTART_OPM_EXTENDED);
     }
 
+    if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
+        ! this->residual_[waterPhaseIdx].empty())
+    {
+        sol.insert("RES_WAT", UnitSystem::measure::liquid_surface_volume,
+                   std::move(this->residual_[waterPhaseIdx]),
+                   data::TargetType::RESTART_OPM_EXTENDED);
+    }
+    if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
+        ! this->residual_[gasPhaseIdx].empty())
+    {
+        sol.insert("RES_GAS", UnitSystem::measure::gas_surface_volume,
+                   std::move(this->residual_[gasPhaseIdx]),
+                   data::TargetType::RESTART_OPM_EXTENDED);
+    }
+    if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
+        ! this->residual_[oilPhaseIdx].empty())
+    {
+        sol.insert("RES_OIL", UnitSystem::measure::liquid_surface_volume,
+                   std::move(this->residual_[oilPhaseIdx]),
+                   data::TargetType::RESTART_OPM_EXTENDED);
+    }
+
     // Fluid in place
     if (this->outputFipRestart_) {
         for (const auto& phase : Inplace::phases()) {
@@ -1135,6 +1157,16 @@ doAllocBuffers(const unsigned bufferSize,
         }
     }
 
+    if (rstKeywords["RESIDUAL"] > 0) {
+        rstKeywords["RESIDUAL"] = 0;
+        for (int phaseIdx = 0; phaseIdx <  numPhases; ++phaseIdx)
+        {
+            if (FluidSystem::phaseIsActive(phaseIdx)) {
+                this->residual_[phaseIdx].resize(bufferSize, 0.0);
+            }
+        }
+    }
+
     // ROCKC
     if (rstKeywords["ROCKC"] > 0) {
         rstKeywords["ROCKC"] = 0;
diff --git a/ebos/eclgenericoutputblackoilmodule.hh b/ebos/eclgenericoutputblackoilmodule.hh
index af65f336c..e86585d17 100644
--- a/ebos/eclgenericoutputblackoilmodule.hh
+++ b/ebos/eclgenericoutputblackoilmodule.hh
@@ -470,6 +470,9 @@ protected:
 
     std::vector<ScalarBuffer> tracerConcentrations_;
 
+    std::array<ScalarBuffer, numPhases> residual_;
+
+
     std::array<ScalarBuffer, numPhases> flowsi_;
     std::array<ScalarBuffer, numPhases> flowsj_;
     std::array<ScalarBuffer, numPhases> flowsk_;
diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh
index 6946c482c..6d803305b 100644
--- a/ebos/ecloutputblackoilmodule.hh
+++ b/ebos/ecloutputblackoilmodule.hh
@@ -278,6 +278,7 @@ public:
             return;
 
         const auto& problem = elemCtx.simulator().problem();
+        const auto& ebosResid = elemCtx.simulator().model().linearizer().residual();
         for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
             const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
             const auto& fs = intQuants.fluidState();
@@ -626,6 +627,15 @@ public:
                         tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
                 }
             }
+
+            // output residual
+            for ( int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx )
+            {
+                if (!this->residual_[phaseIdx].empty()) {
+                    const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
+                    this->residual_[phaseIdx][globalDofIdx] = ebosResid[globalDofIdx][activeCompIdx];
+                }
+            }
         }
     }