From a57d04fde59d4cbe1b5a3d307ecaf3dcb5113a17 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Fri, 8 Dec 2023 16:49:48 +0100
Subject: [PATCH] support source term from deck (SOURCE)

---
 ebos/eclproblem.hh | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index 59e4b34da..5114e6024 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -1407,6 +1407,34 @@ public:
         if (enableAquifers_)
             aquiferModel_.addToSource(rate, globalDofIdx, timeIdx);
 
+        // Add source term from deck
+        const auto& sourceprop = this->simulator().vanguard().schedule()[this->episodeIndex()].sourceprop();
+        if (sourceprop.size() > 0) {
+            std::array<int,3> ijk;
+                this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk);
+            RateVector massRate(0.0);
+            if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                massRate[Indices::canonicalToActiveComponentIndex(oilCompIdx)] = sourceprop.rate({ijk, SourceComponent::OIL}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                massRate[Indices::canonicalToActiveComponentIndex(gasCompIdx)] = sourceprop.rate({ijk, SourceComponent::GAS}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                massRate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = sourceprop.rate({ijk, SourceComponent::WATER}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if constexpr (enableSolvent) {
+                massRate[Indices::solventSaturationIdx] = sourceprop.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx);
+            rate.setMassRate(massRate, pvtRegionIdx);
+            if constexpr (enablePolymer) {
+                rate[Indices::polymerConcentrationIdx] = sourceprop.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
+            }
+            if constexpr (enableEnergy) {
+                rate[Indices::contiEnergyEqIdx] = sourceprop.hrate(ijk) / this->model().dofTotalVolume(globalDofIdx);
+            }
+        }
+
         // if requested, compensate systematic mass loss for cells which were "well
         // behaved" in the last time step
         // Note that we don't allow for drift compensation if there are no active wells.