adding a compositional simulator

flowexp_comp

CMakeLists.txt

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

CMakeLists_files.cmake

@@ -758,8 +758,10 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/flow/FlowMain.hpp
   opm/simulators/flow/FlowProblem.hpp
   opm/simulators/flow/FlowProblemBlackoil.hpp
+  opm/simulators/flow/FlowProblemComp.hpp
   opm/simulators/flow/FlowProblemParameters.hpp
   opm/simulators/flow/FlowProblemBlackoilProperties.hpp
+  opm/simulators/flow/FlowProblemCompProperties.hpp
   opm/simulators/flow/FlowUtils.hpp
   opm/simulators/flow/FlowsData.hpp
   opm/simulators/flow/FlowThresholdPressure.hpp

flowexperimental/comp/flowexp_comp.cpp

@@ -0,0 +1,306 @@
+/*
+  Copyright 2024, SINTEF Digital
+
+  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/models/utils/start.hh>
+#include <opm/material/constraintsolvers/PTFlash.hpp>
+#include "../FlowExpNewtonMethod.hpp"
+#include <opm/models/ptflash/flashmodel.hh>
+#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/models/discretization/common/baseauxiliarymodule.hh>
+
+#include <opm/simulators/flow/FlowProblemComp.hpp>
+#include <opm/simulators/flow/FlowProblemCompProperties.hpp>
+// TODO: not understanding why we need FlowGenericProblem here
+#include <opm/simulators/flow/FlowGenericProblem.hpp>
+#include <opm/simulators/flow/FlowGenericProblem_impl.hpp>
+#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
+
+// // 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 EmptyModel : public BaseAuxiliaryModule<TypeTag>
+    {
+        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        using GridView = GetPropType<TypeTag, Properties::GridView>;
+        using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
+        using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
+    public:
+        using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+        EmptyModel(Simulator& /*simulator*/){
+        };
+        void init(){}
+        template<class Something>
+        void init(Something /*A*/){}
+        void prepareTracerBatches(){};
+        using NeighborSet = std::set<unsigned>;
+        void linearize(SparseMatrixAdapter& /*matrix*/, GlobalEqVector& /*residual*/){};
+        unsigned numDofs() const{return 0;};
+        void addNeighbors(std::vector<NeighborSet>& /*neighbors*/) const{};
+        //void applyInitial(){};
+        void initialSolutionApplied(){};
+        //void initFromRestart(const data::Aquifers& aquiferSoln);
+        template <class Restarter>
+        void serialize(Restarter& /*res*/){};
+
+        template <class Restarter>
+        void deserialize(Restarter& /*res*/){};
+
+        void beginEpisode(){};
+        void beginTimeStep(){};
+        void beginIteration(){};
+        // add the water rate due to aquifers to the source term.
+        template<class RateVector, class Context>
+        void addToSource(RateVector& rates, const Context& context, unsigned spaceIdx, unsigned timeIdx) const{};
+        template<class RateVector>
+        void addToSource(RateVector& rates, unsigned globalSpaceIdx, unsigned timeIdx) const{};
+        void endIteration()const{};
+        void endTimeStep(){};
+        void endEpisode(){};
+        void applyInitial(){};
+        template<class RateType>
+        void computeTotalRatesForDof(RateType& /*rate*/, unsigned /*globalIdx*/) const{};
+    };
+
+}
+
+
+namespace Opm::Properties {
+
+   namespace TTag {
+   struct FlowExpCompProblem {
+       using InheritsFrom = std::tuple<FlowBaseProblemComp, FlashModel>;
+   };
+   }
+
+    template<class TypeTag>
+    struct SparseMatrixAdapter<TypeTag, TTag::FlowExpCompProblem>
+    {
+    private:
+        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+        using Block = MatrixBlock<Scalar, numEq, numEq>;
+
+    public:
+        using type = typename Linear::IstlSparseMatrixAdapter<Block>;
+    };
+
+#if 0
+template<class TypeTag>
+struct SolidEnergyLaw<TypeTag, TTag::FlowExpCompProblem>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    using EclThermalLawManager = ::Opm::EclThermalLawManager<Scalar, FluidSystem>;
+
+    using type = typename EclThermalLawManager::SolidEnergyLaw;
+};
+
+// Set the material law for thermal conduction
+template<class TypeTag>
+struct ThermalConductionLaw<TypeTag, TTag::FlowExpCompProblem>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+public:
+    using EclThermalLawManager = ::Opm::EclThermalLawManager<Scalar, FluidSystem>;
+
+    using type = typename EclThermalLawManager::ThermalConductionLaw;
+};
+
+
+template <class TypeTag>
+struct SpatialDiscretizationSplice<TypeTag, TTag::FlowExpCompProblem>
+{
+    using type = TTag::EcfvDiscretization;
+};
+
+template <class TypeTag>
+struct LocalLinearizerSplice<TypeTag, TTag::FlowExpCompProblem>
+{
+    using type = TTag::AutoDiffLocalLinearizer;
+};
+#endif
+
+// Set the problem property
+template <class TypeTag>
+struct Problem<TypeTag, TTag::FlowExpCompProblem>
+{
+    using type = FlowProblemComp<TypeTag>;
+};
+
+template<class TypeTag>
+struct AquiferModel<TypeTag, TTag::FlowExpCompProblem> {
+    using type = EmptyModel<TypeTag>;
+};
+
+template<class TypeTag>
+struct WellModel<TypeTag, TTag::FlowExpCompProblem> {
+    using type = EmptyModel<TypeTag>;
+};
+
+template<class TypeTag>
+struct TracerModel<TypeTag, TTag::FlowExpCompProblem> {
+    using type = EmptyModel<TypeTag>;
+};
+
+
+template <class TypeTag>
+struct FlashSolver<TypeTag, TTag::FlowExpCompProblem> {
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+
+public:
+    using type = Opm::PTFlash<Scalar, FluidSystem>;
+};
+
+
+template <class TypeTag, class MyTypeTag>
+struct NumComp { using type = UndefinedProperty; };
+
+// TODO: this is unfortunate, have to check why we need to hard-code it
+template <class TypeTag>
+struct NumComp<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr int value = 3;
+};
+#if 0
+struct Temperature { using type = UndefinedProperty; };
+
+ template <class TypeTag>
+ struct Temperature<TypeTag, TTag::FlowExpCompProblem> {
+     using type = GetPropType<TypeTag, Scalar>;
+     static constexpr type value = 423.25;
+ };
+#endif
+
+template <class TypeTag>
+struct FluidSystem<TypeTag, TTag::FlowExpCompProblem>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
+
+public:
+    using type = Opm::GenericOilGasFluidSystem<Scalar, num_comp>;
+};
+template<class TypeTag>
+struct EnableMech<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnableDisgasInWater<TypeTag, TTag::FlowExpCompProblem> { static constexpr bool value = false; };
+
+template<class TypeTag>
+struct Stencil<TypeTag, TTag::FlowExpCompProblem>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+
+public:
+    using type = EcfvStencil<Scalar, GridView>;
+};
+
+template<class TypeTag>
+struct EnableApiTracking<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnableTemperature<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnableSaltPrecipitation<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
+struct EnablePolymerMW<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnablePolymer<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnableBrine<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
+struct EnableVapwat<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+template<class TypeTag>
+struct EnableSolvent<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
+struct EnableEnergy<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
+struct EnableFoam<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
+struct EnableExtbo<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
+struct EnableMICP<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+
+// disable thermal flux boundaries by default
+#if 0
+template<class TypeTag>
+struct EnableThermalFluxBoundaries<TypeTag, TTag::FlowExpCompProblem> {
+    static constexpr bool value = false;
+};
+#endif
+
+} // namespace Opm::Properties
+
+
+int main(int argc, char** argv)
+{
+    using TypeTag = Opm::Properties::TTag::FlowExpCompProblem;
+    Opm::registerEclTimeSteppingParameters<double>();
+    return Opm::start<TypeTag>(argc, argv);
+}

opm/models/ptflash/flashindices.hh

@@ -40,6 +40,7 @@ namespace Opm {
  *
  * \tparam PVOffset The first index in a primary variable vector.
  */
+// TODO: The indices class should handle whether phase is active, not the FluidSystem
 template <class TypeTag, int PVOffset>
 class FlashIndices
     : public EnergyIndices<PVOffset + getPropValue<TypeTag, Properties::NumComponents>(),
@@ -50,6 +51,11 @@ class FlashIndices
     using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
 
 public:
+    static constexpr bool waterEnabled = false;
+    static constexpr bool gasEnabled = true;
+    static constexpr bool oilEnabled = true;
+    static constexpr int waterPhaseIdx = -1;
+    static constexpr int numPhases = 2;
     //! number of equations/primary variables
     static const int numEq = numComponents + EnergyIndices::numEq_;
 

opm/simulators/flow/FlowProblemComp.hpp

@@ -0,0 +1,468 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  Copyright 2024 SINTEF Digital
+
+  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::FlowProblemComp
+ */
+#ifndef OPM_FLOW_PROBLEM_COMP_HPP
+#define OPM_FLOW_PROBLEM_COMP_HPP
+
+
+#include <opm/simulators/flow/FlowProblem.hpp>
+#include <opm/material/fluidstates/CompositionalFluidState.hpp>
+
+#include <opm/material/thermal/EclThermalLawManager.hpp>
+
+
+#include <algorithm>
+#include <functional>
+#include <set>
+#include <string>
+#include <vector>
+
+namespace Opm {
+
+/*!
+ * \ingroup CompositionalSimulator
+ *
+ * \brief This problem simulates an input file given in the data format used by the
+ *        commercial ECLiPSE simulator.
+ */
+template <class TypeTag>
+class FlowProblemComp : public FlowProblem<TypeTag>
+{
+    // TODO: the naming of the Types will be adjusted
+    using FlowProblemType = FlowProblem<TypeTag>;
+
+    using typename FlowProblemType::Scalar;
+    using typename FlowProblemType::Simulator;
+    using typename FlowProblemType::GridView;
+    using typename FlowProblemType::FluidSystem;
+    using typename FlowProblemType::Vanguard;
+
+    // might not be needed
+    using FlowProblemType::dim;
+    using FlowProblemType::dimWorld;
+
+    using FlowProblemType::numPhases;
+    using FlowProblemType::numComponents;
+
+    using FlowProblemType::gasPhaseIdx;
+    using FlowProblemType::oilPhaseIdx;
+    using FlowProblemType::waterPhaseIdx;
+
+    using typename FlowProblemType::Indices;
+    using typename FlowProblemType::PrimaryVariables;
+    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
+    using typename FlowProblemType::Evaluation;
+    using typename FlowProblemType::MaterialLaw;
+    using typename FlowProblemType::RateVector;
+
+    using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
+
+public:
+    using FlowProblemType::porosity;
+    using FlowProblemType::pvtRegionIndex;
+
+    /*!
+     * \copydoc FvBaseProblem::registerParameters
+     */
+    static void registerParameters()
+    {
+        FlowProblemType::registerParameters();
+    }
+
+
+    /*!
+     * \copydoc Doxygen::defaultProblemConstructor
+     */
+    explicit FlowProblemComp(Simulator& simulator)
+        : FlowProblemType(simulator)
+    {
+    }
+
+    /*!
+     * \copydoc FvBaseProblem::finishInit
+     */
+    void finishInit()
+    {
+        // TODO: there should be room to remove duplication for this function,
+        // but there is relatively complicated logic in the function calls in this function
+        // some refactoring is needed for this function
+        FlowProblemType::finishInit();
+
+        auto& simulator = this->simulator();
+
+        auto finishTransmissibilities = [updated = false, this]() mutable {
+            if (updated) {
+                return;
+            }
+            this->transmissibilities_.finishInit(
+                [&vg = this->simulator().vanguard()](const unsigned int it) { return vg.gridIdxToEquilGridIdx(it); });
+            updated = true;
+        };
+
+        finishTransmissibilities();
+
+        const auto& eclState = simulator.vanguard().eclState();
+        const auto& schedule = simulator.vanguard().schedule();
+
+        // Set the start time of the simulation
+        simulator.setStartTime(schedule.getStartTime());
+        simulator.setEndTime(schedule.simTime(schedule.size() - 1));
+
+        // We want the episode index to be the same as the report step index to make
+        // things simpler, so we have to set the episode index to -1 because it is
+        // incremented by endEpisode(). The size of the initial time step and
+        // length of the initial episode is set to zero for the same reason.
+        simulator.setEpisodeIndex(-1);
+        simulator.setEpisodeLength(0.0);
+
+        // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false
+        // disables gravity, else the standard value of the gravity constant at sea level
+        // on earth is used
+        this->gravity_ = 0.0;
+        if (Parameters::Get<Parameters::EnableGravity>())
+            this->gravity_[dim - 1] = 9.80665;
+        if (!eclState.getInitConfig().hasGravity())
+            this->gravity_[dim - 1] = 0.0;
+
+        if (this->enableTuning_) {
+            // if support for the TUNING keyword is enabled, we get the initial time
+            // steping parameters from it instead of from command line parameters
+            const auto& tuning = schedule[0].tuning();
+            this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0;
+            this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
+        }
+
+        this->initFluidSystem_();
+
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
+            && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+            this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0);
+        }
+
+        this->readRockParameters_(simulator.vanguard().cellCenterDepths(), [&simulator](const unsigned idx) {
+            std::array<int, dim> coords;
+            simulator.vanguard().cartesianCoordinate(idx, coords);
+            for (auto& c : coords) {
+                ++c;
+            }
+            return coords;
+        });
+        FlowProblemType::readMaterialParameters_();
+        FlowProblemType::readThermalParameters_();
+
+        const auto& initconfig = eclState.getInitConfig();
+        if (initconfig.restartRequested())
+            readEclRestartSolution_();
+        else
+            this->readInitialCondition_();
+
+        FlowProblemType::updatePffDofData_();
+
+        if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
+            const auto& vanguard = this->simulator().vanguard();
+            const auto& gridView = vanguard.gridView();
+            int numElements = gridView.size(/*codim=*/0);
+            this->polymer_.maxAdsorption.resize(numElements, 0.0);
+        }
+
+        /* readBoundaryConditions_();
+
+        // compute and set eq weights based on initial b values
+        computeAndSetEqWeights_();
+
+        if (enableDriftCompensation_) {
+            drift_.resize(this->model().numGridDof());
+            drift_ = 0.0;
+        } */
+
+        // TODO: check wether the following can work with compostional
+        if (enableVtkOutput_ && eclState.getIOConfig().initOnly()) {
+            simulator.setTimeStepSize(0.0);
+            FlowProblemType::writeOutput(true);
+        }
+
+        // after finishing the initialization and writing the initial solution, we move
+        // to the first "real" episode/report step
+        // for restart the episode index and start is already set
+        if (!initconfig.restartRequested()) {
+            simulator.startNextEpisode(schedule.seconds(1));
+            simulator.setEpisodeIndex(0);
+            simulator.setTimeStepIndex(0);
+        }
+    }
+
+    /*!
+     * \copydoc FvBaseProblem::boundary
+     *
+     * Reservoir simulation uses no-flow conditions as default for all boundaries.
+     */
+    template <class Context>
+    void boundary(BoundaryRateVector& values,
+                  const Context& context,
+                  unsigned spaceIdx,
+                  unsigned /* timeIdx */) const
+    {
+        OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
+        if (!context.intersection(spaceIdx).boundary())
+            return;
+
+        values.setNoFlow();
+
+        if (this->nonTrivialBoundaryConditions()) {
+            throw std::logic_error("boundary condition is not supported by compostional modeling yet");
+        }
+    }
+
+    /*!
+     * \copydoc FvBaseProblem::initial
+     *
+     * The reservoir problem uses a constant boundary condition for
+     * the whole domain.
+     */
+    template <class Context>
+    void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const
+    {
+        const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
+        const auto& initial_fs = initialFluidStates_[globalDofIdx];
+        Opm::CompositionalFluidState<Evaluation, FluidSystem> fs;
+        using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
+        for (unsigned p = 0; p < numPhases; ++p) { // TODO: assuming the phaseidx continuous
+            ComponentVector evals;
+            auto& last_eval = evals[numComponents - 1];
+            last_eval = 1.;
+            for (unsigned c = 0; c < numComponents - 1; ++c) {
+                const auto val = initial_fs.moleFraction(p, c);
+                const Evaluation eval = Evaluation::createVariable(val, c + 1);
+                evals[c] = eval;
+                last_eval -= eval;
+            }
+            for (unsigned c = 0; c < numComponents; ++c) {
+                fs.setMoleFraction(p, c, evals[c]);
+            }
+
+            // pressure
+            const auto p_val = initial_fs.pressure(p);
+            fs.setPressure(p, Evaluation::createVariable(p_val, 0));
+
+            const auto sat_val = initial_fs.saturation(p);
+            fs.setSaturation(p, sat_val);
+
+            const auto temp_val = initial_fs.temperature(p);
+            fs.setTemperature(temp_val);
+        }
+
+        {
+            typename FluidSystem::template ParameterCache<Evaluation> paramCache;
+            paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
+            paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
+            fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
+            fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
+            fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
+            fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
+        }
+
+        // Set initial K and L
+        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+            const Evaluation Ktmp = fs.wilsonK_(compIdx);
+            fs.setKvalue(compIdx, Ktmp);
+        }
+
+        const Evaluation& Ltmp = -1.0;
+        fs.setLvalue(Ltmp);
+
+        values.assignNaive(fs);
+    }
+
+    void addToSourceDense(RateVector&, unsigned, unsigned) const
+    {
+        // we do nothing for now
+    }
+
+    const InitialFluidState& initialFluidState(unsigned globalDofIdx) const
+    { return initialFluidStates_[globalDofIdx]; }
+
+    std::vector<InitialFluidState>& initialFluidStates()
+    { return initialFluidStates_; }
+
+    const std::vector<InitialFluidState>& initialFluidStates() const
+    { return initialFluidStates_; }
+
+    // TODO: do we need this one?
+    template<class Serializer>
+    void serializeOp(Serializer& serializer)
+    {
+        serializer(static_cast<FlowProblemType&>(*this));
+    }
+protected:
+
+    void updateExplicitQuantities_(int /* episodeIdx*/, int /* timeStepSize */, bool /* first_step_after_restart */) override
+    {
+        // we do nothing here for now
+    }
+
+    void readEquilInitialCondition_()
+    {
+        throw std::logic_error("Equilibration is not supported by compositional modeling yet");
+    }
+
+    void readEclRestartSolution_()
+    {
+        throw std::logic_error("Restarting is not supported by compositional modeling yet");
+    }
+
+    void readExplicitInitialCondition_()
+    {
+        readExplicitInitialConditionCompositional_();
+    }
+    
+    void readExplicitInitialConditionCompositional_()
+    {
+        const auto& simulator = this->simulator();
+        const auto& vanguard = simulator.vanguard();
+        const auto& eclState = vanguard.eclState();
+        const auto& fp = eclState.fieldProps();
+        // const bool has_xmf = fp.has_double("XMF");
+        assert(fp.has_double("XMF"));
+        // const bool has_ymf = fp.has_double("YMF");
+        assert(fp.has_double("YMF"));
+        const bool has_temp = fp.has_double("TEMPI");
+        const bool has_pressure = fp.has_double("PRESSURE");
+
+        // const bool has_gas = fp.has_double("SGAS");
+        assert(fp.has_double("SGAS"));
+        if (!has_pressure)
+            throw std::runtime_error("The ECL input file requires the presence of the PRESSURE "
+                                      "keyword if the model is initialized explicitly");
+
+        std::size_t numDof = this->model().numGridDof();
+
+        initialFluidStates_.resize(numDof);
+
+        std::vector<double> xmfData;
+        std::vector<double> ymfData;
+        std::vector<double> waterSaturationData;
+        std::vector<double> gasSaturationData;
+        std::vector<double> soilData;
+        std::vector<double> pressureData;
+        std::vector<double> tempiData;
+
+        if (waterPhaseIdx > 0 && Indices::numPhases > 1)
+            waterSaturationData = fp.get_double("SWAT");
+        else
+            waterSaturationData.resize(numDof);
+
+        pressureData = fp.get_double("PRESSURE");
+
+        assert(waterPhaseIdx < 0);
+
+        xmfData = fp.get_double("XMF");
+        ymfData = fp.get_double("YMF");
+        if (has_temp) {
+            tempiData = fp.get_double("TEMPI");
+        } else {
+            ; // TODO: throw?
+        }
+
+        if (gasPhaseIdx > 0) // && FluidSystem::phaseIsActive(oilPhaseIdx))
+            gasSaturationData = fp.get_double("SGAS");
+        else
+            gasSaturationData.resize(numDof);
+
+        // constexpr std::size_t num_comps = 3;
+
+        for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
+            auto& dofFluidState = initialFluidStates_[dofIdx];
+            // dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx));
+
+            Scalar temperatureLoc = tempiData[dofIdx];
+            assert(std::isfinite(temperatureLoc) && temperatureLoc > 0);
+            dofFluidState.setTemperature(temperatureLoc);
+
+            if (gasPhaseIdx > 0) {
+                    dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
+                                                gasSaturationData[dofIdx]);
+            }
+            if (oilPhaseIdx > 0) {
+                dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
+                                            1.0
+                                            - waterSaturationData[dofIdx]
+                                            - gasSaturationData[dofIdx]);
+            }
+
+            //////
+            // set phase pressures
+            //////
+            const Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase)
+
+            // TODO: zero capillary pressure for now
+            const std::array<Scalar, numPhases> pc = {0};
+            for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                if (!FluidSystem::phaseIsActive(phaseIdx))
+                    continue;
+
+                if (Indices::oilEnabled)
+                    dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
+                else if (Indices::gasEnabled)
+                    dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
+                else if (Indices::waterEnabled)
+                    //single (water) phase
+                    dofFluidState.setPressure(phaseIdx, pressure);
+            }
+
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                const std::size_t data_idx = compIdx * numDof + dofIdx;
+                const Scalar xmf = xmfData[data_idx];
+                const Scalar ymf = xmfData[data_idx];
+
+                dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
+                dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
+            }
+        }
+    }
+
+private:
+
+    void handleSolventBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
+    {
+        throw std::logic_error("solvent is disabled for compositional modeling and you're trying to add solvent to BC");
+    }
+
+    void handlePolymerBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
+    {
+        throw std::logic_error("polymer is disabled for compositional modeling and you're trying to add polymer to BC");
+    }
+
+    std::vector<InitialFluidState> initialFluidStates_;
+
+    bool enableVtkOutput_;
+};
+
+} // namespace Opm
+
+#endif // OPM_FLOW_PROBLEM_COMP_HPP

opm/simulators/flow/FlowProblemCompProperties.hpp

@@ -0,0 +1,96 @@
+// -*- 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::FlowBaseProblemComp
+ */
+#ifndef OPM_FLOW_PROBLEM_COMP_PROPERTIES_HPP
+#define OPM_FLOW_PROBLEM_COMP_PROPERTIES_HPP
+
+
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
+
+#include <opm/models/utils/propertysystem.hh>
+
+#include <opm/simulators/flow/FlowBaseProblemProperties.hpp>
+
+#include <tuple>
+
+namespace Opm {
+template <class TypeTag>
+class FlowProblemComp;
+}
+
+namespace Opm::Properties {
+
+namespace TTag {
+
+struct FlowBaseProblemComp {
+    using InheritsFrom = std::tuple<FlowBaseProblem>;
+};
+
+}
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::FlowBaseProblemComp>
+{ using type = FlowProblemComp<TypeTag>; };
+
+template<class TypeTag>
+struct TracerModel<TypeTag, TTag::FlowBaseProblemComp> {
+    using type = ::Opm::TracerModel<TypeTag>;
+};
+
+// Set the material law for fluid fluxes
+template<class TypeTag>
+struct MaterialLaw<TypeTag, TTag::FlowBaseProblemComp>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+    // using Traits = ThreePhaseMaterialTraits<Scalar,
+    //                                         /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
+    //                                         /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
+    //                                         /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
+
+    // TODO: We should be able to use FluidSystem here and using Indices to handle the active phases
+    // some more development is needed
+    using Traits = ThreePhaseMaterialTraits<Scalar,
+                                            /*wettingPhaseIdx=*/ 0,
+                                            /*nonWettingPhaseIdx=*/ 1,
+                                            /*gasPhaseIdx=*/ 2>;
+
+public:
+    using EclMaterialLawManager = ::Opm::EclMaterialLawManager<Traits>;
+
+    using type = typename EclMaterialLawManager::MaterialLaw;
+};
+
+// Enable diffusion
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::FlowBaseProblemComp>
+{ static constexpr bool value = false; };
+
+} // namespace Opm::Properties
+#endif // OPM_FLOW_PROBLEM_COMP_PROPERTIES_HPP