Merge pull request #6723 from svenn-t/tuningdp

Convergence check on max. solution update (TUNINGDP)

opm/simulators/flow/BlackoilModel.hpp

@@ -75,6 +75,8 @@ public:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using ModelParameters = BlackoilModelParameters<Scalar>;
 
+    static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
+
     static constexpr int numEq = Indices::numEq;
     static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
     static constexpr int contiZfracEqIdx = Indices::contiZfracEqIdx;
@@ -108,6 +110,26 @@ public:
 
     using ComponentName = ::Opm::ComponentName<FluidSystem,Indices>;
 
+    // Helper structs
+    struct CnvPvSplitData {
+        std::pair<std::vector<double>, std::vector<int>> cnvPvSplit;
+        std::vector<unsigned> ixCells;
+    };
+
+    struct MaxSolutionUpdateData {
+        Scalar dPMax = 0.0;
+        Scalar dSMax = 0.0;
+        Scalar dRsMax = 0.0;
+        Scalar dRvMax = 0.0;
+    };
+
+    // Output debug flags for which tolerances used
+    enum class DebugFlags {
+        STRICT = 0,
+        RELAXED = 1,
+        TUNINGDP = 2
+    };
+
     // ---------  Public methods  ---------
 
     /// Construct the model. It will retain references to the
@@ -177,6 +199,11 @@ public:
     /// Apply an update to the primary variables.
     void updateSolution(const BVector& dx);
 
+    /// Get solution update vector as a PrimaryVariable
+    void prepareStoringSolutionUpdate();
+    void storeSolutionUpdate(const BVector& dx);
+    MaxSolutionUpdateData getMaxSolutionUpdate(const std::vector<unsigned>& ixCells);
+
     /// Return true if output to cout is wanted.
     bool terminalOutputEnabled() const
     { return terminal_output_; }
@@ -200,9 +227,9 @@ public:
 
     /// \brief Compute pore-volume/cell count split among "converged",
     /// "relaxed converged", "unconverged" cells based on CNV point
-    /// measures.
-    std::pair<std::vector<double>, std::vector<int>>
-    characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt);
+    /// measures. Also returns list of cells where CNV is greater than
+    /// its strict tolerance
+    CnvPvSplitData characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt);
 
     /// \brief Compute the number of Newtons required by each cell in order to
     /// satisfy the solution change convergence criteria at the last time step.
@@ -213,6 +240,7 @@ public:
                             const int iteration);
 
     void updateTUNING(const Tuning& tuning);
+    void updateTUNINGDP(const TuningDp& tuning_dp);
 
     ConvergenceReport
     getReservoirConvergence(const double reportTime,
@@ -339,6 +367,8 @@ protected:
     Scalar current_relaxation_;
     BVector dx_old_;
 
+    SolutionVector solUpd_;
+
     std::vector<StepReport> convergence_reports_;
     ComponentName compNames_{};
 

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -49,6 +49,10 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     tolerance_cnv_relaxed_ = std::max(tolerance_cnv_, Parameters::Get<Parameters::ToleranceCnvRelaxed<Scalar>>());
     tolerance_cnv_energy_ = Parameters::Get<Parameters::ToleranceCnvEnergy<Scalar>>();
     tolerance_cnv_energy_relaxed_ = std::max(tolerance_cnv_energy_, Parameters::Get<Parameters::ToleranceCnvEnergyRelaxed<Scalar>>());
+    tolerance_max_dp_ = Parameters::Get<Parameters::ToleranceMaxDp<Scalar>>();
+    tolerance_max_ds_ = Parameters::Get<Parameters::ToleranceMaxDs<Scalar>>();
+    tolerance_max_drs_ = Parameters::Get<Parameters::ToleranceMaxDrs<Scalar>>();
+    tolerance_max_drv_ = Parameters::Get<Parameters::ToleranceMaxDrv<Scalar>>();
     tolerance_wells_ = Parameters::Get<Parameters::ToleranceWells<Scalar>>();
     tolerance_well_control_ = Parameters::Get<Parameters::ToleranceWellControl<Scalar>>();
     max_welleq_iter_ = Parameters::Get<Parameters::MaxWelleqIter>();
@@ -159,6 +163,22 @@ void BlackoilModelParameters<Scalar>::registerParameters()
     Parameters::Register<Parameters::ToleranceCnvEnergyRelaxed<Scalar>>
         ("Relaxed local energy convergence tolerance that applies for iterations "
          "after the iterations with the strict tolerance");
+    Parameters::Register<Parameters::ToleranceMaxDp<Scalar>>
+        ("Tolerance for max pressure change during a Newton iteration. "
+         "A value greater than 0.0 allows for convergence regardless "
+         "of residual tolerances. Use with care!");
+    Parameters::Register<Parameters::ToleranceMaxDs<Scalar>>
+        ("Tolerance for max saturation change during a Newton iteration. "
+         "A value greater than 0.0 allows for convergence regardless "
+         "of residual tolerances. Use with care!");
+    Parameters::Register<Parameters::ToleranceMaxDrs<Scalar>>
+        ("Tolerance for max RS change during a Newton iteration. "
+         "A value greater than 0.0 allows for convergence regardless "
+         "of residual tolerances. Use with care!");
+    Parameters::Register<Parameters::ToleranceMaxDrv<Scalar>>
+        ("Tolerance for max RV change during a Newton iteration. "
+         "A value greater than 0.0 allows for convergence regardless "
+         "of residual tolerances. Use with care!");
     Parameters::Register<Parameters::ToleranceWells<Scalar>>
         ("Well convergence tolerance");
     Parameters::Register<Parameters::ToleranceWellControl<Scalar>>

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -56,7 +56,7 @@ template<class Scalar>
 struct ToleranceMbRelaxed { static constexpr Scalar value = 1e-6; };
 
 //TODO change to a simpler number with fewer digits
-//converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water 
+//converting J -> RM3 (entalpy / (cp * deltaK * rho) assuming change of 1e-5K of water
 template<class Scalar>
 struct ToleranceEnergyBalance { static constexpr Scalar value = 1e-7*41.82; };
 
@@ -75,6 +75,18 @@ struct ToleranceCnvEnergy { static constexpr Scalar value = 1e-2*41.82; };
 template<class Scalar>
 struct ToleranceCnvEnergyRelaxed { static constexpr Scalar value = 1.0*41.82; };
 
+template<class Scalar>
+struct ToleranceMaxDp { static constexpr Scalar value = 0.0; };
+
+template<class Scalar>
+struct ToleranceMaxDs { static constexpr Scalar value = 0.0; };
+
+template<class Scalar>
+struct ToleranceMaxDrs { static constexpr Scalar value = 0.0; };
+
+template<class Scalar>
+struct ToleranceMaxDrv { static constexpr Scalar value = 0.0; };
+
 template<class Scalar>
 struct ToleranceWells { static constexpr Scalar value = 1e-4; };
 
@@ -212,6 +224,14 @@ public:
     Scalar tolerance_cnv_energy_;
     /// Relaxed local energy convergence tolerance (can be used when iter >= min_strict_cnv_iter_ && cnvViolatedPV < relaxed_max_pv_fraction_).
     Scalar tolerance_cnv_energy_relaxed_;
+    /// Max pressure change during a Newton iteration (TUNINGDP item = TRGDDP)
+    Scalar tolerance_max_dp_;
+    /// Max saturation change during a Newton iteration (TUNINGDP item = TRGDDS)
+    Scalar tolerance_max_ds_;
+    /// Max RS change during a Newton iteration (TUNINGDP item = TRGDDRS)
+    Scalar tolerance_max_drs_;
+    /// Max RV change during a Newton iteration (TUNINGDP item = TRGDDRV)
+    Scalar tolerance_max_drv_;
     /// Well convergence tolerance.
     Scalar tolerance_wells_;
     /// Tolerance for the well control equations

opm/simulators/flow/BlackoilModel_impl.hpp

@@ -507,6 +507,12 @@ BlackoilModel<TypeTag>::
 updateSolution(const BVector& dx)
 {
     OPM_TIMEBLOCK(updateSolution);
+    // Prepare to store solution update for convergence check
+    if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0
+        || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0) {
+        prepareStoringSolutionUpdate();
+    }
+
     auto& newtonMethod = simulator_.model().newtonMethod();
     SolutionVector& solution = simulator_.model().solution(/*timeIdx=*/0);
 
@@ -522,6 +528,103 @@ updateSolution(const BVector& dx)
         OPM_TIMEBLOCK(invalidateAndUpdateIntensiveQuantities);
         simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
     }
+
+    // Store solution update
+    if (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0
+        || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_  > 0.0) {
+        storeSolutionUpdate(dx);
+    }
+}
+
+template <class TypeTag>
+void
+BlackoilModel<TypeTag>::
+prepareStoringSolutionUpdate()
+{
+    // Init. solution update vector
+    unsigned nc = simulator_.model().numGridDof();
+    solUpd_.resize(nc);
+
+    const auto& elemMapper = simulator_.model().elementMapper();
+    const auto& gridView = simulator_.gridView();
+    for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
+        // Copy solution vector to transfer primary variables meaning
+        unsigned globalElemIdx = elemMapper.index(elem);
+        solUpd_[globalElemIdx] = simulator_.model().solution(/*timeIdx=*/0)[globalElemIdx];
+
+        // Ensure each element is zero
+        std::fill(solUpd_[globalElemIdx].begin(), solUpd_[globalElemIdx].end(), 0.0);
+    }
+}
+
+template <class TypeTag>
+void
+BlackoilModel<TypeTag>::
+storeSolutionUpdate(const BVector& dx)
+{
+    const auto& elemMapper = simulator_.model().elementMapper();
+    const auto& gridView = simulator_.gridView();
+    for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
+        // Get cell vectors
+        unsigned globalElemIdx = elemMapper.index(elem);
+        auto& value = solUpd_[globalElemIdx];
+        const auto& update = dx[globalElemIdx];
+        assert(value.size() == update.size());
+
+        // Transfer update from dx to solution update container (SolutionVector type)
+        std::copy(update.begin(), update.end(), value.begin());
+    }
+}
+
+template <class TypeTag>
+typename BlackoilModel<TypeTag>::MaxSolutionUpdateData
+BlackoilModel<TypeTag>::
+getMaxSolutionUpdate(const std::vector<unsigned>& ixCells)
+{
+    static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
+    static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
+
+    // Init output
+    Scalar dPMax = 0.0;
+    Scalar dSMax = 0.0;
+    Scalar dRsMax = 0.0;
+    Scalar dRvMax = 0.0;
+
+    // Loop over solution update, get the correct variables and calculate max.
+    for (const auto& ix : ixCells) {
+        const auto& value = solUpd_[ix];
+        for (int pvIdx = 0; pvIdx < value.size(); ++pvIdx) {
+            if (pvIdx == Indices::pressureSwitchIdx) {
+                dPMax = std::max(dPMax, std::abs(value[pvIdx]));
+            }
+            else if ( (pvIdx == Indices::waterSwitchIdx
+                       && value.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
+                      || (pvIdx == Indices::compositionSwitchIdx
+                          && value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
+                      || (enableSolvent && pvIdx == Indices::solventSaturationIdx
+                          && value.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss)
+                      || (enableBrine && enableSaltPrecipitation && pvIdx == Indices::saltConcentrationIdx
+                          && value.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) ) {
+                dSMax = std::max(dSMax, std::abs(value[pvIdx]));
+            }
+            else if (pvIdx == Indices::compositionSwitchIdx
+                     && value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
+                dRsMax = std::max(dRsMax, std::abs(value[pvIdx]));
+            }
+            else if (pvIdx == Indices::compositionSwitchIdx
+                     && value.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
+                dRvMax = std::max(dRvMax, std::abs(value[pvIdx]));
+            }
+        }
+    }
+
+    // Communicate max values
+    dPMax = this->grid_.comm().max(dPMax);
+    dSMax = this->grid_.comm().max(dSMax);
+    dRsMax = this->grid_.comm().max(dRsMax);
+    dRvMax = this->grid_.comm().max(dRvMax);
+
+    return { dPMax, dSMax, dRsMax, dRvMax };
 }
 
 template <class TypeTag>
@@ -637,7 +740,7 @@ localConvergenceData(std::vector<Scalar>& R_sum,
 }
 
 template <class TypeTag>
-std::pair<std::vector<double>, std::vector<int>>
+typename BlackoilModel<TypeTag>::CnvPvSplitData
 BlackoilModel<TypeTag>::
 characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
 {
@@ -677,6 +780,8 @@ characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
 
     ElementContext elemCtx(this->simulator());
 
+    std::vector<unsigned> ixCells;
+
     OPM_BEGIN_PARALLEL_TRY_CATCH();
     for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
         // Skip cells of numerical Aquifer
@@ -697,6 +802,13 @@ characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
 
         splitPV[ix] += static_cast<double>(pvValue);
         ++cellCntPV[ix];
+
+        // For dP and dS check, we need cell indices of [1] violations
+        if ( ix > 0 &&
+             (this->param_.tolerance_max_dp_ > 0.0 || this->param_.tolerance_max_ds_ > 0.0
+              || this->param_.tolerance_max_drs_ > 0.0 || this->param_.tolerance_max_drv_ > 0.0 ) ) {
+            ixCells.push_back(cell_idx);
+        }
     }
 
     OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::characteriseCnvPvSplit() failed: ",
@@ -705,7 +817,7 @@ characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt)
     this->grid_.comm().sum(splitPV  .data(), splitPV  .size());
     this->grid_.comm().sum(cellCntPV.data(), cellCntPV.size());
 
-    return cnvPvSplit;
+    return { cnvPvSplit, ixCells };
 }
 
 template <class TypeTag>
@@ -721,6 +833,18 @@ updateTUNING(const Tuning& tuning)
     this->param_.newton_min_iter_ = tuning.NEWTMN;
 }
 
+template <class TypeTag>
+void
+BlackoilModel<TypeTag>::
+updateTUNINGDP(const TuningDp& tuning_dp)
+{
+    // NOTE: If TUNINGDP item is _not_ set it should be 0.0
+    this->param_.tolerance_max_dp_ = tuning_dp.TRGDDP;
+    this->param_.tolerance_max_ds_ = tuning_dp.TRGDDS;
+    this->param_.tolerance_max_drs_ = tuning_dp.TRGDDRS;
+    this->param_.tolerance_max_drv_ = tuning_dp.TRGDDRV;
+}
+
 template <class TypeTag>
 ConvergenceReport
 BlackoilModel<TypeTag>::
@@ -752,7 +876,8 @@ getReservoirConvergence(const double reportTime,
                                    numAquiferPvSumLocal,
                                    R_sum, maxCoeff, B_avg);
 
-    report.setCnvPoreVolSplit(this->characteriseCnvPvSplit(B_avg, dt),
+    auto cnvSplitData = this->characteriseCnvPvSplit(B_avg, dt);
+    report.setCnvPoreVolSplit(cnvSplitData.cnvPvSplit,
                               pvSum - numAquiferPvSum);
 
     // For each iteration, we need to determine whether to use the
@@ -765,9 +890,12 @@ getReservoirConvergence(const double reportTime,
     const bool relax_final_iteration_mb =
         this->param_.min_strict_mb_iter_ < 0 && iteration == maxIter;
 
-    const bool use_relaxed_mb = relax_final_iteration_mb ||
-        (this->param_.min_strict_mb_iter_ >= 0 &&
-         iteration >= this->param_.min_strict_mb_iter_);
+    const bool relax_iter_mb =
+        this->param_.min_strict_mb_iter_ >= 0 &&
+        iteration >= this->param_.min_strict_mb_iter_;
+
+    const bool use_relaxed_mb = relax_final_iteration_mb
+        || relax_iter_mb;
 
     // If min_strict_cnv_iter = -1 we use a relaxed tolerance for
     // the cnv for the last iterations.  For positive values we use
@@ -792,33 +920,41 @@ getReservoirConvergence(const double reportTime,
 
         // [1]: tol < cnv <= relaxed
         // [2]: relaxed < cnv
-        return static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]) <
-            this->param_.relaxed_max_pv_fraction_ * eligible;
+        Scalar cnvPvSum = static_cast<Scalar>(cnvPvSplit[1] + cnvPvSplit[2]);
+        return cnvPvSum < this->param_.relaxed_max_pv_fraction_ * eligible &&
+            cnvPvSum > 0.0;
     }();
 
-    const bool use_relaxed_cnv = relax_final_iteration_cnv
-        || relax_pv_fraction_cnv
-        || relax_iter_cnv;
-
-    if ((relax_final_iteration_mb || relax_final_iteration_cnv) &&
-        this->terminal_output_)
-    {
-        std::string message =
-            "Number of newton iterations reached its maximum "
-            "try to continue with relaxed tolerances:";
-
-        if (relax_final_iteration_mb) {
-            message += fmt::format("  MB: {:.1e}", param_.tolerance_mb_relaxed_);
-        }
-
-        if (relax_final_iteration_cnv) {
-            message += fmt::format(" CNV: {:.1e}", param_.tolerance_cnv_relaxed_);
-        }
-
-        OpmLog::debug(message);
+    // If tolerances for solution changes are met, we use the
+    // relaxed cnv tolerance. Note that all tolerances > 0.0
+    // must be met to enable use of relaxed cnv tolerance.
+    MaxSolutionUpdateData maxSolUpd;
+    const bool use_dp_tol = this->param_.tolerance_max_dp_ > 0.0;
+    const bool use_ds_tol = this->param_.tolerance_max_ds_ > 0.0;
+    const bool use_drs_tol = this->param_.tolerance_max_drs_ > 0.0;
+    const bool use_drv_tol = this->param_.tolerance_max_drv_ > 0.0;
+    const bool use_dsol_tol = use_dp_tol || use_ds_tol || use_drs_tol || use_drv_tol;
+    bool relax_dsol_cnv = false;
+    if (iteration > 0 && use_dsol_tol) {
+        maxSolUpd = getMaxSolutionUpdate(cnvSplitData.ixCells);
+        relax_dsol_cnv =
+            (!use_dp_tol || (maxSolUpd.dPMax > 0.0 && maxSolUpd.dPMax < this->param_.tolerance_max_dp_)) &&
+            (!use_ds_tol || (maxSolUpd.dSMax > 0.0 && maxSolUpd.dSMax < this->param_.tolerance_max_ds_)) &&
+            (!use_drs_tol || (maxSolUpd.dRsMax > 0.0 && maxSolUpd.dRsMax < this->param_.tolerance_max_drs_)) &&
+            (!use_drv_tol || (maxSolUpd.dRvMax > 0.0 && maxSolUpd.dRvMax < this->param_.tolerance_max_drv_));
     }
 
-    const auto tol_cnv = use_relaxed_cnv ? param_.tolerance_cnv_relaxed_ : param_.tolerance_cnv_;
+    // Determine if relaxed CNV tolerances should be used
+    const bool use_relaxed_cnv = relax_final_iteration_cnv
+        || relax_pv_fraction_cnv
+        || relax_iter_cnv
+        || relax_dsol_cnv;
+
+    // Ensure that CNV convergence criteria is met when max.
+    // solution change tolerances have been fulfilled
+    Scalar tolerance_cnv_relaxed = relax_dsol_cnv ? 1e20 : param_.tolerance_cnv_relaxed_;
+
+    const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : param_.tolerance_cnv_;
     const auto tol_mb  = use_relaxed_mb ? param_.tolerance_mb_relaxed_ : param_.tolerance_mb_;
     const auto tol_cnv_energy = use_relaxed_cnv ? param_.tolerance_cnv_energy_relaxed_ : param_.tolerance_cnv_energy_;
     const auto tol_eb = use_relaxed_mb ? param_.tolerance_energy_balance_relaxed_ : param_.tolerance_energy_balance_;
@@ -897,6 +1033,16 @@ getReservoirConvergence(const double reportTime,
                 msg += ") ";
             }
 
+            if (use_dsol_tol) {
+                msg += use_dp_tol ? "    DP     " : "";
+                msg += use_ds_tol ? "    DS     " : "";
+                msg += use_drs_tol ? "    DRS    " : "";
+                msg += use_drv_tol ? "    DRV    " : "";
+            }
+
+            msg += " MBFLAG";
+            msg += " CNVFLAG";
+
             OpmLog::debug(msg);
         }
 
@@ -913,6 +1059,31 @@ getReservoirConvergence(const double reportTime,
             ss << std::setw(11) << CNV[compIdx];
         }
 
+        if (use_dsol_tol) {
+            auto print_dsol =
+            [&] (bool use_tol, Scalar dsol) {
+                if (!use_tol) {
+                    return;
+                }
+                if (iteration == 0 || dsol <= 0.0) {
+                    ss << std::string(5, ' ') << "-" << std::string(5, ' ');
+                }
+                else {
+                    ss << std::setw(11) << dsol;
+                }
+            };
+            print_dsol(use_dp_tol, maxSolUpd.dPMax);
+            print_dsol(use_ds_tol, maxSolUpd.dSMax);
+            print_dsol(use_drs_tol, maxSolUpd.dRsMax);
+            print_dsol(use_drv_tol, maxSolUpd.dRvMax);
+        }
+
+        int mb_flag = use_relaxed_mb ? static_cast<int>(DebugFlags::RELAXED) : static_cast<int>(DebugFlags::STRICT);
+        int cnv_flag = relax_dsol_cnv ? static_cast<int>(DebugFlags::TUNINGDP) :
+            (use_relaxed_cnv ? static_cast<int>(DebugFlags::RELAXED) : static_cast<int>(DebugFlags::STRICT));
+        ss << std::setw(7) << mb_flag;
+        ss << std::setw(8) << cnv_flag;
+
         ss.precision(oprec);
         ss.flags(oflags);
 

opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp

@@ -355,6 +355,26 @@ public:
         }
     }
 
+    void updateTUNINGDP(const TuningDp& tuning_dp)
+    {
+        // NOTE: If TUNINGDP item is _not_ set it should be 0.0
+        modelParam_.tolerance_max_dp_ = tuning_dp.TRGDDP;
+        modelParam_.tolerance_max_ds_ = tuning_dp.TRGDDS;
+        modelParam_.tolerance_max_drs_ = tuning_dp.TRGDDRS;
+        modelParam_.tolerance_max_drv_ = tuning_dp.TRGDDRV;
+
+        // Terminal warnings
+        if (terminalOutput_) {
+            // Warnings unsupported items
+            if (tuning_dp.TRGLCV_has_value) {
+                OpmLog::warning("TUNINGDP item 1 (TRGLCV) is not supported.");
+            }
+            if (tuning_dp.XXXLCV_has_value) {
+                OpmLog::warning("TUNINGDP item 2 (XXXLCV) is not supported.");
+            }
+        }
+    }
+
     bool runStep(SimulatorTimer& timer)
     {
         if (schedule().exitStatus().has_value()) {
@@ -449,6 +469,19 @@ public:
                     result = max_next_tstep > 0;
                 }
 
+                if (events.hasEvent(ScheduleEvents::TUNINGDP_CHANGE)) {
+                    // Unset the event to not trigger it again on the next sub step
+                    schedule.clear_event(ScheduleEvents::TUNINGDP_CHANGE, reportStep);
+
+                    // Update TUNINGDP parameters
+                    // NOTE: Need to update both solver (model) and simulator since solver is re-created each report
+                    // step.
+                    const auto& sched_state = schedule[reportStep];
+                    const auto& tuning_dp = sched_state.tuning_dp();
+                    solver_->model().updateTUNINGDP(tuning_dp);
+                    this->updateTUNINGDP(tuning_dp);
+                }
+
                 const auto& wcycle = schedule[reportStep].wcycle.get();
                 if (wcycle.empty()) {
                     return result;

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -631,7 +631,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"TRNHD", {true, std::nullopt}},
         {"TRPLPORO", {true, std::nullopt}},
         {"TRROCK", {true, std::nullopt}},
-        {"TUNINGDP", {false, std::nullopt}},
         {"TUNINGH", {false, std::nullopt}},
         {"TUNINGL", {false, std::nullopt}},
         {"TUNINGS", {false, std::nullopt}},