copying the ebos from opm-flowexperimental

CMakeLists.txt

@@ -672,6 +672,16 @@ opm_add_test(flow_distribute_z
   $<TARGET_OBJECTS:moduleVersion>
   )
 
+opm_add_test(ebos
+  ONLY_COMPILE
+  ALWAYS_ENABLE
+  DEPENDS opmsimulators
+  LIBRARIES opmsimulators
+  SOURCES
+  flowexperimental/ebos.cpp
+  $<TARGET_OBJECTS:moduleVersion>
+)
+
 if(dune-alugrid_FOUND)
   if (NOT BUILD_FLOW_ALU_GRID)
     set(FLOW_ALUGRID_ONLY_DEFAULT_ENABLE_IF "FALSE")

flowexperimental/BlackOilModelFvNoCache.hpp

@@ -0,0 +1,34 @@
+#ifndef BLACK_OIL_MODEL_FV_NOCACHE_HPP
+#define BLACK_OIL_MODEL_FV_NOCACHE_HPP
+#include <opm/simulators/flow/FIBlackoilModel.hpp>
+//#include "BlackOilModelFvFast.hpp"
+namespace Opm{
+    template<typename TypeTag>
+    class BlackOilModelFvNoCache: public FIBlackOilModel<TypeTag>{
+        using Parent = FIBlackOilModel<TypeTag>;
+        using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+        using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    public:
+        BlackOilModelFvNoCache(Simulator& simulator): FIBlackOilModel<TypeTag>(simulator){
+        }
+
+        IntensiveQuantities intensiveQuantities(unsigned globalIdx, unsigned timeIdx) const{
+            OPM_TIMEBLOCK_LOCAL(intensiveQuantitiesNoCache);
+            const auto& primaryVar = this->solution(timeIdx)[globalIdx];
+            const auto& problem = this->simulator_.problem();
+            //IntensiveQuantities* intQuants = &(this->intensiveQuantityCache_[timeIdx][globalIdx]);
+            if (!(this->enableIntensiveQuantityCache_) ||
+                !(this->intensiveQuantityCacheUpToDate_[timeIdx][globalIdx])){
+                IntensiveQuantities intQuants;
+                intQuants.update(problem,primaryVar, globalIdx, timeIdx);
+                return intQuants;// reqiored for updating extrution factor
+            }else{
+                IntensiveQuantities intQuants = (this->intensiveQuantityCache_[timeIdx][globalIdx]);
+                return intQuants;
+            }
+
+        }
+
+    };
+}
+#endif

flowexperimental/blackoilenergymodulesfv.hh

@@ -0,0 +1,141 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief Contains the classes required to extend the black-oil model by energy.
+ */
+#ifndef EWOMS_BLACK_OIL_ENERGY_MODULE_TPFA_HH
+#define EWOMS_BLACK_OIL_ENERGY_MODULE_TPFA_HH
+#include <opm/models/blackoil/blackoilenergymodules.hh>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOil
+ * \brief Contains the high level supplements required to extend the black oil
+ *        model by energy.
+ */
+// template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
+// class BlackOilEnergyModuleFV: public  BlackOilEnergyModule<TypeTag,enableEnergyV>
+// {
+
+// }
+template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
+class BlackOilEnergyIntensiveQuantitiesFV: public BlackOilEnergyIntensiveQuantities<TypeTag,enableEnergyV>
+{
+    using Parent =  BlackOilEnergyIntensiveQuantities<TypeTag, enableEnergyV>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
+    using ThermalConductionLaw  = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
+    static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
+    
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
+    static constexpr unsigned numPhases = FluidSystem::numPhases;
+public:
+    void updateTemperature_(const Problem& problem,
+                            const PrimaryVariables& priVars,
+                            unsigned globalSpaceIndex,
+                            unsigned timeIdx)
+    {
+        auto& fs = Parent::asImp_().fluidState_;
+        //const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        // set temperature
+        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx));
+    }
+    void updateEnergyQuantities_(const Problem& problem,
+                                 const PrimaryVariables& priVars,
+                                 unsigned globalSpaceIndex,
+                                 unsigned timeIdx,
+                                 const typename FluidSystem::template ParameterCache<Evaluation>& paramCache)
+    {
+        auto& fs = Parent::asImp_().fluidState_;
+        
+        // compute the specific enthalpy of the fluids, the specific enthalpy of the rock
+        // and the thermal condictivity coefficients
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+
+            const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
+            fs.setEnthalpy(phaseIdx, h);
+        }
+
+        const auto& solidEnergyLawParams = problem().solidEnergyLawParams(globalSpaceIndex, timeIdx);
+        this->rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
+
+        const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIndex, timeIdx);
+        this->totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
+
+        // Retrieve the rock fraction from the problem
+        // Usually 1 - porosity, but if pvmult is used to modify porosity
+        // we will apply the same multiplier to the rock fraction
+        // i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
+        // multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
+        this->rockFraction_ = problem.rockFraction(globalSpaceIndex, timeIdx);
+    }
+
+};
+template <class TypeTag>
+class BlackOilEnergyIntensiveQuantitiesFV<TypeTag, false>: public BlackOilEnergyIntensiveQuantities<TypeTag, false>
+    {
+        using Parent =  BlackOilEnergyIntensiveQuantities<TypeTag, false>;
+        using Problem = GetPropType<TypeTag, Properties::Problem>;
+        using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+        using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+        using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
+    public:
+        void updateTemperature_([[maybe_unused]] const Problem& problem,
+                                [[maybe_unused]] const PrimaryVariables& priVars,
+                                [[maybe_unused]] unsigned globalSpaceIdx,
+                                [[maybe_unused]] unsigned timeIdx)
+        {
+            if constexpr (enableTemperature) {
+                // even if energy is conserved, the temperature can vary over the spatial
+                // domain if the EnableTemperature property is set to true
+                //auto& fs = Parent::asImp_().fluidState_;
+                auto& fs = this->asImp_().fluidState_;
+                Scalar T = problem.temperature(globalSpaceIdx, timeIdx);
+                fs.setTemperature(T);
+            }
+        }
+        void updateEnergyQuantities_([[maybe_unused]] const Problem& problem,
+                                     [[maybe_unused]] const PrimaryVariables& priVars,
+                                     [[maybe_unused]] unsigned globalSpaceIdx,
+                                     [[maybe_unused]] unsigned timeIdx,
+                                     const typename FluidSystem::template ParameterCache<Evaluation>&)
+        { }
+
+
+};
+    
+}
+
+#endif    

flowexperimental/blackoilintensivequantitiessimple.hh

@@ -0,0 +1,634 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::BlackOilIntensiveQuantities
+ */
+#ifndef EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_SIMPLE_HH
+#define EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_SIMPLE_HH
+#include <opm/models/blackoil/blackoilproperties.hh>
+#include <opm/models/blackoil/blackoilsolventmodules.hh>
+#include <opm/models/blackoil/blackoilextbomodules.hh>
+#include <opm/models/blackoil/blackoilpolymermodules.hh>
+#include <opm/models/blackoil/blackoilfoammodules.hh>
+#include <opm/models/blackoil/blackoilbrinemodules.hh>
+#include <opm/models/blackoil/blackoilenergymodules.hh>
+#include "blackoilenergymodulesfv.hh"
+#include <opm/models/blackoil/blackoildiffusionmodule.hh>
+#include <opm/models/blackoil/blackoilmicpmodules.hh>
+#include <opm/material/fluidstates/BlackOilFluidState.hpp>
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/utility/CopyablePtr.hpp>
+
+#include <dune/common/fmatrix.hh>
+
+#include <cstring>
+#include <utility>
+
+#include <fmt/format.h>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOilModel
+ * \ingroup IntensiveQuantities
+ *
+ * \brief Contains the quantities which are are constant within a
+ *        finite volume in the black-oil model.
+ */
+template <class TypeTag>
+class BlackOilIntensiveQuantitiesSimple
+    : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
+    , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
+    , public BlackOilDiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
+    , public BlackOilSolventIntensiveQuantities<TypeTag>
+    , public BlackOilExtboIntensiveQuantities<TypeTag>
+    , public BlackOilPolymerIntensiveQuantities<TypeTag>
+    , public BlackOilFoamIntensiveQuantities<TypeTag>
+    , public BlackOilBrineIntensiveQuantities<TypeTag>
+    , public BlackOilEnergyIntensiveQuantitiesFV<TypeTag>
+    , public BlackOilMICPIntensiveQuantities<TypeTag>
+{
+    using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
+    using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using FluxModule = GetPropType<TypeTag, Properties::FluxModule>;
+
+    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+    enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+    enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
+    enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
+    enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
+    enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
+    enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
+    enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
+    enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
+    enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
+    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { dimWorld = GridView::dimensionworld };
+    enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
+
+    static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
+    static constexpr bool waterEnabled = Indices::waterEnabled;
+    static constexpr bool gasEnabled = Indices::gasEnabled;
+    static constexpr bool oilEnabled = Indices::oilEnabled;
+
+    using Toolbox = MathToolbox<Evaluation>;
+    using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
+    using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
+    using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
+
+    using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
+
+
+public:
+    using FluidState = BlackOilFluidState<Evaluation,
+                                          FluidSystem,
+                                          enableTemperature,
+                                          enableEnergy,
+                                          compositionSwitchEnabled,
+                                          enableVapwat,
+                                          enableBrine,
+                                          enableSaltPrecipitation,
+                                          false,
+                                          Indices::numPhases>;
+    using ScalarFluidState = BlackOilFluidState<Scalar,
+                                                FluidSystem,
+                                                enableTemperature,
+                                                enableEnergy,
+                                                compositionSwitchEnabled,
+                                                enableVapwat,
+                                                enableBrine,
+                                                enableSaltPrecipitation,
+                                                false,// no gas in water
+                                                Indices::numPhases>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+
+    BlackOilIntensiveQuantitiesSimple()
+    {
+        if (compositionSwitchEnabled) {
+            fluidState_.setRs(0.0);
+            fluidState_.setRv(0.0);
+        }
+        if (enableVapwat) {
+            fluidState_.setRvw(0.0);
+        }
+    }
+    BlackOilIntensiveQuantitiesSimple(const BlackOilIntensiveQuantitiesSimple& other) = default;
+
+    BlackOilIntensiveQuantitiesSimple& operator=(const BlackOilIntensiveQuantitiesSimple& other) = default;
+
+    /*!
+     * \copydoc IntensiveQuantities::update
+     */
+    void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
+    {
+        ParentType::update(elemCtx, dofIdx, timeIdx);
+
+        const auto& problem = elemCtx.problem();
+        const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        //const auto& linearizationType = problem.model().linearizer().getLinearizationType();
+        unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
+        this->update(problem,priVars,globalSpaceIdx,timeIdx);
+    }
+
+    void update(const Problem& problem,
+                const PrimaryVariables& priVars,
+                unsigned globalSpaceIdx,
+                unsigned timeIdx)
+    {
+
+        this->extrusionFactor_ = 1.0;// to avoid fixing parent update
+        OPM_TIMEBLOCK_LOCAL(UpdateIntensiveQuantities);
+        //const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        //const auto& linearizationType = problem.model().linearizer().getLinearizationType();
+        Scalar RvMax = FluidSystem::enableVaporizedOil()
+            ? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
+            : 0.0;
+        Scalar RsMax = FluidSystem::enableDissolvedGas()
+            ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
+            : 0.0;
+
+        asImp_().updateTemperature_(problem, priVars, globalSpaceIdx, timeIdx);
+
+        unsigned pvtRegionIdx = priVars.pvtRegionIndex();
+        fluidState_.setPvtRegionIndex(pvtRegionIdx);
+
+        //asImp_().updateSaltConcentration_(elemCtx, dofIdx, timeIdx);
+
+        // extract the water and the gas saturations for convenience
+        Evaluation Sw = 0.0;
+        if constexpr (waterEnabled) {
+            if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
+                Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
+            } else if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled){
+                // water is enabled but is not a primary variable i.e. one phase case
+                Sw = 1.0;
+            }
+        }
+        Evaluation Sg = 0.0;
+        if constexpr (gasEnabled) {
+            if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
+                Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            } else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
+                Sg = 1.0 - Sw;
+            } else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Disabled) {
+                if constexpr (waterEnabled) {
+                    Sg = 1.0 - Sw; // two phase water + gas
+                } else {
+                    // one phase case
+                    Sg = 1.0;
+                }
+            }
+        }
+        Valgrind::CheckDefined(Sg);
+        Valgrind::CheckDefined(Sw);
+
+        Evaluation So = 1.0 - Sw - Sg;
+
+        // deal with solvent
+        if constexpr (enableSolvent)
+            So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+
+        if (FluidSystem::phaseIsActive(waterPhaseIdx))
+            fluidState_.setSaturation(waterPhaseIdx, Sw);
+
+        if (FluidSystem::phaseIsActive(gasPhaseIdx))
+            fluidState_.setSaturation(gasPhaseIdx, Sg);
+
+        if (FluidSystem::phaseIsActive(oilPhaseIdx))
+            fluidState_.setSaturation(oilPhaseIdx, So);
+
+        //asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
+
+        // now we compute all phase pressures
+
+        std::array<Evaluation, numPhases> pC;// = {0, 0, 0};
+        computeRelpermAndPC(mobility_, pC, problem, fluidState_, globalSpaceIdx);
+        // oil is the reference phase for pressure
+        if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
+            const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                if (FluidSystem::phaseIsActive(phaseIdx))
+                    fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
+        //} else if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::Pw) {
+        //    const Evaluation& pw = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+        //   for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+        //        if (FluidSystem::phaseIsActive(phaseIdx))
+        //            fluidState_.setPressure(phaseIdx, pw + (pC[phaseIdx] - pC[waterPhaseIdx]));
+        } else {
+            //assert(FluidSystem::phaseIsActive(oilPhaseIdx));
+            const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                if (FluidSystem::phaseIsActive(phaseIdx))
+                    fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
+        }
+
+        // update the Saturation functions for the blackoil solvent module.
+        //asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
+
+        // update extBO parameters
+        //asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
+
+        Evaluation SoMax = 0.0;
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+            SoMax = max(fluidState_.saturation(oilPhaseIdx),
+                        problem.maxOilSaturation(globalSpaceIdx));
+        }
+
+        // take the meaning of the switching primary variable into account for the gas
+        // and oil phase compositions
+        if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
+            const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            fluidState_.setRs(Rs);
+        } else {
+            if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
+                OPM_TIMEBLOCK_LOCAL(UpdateSaturatedRs);
+                const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            oilPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+                fluidState_.setRs(min(RsMax, RsSat));
+            }
+            else if constexpr (compositionSwitchEnabled)
+                fluidState_.setRs(0.0);
+        }
+
+        if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
+            const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            fluidState_.setRv(Rv);
+        } else {
+            if (FluidSystem::enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
+                OPM_TIMEBLOCK_LOCAL(UpdateSaturatedRv);
+                //NB! should save the indexing for later evalustion
+                const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+                fluidState_.setRv(min(RvMax, RvSat));
+            }
+            else if constexpr (compositionSwitchEnabled)
+                fluidState_.setRv(0.0);
+        }
+
+        if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) {
+            const auto& Rvw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
+            fluidState_.setRvw(Rvw);
+        } else {
+            //NB! should save the indexing for later evaluation
+            if (FluidSystem::enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
+                OPM_TIMEBLOCK_LOCAL(UpdateSaturatedRv);
+                const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx);
+                fluidState_.setRvw(RvwSat);
+            }
+        }
+
+        // typename FluidSystem::template ParameterCache<Evaluation> paramCache;
+        // paramCache.setRegionIndex(pvtRegionIdx);
+        // if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+        //     paramCache.setMaxOilSat(SoMax);
+        // }
+        // paramCache.updateAll(fluidState_);
+
+        // compute the phase densities and transform the phase permeabilities into mobilities
+        // int nmobilities = 1;
+        // std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
+        // if (dirMob_) {
+        //     for (int i=0; i<3; i++) {
+        //         nmobilities += 1;
+        //         mobilities.push_back(&(dirMob_->getArray(i)));
+        //     }
+        // }
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+            computeInverseFormationVolumeFactorAndViscosity(fluidState_, phaseIdx, pvtRegionIdx, SoMax);
+        }
+        Valgrind::CheckDefined(mobility_);
+
+        // calculate the phase densities
+        Evaluation rho;
+        if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
+            OPM_TIMEBLOCK_LOCAL(UpdateWDensity);
+            rho = fluidState_.invB(waterPhaseIdx);
+            rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+            fluidState_.setDensity(waterPhaseIdx, rho);
+        }
+
+        if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+            OPM_TIMEBLOCK_LOCAL(UpdateGDensity);
+            rho = fluidState_.invB(gasPhaseIdx);
+            rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+            if (FluidSystem::enableVaporizedOil()) {
+                rho +=
+                    fluidState_.invB(gasPhaseIdx) *
+                    fluidState_.Rv() *
+                    FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+            }
+            if (FluidSystem::enableVaporizedWater()) {
+                rho +=
+                    fluidState_.invB(gasPhaseIdx) *
+                    fluidState_.Rvw() *
+                    FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+            }
+            fluidState_.setDensity(gasPhaseIdx, rho);
+        }
+
+        if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+            OPM_TIMEBLOCK_LOCAL(UpdateODensity);
+            rho = fluidState_.invB(oilPhaseIdx);
+            rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+            if (FluidSystem::enableDissolvedGas()) {
+                rho +=
+                    fluidState_.invB(oilPhaseIdx) *
+                    fluidState_.Rs() *
+                    FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+            }
+            fluidState_.setDensity(oilPhaseIdx, rho);
+        }
+
+        // retrieve the porosity from the problem
+        referencePorosity_ = problem.porosity(globalSpaceIdx,timeIdx);
+        porosity_ = referencePorosity_;
+        // the porosity must be modified by the compressibility of the
+        // rock...
+
+        Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
+        if (rockCompressibility > 0.0) {
+            OPM_TIMEBLOCK_LOCAL(UpdateRockCompressibility);
+            Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
+            Evaluation x;
+            if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
+            } else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
+                x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
+            } else {
+                x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
+            }
+            porosity_ *= 1.0 + x + 0.5*x*x;
+        }
+
+        // deal with water induced rock compaction
+        porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
+
+        // the MICP processes change the porosity
+        // if constexpr (enableMICP){
+        //   Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
+        //   Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
+        //   porosity_ += - biofilm_ - calcite_;
+        // }
+
+        // // deal with salt-precipitation
+        // if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
+        //     Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
+        //     porosity_ *= (1.0 - Sp);
+        // }
+
+        rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
+
+        // asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
+        // asImp_().zPvtUpdate_();
+        // asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        // asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
+        // asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        // asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        // asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+
+        // update the quantities which are required by the chosen
+        // velocity model
+//        FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
+
+        // update the diffusion specific quantities of the intensive quantities
+//        DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
+
+#ifndef NDEBUG
+        // some safety checks in debug mode
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+
+            assert(isfinite(fluidState_.density(phaseIdx)));
+            assert(isfinite(fluidState_.saturation(phaseIdx)));
+            assert(isfinite(fluidState_.temperature(phaseIdx)));
+            assert(isfinite(fluidState_.pressure(phaseIdx)));
+            assert(isfinite(fluidState_.invB(phaseIdx)));
+        }
+        assert(isfinite(fluidState_.Rs()));
+        assert(isfinite(fluidState_.Rv()));
+#endif
+    }
+
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::fluidState
+     */
+    const FluidState& fluidState() const
+    { return fluidState_; }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::mobility
+     */
+    const Evaluation& mobility(unsigned phaseIdx) const
+    { return mobility_[phaseIdx]; }
+
+    const Evaluation& mobility(unsigned phaseIdx, FaceDir::DirEnum facedir) const
+    {
+        using Dir = FaceDir::DirEnum;
+        if (dirMob_) {
+            switch(facedir) {
+                case Dir::XPlus:
+                    return dirMob_->mobilityX_[phaseIdx];
+                case Dir::YPlus:
+                    return dirMob_->mobilityY_[phaseIdx];
+                case Dir::ZPlus:
+                    return dirMob_->mobilityZ_[phaseIdx];
+                default:
+                    throw std::runtime_error("Unexpected face direction");
+            }
+        }
+        else {
+            return mobility_[phaseIdx];
+        }
+
+    }
+
+    void computeInverseFormationVolumeFactorAndViscosity(FluidState& fluidState,
+                                                         unsigned phaseIdx,
+                                                         unsigned pvtRegionIdx,
+                                                         const Evaluation& SoMax){
+        OPM_TIMEBLOCK_LOCAL(UpdateInverseFormationFactorAndViscosity);
+        {
+        OPM_TIMEBLOCK_LOCAL(UpdateFormationFactor);
+        const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx);
+        fluidState_.setInvB(phaseIdx, b);
+        }
+        {
+            OPM_TIMEBLOCK_LOCAL(UpdateViscosity);
+        typename FluidSystem::template ParameterCache<Evaluation> paramCache;
+        paramCache.setRegionIndex(pvtRegionIdx);
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+            paramCache.setMaxOilSat(SoMax);
+        }
+        paramCache.updateAll(fluidState_);
+
+        const auto& mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
+        // for (int i = 0; i<nmobilities; i++) {
+        //     if (enableExtbo && phaseIdx == oilPhaseIdx) {
+        //         (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
+        //     }
+        //     else if (enableExtbo && phaseIdx == gasPhaseIdx) {
+        //         (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
+        //     }
+        //     else {
+        //         (*mobilities[i])[phaseIdx] /= mu;
+        //     }
+        // }
+        mobility_[phaseIdx] /=mu;
+        }
+        //mobility_[phaseIdx] /= 1e-3;
+    }
+
+    void computeRelpermAndPC(std::array<Evaluation,numPhases>& mobility,
+                             std::array<Evaluation, numPhases>& pC,
+                             const Problem& problem,
+                             const FluidState& fluidState,
+                             unsigned globalSpaceIdx){
+        OPM_TIMEBLOCK_LOCAL(UpdateRelperm);
+        const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
+        MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
+        problem.updateRelperms(mobility, dirMob_, fluidState, globalSpaceIdx);
+
+        //mobility_[waterPhaseIdx] = Sw;
+        //mobility_[oilPhaseIdx] = So;
+        //mobility_[gasPhaseIdx] = Sg;
+    }
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::porosity
+     */
+    const Evaluation& porosity() const
+    { return porosity_; }
+
+    /*!
+     * The pressure-dependent transmissibility multiplier due to rock compressibility.
+     */
+    const Evaluation& rockCompTransMultiplier() const
+    { return rockCompTransMultiplier_; }
+
+    /*!
+     * \brief Returns the index of the PVT region used to calculate the thermodynamic
+     *        quantities.
+     *
+     * This allows to specify different Pressure-Volume-Temperature (PVT) relations in
+     * different parts of the spatial domain. Note that this concept should be seen as a
+     * work-around of the fact that the black-oil model does not capture the
+     * thermodynamics well enough. (Because there is, err, only a single real world with
+     * in which all substances follow the same physical laws and hence the same
+     * thermodynamics.) Anyway: Since the ECL file format uses multiple PVT regions, we
+     * support it as well in our black-oil model. (Note that, if it is not explicitly
+     * specified, the PVT region index is 0.)
+     */
+    auto pvtRegionIndex() const
+        -> decltype(std::declval<FluidState>().pvtRegionIndex())
+    { return fluidState_.pvtRegionIndex(); }
+
+    /*!
+     * \copydoc ImmiscibleIntensiveQuantities::relativePermeability
+     */
+    Evaluation relativePermeability(unsigned phaseIdx) const
+    {
+        // warning: slow
+        return fluidState_.viscosity(phaseIdx)*mobility(phaseIdx);
+    }
+
+    /*!
+     * \brief Returns the porosity of the rock at reference conditions.
+     *
+     * I.e., the porosity of rock which is not perturbed by pressure and temperature
+     * changes.
+     */
+    Scalar referencePorosity() const
+    { return referencePorosity_; }
+
+private:
+    friend BlackOilSolventIntensiveQuantities<TypeTag>;
+    friend BlackOilExtboIntensiveQuantities<TypeTag>;
+    friend BlackOilPolymerIntensiveQuantities<TypeTag>;
+    friend BlackOilEnergyIntensiveQuantitiesFV<TypeTag>;
+    friend BlackOilFoamIntensiveQuantities<TypeTag>;
+    friend BlackOilBrineIntensiveQuantities<TypeTag>;
+    friend BlackOilMICPIntensiveQuantities<TypeTag>;
+
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    FluidState fluidState_;
+    Scalar referencePorosity_;
+    Evaluation porosity_;
+    Evaluation rockCompTransMultiplier_;
+    std::array<Evaluation,numPhases> mobility_;
+
+    // Instead of writing a custom copy constructor and a custom assignment operator just to handle
+    // the dirMob_ unique ptr member variable when copying BlackOilIntensiveQuantites (see for example
+    // updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we write the below
+    // custom wrapper class CopyablePtr which wraps the unique ptr and makes it copyable.
+    //
+    // The advantage of this approach is that we avoid having to call all the base class copy constructors and
+    // assignment operators explicitly (which is needed when writing the custom copy constructor and assignment
+    // operators) which could become a maintenance burden. For example, when adding a new base class (if that should
+    // be needed sometime in the future) to BlackOilIntensiveQuantites we could forget to update the copy
+    // constructor and assignment operators.
+    //
+    // We want each copy of the BlackOilIntensiveQuantites to be unique, (TODO: why?) so we have to make a copy
+    // of the unique_ptr each time we copy construct or assign to it from another BlackOilIntensiveQuantites.
+    // (On the other hand, if a copy could share the ptr with the original, a shared_ptr could be used instead and the
+    // wrapper would not be needed)
+    DirectionalMobilityPtr dirMob_;
+};
+
+} // namespace Opm
+
+#endif

flowexperimental/ebos.cpp

@@ -0,0 +1,145 @@
+/*
+  Copyright 2020, NORCE AS
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include "config.h"
+#include <opm/simulators/flow/FlowProblem.hpp>
+#include "eclnewtonmethod.hh"
+#include "ebos.hh"
+#include <opm/simulators/flow/Main.hpp>
+
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
+#include "blackoilintensivequantitiessimple.hh"
+#include "BlackOilModelFvNoCache.hpp"
+// the current code use eclnewtonmethod adding other conditions to proceed_ should do the trick for KA
+// adding linearshe sould be chaning the update_ function in the same class with condition that the error is reduced.
+// the trick is to be able to recalculate the residual from here.
+// unsure where the timestepping is done from suggestedtime??
+// suggestTimeStep is taken from newton solver in problem.limitTimestep
+namespace Opm{
+    template<typename TypeTag>
+    class OutputAuxModule : public BaseAuxiliaryModule<TypeTag>
+    {
+
+    };
+
+}
+
+
+namespace Opm {
+namespace Properties {
+namespace TTag {
+struct EclFlowProblemEbos {
+    using InheritsFrom = std::tuple<EbosTypeTag>;
+};
+}
+
+template<class TypeTag>
+struct Model<TypeTag, TTag::EclFlowProblemEbos> {
+    using type = BlackOilModelFvNoCache<TypeTag>;
+};
+template<class TypeTag>
+struct IntensiveQuantities<TypeTag, TTag::EclFlowProblemEbos> {
+     using type = BlackOilIntensiveQuantitiesSimple<TypeTag>;
+};
+// Set the problem class
+template<class TypeTag>
+struct Problem<TypeTag, TTag::EclFlowProblemEbos> {
+    using type = EbosProblem<TypeTag>;
+};
+
+
+template<class TypeTag>
+struct ThreadsPerProcess<TypeTag, TTag::EclFlowProblemEbos> {
+    static constexpr int value = 1;
+};
+
+template<class TypeTag>
+struct ContinueOnConvergenceError<TypeTag, TTag::EclFlowProblemEbos> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EclNewtonSumTolerance<TypeTag, TTag::EclFlowProblemEbos> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1e-5;
+};
+
+// the default for the allowed volumetric error for oil per second
+template<class TypeTag>
+struct NewtonTolerance<TypeTag, TTag::EclFlowProblemEbos> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1e-2;
+};
+
+// set fraction of the pore volume where the volumetric residual may be violated during
+// strict Newton iterations
+template<class TypeTag>
+struct EclNewtonRelaxedVolumeFraction<TypeTag, TTag::EclFlowProblemEbos> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 0.0;
+};
+
+template<class TypeTag>
+struct EclNewtonRelaxedTolerance<TypeTag, TTag::EclFlowProblemEbos> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 10*getPropValue<TypeTag, Properties::NewtonTolerance>();
+};
+
+//template<class TypeTag>
+//struct Linearizer<TypeTag, TTag::EclFlowProblemEbos> { using type = TpfaLinearizer<TypeTag>; };
+
+// template<class TypeTag>
+// struct LocalResidual<TypeTag, TTag::EclFlowProblemEbos> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
+
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::EclFlowProblemEbos> { static constexpr bool value = false; };
+
+template<class TypeTag>
+struct EnableDisgasInWater<TypeTag, TTag::EclFlowProblemEbos> { static constexpr bool value = false; };
+
+//static constexpr bool has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>();
+
+template<class TypeTag>
+struct Simulator<TypeTag, TTag::EclFlowProblemEbos> { using type = Opm::Simulator<TypeTag>; };
+
+// template<class TypeTag>
+// struct LinearSolverBackend<TypeTag, TTag::EclFlowProblemEbos> {
+//     using type = ISTLSolver<TypeTag>;
+// };
+
+// // Set the problem class
+// template<class TypeTag>
+// struct Problem<TypeTag, TTag::EbosTypeTag> {
+//     using type = EbosProblem<TypeTag>;
+// };
+
+// template<class TypeTag>
+// struct EclEnableAquifers<TypeTag, TTag::EclFlowProblemTest> {
+//     static constexpr bool value = false;
+// };
+}
+}
+
+
+int main(int argc, char** argv)
+{
+    using TypeTag = Opm::Properties::TTag::EclFlowProblemEbos;
+    Opm::registerEclTimeSteppingParameters<TypeTag>();
+    return Opm::start<TypeTag>(argc, argv);
+}

flowexperimental/ebos.hh

@@ -0,0 +1,245 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief The common settings for all ebos variants.
+ */
+#ifndef EBOS_HH
+#define EBOS_HH
+
+#include <opm/simulators/flow/FlowProblem.hpp>
+#include <opm/simulators/flow/FlowProblemProperties.hpp>
+
+#include <opm/models/utils/start.hh>
+#include <opm/models/discretization/common/fvbaseproblem.hh>
+
+#include <opm/simulators/aquifers/BlackoilAquiferModel.hpp>
+#include <opm/simulators/linalg/ISTLSolver.hpp>
+#include <opm/simulators/timestepping/EclTimeSteppingParams.hpp>
+#include <opm/simulators/wells/BlackoilWellModel.hpp>
+
+namespace Opm {
+template <class TypeTag>
+class EbosProblem;
+}
+
+namespace Opm::Properties {
+
+namespace TTag {
+struct EbosTypeTag {
+    using InheritsFrom = std::tuple<FlowModelParameters, FlowBaseProblem, BlackOilModel, EclTimeSteppingParameters>;
+};
+}
+
+// Set the problem class
+template<class TypeTag>
+struct Problem<TypeTag, TTag::EbosTypeTag> {
+    using type = EbosProblem<TypeTag>;
+};
+
+// Enable experimental features for ebos: ebos is the research simulator of the OPM
+// project. If you're looking for a more stable "production quality" simulator, consider
+// using `flow`
+template<class TypeTag>
+struct EnableExperiments<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = true;
+};
+
+// use flow's well model for now
+template<class TypeTag>
+struct WellModel<TypeTag, TTag::EbosTypeTag> {
+    using type = BlackoilWellModel<TypeTag>;
+};
+
+// currently, ebos uses the non-multisegment well model by default to avoid
+// regressions. the --use-multisegment-well=true|false command line parameter is still
+// available in ebos, but hidden from view.
+template<class TypeTag>
+struct UseMultisegmentWell<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = false;
+};
+
+// set some properties that are only required by the well model
+template<class TypeTag>
+struct MatrixAddWellContributions<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = true;
+};
+
+template<class TypeTag>
+struct EnableTerminalOutput<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = false;
+};
+
+// flow's well model only works with surface volumes
+template<class TypeTag>
+struct BlackoilConserveSurfaceVolume<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = true;
+};
+
+// the values for the residual are for the whole cell instead of for a cubic meter of the cell
+template<class TypeTag>
+struct UseVolumetricResidual<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = false;
+};
+
+// by default use flow's aquifer model for now
+template<class TypeTag>
+struct AquiferModel<TypeTag, TTag::EbosTypeTag> {
+    using type = BlackoilAquiferModel<TypeTag>;
+};
+
+// use flow's linear solver backend for now
+template<class TypeTag>
+struct LinearSolverSplice<TypeTag, TTag::EbosTypeTag> {
+    using type = TTag::FlowIstlSolver;
+};
+
+template<>
+struct LinearSolverBackend<TTag::EbosTypeTag, TTag::FlowIstlSolverParams> {
+    using type = ISTLSolver<TTag::EbosTypeTag>;
+};
+
+// the default for the allowed volumetric error for oil per second
+template<class TypeTag>
+struct NewtonTolerance<TypeTag, TTag::EbosTypeTag> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1e-1;
+};
+
+// set fraction of the pore volume where the volumetric residual may be violated during
+// strict Newton iterations
+template<class TypeTag>
+struct EclNewtonRelaxedVolumeFraction<TypeTag, TTag::EbosTypeTag> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 0.05;
+};
+
+// the maximum volumetric error of a cell in the relaxed region
+template<class TypeTag>
+struct EclNewtonRelaxedTolerance<TypeTag, TTag::EbosTypeTag> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1e6*getPropValue<TypeTag, Properties::NewtonTolerance>();
+};
+
+// the tolerated amount of "incorrect" amount of oil per time step for the complete
+// reservoir. this is scaled by the pore volume of the reservoir, i.e., larger reservoirs
+// will tolerate larger residuals.
+template<class TypeTag>
+struct EclNewtonSumTolerance<TypeTag, TTag::EbosTypeTag> {
+    using type = GetPropType<TypeTag, Scalar>;
+    static constexpr type value = 1e-5;
+};
+
+// make all Newton iterations strict, i.e., the volumetric Newton tolerance must be
+// always be upheld in the majority of the spatial domain. In this context, "majority"
+// means 1 - EclNewtonRelaxedVolumeFraction.
+template<class TypeTag>
+struct EclNewtonStrictIterations<TypeTag, TTag::EbosTypeTag> {
+    static constexpr int value = 100;
+};
+
+// set the maximum number of Newton iterations to 8 so that we fail quickly (albeit
+// relatively often)
+template<class TypeTag>
+struct NewtonMaxIterations<TypeTag, TTag::EbosTypeTag> {
+    static constexpr int value = 8;
+};
+
+// if openMP is available, set the default the number of threads per process for the main
+// simulation to 2 (instead of grabbing everything that is available).
+#if _OPENMP
+template<class TypeTag>
+struct ThreadsPerProcess<TypeTag, TTag::EbosTypeTag> {
+    static constexpr int value = 2;
+};
+#endif
+
+// By default, ebos accepts the result of the time integration unconditionally if the
+// smallest time step size is reached.
+template<class TypeTag>
+struct ContinueOnConvergenceError<TypeTag, TTag::EbosTypeTag> {
+    static constexpr bool value = true;
+};
+template<class TypeTag>
+struct LinearSolverBackend<TypeTag, TTag::EbosTypeTag> {
+    using type = ISTLSolver<TypeTag>;
+};
+
+} // namespace Opm::Properties
+
+namespace Opm {
+template <class TypeTag>
+class EbosProblem : public FlowProblem<TypeTag> //, public FvBaseProblem<TypeTag>
+{
+    typedef FlowProblem<TypeTag> ParentType;
+    using BaseType = GetPropType<TypeTag, Properties::BaseProblem>;
+public:
+    void writeOutput(bool verbose = true)
+    {
+        OPM_TIMEBLOCK(problemWriteOutput);
+        // use the generic code to prepare the output fields and to
+        // write the desired VTK files.
+        if (EWOMS_GET_PARAM(TypeTag, bool, EnableWriteAllSolutions) || this->simulator().episodeWillBeOver()){
+            BaseType::writeOutput(verbose);
+        }
+    }
+
+    static void registerParameters()
+    {
+        ParentType::registerParameters();
+
+        BlackoilModelParameters<TypeTag>::registerParameters();
+        Parameters::registerParam<TypeTag, Properties::EnableTerminalOutput>("Do *NOT* use!");
+        Parameters::hideParam<TypeTag, Properties::DbhpMaxRel>();
+        Parameters::hideParam<TypeTag, Properties::DwellFractionMax>();
+        Parameters::hideParam<TypeTag, Properties::MaxResidualAllowed>();
+        Parameters::hideParam<TypeTag, Properties::ToleranceMb>();
+        Parameters::hideParam<TypeTag, Properties::ToleranceMbRelaxed>();
+        Parameters::hideParam<TypeTag, Properties::ToleranceCnv>();
+        Parameters::hideParam<TypeTag, Properties::ToleranceCnvRelaxed>();
+        Parameters::hideParam<TypeTag, Properties::ToleranceWells>();
+        Parameters::hideParam<TypeTag, Properties::ToleranceWellControl>();
+        Parameters::hideParam<TypeTag, Properties::MaxWelleqIter>();
+        Parameters::hideParam<TypeTag, Properties::UseMultisegmentWell>();
+        Parameters::hideParam<TypeTag, Properties::TolerancePressureMsWells>();
+        Parameters::hideParam<TypeTag, Properties::MaxPressureChangeMsWells>();
+        Parameters::hideParam<TypeTag, Properties::MaxInnerIterMsWells>();
+        Parameters::hideParam<TypeTag, Properties::MaxNewtonIterationsWithInnerWellIterations>();
+        Parameters::hideParam<TypeTag, Properties::MaxInnerIterWells>();
+        Parameters::hideParam<TypeTag, Properties::MaxSinglePrecisionDays>();
+        Parameters::hideParam<TypeTag, Properties::MinStrictCnvIter>();
+        Parameters::hideParam<TypeTag, Properties::MinStrictMbIter>();
+        Parameters::hideParam<TypeTag, Properties::SolveWelleqInitially>();
+        Parameters::hideParam<TypeTag, Properties::UpdateEquationsScaling>();
+        Parameters::hideParam<TypeTag, Properties::UseUpdateStabilization>();
+        Parameters::hideParam<TypeTag, Properties::MatrixAddWellContributions>();
+        Parameters::hideParam<TypeTag, Properties::EnableTerminalOutput>();
+    }
+
+    // inherit the constructors
+    using ParentType::FlowProblem;
+};
+}
+
+#endif // EBOS_HH

flowexperimental/eclnewtonmethod.hh

@@ -0,0 +1,266 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::EclNewtonMethod
+ */
+#ifndef EWOMS_ECL_NEWTON_METHOD_HH
+#define EWOMS_ECL_NEWTON_METHOD_HH
+
+#include <opm/common/Exceptions.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+
+#include <opm/models/blackoil/blackoilnewtonmethod.hh>
+#include <opm/models/utils/signum.hh>
+
+namespace Opm::Properties {
+
+template<class TypeTag, class MyTypeTag>
+struct EclNewtonSumTolerance {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct EclNewtonStrictIterations {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct EclNewtonRelaxedVolumeFraction {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct EclNewtonSumToleranceExponent {
+    using type = UndefinedProperty;
+};
+template<class TypeTag, class MyTypeTag>
+struct EclNewtonRelaxedTolerance {
+    using type = UndefinedProperty;
+};
+
+} // namespace Opm::Properties
+
+namespace Opm {
+
+/*!
+ * \brief A newton solver which is ebos specific.
+ */
+template <class TypeTag>
+class EclNewtonMethod : public BlackOilNewtonMethod<TypeTag>
+{
+    using ParentType = BlackOilNewtonMethod<TypeTag>;
+    using DiscNewtonMethod = GetPropType<TypeTag, Properties::DiscNewtonMethod>;
+
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Linearizer = GetPropType<TypeTag, Properties::Linearizer>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+
+    static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
+
+    static constexpr int contiSolventEqIdx = Indices::contiSolventEqIdx;
+    static constexpr int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
+    static constexpr int contiEnergyEqIdx = Indices::contiEnergyEqIdx;
+
+    friend NewtonMethod<TypeTag>;
+    friend DiscNewtonMethod;
+    friend ParentType;
+
+public:
+    EclNewtonMethod(Simulator& simulator) : ParentType(simulator)
+    {
+        errorPvFraction_ = 1.0;
+        relaxedMaxPvFraction_ = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonRelaxedVolumeFraction);
+
+        sumTolerance_ = 0.0; // this gets determined in the error calculation proceedure
+        relaxedTolerance_ = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonRelaxedTolerance);
+
+        numStrictIterations_ = EWOMS_GET_PARAM(TypeTag, int, EclNewtonStrictIterations);
+    }
+
+    /*!
+     * \brief Register all run-time parameters for the Newton method.
+     */
+    static void registerParameters()
+    {
+        ParentType::registerParameters();
+
+        Parameters::registerParam<TypeTag, Properties::EclNewtonSumTolerance>
+                ("The maximum error tolerated by the Newton "
+                 "method for considering a solution to be converged");
+        Parameters::registerParam<TypeTag, Properties::EclNewtonStrictIterations>
+                 ("The number of Newton iterations where the "
+                  "volumetric error is considered.");
+        Parameters::registerParam<TypeTag, Properties::EclNewtonRelaxedVolumeFraction>
+                  ("The fraction of the pore volume of the reservoir "
+                   "where the volumetric error may be violated during strict Newton iterations.");
+        Parameters::registerParam<TypeTag, Properties::EclNewtonSumToleranceExponent>
+                  ("The the exponent used to scale the sum tolerance by "
+                   "the total pore volume of the reservoir.");
+        Parameters::registerParam<TypeTag, Properties::EclNewtonRelaxedTolerance>
+                   ("The maximum error which the volumetric residual "
+                     "may exhibit if it is in a 'relaxed' region during a strict iteration.");
+    }
+
+    /*!
+     * \brief Returns true if the error of the solution is below the
+     *        tolerance.
+     */
+    bool converged() const
+    {
+        if (errorPvFraction_ < relaxedMaxPvFraction_)
+            return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ;
+        else if (this->numIterations() > numStrictIterations_)
+            return (this->error_ < relaxedTolerance_ && errorSum_ < sumTolerance_) ;
+
+        return this->error_ <= this->tolerance() && errorSum_ <= sumTolerance_;
+    }
+
+    void preSolve_(const SolutionVector&,
+                   const GlobalEqVector& currentResidual)
+    {
+        const auto& constraintsMap = this->model().linearizer().constraintsMap();
+        this->lastError_ = this->error_;
+        Scalar newtonMaxError = EWOMS_GET_PARAM(TypeTag, Scalar, NewtonMaxError);
+
+        // calculate the error as the maximum weighted tolerance of
+        // the solution's residual
+        this->error_ = 0.0;
+        Dune::FieldVector<Scalar, numEq> componentSumError;
+        std::fill(componentSumError.begin(), componentSumError.end(), 0.0);
+        Scalar sumPv = 0.0;
+        errorPvFraction_ = 0.0;
+        const Scalar dt = this->simulator_.timeStepSize();
+        for (unsigned dofIdx = 0; dofIdx < currentResidual.size(); ++dofIdx) {
+            // do not consider auxiliary DOFs for the error
+            if (dofIdx >= this->model().numGridDof()
+                || this->model().dofTotalVolume(dofIdx) <= 0.0)
+                continue;
+
+            if (!this->model().isLocalDof(dofIdx))
+                continue;
+
+            // also do not consider DOFs which are constraint
+            if (this->enableConstraints_()) {
+                if (constraintsMap.count(dofIdx) > 0)
+                    continue;
+            }
+
+            const auto& r = currentResidual[dofIdx];
+            Scalar pvValue =
+                this->simulator_.problem().referencePorosity(dofIdx, /*timeIdx=*/0)
+                * this->model().dofTotalVolume(dofIdx);
+            sumPv += pvValue;
+            bool cnvViolated = false;
+
+            Scalar dofVolume = this->model().dofTotalVolume(dofIdx);
+
+            for (unsigned eqIdx = 0; eqIdx < r.size(); ++eqIdx) {
+                Scalar tmpError = r[eqIdx] * dt * this->model().eqWeight(dofIdx, eqIdx) / pvValue;
+                Scalar tmpError2 = r[eqIdx] * this->model().eqWeight(dofIdx, eqIdx);
+
+                // in the case of a volumetric formulation, the residual in the above is
+                // per cubic meter
+                if (getPropValue<TypeTag, Properties::UseVolumetricResidual>()) {
+                    tmpError *= dofVolume;
+                    tmpError2 *= dofVolume;
+                }
+
+                this->error_ = max(std::abs(tmpError), this->error_);
+
+                if (std::abs(tmpError) > this->tolerance_)
+                    cnvViolated = true;
+
+                componentSumError[eqIdx] += std::abs(tmpError2);
+            }
+            if (cnvViolated)
+                errorPvFraction_ += pvValue;
+        }
+
+        // take the other processes into account
+        this->error_ = this->comm_.max(this->error_);
+        componentSumError = this->comm_.sum(componentSumError);
+        sumPv = this->comm_.sum(sumPv);
+        errorPvFraction_ = this->comm_.sum(errorPvFraction_);
+
+        componentSumError /= sumPv;
+        componentSumError *= dt;
+
+        errorPvFraction_ /= sumPv;
+
+        errorSum_ = 0;
+        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+            errorSum_ = std::max(std::abs(componentSumError[eqIdx]), errorSum_);
+
+        // scale the tolerance for the total error with the pore volume. by default, the
+        // exponent is 1/3, i.e., cubic root.
+        Scalar x = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumTolerance);
+        Scalar y = EWOMS_GET_PARAM(TypeTag, Scalar, EclNewtonSumToleranceExponent);
+        sumTolerance_ = x*std::pow(sumPv, y);
+
+        this->endIterMsg() << " (max: " << this->tolerance_ << ", violated for " << errorPvFraction_*100 << "% of the pore volume), aggegate error: " << errorSum_ << " (max: " << sumTolerance_ << ")";
+
+        // make sure that the error never grows beyond the maximum
+        // allowed one
+        if (this->error_ > newtonMaxError)
+            throw NumericalProblem("Newton: Error "+std::to_string(double(this->error_))
+                                   + " is larger than maximum allowed error of "
+                                   + std::to_string(double(newtonMaxError)));
+
+        // make sure that the error never grows beyond the maximum
+        // allowed one
+        if (errorSum_ > newtonMaxError)
+            throw NumericalProblem("Newton: Sum of the error "+std::to_string(double(errorSum_))
+                                   + " is larger than maximum allowed error of "
+                                   + std::to_string(double(newtonMaxError)));
+    }
+
+    void endIteration_(SolutionVector& nextSolution,
+                       const SolutionVector& currentSolution)
+    {
+        ParentType::endIteration_(nextSolution, currentSolution);
+        OpmLog::debug( "Newton iteration " + std::to_string(this->numIterations_) + ""
+                  + " error: " + std::to_string(double(this->error_))
+                  + this->endIterMsg().str());
+        this->endIterMsg().str("");
+    }
+
+private:
+    Scalar errorPvFraction_;
+    Scalar errorSum_;
+
+    Scalar relaxedTolerance_;
+    Scalar relaxedMaxPvFraction_;
+
+    Scalar sumTolerance_;
+
+    int numStrictIterations_;
+};
+} // namespace Opm
+
+#endif