moved cd to xxx_impl.hh to avoid

this way we do not have to add explicit initiation of all
variants of grids and gridView iterators in the library

ebos/eclgenericproblem.cc

@@ -21,9 +21,40 @@
   copyright holders.
 */
 #include <config.h>
+
 #include <ebos/eclgenericproblem.hh>
 #include <ebos/eclgenericproblem_impl.hh>
 
-template class Opm::EclGenericProblem<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
+#include <dune/grid/common/defaultgridview.hh>
-                                      Opm::BlackOilFluidSystem<double,Opm::BlackOilDefaultIndexTraits>,
+#include <dune/grid/common/gridview.hh>
-                                      double>;
+
+#include <opm/grid/CpGrid.hpp>
+
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+#include <dune/fem/gridpart/common/gridpart2gridview.hh>
+#include <ebos/femcpgridcompat.hh>
+#endif // HAVE_DUNE_FEM
+
+namespace Opm {
+
+#if HAVE_DUNE_FEM
+template class EclGenericProblem<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
+                                 BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
+                                 double>;
+template class EclGenericProblem<Dune::Fem::GridPart2GridViewImpl<
+                                    Dune::Fem::AdaptiveLeafGridPart<
+                                        Dune::CpGrid,
+                                        Dune::PartitionIteratorType(4),
+                                        false> >,
+                                 BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
+                                 double>;
+#else
+template class EclGenericProblem<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
+                                 BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
+                                 double>;
+#endif
+
+} // end namespace Opm

ebos/eclgenericproblem_impl.hh

@@ -27,14 +27,14 @@
 
 #include <dune/common/parametertree.hh>
 
+#include <ebos/eclsolutioncontainers.hh>
+
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/OverburdTable.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/RockwnodTable.hpp>
 #include <opm/input/eclipse/Schedule/Schedule.hpp>
 #include <opm/input/eclipse/Units/Units.hpp>
 
-#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
-
 #include <boost/date_time.hpp>
 
 #include <algorithm>
@@ -97,18 +97,12 @@ serializationTestObject(const EclipseState& eclState,
 {
     EclGenericProblem result(eclState, schedule, gridView);
     result.maxOilSaturation_ = {1.0, 2.0};
-    result.maxPolymerAdsorption_ = {3.0, 4.0, 5.0};
     result.maxWaterSaturation_ = {6.0};
     result.minOilPressure_ = {7.0, 8.0, 9.0, 10.0};
     result.overburdenPressure_ = {11.0};
-    result.polymerConcentration_ = {12.0};
-    result.polymerMoleWeight_ = {13.0, 14.0};
     result.solventSaturation_ = {15.0};
-    result.microbialConcentration_ = {16.0};
+    result.polymer_ = PolymerSolutionContainer<Scalar>::serializationTestObject();
-    result.oxygenConcentration_ = {17.0};
+    result.micp_ = MICPSolutionContainer<Scalar>::serializationTestObject();
-    result.ureaConcentration_ = {18.0};
-    result.biofilmConcentration_ = {19.0};
-    result.calciteConcentration_ = {20.0};
     result.lastRv_ = {21.0};
     result.maxDRv_ = {22.0, 23.0};
     result.convectiveDrs_ = {24.0, 25.0, 26.0};
@@ -558,40 +552,47 @@ readBlackoilExtentionsInitialConditions_(size_t numDof,
     }
 
     if (enablePolymer) {
-        if (eclState_.fieldProps().has_double("SPOLY"))
+        if (eclState_.fieldProps().has_double("SPOLY")) {
-            polymerConcentration_ = eclState_.fieldProps().get_double("SPOLY");
+            polymer_.concentration = eclState_.fieldProps().get_double("SPOLY");
-        else
+        } else {
-            polymerConcentration_.resize(numDof, 0.0);
+            polymer_.concentration.resize(numDof, 0.0);
+        }
     }
 
     if (enablePolymerMolarWeight) {
-        if (eclState_.fieldProps().has_double("SPOLYMW"))
+        if (eclState_.fieldProps().has_double("SPOLYMW")) {
-            polymerMoleWeight_ = eclState_.fieldProps().get_double("SPOLYMW");
+            polymer_.moleWeight = eclState_.fieldProps().get_double("SPOLYMW");
-        else
+        } else {
-            polymerMoleWeight_.resize(numDof, 0.0);
+            polymer_.moleWeight.resize(numDof, 0.0);
+        }
     }
 
     if (enableMICP) {
-        if (eclState_.fieldProps().has_double("SMICR"))
+        if (eclState_.fieldProps().has_double("SMICR")) {
-            microbialConcentration_ = eclState_.fieldProps().get_double("SMICR");
+            micp_.microbialConcentration = eclState_.fieldProps().get_double("SMICR");
-        else
+        } else {
-            microbialConcentration_.resize(numDof, 0.0);
+            micp_.microbialConcentration.resize(numDof, 0.0);
-        if (eclState_.fieldProps().has_double("SOXYG"))
+        }
-            oxygenConcentration_ = eclState_.fieldProps().get_double("SOXYG");
+        if (eclState_.fieldProps().has_double("SOXYG")) {
-        else
+            micp_.oxygenConcentration = eclState_.fieldProps().get_double("SOXYG");
-            oxygenConcentration_.resize(numDof, 0.0);
+        } else {
-        if (eclState_.fieldProps().has_double("SUREA"))
+            micp_.oxygenConcentration.resize(numDof, 0.0);
-            ureaConcentration_ = eclState_.fieldProps().get_double("SUREA");
+        }
-        else
+        if (eclState_.fieldProps().has_double("SUREA")) {
-            ureaConcentration_.resize(numDof, 0.0);
+            micp_.ureaConcentration = eclState_.fieldProps().get_double("SUREA");
-        if (eclState_.fieldProps().has_double("SBIOF"))
+        } else {
-            biofilmConcentration_ = eclState_.fieldProps().get_double("SBIOF");
+            micp_.ureaConcentration.resize(numDof, 0.0);
-        else
+        }
-            biofilmConcentration_.resize(numDof, 0.0);
+        if (eclState_.fieldProps().has_double("SBIOF")) {
-        if (eclState_.fieldProps().has_double("SCALC"))
+            micp_.biofilmConcentration = eclState_.fieldProps().get_double("SBIOF");
-            calciteConcentration_ = eclState_.fieldProps().get_double("SCALC");
+        } else {
-        else
+            micp_.biofilmConcentration.resize(numDof, 0.0);
-            calciteConcentration_.resize(numDof, 0.0);
+        }
+        if (eclState_.fieldProps().has_double("SCALC")) {
+            micp_.calciteConcentration = eclState_.fieldProps().get_double("SCALC");
+        } else {
+            micp_.calciteConcentration.resize(numDof, 0.0);
+        }
 }
 }
 
@@ -655,70 +656,77 @@ template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 polymerConcentration(unsigned elemIdx) const
 {
-    if (polymerConcentration_.empty())
+    if (polymer_.concentration.empty()) {
         return 0;
+    }
 
-    return polymerConcentration_[elemIdx];
+    return polymer_.concentration[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 polymerMolecularWeight(const unsigned elemIdx) const
 {
-    if (polymerMoleWeight_.empty())
+    if (polymer_.moleWeight.empty()) {
         return 0.0;
+    }
 
-    return polymerMoleWeight_[elemIdx];
+    return polymer_.moleWeight[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 microbialConcentration(unsigned elemIdx) const
 {
-    if (microbialConcentration_.empty())
+    if (micp_.microbialConcentration.empty()) {
         return 0;
+    }
 
-    return microbialConcentration_[elemIdx];
+    return micp_.microbialConcentration[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 oxygenConcentration(unsigned elemIdx) const
 {
-    if (oxygenConcentration_.empty())
+    if (micp_.oxygenConcentration.empty()) {
         return 0;
+    }
 
-    return oxygenConcentration_[elemIdx];
+    return micp_.oxygenConcentration[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 ureaConcentration(unsigned elemIdx) const
 {
-    if (ureaConcentration_.empty())
+    if (micp_.ureaConcentration.empty()) {
         return 0;
+    }
 
-    return ureaConcentration_[elemIdx];
+    return micp_.ureaConcentration[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 biofilmConcentration(unsigned elemIdx) const
 {
-    if (biofilmConcentration_.empty())
+    if (micp_.biofilmConcentration.empty()) {
         return 0;
+    }
 
-    return biofilmConcentration_[elemIdx];
+    return micp_.biofilmConcentration[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 calciteConcentration(unsigned elemIdx) const
 {
-    if (calciteConcentration_.empty())
+    if (micp_.calciteConcentration.empty()) {
         return 0;
+    }
 
-    return calciteConcentration_[elemIdx];
+    return micp_.calciteConcentration[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
@@ -765,10 +773,11 @@ template<class GridView, class FluidSystem, class Scalar>
 Scalar EclGenericProblem<GridView,FluidSystem,Scalar>::
 maxPolymerAdsorption(unsigned elemIdx) const
 {
-    if (maxPolymerAdsorption_.empty())
+    if (polymer_.maxAdsorption.empty()) {
         return 0;
+    }
 
-    return maxPolymerAdsorption_[elemIdx];
+    return polymer_.maxAdsorption[elemIdx];
 }
 
 template<class GridView, class FluidSystem, class Scalar>
@@ -801,18 +810,12 @@ bool EclGenericProblem<GridView,FluidSystem,Scalar>::
 operator==(const EclGenericProblem& rhs) const
 {
     return this->maxOilSaturation_ == rhs.maxOilSaturation_ &&
-           this->maxPolymerAdsorption_ == rhs.maxPolymerAdsorption_ &&
            this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
            this->minOilPressure_ == rhs.minOilPressure_ &&
            this->overburdenPressure_ == rhs.overburdenPressure_ &&
-           this->polymerConcentration_ == rhs.polymerConcentration_ &&
-           this->polymerMoleWeight_ == rhs.polymerMoleWeight_ &&
            this->solventSaturation_ == rhs.solventSaturation_ &&
-           this->microbialConcentration_ == rhs.microbialConcentration_ &&
+           this->polymer_ == rhs.polymer_ &&
-           this->oxygenConcentration_ == rhs.oxygenConcentration_ &&
+           this->micp_ == rhs.micp_ &&
-           this->ureaConcentration_ == rhs.ureaConcentration_ &&
-           this->biofilmConcentration_ == rhs.biofilmConcentration_ &&
-           this->calciteConcentration_ == rhs.calciteConcentration_ &&
            this->lastRv_ == rhs.lastRv_ &&
            this->maxDRv_ == rhs.maxDRv_ &&
            this->convectiveDrs_ == rhs.convectiveDrs_ &&

ebos/eclgenericthresholdpressure.cc

@@ -24,22 +24,9 @@
 #include <config.h>
 #include <limits>
 #include <ebos/eclgenericthresholdpressure.hh>
-
+#include <ebos/eclgenericthresholdpressure_impl.hh>
-#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
-#include <opm/input/eclipse/EclipseState/Tables/Eqldims.hpp>
-#include <opm/input/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
-#include <opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
-
-#include <dune/grid/common/mcmgmapper.hh>
 
 #include <opm/grid/CpGrid.hpp>
-#include <opm/grid/polyhedralgrid.hh>
-#if HAVE_DUNE_ALUGRID
-#include <dune/alugrid/grid.hh>
-#include <dune/alugrid/3d/gridview.hh>
-#include "alucartesianindexmapper.hh"
-#endif // HAVE_DUNE_ALUGRID
 
 #if HAVE_DUNE_FEM
 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -47,272 +34,8 @@
 #include <ebos/femcpgridcompat.hh>
 #endif // HAVE_DUNE_FEM
 
-#include <fmt/format.h>
-
-#include <algorithm>
-#include <cassert>
-#include <stdexcept>
-
 namespace Opm {
 
-template<class Grid, class GridView, class ElementMapper, class Scalar>
-EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper,
-                            const GridView& gridView,
-                            const ElementMapper& elementMapper,
-                            const EclipseState& eclState)
-    : cartMapper_(cartMapper)
-    , gridView_(gridView)
-    , elementMapper_(elementMapper)
-    , eclState_(eclState)
-{
-}
-
-template<class Grid, class GridView, class ElementMapper,class Scalar>
-Scalar EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-thresholdPressure(int elem1Idx, int elem2Idx) const
-{
-    if (!enableThresholdPressure_)
-        return 0.0;
-
-    // threshold pressure accross faults
-    if (!thpresftValues_.empty()) {
-        int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx);
-        int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx);
-
-        assert(0 <= cartElem1Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem1Idx);
-        assert(0 <= cartElem2Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem2Idx);
-
-        int fault1Idx = cartElemFaultIdx_[cartElem1Idx];
-        int fault2Idx = cartElemFaultIdx_[cartElem2Idx];
-        if (fault1Idx != -1 && fault1Idx == fault2Idx)
-            // inside a fault there's no threshold pressure, even accross EQUIL
-            // regions.
-            return 0.0;
-        if (fault1Idx != fault2Idx) {
-            // TODO: which value if a cell is part of multiple faults? we take
-            // the maximum here.
-            Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0;
-            Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0;
-            return std::max(val1, val2);
-        }
-    }
-
-    // threshold pressure accross EQUIL regions
-    auto equilRegion1Idx = elemEquilRegion_[elem1Idx];
-    auto equilRegion2Idx = elemEquilRegion_[elem2Idx];
-
-    if (equilRegion1Idx == equilRegion2Idx)
-        return 0.0;
-
-    return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
-}
-
-template<class Grid, class GridView, class ElementMapper, class Scalar>
-void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-finishInit()
-{
-    unsigned numElements = gridView_.size(/*codim=*/0);
-
-    const auto& simConfig = eclState_.getSimulationConfig();
-
-    enableThresholdPressure_ = simConfig.useThresholdPressure();
-    if (!enableThresholdPressure_)
-        return;
-
-    numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions();
-    const decltype(numEquilRegions_) maxRegions =
-        std::numeric_limits<std::decay_t<decltype(elemEquilRegion_[0])>>::max();
-
-    if (numEquilRegions_ > maxRegions) {
-        // make sure that the index of an equilibration region can be stored
-        // in the vector
-        OPM_THROW(std::invalid_argument,
-                  (fmt::format("The maximum number of supported "
-                               "equilibration regions by OPM flow is {}, but "
-                               "{} are used!",
-                               maxRegions, numEquilRegions_)));
-    }
-
-    if (numEquilRegions_ > 2048) {
-        // warn about performance
-        OpmLog::warning(fmt::format("Number of equilibration regions is {}, which is "
-                                 "rather large. Note, that this might "
-                                 "have a negative impact on performance "
-                                 "and memory consumption.", numEquilRegions_));
-    }
-
-    // internalize the data specified using the EQLNUM keyword
-    const auto& fp = eclState_.fieldProps();
-    const auto& equilRegionData = fp.get_int("EQLNUM");
-    elemEquilRegion_.resize(numElements, 0);
-    for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
-        elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1;
-    }
-
-    /*
-      If this is a restart run the ThresholdPressure object will be active,
-      but it will *not* be properly initialized with numerical values. The
-      values must instead come from the THPRES vector in the restart file.
-    */
-    if (simConfig.getThresholdPressure().restart())
-        return;
-
-    // allocate the array which specifies the threshold pressures
-    thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
-    thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
-}
-
-template<class Grid, class GridView, class ElementMapper, class Scalar>
-void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-applyExplicitThresholdPressures_()
-{
-    const SimulationConfig& simConfig = eclState_.getSimulationConfig();
-    const auto& thpres = simConfig.getThresholdPressure();
-
-    // set the threshold pressures for all EQUIL region boundaries which have a
-    // intersection in the grid
-    for (const auto& elem : elements(gridView_, Dune::Partitions::interior)) {
-        for (const auto& intersection : intersections(gridView_, elem)) {
-            if (intersection.boundary())
-                continue; // ignore boundary intersections for now (TODO?)
-            else if (!intersection.neighbor()) //processor boundary but not domain boundary
-                continue;
-
-            const auto& inside = intersection.inside();
-            const auto& outside = intersection.outside();
-
-            unsigned insideElemIdx = elementMapper_.index(inside);
-            unsigned outsideElemIdx = elementMapper_.index(outside);
-
-            auto equilRegionInside = elemEquilRegion_[insideElemIdx];
-            auto equilRegionOutside = elemEquilRegion_[outsideElemIdx];
-            if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
-                Scalar pth = 0.0;
-                if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
-                    // threshold pressure explicitly specified
-                    pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
-                }
-                else {
-                    // take the threshold pressure from the initial condition
-                    unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
-                    pth = thpresDefault_[offset];
-                }
-
-                unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
-                unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
-
-                thpres_[offset1] = pth;
-                thpres_[offset2] = pth;
-            }
-        }
-    }
-
-    // apply threshold pressures across faults
-    if (thpres.ftSize() > 0)
-        configureThpresft_();
-}
-
-template<class Grid, class GridView, class ElementMapper, class Scalar>
-void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-configureThpresft_()
-{
-    // retrieve the faults collection.
-    const FaultCollection& faults = eclState_.getFaults();
-    const SimulationConfig& simConfig = eclState_.getSimulationConfig();
-    const auto& thpres = simConfig.getThresholdPressure();
-
-    // extract the multipliers
-    int numFaults = faults.size();
-    int numCartesianElem = eclState_.getInputGrid().getCartesianSize();
-    thpresftValues_.resize(numFaults, -1.0);
-    cartElemFaultIdx_.resize(numCartesianElem, -1);
-    for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) {
-        auto& fault = faults.getFault(faultIdx);
-        thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx);
-        for (const FaultFace& face : fault)
-            // "face" is a misnomer because the object describes a set of cell
-            // indices, but we go with the conventions of the parser here...
-            for (size_t cartElemIdx : face)
-                cartElemFaultIdx_[cartElemIdx] = faultIdx;
-    }
-}
-
-template<class Grid, class GridView, class ElementMapper, class Scalar>
-std::vector<Scalar>
-EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-getRestartVector() const
-{
-    if (!enableThresholdPressure_)
-        return {};
-
-    std::vector<Scalar> result(numEquilRegions_ * numEquilRegions_, 0.0);
-    const auto& simConfig = eclState_.getSimulationConfig();
-    const auto& thpres = simConfig.getThresholdPressure();
-
-    std::size_t idx = 0;
-    for (unsigned j = 1; j <= numEquilRegions_; ++j) {
-        for (unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) {
-            if (thpres.hasRegionBarrier(i, j)) {
-                if (thpres.hasThresholdPressure(i, j)) {
-                    result[idx] = thpres.getThresholdPressure(i, j);
-                } else {
-                    result[idx] = this->thpresDefault_[idx];
-                }
-            }
-        }
-    }
-
-    return result;
-}
-
-template<class Grid, class GridView, class ElementMapper, class Scalar>
-void
-EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
-logPressures()
-{
-    if (!enableThresholdPressure_)
-        return;
-
-    auto lineFormat = [this](unsigned i, unsigned j, double val)
-    {
-        const auto& units = eclState_.getUnits();
-        return fmt::format("{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n",
-                           " ", i,
-                           " ", j,
-                           " ", units.from_si(UnitSystem::measure::pressure, val),
-                           " ", units.name(UnitSystem::measure::pressure));
-    };
-    auto lineFormatS = [](auto s1, auto s2, auto s3)
-    {
-        return fmt::format("{:4}{:^16}{:13}{:^9}{:21}{:^18}\n",
-                           " ", s1, " ", s2, " ", s3);
-    };
-
-    std::string str = "\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n"
-                      "----------------------------------------\n"
-                      "\n";
-    str += lineFormatS("FLOW FROM REGION", "TO REGION", "THRESHOLD PRESSURE");
-    str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
-    const auto& simConfig = eclState_.getSimulationConfig();
-    const auto& thpres = simConfig.getThresholdPressure();
-
-    for (unsigned i = 1; i <= numEquilRegions_; ++i) {
-        for (unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) {
-            if (thpres.hasRegionBarrier(i, j)) {
-                if (thpres.hasThresholdPressure(i, j)) {
-                    str += lineFormat(i, j, thpres.getThresholdPressure(j, i));
-                } else {
-                    std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1);
-                    str += lineFormat(i, j, this->thpresDefault_[idx]);
-                }
-            }
-        }
-    }
-    str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
-    OpmLog::note(str);
-}
-
 #if HAVE_DUNE_FEM
 template class EclGenericThresholdPressure<Dune::CpGrid,
                                            Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
@@ -331,57 +54,12 @@ template class EclGenericThresholdPressure<Dune::CpGrid,
                                                         Dune::PartitionIteratorType(4),
                                                         false>>>,
                                             double>;
-#if HAVE_DUNE_ALUGRID
-#if HAVE_MPI
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
-#else
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
-#endif //HAVE_MPI
-
-template class EclGenericThresholdPressure<ALUGrid3CN,
-                                           Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,
-                                           Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
-                                           double>;
-template class EclGenericThresholdPressure<ALUGrid3CN,
-                                            Dune::Fem::GridPart2GridViewImpl<
-                                                Dune::Fem::AdaptiveLeafGridPart<
-                                                    ALUGrid3CN,
-                                                    Dune::PartitionIteratorType(4),
-                                                    false> >,
-                                            Dune::MultipleCodimMultipleGeomTypeMapper<
-                                                Dune::Fem::GridPart2GridViewImpl<
-                                                    Dune::Fem::AdaptiveLeafGridPart<
-                                                        ALUGrid3CN,
-                                                        Dune::PartitionIteratorType(4),
-                                                        false>>>,
-                                            double>;                                           
-                                           
-                                           
-
-#endif //HAVE_DUNE_ALUGRID
 #else
 template class EclGenericThresholdPressure<Dune::CpGrid,
                                            Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
                                            Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
                                            double>;
-#if HAVE_DUNE_ALUGRID
-
-#if HAVE_MPI
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
-#else
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
-#endif //HAVE_MPI
-template class EclGenericThresholdPressure<ALUGrid3CN,
-                                           Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>,
-                                           Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>>,
-                                           double>;
-#endif //HAVE_DUNE_ALUGRID
-
 #endif
 
-template class EclGenericThresholdPressure<Dune::PolyhedralGrid<3,3,double>,
-                                           Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
-                                           Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
-                                           double>;
 
 } // namespace Opm

ebos/eclgenericthresholdpressure_impl.hh

@@ -0,0 +1,307 @@
+// -*- 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.
+*/
+#ifndef EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_IMPL_HH
+#define EWOMS_ECL_GENERIC_THRESHOLD_PRESSURE_IMPL_HH
+
+#include <ebos/eclgenericthresholdpressure.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/common/rangegenerators.hh>
+
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/Eqldims.hpp>
+#include <opm/input/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
+#include <opm/input/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
+
+#include <fmt/format.h>
+
+#include <algorithm>
+#include <cassert>
+#include <limits>
+#include <stdexcept>
+
+namespace Opm {
+
+template<class Grid, class GridView, class ElementMapper, class Scalar>
+EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+EclGenericThresholdPressure(const CartesianIndexMapper& cartMapper,
+                            const GridView& gridView,
+                            const ElementMapper& elementMapper,
+                            const EclipseState& eclState)
+    : cartMapper_(cartMapper)
+    , gridView_(gridView)
+    , elementMapper_(elementMapper)
+    , eclState_(eclState)
+{
+}
+
+template<class Grid, class GridView, class ElementMapper,class Scalar>
+Scalar EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+thresholdPressure(int elem1Idx, int elem2Idx) const
+{
+    if (!enableThresholdPressure_)
+        return 0.0;
+
+    // threshold pressure accross faults
+    if (!thpresftValues_.empty()) {
+        int cartElem1Idx = cartMapper_.cartesianIndex(elem1Idx);
+        int cartElem2Idx = cartMapper_.cartesianIndex(elem2Idx);
+
+        assert(0 <= cartElem1Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem1Idx);
+        assert(0 <= cartElem2Idx && static_cast<int>(cartElemFaultIdx_.size()) > cartElem2Idx);
+
+        int fault1Idx = cartElemFaultIdx_[cartElem1Idx];
+        int fault2Idx = cartElemFaultIdx_[cartElem2Idx];
+        if (fault1Idx != -1 && fault1Idx == fault2Idx)
+            // inside a fault there's no threshold pressure, even accross EQUIL
+            // regions.
+            return 0.0;
+        if (fault1Idx != fault2Idx) {
+            // TODO: which value if a cell is part of multiple faults? we take
+            // the maximum here.
+            Scalar val1 = (fault1Idx >= 0) ? thpresftValues_[fault1Idx] : 0.0;
+            Scalar val2 = (fault2Idx >= 0) ? thpresftValues_[fault2Idx] : 0.0;
+            return std::max(val1, val2);
+        }
+    }
+
+    // threshold pressure accross EQUIL regions
+    auto equilRegion1Idx = elemEquilRegion_[elem1Idx];
+    auto equilRegion2Idx = elemEquilRegion_[elem2Idx];
+
+    if (equilRegion1Idx == equilRegion2Idx)
+        return 0.0;
+
+    return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
+}
+
+template<class Grid, class GridView, class ElementMapper, class Scalar>
+void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+finishInit()
+{
+    unsigned numElements = gridView_.size(/*codim=*/0);
+
+    const auto& simConfig = eclState_.getSimulationConfig();
+
+    enableThresholdPressure_ = simConfig.useThresholdPressure();
+    if (!enableThresholdPressure_)
+        return;
+
+    numEquilRegions_ = eclState_.getTableManager().getEqldims().getNumEquilRegions();
+    const decltype(numEquilRegions_) maxRegions =
+        std::numeric_limits<std::decay_t<decltype(elemEquilRegion_[0])>>::max();
+
+    if (numEquilRegions_ > maxRegions) {
+        // make sure that the index of an equilibration region can be stored
+        // in the vector
+        OPM_THROW(std::invalid_argument,
+                  (fmt::format("The maximum number of supported "
+                               "equilibration regions by OPM flow is {}, but "
+                               "{} are used!",
+                               maxRegions, numEquilRegions_)));
+    }
+
+    if (numEquilRegions_ > 2048) {
+        // warn about performance
+        OpmLog::warning(fmt::format("Number of equilibration regions is {}, which is "
+                                 "rather large. Note, that this might "
+                                 "have a negative impact on performance "
+                                 "and memory consumption.", numEquilRegions_));
+    }
+
+    // internalize the data specified using the EQLNUM keyword
+    const auto& fp = eclState_.fieldProps();
+    const auto& equilRegionData = fp.get_int("EQLNUM");
+    elemEquilRegion_.resize(numElements, 0);
+    for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
+        elemEquilRegion_[elemIdx] = equilRegionData[elemIdx] - 1;
+    }
+
+    /*
+      If this is a restart run the ThresholdPressure object will be active,
+      but it will *not* be properly initialized with numerical values. The
+      values must instead come from the THPRES vector in the restart file.
+    */
+    if (simConfig.getThresholdPressure().restart())
+        return;
+
+    // allocate the array which specifies the threshold pressures
+    thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
+    thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
+}
+
+template<class Grid, class GridView, class ElementMapper, class Scalar>
+void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+applyExplicitThresholdPressures_()
+{
+    const SimulationConfig& simConfig = eclState_.getSimulationConfig();
+    const auto& thpres = simConfig.getThresholdPressure();
+
+    // set the threshold pressures for all EQUIL region boundaries which have a
+    // intersection in the grid
+    for (const auto& elem : elements(gridView_, Dune::Partitions::interior)) {
+        for (const auto& intersection : intersections(gridView_, elem)) {
+            if (intersection.boundary())
+                continue; // ignore boundary intersections for now (TODO?)
+            else if (!intersection.neighbor()) //processor boundary but not domain boundary
+                continue;
+
+            const auto& inside = intersection.inside();
+            const auto& outside = intersection.outside();
+
+            unsigned insideElemIdx = elementMapper_.index(inside);
+            unsigned outsideElemIdx = elementMapper_.index(outside);
+
+            auto equilRegionInside = elemEquilRegion_[insideElemIdx];
+            auto equilRegionOutside = elemEquilRegion_[outsideElemIdx];
+            if (thpres.hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1)) {
+                Scalar pth = 0.0;
+                if (thpres.hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
+                    // threshold pressure explicitly specified
+                    pth = thpres.getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
+                }
+                else {
+                    // take the threshold pressure from the initial condition
+                    unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
+                    pth = thpresDefault_[offset];
+                }
+
+                unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
+                unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
+
+                thpres_[offset1] = pth;
+                thpres_[offset2] = pth;
+            }
+        }
+    }
+
+    // apply threshold pressures across faults
+    if (thpres.ftSize() > 0)
+        configureThpresft_();
+}
+
+template<class Grid, class GridView, class ElementMapper, class Scalar>
+void EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+configureThpresft_()
+{
+    // retrieve the faults collection.
+    const FaultCollection& faults = eclState_.getFaults();
+    const SimulationConfig& simConfig = eclState_.getSimulationConfig();
+    const auto& thpres = simConfig.getThresholdPressure();
+
+    // extract the multipliers
+    int numFaults = faults.size();
+    int numCartesianElem = eclState_.getInputGrid().getCartesianSize();
+    thpresftValues_.resize(numFaults, -1.0);
+    cartElemFaultIdx_.resize(numCartesianElem, -1);
+    for (size_t faultIdx = 0; faultIdx < faults.size(); faultIdx++) {
+        auto& fault = faults.getFault(faultIdx);
+        thpresftValues_[faultIdx] = thpres.getThresholdPressureFault(faultIdx);
+        for (const FaultFace& face : fault)
+            // "face" is a misnomer because the object describes a set of cell
+            // indices, but we go with the conventions of the parser here...
+            for (size_t cartElemIdx : face)
+                cartElemFaultIdx_[cartElemIdx] = faultIdx;
+    }
+}
+
+template<class Grid, class GridView, class ElementMapper, class Scalar>
+std::vector<Scalar>
+EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+getRestartVector() const
+{
+    if (!enableThresholdPressure_)
+        return {};
+
+    std::vector<Scalar> result(numEquilRegions_ * numEquilRegions_, 0.0);
+    const auto& simConfig = eclState_.getSimulationConfig();
+    const auto& thpres = simConfig.getThresholdPressure();
+
+    std::size_t idx = 0;
+    for (unsigned j = 1; j <= numEquilRegions_; ++j) {
+        for (unsigned i = 1; i <= numEquilRegions_; ++i, ++idx) {
+            if (thpres.hasRegionBarrier(i, j)) {
+                if (thpres.hasThresholdPressure(i, j)) {
+                    result[idx] = thpres.getThresholdPressure(i, j);
+                } else {
+                    result[idx] = this->thpresDefault_[idx];
+                }
+            }
+        }
+    }
+
+    return result;
+}
+
+template<class Grid, class GridView, class ElementMapper, class Scalar>
+void
+EclGenericThresholdPressure<Grid,GridView,ElementMapper,Scalar>::
+logPressures()
+{
+    if (!enableThresholdPressure_)
+        return;
+
+    auto lineFormat = [this](unsigned i, unsigned j, double val)
+    {
+        const auto& units = eclState_.getUnits();
+        return fmt::format("{:4}{:>6}{:23}{:>6}{:24}{:>11.07}{:7}{}\n",
+                           " ", i,
+                           " ", j,
+                           " ", units.from_si(UnitSystem::measure::pressure, val),
+                           " ", units.name(UnitSystem::measure::pressure));
+    };
+    auto lineFormatS = [](auto s1, auto s2, auto s3)
+    {
+        return fmt::format("{:4}{:^16}{:13}{:^9}{:21}{:^18}\n",
+                           " ", s1, " ", s2, " ", s3);
+    };
+
+    std::string str = "\nLIST OF ALL NON-ZERO THRESHOLD PRESSURES\n"
+                      "----------------------------------------\n"
+                      "\n";
+    str += lineFormatS("FLOW FROM REGION", "TO REGION", "THRESHOLD PRESSURE");
+    str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
+    const auto& simConfig = eclState_.getSimulationConfig();
+    const auto& thpres = simConfig.getThresholdPressure();
+
+    for (unsigned i = 1; i <= numEquilRegions_; ++i) {
+        for (unsigned j = (thpres.irreversible() ? 1 : i); j <= numEquilRegions_; ++j) {
+            if (thpres.hasRegionBarrier(i, j)) {
+                if (thpres.hasThresholdPressure(i, j)) {
+                    str += lineFormat(i, j, thpres.getThresholdPressure(j, i));
+                } else {
+                    std::size_t idx = (j - 1) * numEquilRegions_ + (i - 1);
+                    str += lineFormat(i, j, this->thpresDefault_[idx]);
+                }
+            }
+        }
+    }
+    str += lineFormatS(std::string(16, '-'), std::string(9, '-'), std::string(18, '-'));
+    OpmLog::note(str);
+}
+
+} // namespace Opm
+#endif

ebos/eclgenerictracermodel.cc

@@ -23,7 +23,14 @@
 #include <config.h>
 #include "eclgenerictracermodel_impl.hh"
 
+#if HAVE_DUNE_FEM
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+#include <dune/fem/gridpart/common/gridpart2gridview.hh>
+#include <ebos/femcpgridcompat.hh>
+#endif // HAVE_DUNE_FEM
+
 namespace Opm {
+
 #if HAVE_DUNE_FEM
 template class EclGenericTracerModel<Dune::CpGrid,
                                      Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
@@ -39,54 +46,12 @@ template class EclGenericTracerModel<Dune::CpGrid,
                                                                   Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false> >,
                                                       false, false>,
                                      double>;
-#if HAVE_DUNE_ALUGRID
-
-#if HAVE_MPI
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
-#else
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
-#endif //HAVE_MPI
-                                    
-template class EclGenericTracerModel<ALUGrid3CN, Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>, Opm::EcfvStencil<double,Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>,false,false>,
-                                     double>;
-
-template class EclGenericTracerModel<ALUGrid3CN,
-Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> >,
-                                     Dune::MultipleCodimMultipleGeomTypeMapper<
-                                         Dune::Fem::GridPart2GridViewImpl<
-Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> > >,
-                                     Opm::EcfvStencil<double, Dune::Fem::GridPart2GridViewImpl<
-                                                                  Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false> >,
-                                                      false, false>,
-                                     double>;                                     
-#endif //HAVE_DUNE_ALUGRID                                     
 #else
 template class EclGenericTracerModel<Dune::CpGrid,
                                      Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
                                      Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
                                      Opm::EcfvStencil<double,Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,false,false>,
                                      double>;
-#if HAVE_DUNE_ALUGRID
-#if HAVE_MPI
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
-#else
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
-#endif //HAVE_MPI
-
-template class EclGenericTracerModel<ALUGrid3CN,
-                                     Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>,
-                                     Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,
-                                         Dune::PartitionIteratorType(4)>>>,
-                                     Opm::EcfvStencil<double,Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,
-                                        Dune::PartitionIteratorType(4)>>,false,false>,
-                                     double>;
-#endif //HAVE_DUNE_ALUGRID
 #endif //HAVE_DUNE_FEM
 
-template class EclGenericTracerModel<Dune::PolyhedralGrid<3,3,double>,
-                                     Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,
-                                     Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
-                                     Opm::EcfvStencil<double, Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>,false,false>,
-                                     double>;
-
 } // namespace Opm

ebos/eclgenerictracermodel_impl.hh

@@ -39,7 +39,6 @@
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/grid/CpGrid.hpp>
-#include <opm/grid/polyhedralgrid.hh>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/input/eclipse/EclipseState/Phase.hpp>
@@ -53,18 +52,6 @@
 #include <opm/simulators/linalg/PropertyTree.hpp>
 #include <opm/simulators/linalg/FlexibleSolver.hpp>
 
-#if HAVE_DUNE_FEM
-#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
-#include <dune/fem/gridpart/common/gridpart2gridview.hh>
-#include <ebos/femcpgridcompat.hh>
-#endif // HAVE_DUNE_FEM
-
-#if HAVE_DUNE_ALUGRID
-#include <dune/alugrid/grid.hh>
-#include <dune/alugrid/3d/gridview.hh>
-#include "alucartesianindexmapper.hh"
-#endif // HAVE_DUNE_ALUGRID
-
 #include <fmt/format.h>
 
 #include <array>

ebos/eclgenericwriter.cc

@@ -23,40 +23,12 @@
 
 #include <config.h>
 
-#include <ebos/eclgenericwriter.hh>
+// need to include these before impl because of specializations.
-
 #include <opm/grid/CpGrid.hpp>
 #include <opm/grid/cpgrid/GridHelpers.hpp>
-#include <opm/grid/polyhedralgrid.hh>
-#include <opm/grid/utility/cartesianToCompressed.hpp>
 
-#if HAVE_DUNE_ALUGRID
+#include <ebos/eclgenericwriter.hh>
-#include "eclalugridvanguard.hh"
+#include <ebos/eclgenericwriter_impl.hh>
-#include <dune/alugrid/grid.hh>
-#include <dune/alugrid/3d/gridview.hh>
-#endif // HAVE_DUNE_ALUGRID
-
-#include <opm/output/eclipse/EclipseIO.hpp>
-#include <opm/output/eclipse/RestartValue.hpp>
-#include <opm/output/eclipse/Summary.hpp>
-
-#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
-
-#include <opm/input/eclipse/Schedule/Action/State.hpp>
-#include <opm/input/eclipse/Schedule/Schedule.hpp>
-#include <opm/input/eclipse/Schedule/SummaryState.hpp>
-#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
-#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
-#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
-
-#include <opm/input/eclipse/Units/UnitSystem.hpp>
-
-#include <dune/grid/common/mcmgmapper.hh>
-
-#if HAVE_MPI
-#include <ebos/eclmpiserializer.hh>
-#endif
 
 #if HAVE_DUNE_FEM
 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -64,543 +36,8 @@
 #include <ebos/femcpgridcompat.hh>
 #endif // HAVE_DUNE_FEM
 
-#if HAVE_MPI
-#include <mpi.h>
-#endif
-
-#include <algorithm>
-#include <array>
-#include <cassert>
-#include <functional>
-#include <map>
-#include <memory>
-#include <string>
-#include <unordered_map>
-#include <utility>
-#include <vector>
-
-namespace {
-
-/*!
- * \brief Detect whether two cells are direct vertical neighbours.
- *
- * I.e. have the same i and j index and all cartesian cells between them
- * along the vertical column are inactive.
- *
- * \tparam CM The type of the cartesian index mapper.
- * \param cartMapper The mapper onto cartesian indices.
- * \param cartesianToActive The mapping of cartesian indices to active indices.
- * \param smallGlobalIndex The cartesian cell index of the cell with smaller index
- * \param largeGlobalIndex The cartesian cell index of the cell with larger index
- * \return True if the cells have the same i and j indices and all cartesian cells
- *         between them are inactive.
- */
-bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
-                             const std::unordered_map<int,int>& cartesianToActive,
-                             int smallGlobalIndex, int largeGlobalIndex)
-{
-    assert(smallGlobalIndex <= largeGlobalIndex);
-    std::array<int, 3> ijk1, ijk2;
-    auto globalToIjk = [cartDims](int gc) {
-                           std::array<int, 3> ijk;
-                           ijk[0] = gc % cartDims[0];
-                           gc /= cartDims[0];
-                           ijk[1] = gc % cartDims[1];
-                           ijk[2] = gc / cartDims[1];
-                           return ijk;
-                       };
-    ijk1 = globalToIjk(smallGlobalIndex);
-    ijk2 = globalToIjk(largeGlobalIndex);
-    assert(ijk2[2]>=ijk1[2]);
-
-    if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
-    {
-        assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
-        for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
-              gi += cartDims[0] * cartDims[1] )
-        {
-            if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
-            {
-                return false;
-            }
-        }
-        return true;
-    } else
-        return false;
-}
-
-std::unordered_map<std::string, Opm::data::InterRegFlowMap>
-getInterRegFlowsAsMap(const Opm::EclInterRegFlowMap& map)
-{
-    auto maps = std::unordered_map<std::string, Opm::data::InterRegFlowMap>{};
-
-    const auto& regionNames = map.names();
-    auto flows = map.getInterRegFlows();
-    const auto nmap = regionNames.size();
-
-    maps.reserve(nmap);
-    for (auto mapID = 0*nmap; mapID < nmap; ++mapID) {
-        maps.emplace(regionNames[mapID], std::move(flows[mapID]));
-    }
-
-    return maps;
-}
-
-struct EclWriteTasklet : public Opm::TaskletInterface
-{
-    Opm::Action::State actionState_;
-    Opm::WellTestState wtestState_;
-    Opm::SummaryState summaryState_;
-    Opm::UDQState udqState_;
-    Opm::EclipseIO& eclIO_;
-    int reportStepNum_;
-    bool isSubStep_;
-    double secondsElapsed_;
-    Opm::RestartValue restartValue_;
-    bool writeDoublePrecision_;
-
-    explicit EclWriteTasklet(const Opm::Action::State& actionState,
-                             const Opm::WellTestState& wtestState,
-                             const Opm::SummaryState& summaryState,
-                             const Opm::UDQState& udqState,
-                             Opm::EclipseIO& eclIO,
-                             int reportStepNum,
-                             bool isSubStep,
-                             double secondsElapsed,
-                             Opm::RestartValue restartValue,
-                             bool writeDoublePrecision)
-        : actionState_(actionState)
-        , wtestState_(wtestState)
-        , summaryState_(summaryState)
-        , udqState_(udqState)
-        , eclIO_(eclIO)
-        , reportStepNum_(reportStepNum)
-        , isSubStep_(isSubStep)
-        , secondsElapsed_(secondsElapsed)
-        , restartValue_(std::move(restartValue))
-        , writeDoublePrecision_(writeDoublePrecision)
-    {}
-
-    // callback to eclIO serial writeTimeStep method
-    void run()
-    {
-        this->eclIO_.writeTimeStep(this->actionState_,
-                                   this->wtestState_,
-                                   this->summaryState_,
-                                   this->udqState_,
-                                   this->reportStepNum_,
-                                   this->isSubStep_,
-                                   this->secondsElapsed_,
-                                   std::move(this->restartValue_),
-                                   this->writeDoublePrecision_);
-    }
-};
-
-}
-
 namespace Opm {
 
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-EclGenericWriter(const Schedule& schedule,
-                 const EclipseState& eclState,
-                 const SummaryConfig& summaryConfig,
-                 const Grid& grid,
-                 const EquilGrid* equilGrid,
-                 const GridView& gridView,
-                 const Dune::CartesianIndexMapper<Grid>& cartMapper,
-                 const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper,
-                 bool enableAsyncOutput,
-                 bool enableEsmry )
-    : collectToIORank_(grid,
-                       equilGrid,
-                       gridView,
-                       cartMapper,
-                       equilCartMapper,
-                       summaryConfig.fip_regions_interreg_flow())
-    , grid_           (grid)
-    , gridView_       (gridView)
-    , schedule_       (schedule)
-    , eclState_       (eclState)
-    , summaryConfig_  (summaryConfig)
-    , cartMapper_     (cartMapper)
-    , equilCartMapper_(equilCartMapper)
-    , equilGrid_      (equilGrid)
-{
-    if (this->collectToIORank_.isIORank()) {
-        this->eclIO_ = std::make_unique<EclipseIO>
-            (this->eclState_,
-             UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
-             this->schedule_, this->summaryConfig_, "", enableEsmry);
-    }
-
-    // create output thread if enabled and rank is I/O rank
-    // async output is enabled by default if pthread are enabled
-    int numWorkerThreads = 0;
-    if (enableAsyncOutput && collectToIORank_.isIORank()) {
-        numWorkerThreads = 1;
-    }
-
-    this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-const EclipseIO& EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-eclIO() const
-{
-    assert(eclIO_);
-    return *eclIO_;
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-writeInit(const std::function<unsigned int(unsigned int)>& map)
-{
-    if (collectToIORank_.isIORank()) {
-        std::map<std::string, std::vector<int>> integerVectors;
-        if (collectToIORank_.isParallel()) {
-            integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
-        }
-
-        auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
-
-        eclIO_->writeInitial(computeTrans_(cartMap, map),
-                             integerVectors,
-                             exportNncStructure_(cartMap, map));
-    }
-
-#if HAVE_MPI
-    if (collectToIORank_.isParallel()) {
-        const auto& comm = grid_.comm();
-        Opm::EclMpiSerializer ser(comm);
-        ser.broadcast(outputNnc_);
-    }
-#endif
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-computeTrans_(const std::unordered_map<int,int>& cartesianToActive,
-              const std::function<unsigned int(unsigned int)>& map) const
-{
-    const auto& cartMapper = *equilCartMapper_;
-    const auto& cartDims = cartMapper.cartesianDimensions();
-    const int globalSize = cartDims[0] * cartDims[1] * cartDims[2];
-
-    auto tranx = data::CellData {
-        UnitSystem::measure::transmissibility,
-        std::vector<double>(globalSize, 0.0),
-        data::TargetType::INIT
-    };
-
-    auto trany = tranx;
-    auto tranz = tranx;
-
-    using GlobalGridView = typename EquilGrid::LeafGridView;
-    const GlobalGridView& globalGridView = equilGrid_->leafGridView();
-    using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
-    GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
-
-    for (const auto& elem : elements(globalGridView)) {
-        for (const auto& is : intersections(globalGridView, elem)) {
-            if (!is.neighbor())
-                continue; // intersection is on the domain boundary
-
-            unsigned c1 = globalElemMapper.index(is.inside());
-            unsigned c2 = globalElemMapper.index(is.outside());
-
-            if (c1 > c2)
-                continue; // we only need to handle each connection once, thank you.
-
-            // Ordering of compressed and uncompressed index should be the same
-            const int cartIdx1 = cartMapper.cartesianIndex( c1 );
-            const int cartIdx2 = cartMapper.cartesianIndex( c2 );
-            // Ordering of compressed and uncompressed index should be the same
-            assert(cartIdx1 <= cartIdx2);
-            int gc1 = std::min(cartIdx1, cartIdx2);
-            int gc2 = std::max(cartIdx1, cartIdx2);
-
-            // Re-ordering in case of non-empty mapping between equilGrid to grid
-            if (map) {
-                c1 = map(c1); // equilGridToGrid map
-                c2 = map(c2);
-            }
-
-            if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
-                tranx.data[gc1] = globalTrans().transmissibility(c1, c2);
-                continue; // skip other if clauses as they are false, last one needs some computation
-            }
-
-            if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
-                trany.data[gc1] = globalTrans().transmissibility(c1, c2);
-                continue; // skipt next if clause as it needs some computation
-            }
-
-            if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
-                 directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
-                tranz.data[gc1] = globalTrans().transmissibility(c1, c2);
-        }
-    }
-
-    return {
-        {"TRANX", tranx},
-        {"TRANY", trany},
-        {"TRANZ", tranz},
-    };
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
-{
-    std::size_t nx = eclState_.getInputGrid().getNX();
-    std::size_t ny = eclState_.getInputGrid().getNY();
-    auto nncData = eclState_.getInputNNC().input();
-    const auto& unitSystem = eclState_.getDeckUnitSystem();
-
-    for( const auto& entry : nncData ) {
-        // test whether NNC is not a neighboring connection
-        // cell2>=cell1 holds due to sortNncAndApplyEditnnc
-        assert( entry.cell2 >= entry.cell1 );
-        auto cellDiff = entry.cell2 - entry.cell1;
-
-        if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) {
-            auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans);
-            // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6
-            if ( tt >= 1.0e-6 )
-                outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans);
-        }
-    }
-
-    using GlobalGridView = typename EquilGrid::LeafGridView;
-    const GlobalGridView& globalGridView = equilGrid_->leafGridView();
-    using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
-    GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
-
-    // Cartesian index mapper for the serial I/O grid
-    const auto& equilCartMapper =  *equilCartMapper_;
-    const auto& cartDims = cartMapper_.cartesianDimensions();
-    for (const auto& elem : elements(globalGridView)) {
-        for (const auto& is : intersections(globalGridView, elem)) {
-            if (!is.neighbor())
-                continue; // intersection is on the domain boundary
-
-            unsigned c1 = globalElemMapper.index(is.inside());
-            unsigned c2 = globalElemMapper.index(is.outside());
-
-            if (c1 > c2)
-                continue; // we only need to handle each connection once, thank you.
-
-            std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
-            std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
-
-            if ( cc2 < cc1 )
-                std::swap(cc1, cc2);
-
-            auto cellDiff = cc2 - cc1;
-
-            // Re-ordering in case of non-empty mapping between equilGrid to grid
-            if (map) {
-                c1 = map(c1); // equilGridToGrid map
-                c2 = map(c2);
-            }
-
-            if (cellDiff != 1 &&
-                cellDiff != nx &&
-                cellDiff != nx*ny &&
-                !directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) {
-                // We need to check whether an NNC for this face was also specified
-                // via the NNC keyword in the deck (i.e. in the first origNncSize entries.
-                auto t = globalTrans().transmissibility(c1, c2);
-                auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0));
-
-                while ( candidate != nncData.end() && candidate->cell1 == cc1
-                     && candidate->cell2 == cc2) {
-                    t -= candidate->trans;
-                    ++candidate;
-                }
-                // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC
-                // with corresponding EDITNNC above). In addition we do set small transmissibilties
-                // to zero when setting up the simulator. These will be ignored here, too.
-                auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t));
-                if ( tt > 1e-12 )
-                    outputNnc_.push_back({cc1, cc2, t});
-            }
-        }
-    }
-    return outputNnc_;
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-doWriteOutput(const int                     reportStepNum,
-              const bool                    isSubStep,
-              data::Solution&&              localCellData,
-              data::Wells&&                 localWellData,
-              data::GroupAndNetworkValues&& localGroupAndNetworkData,
-              data::Aquifers&&              localAquiferData,
-              WellTestState&&               localWTestState,
-              const Action::State& actionState,
-              const UDQState& udqState,
-              const SummaryState& summaryState,
-              const std::vector<Scalar>& thresholdPressure,
-              Scalar curTime,
-              Scalar nextStepSize,
-              bool doublePrecision,
-              bool isFlowsn,
-              std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& flowsn,
-              bool isFloresn,
-              std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& floresn)
-{
-    const auto isParallel = this->collectToIORank_.isParallel();
-    const bool needsReordering = this->collectToIORank_.doesNeedReordering();
-
-    RestartValue restartValue {
-        (isParallel || needsReordering)
-        ? this->collectToIORank_.globalCellData()
-        : std::move(localCellData),
-
-        isParallel ? this->collectToIORank_.globalWellData()
-                   : std::move(localWellData),
-
-        isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
-                   : std::move(localGroupAndNetworkData),
-
-        isParallel ? this->collectToIORank_.globalAquiferData()
-                   : std::move(localAquiferData)
-    };
-
-    if (eclState_.getSimulationConfig().useThresholdPressure()) {
-        restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure,
-                              thresholdPressure);
-    }
-
-    // Add suggested next timestep to extra data.
-    if (! isSubStep) {
-        restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
-    }
-
-    // Add nnc flows and flores.
-    if (isFlowsn) {
-        const auto flowsn_global = isParallel ? this->collectToIORank_.globalFlowsn() : std::move(flowsn);
-        for (const auto& flows : flowsn_global) { 
-            if (flows.first.empty())
-                continue;
-            if (flows.first == "FLOGASN+") {
-                restartValue.addExtra(flows.first, UnitSystem::measure::gas_surface_rate, flows.second.second);
-            }
-            else {
-                restartValue.addExtra(flows.first, UnitSystem::measure::liquid_surface_rate, flows.second.second);
-            }
-        }
-    }
-    if (isFloresn) {
-        const auto floresn_global = isParallel ? this->collectToIORank_.globalFloresn() : std::move(floresn);
-        for (const auto& flores : floresn_global) {
-            if (flores.first.empty())
-                continue;
-            restartValue.addExtra(flores.first, UnitSystem::measure::rate, flores.second.second);
-        }
-    }
-
-    // first, create a tasklet to write the data for the current time
-    // step to disk
-    auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
-        actionState,
-        isParallel ? this->collectToIORank_.globalWellTestState() : std::move(localWTestState),
-        summaryState, udqState, *this->eclIO_,
-        reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision);
-
-    // then, make sure that the previous I/O request has been completed
-    // and the number of incomplete tasklets does not increase between
-    // time steps
-    this->taskletRunner_->barrier();
-
-    // finally, start a new output writing job
-    this->taskletRunner_->dispatch(std::move(eclWriteTasklet));
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-evalSummary(const int                                            reportStepNum,
-            const Scalar                                         curTime,
-            const data::Wells&                                   localWellData,
-            const data::WellBlockAveragePressures&               localWBPData,
-            const data::GroupAndNetworkValues&                   localGroupAndNetworkData,
-            const std::map<int,data::AquiferData>&               localAquiferData,
-            const std::map<std::pair<std::string, int>, double>& blockData,
-            const std::map<std::string, double>&                 miscSummaryData,
-            const std::map<std::string, std::vector<double>>&    regionData,
-            const Inplace&                                       inplace,
-            const Inplace&                                       initialInPlace,
-            const EclInterRegFlowMap&                            interRegFlows,
-            SummaryState&                                        summaryState,
-            UDQState&                                            udqState)
-{
-    if (collectToIORank_.isIORank()) {
-        const auto& summary = eclIO_->summary();
-
-        const auto& wellData = this->collectToIORank_.isParallel()
-            ? this->collectToIORank_.globalWellData()
-            : localWellData;
-
-        const auto& wbpData = this->collectToIORank_.isParallel()
-            ? this->collectToIORank_.globalWBPData()
-            : localWBPData;
-
-        const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
-            ? this->collectToIORank_.globalGroupAndNetworkData()
-            : localGroupAndNetworkData;
-
-        const auto& aquiferData = this->collectToIORank_.isParallel()
-            ? this->collectToIORank_.globalAquiferData()
-            : localAquiferData;
-
-        summary.eval(summaryState,
-                     reportStepNum,
-                     curTime,
-                     wellData,
-                     wbpData,
-                     groupAndNetworkData,
-                     miscSummaryData,
-                     initialInPlace,
-                     inplace,
-                     regionData,
-                     blockData,
-                     aquiferData,
-                     getInterRegFlowsAsMap(interRegFlows));
-
-        // Off-by-one-fun: The reportStepNum argument corresponds to the
-        // report step these results will be written to, whereas the
-        // argument to UDQ function evaluation corresponds to the report
-        // step we are currently on.
-        const auto udq_step = reportStepNum - 1;
-
-        this->schedule_.getUDQConfig(udq_step)
-            .eval(udq_step,
-                  this->schedule_,
-                  this->schedule_.wellMatcher(udq_step),
-                  summaryState,
-                  udqState);
-    }
-
-#if HAVE_MPI
-    if (collectToIORank_.isParallel()) {
-        EclMpiSerializer ser(grid_.comm());
-        ser.append(summaryState);
-    }
-#endif
-}
-
-template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
-const typename EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::TransmissibilityType&
-EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
-globalTrans() const
-{
-    assert (globalTrans_);
-    return *globalTrans_;
-}
-
 #if HAVE_DUNE_FEM
 template class EclGenericWriter<Dune::CpGrid,
                                 Dune::CpGrid,
@@ -620,59 +57,12 @@ template class EclGenericWriter<Dune::CpGrid,
                                             Dune::PartitionIteratorType(4),
                                             false>>>,
                                 double>;
-
-
-#ifdef HAVE_DUNE_ALUGRID
-#if HAVE_MPI
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
-#else
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
-#endif //HAVE_MPI
-                                                               
-template class EclGenericWriter<ALUGrid3CN,
-                                Dune::CpGrid,
-                                Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>, Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<ALUGrid3CN, Dune::PartitionIteratorType(4), false>>>>,
-                                double>;
-                                
-template class EclGenericWriter<ALUGrid3CN,
-                                Dune::CpGrid,
-                                Dune::Fem::GridPart2GridViewImpl<
-                                    Dune::Fem::AdaptiveLeafGridPart<
-                                        ALUGrid3CN,
-                                        Dune::PartitionIteratorType(4),
-                                        false>>,
-                                Dune::MultipleCodimMultipleGeomTypeMapper<
-                                    Dune::Fem::GridPart2GridViewImpl<
-                                        Dune::Fem::AdaptiveLeafGridPart<
-                                            ALUGrid3CN,
-                                            Dune::PartitionIteratorType(4),
-                                            false>>>,
-                                double>;                                
-#endif // HAVE_DUNE_ALUGRID
-
 #else // !HAVE_DUNE_FEM
 template class EclGenericWriter<Dune::CpGrid,
                                 Dune::CpGrid,
                                 Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
                                 Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
                                 double>;
-#if HAVE_DUNE_ALUGRID
+#endif // HAVE_DUNE_FEM
-#if HAVE_MPI
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
-#else
-    using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
-#endif //HAVE_MPI
-template class EclGenericWriter<ALUGrid3CN,
-                                Dune::CpGrid,
-                                Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN,Dune::PartitionIteratorType(4)>>,  Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>>,
-                                double>;
 
-#endif // HAVE_DUNE_ALUGRID
-#endif // !HAVE_DUNE_FEM
-
-template class EclGenericWriter<Dune::PolyhedralGrid<3,3,double>,
-                                Dune::PolyhedralGrid<3,3,double>,
-                                Dune::GridView<Dune::PolyhedralGridViewTraits<3, 3, double, Dune::PartitionIteratorType(4)>>,                              Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::PolyhedralGridViewTraits<3,3,double,Dune::PartitionIteratorType(4)>>>,
-                                double>;
-                                                                  
 } // namespace Opm

ebos/eclgenericwriter_impl.hh

@@ -0,0 +1,590 @@
+// -*- 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.
+*/
+#ifndef EWOMS_ECL_GENERIC_WRITER_IMPL_HH
+#define EWOMS_ECL_GENERIC_WRITER_IMPL_HH
+
+#include <ebos/eclgenericwriter.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <opm/grid/GridHelpers.hpp>
+#include <opm/grid/utility/cartesianToCompressed.hpp>
+
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
+#include <opm/input/eclipse/Schedule/Action/State.hpp>
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/input/eclipse/Schedule/SummaryState.hpp>
+#include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
+#include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
+#include <opm/input/eclipse/Schedule/Well/WellMatcher.hpp>
+#include <opm/input/eclipse/Units/UnitSystem.hpp>
+
+#include <opm/output/eclipse/EclipseIO.hpp>
+#include <opm/output/eclipse/RestartValue.hpp>
+#include <opm/output/eclipse/Summary.hpp>
+
+#if HAVE_MPI
+#include <ebos/eclmpiserializer.hh>
+#endif
+
+#if HAVE_MPI
+#include <mpi.h>
+#endif
+
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <functional>
+#include <map>
+#include <memory>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace {
+
+/*!
+ * \brief Detect whether two cells are direct vertical neighbours.
+ *
+ * I.e. have the same i and j index and all cartesian cells between them
+ * along the vertical column are inactive.
+ *
+ * \tparam CM The type of the cartesian index mapper.
+ * \param cartMapper The mapper onto cartesian indices.
+ * \param cartesianToActive The mapping of cartesian indices to active indices.
+ * \param smallGlobalIndex The cartesian cell index of the cell with smaller index
+ * \param largeGlobalIndex The cartesian cell index of the cell with larger index
+ * \return True if the cells have the same i and j indices and all cartesian cells
+ *         between them are inactive.
+ */
+bool directVerticalNeighbors(const std::array<int, 3>& cartDims,
+                             const std::unordered_map<int,int>& cartesianToActive,
+                             int smallGlobalIndex, int largeGlobalIndex)
+{
+    assert(smallGlobalIndex <= largeGlobalIndex);
+    std::array<int, 3> ijk1, ijk2;
+    auto globalToIjk = [cartDims](int gc) {
+                           std::array<int, 3> ijk;
+                           ijk[0] = gc % cartDims[0];
+                           gc /= cartDims[0];
+                           ijk[1] = gc % cartDims[1];
+                           ijk[2] = gc / cartDims[1];
+                           return ijk;
+                       };
+    ijk1 = globalToIjk(smallGlobalIndex);
+    ijk2 = globalToIjk(largeGlobalIndex);
+    assert(ijk2[2]>=ijk1[2]);
+
+    if ( ijk1[0] == ijk2[0] && ijk1[1] == ijk2[1] && (ijk2[2] - ijk1[2]) > 1)
+    {
+        assert((largeGlobalIndex-smallGlobalIndex)%(cartDims[0]*cartDims[1])==0);
+        for ( int gi = smallGlobalIndex + cartDims[0] * cartDims[1]; gi < largeGlobalIndex;
+              gi += cartDims[0] * cartDims[1] )
+        {
+            if ( cartesianToActive.find( gi ) != cartesianToActive.end() )
+            {
+                return false;
+            }
+        }
+        return true;
+    } else
+        return false;
+}
+
+std::unordered_map<std::string, Opm::data::InterRegFlowMap>
+getInterRegFlowsAsMap(const Opm::EclInterRegFlowMap& map)
+{
+    auto maps = std::unordered_map<std::string, Opm::data::InterRegFlowMap>{};
+
+    const auto& regionNames = map.names();
+    auto flows = map.getInterRegFlows();
+    const auto nmap = regionNames.size();
+
+    maps.reserve(nmap);
+    for (auto mapID = 0*nmap; mapID < nmap; ++mapID) {
+        maps.emplace(regionNames[mapID], std::move(flows[mapID]));
+    }
+
+    return maps;
+}
+
+struct EclWriteTasklet : public Opm::TaskletInterface
+{
+    Opm::Action::State actionState_;
+    Opm::WellTestState wtestState_;
+    Opm::SummaryState summaryState_;
+    Opm::UDQState udqState_;
+    Opm::EclipseIO& eclIO_;
+    int reportStepNum_;
+    bool isSubStep_;
+    double secondsElapsed_;
+    Opm::RestartValue restartValue_;
+    bool writeDoublePrecision_;
+
+    explicit EclWriteTasklet(const Opm::Action::State& actionState,
+                             const Opm::WellTestState& wtestState,
+                             const Opm::SummaryState& summaryState,
+                             const Opm::UDQState& udqState,
+                             Opm::EclipseIO& eclIO,
+                             int reportStepNum,
+                             bool isSubStep,
+                             double secondsElapsed,
+                             Opm::RestartValue restartValue,
+                             bool writeDoublePrecision)
+        : actionState_(actionState)
+        , wtestState_(wtestState)
+        , summaryState_(summaryState)
+        , udqState_(udqState)
+        , eclIO_(eclIO)
+        , reportStepNum_(reportStepNum)
+        , isSubStep_(isSubStep)
+        , secondsElapsed_(secondsElapsed)
+        , restartValue_(std::move(restartValue))
+        , writeDoublePrecision_(writeDoublePrecision)
+    {}
+
+    // callback to eclIO serial writeTimeStep method
+    void run()
+    {
+        this->eclIO_.writeTimeStep(this->actionState_,
+                                   this->wtestState_,
+                                   this->summaryState_,
+                                   this->udqState_,
+                                   this->reportStepNum_,
+                                   this->isSubStep_,
+                                   this->secondsElapsed_,
+                                   std::move(this->restartValue_),
+                                   this->writeDoublePrecision_);
+    }
+};
+
+}
+
+namespace Opm {
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+EclGenericWriter(const Schedule& schedule,
+                 const EclipseState& eclState,
+                 const SummaryConfig& summaryConfig,
+                 const Grid& grid,
+                 const EquilGrid* equilGrid,
+                 const GridView& gridView,
+                 const Dune::CartesianIndexMapper<Grid>& cartMapper,
+                 const Dune::CartesianIndexMapper<EquilGrid>* equilCartMapper,
+                 bool enableAsyncOutput,
+                 bool enableEsmry )
+    : collectToIORank_(grid,
+                       equilGrid,
+                       gridView,
+                       cartMapper,
+                       equilCartMapper,
+                       summaryConfig.fip_regions_interreg_flow())
+    , grid_           (grid)
+    , gridView_       (gridView)
+    , schedule_       (schedule)
+    , eclState_       (eclState)
+    , summaryConfig_  (summaryConfig)
+    , cartMapper_     (cartMapper)
+    , equilCartMapper_(equilCartMapper)
+    , equilGrid_      (equilGrid)
+{
+    if (this->collectToIORank_.isIORank()) {
+        this->eclIO_ = std::make_unique<EclipseIO>
+            (this->eclState_,
+             UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
+             this->schedule_, this->summaryConfig_, "", enableEsmry);
+    }
+
+    // create output thread if enabled and rank is I/O rank
+    // async output is enabled by default if pthread are enabled
+    int numWorkerThreads = 0;
+    if (enableAsyncOutput && collectToIORank_.isIORank()) {
+        numWorkerThreads = 1;
+    }
+
+    this->taskletRunner_.reset(new TaskletRunner(numWorkerThreads));
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+const EclipseIO& EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+eclIO() const
+{
+    assert(eclIO_);
+    return *eclIO_;
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+writeInit(const std::function<unsigned int(unsigned int)>& map)
+{
+    if (collectToIORank_.isIORank()) {
+        std::map<std::string, std::vector<int>> integerVectors;
+        if (collectToIORank_.isParallel()) {
+            integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks());
+        }
+
+        auto cartMap = cartesianToCompressed(equilGrid_->size(0), UgGridHelpers::globalCell(*equilGrid_));
+
+        eclIO_->writeInitial(computeTrans_(cartMap, map),
+                             integerVectors,
+                             exportNncStructure_(cartMap, map));
+    }
+
+#if HAVE_MPI
+    if (collectToIORank_.isParallel()) {
+        const auto& comm = grid_.comm();
+        Opm::EclMpiSerializer ser(comm);
+        ser.broadcast(outputNnc_);
+    }
+#endif
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+data::Solution EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+computeTrans_(const std::unordered_map<int,int>& cartesianToActive,
+              const std::function<unsigned int(unsigned int)>& map) const
+{
+    const auto& cartMapper = *equilCartMapper_;
+    const auto& cartDims = cartMapper.cartesianDimensions();
+    const int globalSize = cartDims[0] * cartDims[1] * cartDims[2];
+
+    auto tranx = data::CellData {
+        UnitSystem::measure::transmissibility,
+        std::vector<double>(globalSize, 0.0),
+        data::TargetType::INIT
+    };
+
+    auto trany = tranx;
+    auto tranz = tranx;
+
+    using GlobalGridView = typename EquilGrid::LeafGridView;
+    const GlobalGridView& globalGridView = equilGrid_->leafGridView();
+    using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
+    GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
+
+    for (const auto& elem : elements(globalGridView)) {
+        for (const auto& is : intersections(globalGridView, elem)) {
+            if (!is.neighbor())
+                continue; // intersection is on the domain boundary
+
+            unsigned c1 = globalElemMapper.index(is.inside());
+            unsigned c2 = globalElemMapper.index(is.outside());
+
+            if (c1 > c2)
+                continue; // we only need to handle each connection once, thank you.
+
+            // Ordering of compressed and uncompressed index should be the same
+            const int cartIdx1 = cartMapper.cartesianIndex( c1 );
+            const int cartIdx2 = cartMapper.cartesianIndex( c2 );
+            // Ordering of compressed and uncompressed index should be the same
+            assert(cartIdx1 <= cartIdx2);
+            int gc1 = std::min(cartIdx1, cartIdx2);
+            int gc2 = std::max(cartIdx1, cartIdx2);
+
+            // Re-ordering in case of non-empty mapping between equilGrid to grid
+            if (map) {
+                c1 = map(c1); // equilGridToGrid map
+                c2 = map(c2);
+            }
+
+            if (gc2 - gc1 == 1 && cartDims[0] > 1 ) {
+                tranx.data[gc1] = globalTrans().transmissibility(c1, c2);
+                continue; // skip other if clauses as they are false, last one needs some computation
+            }
+
+            if (gc2 - gc1 == cartDims[0] && cartDims[1] > 1) {
+                trany.data[gc1] = globalTrans().transmissibility(c1, c2);
+                continue; // skipt next if clause as it needs some computation
+            }
+
+            if ( gc2 - gc1 == cartDims[0]*cartDims[1] ||
+                 directVerticalNeighbors(cartDims, cartesianToActive, gc1, gc2))
+                tranz.data[gc1] = globalTrans().transmissibility(c1, c2);
+        }
+    }
+
+    return {
+        {"TRANX", tranx},
+        {"TRANY", trany},
+        {"TRANZ", tranz},
+    };
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+std::vector<NNCdata> EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+exportNncStructure_(const std::unordered_map<int,int>& cartesianToActive, const std::function<unsigned int(unsigned int)>& map) const
+{
+    std::size_t nx = eclState_.getInputGrid().getNX();
+    std::size_t ny = eclState_.getInputGrid().getNY();
+    auto nncData = eclState_.getInputNNC().input();
+    const auto& unitSystem = eclState_.getDeckUnitSystem();
+
+    for( const auto& entry : nncData ) {
+        // test whether NNC is not a neighboring connection
+        // cell2>=cell1 holds due to sortNncAndApplyEditnnc
+        assert( entry.cell2 >= entry.cell1 );
+        auto cellDiff = entry.cell2 - entry.cell1;
+
+        if (cellDiff != 1 && cellDiff != nx && cellDiff != nx*ny) {
+            auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, entry.trans);
+            // Eclipse ignores NNCs (with EDITNNC applied) that are small. Seems like the threshold is 1.0e-6
+            if ( tt >= 1.0e-6 )
+                outputNnc_.emplace_back(entry.cell1, entry.cell2, entry.trans);
+        }
+    }
+
+    using GlobalGridView = typename EquilGrid::LeafGridView;
+    const GlobalGridView& globalGridView = equilGrid_->leafGridView();
+    using GlobElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GlobalGridView>;
+    GlobElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout());
+
+    // Cartesian index mapper for the serial I/O grid
+    const auto& equilCartMapper =  *equilCartMapper_;
+    const auto& cartDims = cartMapper_.cartesianDimensions();
+    for (const auto& elem : elements(globalGridView)) {
+        for (const auto& is : intersections(globalGridView, elem)) {
+            if (!is.neighbor())
+                continue; // intersection is on the domain boundary
+
+            unsigned c1 = globalElemMapper.index(is.inside());
+            unsigned c2 = globalElemMapper.index(is.outside());
+
+            if (c1 > c2)
+                continue; // we only need to handle each connection once, thank you.
+
+            std::size_t cc1 = equilCartMapper.cartesianIndex( c1 );
+            std::size_t cc2 = equilCartMapper.cartesianIndex( c2 );
+
+            if ( cc2 < cc1 )
+                std::swap(cc1, cc2);
+
+            auto cellDiff = cc2 - cc1;
+
+            // Re-ordering in case of non-empty mapping between equilGrid to grid
+            if (map) {
+                c1 = map(c1); // equilGridToGrid map
+                c2 = map(c2);
+            }
+
+            if (cellDiff != 1 &&
+                cellDiff != nx &&
+                cellDiff != nx*ny &&
+                !directVerticalNeighbors(cartDims, cartesianToActive, cc1, cc2)) {
+                // We need to check whether an NNC for this face was also specified
+                // via the NNC keyword in the deck (i.e. in the first origNncSize entries.
+                auto t = globalTrans().transmissibility(c1, c2);
+                auto candidate = std::lower_bound(nncData.begin(), nncData.end(), NNCdata(cc1, cc2, 0.0));
+
+                while ( candidate != nncData.end() && candidate->cell1 == cc1
+                     && candidate->cell2 == cc2) {
+                    t -= candidate->trans;
+                    ++candidate;
+                }
+                // eclipse ignores NNCs with zero transmissibility (different threshold than for NNC
+                // with corresponding EDITNNC above). In addition we do set small transmissibilties
+                // to zero when setting up the simulator. These will be ignored here, too.
+                auto tt = unitSystem.from_si(UnitSystem::measure::transmissibility, std::abs(t));
+                if ( tt > 1e-12 )
+                    outputNnc_.push_back({cc1, cc2, t});
+            }
+        }
+    }
+    return outputNnc_;
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+doWriteOutput(const int                     reportStepNum,
+              const bool                    isSubStep,
+              data::Solution&&              localCellData,
+              data::Wells&&                 localWellData,
+              data::GroupAndNetworkValues&& localGroupAndNetworkData,
+              data::Aquifers&&              localAquiferData,
+              WellTestState&&               localWTestState,
+              const Action::State& actionState,
+              const UDQState& udqState,
+              const SummaryState& summaryState,
+              const std::vector<Scalar>& thresholdPressure,
+              Scalar curTime,
+              Scalar nextStepSize,
+              bool doublePrecision,
+              bool isFlowsn,
+              std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& flowsn,
+              bool isFloresn,
+              std::array<std::pair<std::string, std::pair<std::vector<int>, std::vector<double>>>, 3>&& floresn)
+{
+    const auto isParallel = this->collectToIORank_.isParallel();
+    const bool needsReordering = this->collectToIORank_.doesNeedReordering();
+
+    RestartValue restartValue {
+        (isParallel || needsReordering)
+        ? this->collectToIORank_.globalCellData()
+        : std::move(localCellData),
+
+        isParallel ? this->collectToIORank_.globalWellData()
+                   : std::move(localWellData),
+
+        isParallel ? this->collectToIORank_.globalGroupAndNetworkData()
+                   : std::move(localGroupAndNetworkData),
+
+        isParallel ? this->collectToIORank_.globalAquiferData()
+                   : std::move(localAquiferData)
+    };
+
+    if (eclState_.getSimulationConfig().useThresholdPressure()) {
+        restartValue.addExtra("THRESHPR", UnitSystem::measure::pressure,
+                              thresholdPressure);
+    }
+
+    // Add suggested next timestep to extra data.
+    if (! isSubStep) {
+        restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
+    }
+
+    // Add nnc flows and flores.
+    if (isFlowsn) {
+        const auto flowsn_global = isParallel ? this->collectToIORank_.globalFlowsn() : std::move(flowsn);
+        for (const auto& flows : flowsn_global) {
+            if (flows.first.empty())
+                continue;
+            if (flows.first == "FLOGASN+") {
+                restartValue.addExtra(flows.first, UnitSystem::measure::gas_surface_rate, flows.second.second);
+            }
+            else {
+                restartValue.addExtra(flows.first, UnitSystem::measure::liquid_surface_rate, flows.second.second);
+            }
+        }
+    }
+    if (isFloresn) {
+        const auto floresn_global = isParallel ? this->collectToIORank_.globalFloresn() : std::move(floresn);
+        for (const auto& flores : floresn_global) {
+            if (flores.first.empty())
+                continue;
+            restartValue.addExtra(flores.first, UnitSystem::measure::rate, flores.second.second);
+        }
+    }
+
+    // first, create a tasklet to write the data for the current time
+    // step to disk
+    auto eclWriteTasklet = std::make_shared<EclWriteTasklet>(
+        actionState,
+        isParallel ? this->collectToIORank_.globalWellTestState() : std::move(localWTestState),
+        summaryState, udqState, *this->eclIO_,
+        reportStepNum, isSubStep, curTime, std::move(restartValue), doublePrecision);
+
+    // then, make sure that the previous I/O request has been completed
+    // and the number of incomplete tasklets does not increase between
+    // time steps
+    this->taskletRunner_->barrier();
+
+    // finally, start a new output writing job
+    this->taskletRunner_->dispatch(std::move(eclWriteTasklet));
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+void EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+evalSummary(const int                                            reportStepNum,
+            const Scalar                                         curTime,
+            const data::Wells&                                   localWellData,
+            const data::WellBlockAveragePressures&               localWBPData,
+            const data::GroupAndNetworkValues&                   localGroupAndNetworkData,
+            const std::map<int,data::AquiferData>&               localAquiferData,
+            const std::map<std::pair<std::string, int>, double>& blockData,
+            const std::map<std::string, double>&                 miscSummaryData,
+            const std::map<std::string, std::vector<double>>&    regionData,
+            const Inplace&                                       inplace,
+            const Inplace&                                       initialInPlace,
+            const EclInterRegFlowMap&                            interRegFlows,
+            SummaryState&                                        summaryState,
+            UDQState&                                            udqState)
+{
+    if (collectToIORank_.isIORank()) {
+        const auto& summary = eclIO_->summary();
+
+        const auto& wellData = this->collectToIORank_.isParallel()
+            ? this->collectToIORank_.globalWellData()
+            : localWellData;
+
+        const auto& wbpData = this->collectToIORank_.isParallel()
+            ? this->collectToIORank_.globalWBPData()
+            : localWBPData;
+
+        const auto& groupAndNetworkData = this->collectToIORank_.isParallel()
+            ? this->collectToIORank_.globalGroupAndNetworkData()
+            : localGroupAndNetworkData;
+
+        const auto& aquiferData = this->collectToIORank_.isParallel()
+            ? this->collectToIORank_.globalAquiferData()
+            : localAquiferData;
+
+        summary.eval(summaryState,
+                     reportStepNum,
+                     curTime,
+                     wellData,
+                     wbpData,
+                     groupAndNetworkData,
+                     miscSummaryData,
+                     initialInPlace,
+                     inplace,
+                     regionData,
+                     blockData,
+                     aquiferData,
+                     getInterRegFlowsAsMap(interRegFlows));
+
+        // Off-by-one-fun: The reportStepNum argument corresponds to the
+        // report step these results will be written to, whereas the
+        // argument to UDQ function evaluation corresponds to the report
+        // step we are currently on.
+        const auto udq_step = reportStepNum - 1;
+
+        this->schedule_.getUDQConfig(udq_step)
+            .eval(udq_step,
+                  this->schedule_,
+                  this->schedule_.wellMatcher(udq_step),
+                  summaryState,
+                  udqState);
+    }
+
+#if HAVE_MPI
+    if (collectToIORank_.isParallel()) {
+        EclMpiSerializer ser(grid_.comm());
+        ser.append(summaryState);
+    }
+#endif
+}
+
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+const typename EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::TransmissibilityType&
+EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+globalTrans() const
+{
+    assert (globalTrans_);
+    return *globalTrans_;
+}
+
+
+} // namespace Opm
+#endif

ebos/eclpolyhedralgridvanguard.hh

@@ -94,15 +94,12 @@ public:
     using GridView = GetPropType<TypeTag, Properties::GridView>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     using EquilCartesianIndexMapper = Dune::CartesianIndexMapper<EquilGrid>;
-    using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper, CartesianIndexMapper, Scalar>;
     static constexpr int dimension = Grid::dimension;
     static constexpr int dimensionworld = Grid::dimensionworld;
 
 private:
     using GridPointer = Grid*;
     using EquilGridPointer = EquilGrid*;
-    //using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
-    //using CartesianIndexMapperPointer = std::unique_ptr<CartesianIndexMapper>;
 
 public:
     using TransmissibilityType = EclTransmissibility<Grid, GridView, ElementMapper,

ebos/equil/initstateequil.cc

@@ -21,13 +21,22 @@
   copyright holders.
 */
 
-#include "config.h"
+#include <config.h>
-#include "initstateequil.hh"
+#include <ebos/equil/initstateequil.hh>
-#include "initstateequil_impl.hh"
+#include <ebos/equil/initstateequil_impl.hh>
+
+#include <opm/grid/CpGrid.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
+#include <dune/fem/gridpart/common/gridpart2gridview.hh>
+#include <ebos/femcpgridcompat.hh>
+#endif
 
 namespace Opm {
 namespace EQUIL {
 namespace DeckDependent {
+
 #if HAVE_DUNE_FEM
 using GridView = Dune::Fem::GridPart2GridViewImpl<
                                      Dune::Fem::AdaptiveLeafGridPart<
@@ -59,37 +68,6 @@ template InitialStateComputer<BlackOilFluidSystem<double>,
                          const double,
                          const int,
                          const bool);
-#if HAVE_DUNE_ALUGRID
-#if HAVE_MPI
-using ALUGridComm = Dune::ALUGridMPIComm;
-#else
-using ALUGridComm = Dune::ALUGridNoComm;
-#endif //HAVE_MPI
-using ALUGrid3CN = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, ALUGridComm>;
-using ALUGridView = Dune::GridView<Dune::ALU3dLeafGridViewTraits<const ALUGrid3CN, Dune::PartitionIteratorType(4)>>;
-using ALUGridMapper = Dune::MultipleCodimMultipleGeomTypeMapper<ALUGridView>;
-template class InitialStateComputer<BlackOilFluidSystem<double>,
-                                    ALUGrid3CN,
-                                    ALUGridView,
-                                    ALUGridMapper,
-                                    Dune::CartesianIndexMapper<ALUGrid3CN>>;
-
-template InitialStateComputer<BlackOilFluidSystem<double>,
-                              ALUGrid3CN,
-                              ALUGridView,
-                              ALUGridMapper,
-                              Dune::CartesianIndexMapper<ALUGrid3CN>>::
-    InitialStateComputer(MatLaw&,
-                         const EclipseState&,
-                         const ALUGrid3CN&,
-                         const ALUGridView&,
-                         const Dune::CartesianIndexMapper<ALUGrid3CN>&,
-                         const double,
-                         const int,
-                         const bool);
-#endif //HAVE_DUNE_ALUGRID
-
-
 } // namespace DeckDependent
 
 namespace Details {

ebos/equil/initstateequil_impl.hh

@@ -29,7 +29,6 @@
 
 #include <opm/common/OpmLog/OpmLog.hpp>
 
-#include <opm/grid/CpGrid.hpp>
 #include <opm/grid/utility/RegionMapping.hpp>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
@@ -54,17 +53,6 @@
 #include <limits>
 #include <stdexcept>
 
-#if HAVE_DUNE_FEM
-#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
-#include <dune/fem/gridpart/common/gridpart2gridview.hh>
-#include <ebos/femcpgridcompat.hh>
-#endif
-
-#if HAVE_DUNE_ALUGRID
-#include <dune/alugrid/grid.hh>
-#include <dune/alugrid/3d/gridview.hh>
-#endif // HAVE_DUNE_ALUGRID
-
 namespace Opm {
 namespace EQUIL {
 

flow/flow_blackoil_alugrid.cpp

@@ -16,11 +16,26 @@
   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 <config.h>
-#include <opm/simulators/flow/Main.hpp>
+
 #include <dune/alugrid/grid.hh>
 #include <ebos/eclalugridvanguard.hh>
+#include <opm/simulators/flow/Main.hpp>
+
+// for equilgrid in writer
+// need to include this before eclgenericwriter_impl.hh due to specializations.
+#include <opm/grid/CpGrid.hpp>
+#include <opm/grid/cpgrid/GridHelpers.hpp>
+
+// these are not explicitly instanced in library
+#include <ebos/collecttoiorank_impl.hh>
 #include <ebos/eclgenericproblem_impl.hh>
+#include <ebos/eclgenericthresholdpressure_impl.hh>
+#include <ebos/eclgenerictracermodel_impl.hh>
+#include <ebos/eclgenericwriter_impl.hh>
+#include <ebos/ecltransmissibility_impl.hh>
+#include <ebos/equil/initstateequil_impl.hh>
+
 namespace Opm {
 namespace Properties {
 namespace TTag {

flow/flow_blackoil_polyhedralgrid.cpp

@@ -16,14 +16,24 @@
   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 <config.h>
-#include <opm/simulators/flow/Main.hpp>
+
-#include <opm/grid/polyhedralgrid.hh>
 #include <ebos/eclpolyhedralgridvanguard.hh>
-#include <ebos/equil/initstateequil_impl.hh>
+
+#include <opm/grid/polyhedralgrid.hh>
 #include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
 #include <opm/models/discretization/common/tpfalinearizer.hh>
+#include <opm/simulators/flow/Main.hpp>
+
+// these are not explicitly instanced in library
+#include <ebos/collecttoiorank_impl.hh>
 #include <ebos/eclgenericproblem_impl.hh>
+#include <ebos/eclgenericthresholdpressure_impl.hh>
+#include <ebos/eclgenerictracermodel_impl.hh>
+#include <ebos/ecltransmissibility_impl.hh>
+#include <ebos/eclgenericwriter_impl.hh>
+#include <ebos/equil/initstateequil_impl.hh>
+
 namespace Opm {
 namespace Properties {
     namespace TTag {