Add option for drift compensation for TEMP

opm/simulators/flow/FlowGenericProblem.hpp

@@ -351,6 +351,7 @@ class FlowGenericProblem
     int numPressurePointsEquil_;
 
     bool enableDriftCompensation_;
+    bool enableDriftCompensationTemp_{false};
     bool explicitRockCompaction_;
 
     // To lookup origin cell indices

opm/simulators/flow/FlowProblem.hpp

@@ -241,6 +241,11 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
             relpermDiagnostics.diagnosis(simulator.vanguard().eclState(),
                                          simulator.vanguard().levelCartesianIndexMapper());
         }
+
+        if (energyModuleType == EnergyModules::SequentialImplicitThermal) {
+            this->enableDriftCompensationTemp_ = Parameters::Get<Parameters::EnableDriftCompensationTemp>();
+        }
+
     }
 
     virtual ~FlowProblem() = default;
@@ -434,14 +439,11 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
         this->wellModel_.endTimeStep();
         this->aquiferModel_.endTimeStep();
         this->tracerModel_.endTimeStep();
-        if constexpr(energyModuleType == EnergyModules::SequentialImplicitThermal) {
-            this->temperatureModel_.endTimeStep(wellModel_.wellState());
-        }
 
         // Compute flux for output
         this->model().linearizer().updateFlowsInfo();
 
-        if (this->enableDriftCompensation_) {
+        if (this->enableDriftCompensation_ || this->enableDriftCompensationTemp_) {
             OPM_TIMEBLOCK(driftCompansation);
 
             const auto& residual = this->model().linearizer().residual();
@@ -455,6 +457,11 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
                 }
             }
         }
+
+        // Drift compensation needs to be updated before calling the temperature equation
+        if constexpr(energyModuleType == EnergyModules::SequentialImplicitThermal) {
+            this->temperatureModel_.endTimeStep(wellModel_.wellState());
+        }
     }
 
     /*!
@@ -1216,6 +1223,11 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
         serializer(*materialLawManager_);
     }
 
+    const GlobalEqVector& drift() const
+    {
+        return drift_;
+    }
+
 private:
     Implementation& asImp_()
     { return *static_cast<Implementation *>(this); }

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -410,7 +410,7 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         // compute and set eq weights based on initial b values
         this->computeAndSetEqWeights_();
 
-        if (this->enableDriftCompensation_) {
+        if (this->enableDriftCompensation_ || this->enableDriftCompensationTemp_) {
             this->drift_.resize(this->model().numGridDof());
             this->drift_ = 0.0;
         }

opm/simulators/flow/FlowProblemParameters.cpp

@@ -53,6 +53,9 @@ void registerFlowProblemParameters()
     Parameters::Register<Parameters::EnableDriftCompensation>
         ("Enable partial compensation of systematic mass losses via "
          "the source term of the next time step");
+    Parameters::Register<Parameters::EnableDriftCompensationTemp>
+        ("Enable compensation of systematic mass losses "
+         "in the energy equation (only TEMP option)");
     Parameters::Register<Parameters::OutputMode>
         ("Specify which messages are going to be printed. "
          "Valid values are: none, log, all (default)");

opm/simulators/flow/FlowProblemParameters.hpp

@@ -36,6 +36,10 @@ namespace Opm::Parameters {
 // the source term of the next time step
 struct EnableDriftCompensation { static constexpr bool value = false; };
 
+// Enable compensation of systematic mass losses in
+// the sequential energy equation
+struct EnableDriftCompensationTemp { static constexpr bool value = true; };
+
 // implicit or explicit pressure in rock compaction
 struct ExplicitRockCompaction { static constexpr bool value = false; };
 

opm/simulators/flow/TemperatureModel.hpp

@@ -408,6 +408,24 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
             this->assembleEquationWell(*wellPtr);
         }
 
+        const auto& problem = simulator_.problem();
+
+        bool enableDriftCompensation = Parameters::Get<Parameters::EnableDriftCompensationTemp>();
+        if (enableDriftCompensation) {
+            const auto& model = simulator_.model();
+            for (unsigned globalDofIdx = 0; globalDofIdx < numCells; ++globalDofIdx) {
+                Scalar dt = simulator_.timeStepSize();
+                auto dofDriftRate = problem.drift()[globalDofIdx]/dt;
+                const auto& fs = intQuants_[globalDofIdx].fluidState();
+                for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+                   const unsigned activeCompIdx =
+                        FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
+                   auto drift_hrate = dofDriftRate[activeCompIdx]*getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx)) /  getValue(fs.invB(phaseIdx));
+                   this->energyVector_[globalDofIdx] -= drift_hrate*getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+                }
+            }
+        }
+
         if (simulator_.gridView().comm().size() > 1) {
             // Set dirichlet conditions for overlapping cells
             // loop over precalculated overlap rows and columns