Merge pull request #2302 from joakim-hove/use-bc

Use bc

ebos/eclproblem.hh

@@ -2996,13 +2996,13 @@ private:
         nonTrivialBoundaryConditions_ = false;
         const auto& simulator = this->simulator();
         const auto& vanguard = simulator.vanguard();
-
+        const auto& bcconfig = vanguard.eclState().getSimulationConfig().bcconfig();
-        if (vanguard.deck().hasKeyword("BC")) {
+        if (bcconfig.size() > 0) {
             nonTrivialBoundaryConditions_ = true;
 
             size_t numCartDof = vanguard.cartesianSize();
             unsigned numElems = vanguard.gridView().size(/*codim=*/0);
-            std::vector<int> cartesianToCompressedElemIdx(numCartDof);
+            std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
 
             for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
                 cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
@@ -3020,115 +3020,114 @@ private:
             freebcZ_.resize(numElems, false);
             freebcZMinus_.resize(numElems, false);
 
-            const auto& bcs = vanguard.deck().getKeywordList("BC");
+            for (const auto& bcface : bcconfig) {
-            for (size_t listIdx = 0; listIdx < bcs.size(); ++listIdx) {
+                const auto& type = bcface.bctype;
-                const auto& bc = *bcs[listIdx];
+                if (type == BCType::RATE) {
+                    int compIdx;
 
-                for (size_t record = 0; record < bc.size(); ++record) {
+                    switch (bcface.component) {
-
+                    case BCComponent::OIL:
-                    std::string type = bc.getRecord(record).getItem("TYPE").getTrimmedString(0);
-                    std::string compName = bc.getRecord(record).getItem("COMPONENT").getTrimmedString(0);
-                    int compIdx = -999;
-
-                    if (compName == "OIL")
                         compIdx = oilCompIdx;
-                    else if (compName == "GAS")
+                        break;
+                    case BCComponent::GAS:
                         compIdx = gasCompIdx;
-                    else if (compName == "WATER")
+                        break;
+                    case BCComponent::WATER:
                         compIdx = waterCompIdx;
-                    else if (compName == "SOLVENT")
+                        break;
-                    {
+                    case BCComponent::SOLVENT:
                         if (!enableSolvent)
                             throw std::logic_error("solvent is disabled and you're trying to add solvent to BC");
 
                         compIdx = Indices::solventSaturationIdx;
-                    }
+                        break;
-                    else if (compName == "POLYMER")
+                    case BCComponent::POLYMER:
-                    {
                         if (!enablePolymer)
                             throw std::logic_error("polymer is disabled and you're trying to add polymer to BC");
 
                         compIdx = Indices::polymerConcentrationIdx;
-                    }
+                        break;
-                    else if (compName == "NONE")
+                    case BCComponent::NONE:
-                    {
+                        if (type == BCType::RATE)
-                        if ( type == "RATE")
                             throw std::logic_error("you need to specify the component when RATE type is set in BC");
+                        break;
                     }
-                    else
-                        throw std::logic_error("invalid component name for BC");
 
-                    int i1 = bc.getRecord(record).getItem("I1").template get< int >(0) - 1;
+                    std::vector<RateVector>* data = nullptr;
-                    int i2 = bc.getRecord(record).getItem("I2").template get< int >(0) - 1;
+                    switch (bcface.dir) {
-                    int j1 = bc.getRecord(record).getItem("J1").template get< int >(0) - 1;
+                    case FaceDir::XMinus:
-                    int j2 = bc.getRecord(record).getItem("J2").template get< int >(0) - 1;
+                        data = &massratebcXMinus_;
-                    int k1 = bc.getRecord(record).getItem("K1").template get< int >(0) - 1;
+                        break;
-                    int k2 = bc.getRecord(record).getItem("K2").template get< int >(0) - 1;
+                    case FaceDir::XPlus:
-                    std::string direction = bc.getRecord(record).getItem("DIRECTION").getTrimmedString(0);
+                        data = &massratebcX_;
+                        break;
+                    case FaceDir::YMinus:
+                        data = &massratebcYMinus_;
+                        break;
+                    case FaceDir::YPlus:
+                        data = &massratebcY_;
+                        break;
+                    case FaceDir::ZMinus:
+                        data = &massratebcZMinus_;
+                        break;
+                    case FaceDir::ZPlus:
+                        data = &massratebcZ_;
+                        break;
+                    }
 
-                    if (type == "RATE") {
+                    const Evaluation rate = bcface.rate;
-                        assert(compIdx >= 0);
+                    for (int i = bcface.i1; i <= bcface.i2; ++i) {
-                        std::vector<RateVector>* data = 0;
+                        for (int j = bcface.j1; j <= bcface.j2; ++j) {
-                        if (direction == "X-")
+                            for (int k = bcface.k1; k <= bcface.k2; ++k) {
-                            data = &massratebcXMinus_;
+                                std::array<int, 3> tmp = {i,j,k};
-                        else if (direction == "X")
+                                auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
-                            data = &massratebcX_;
+                                if (elemIdx >= 0)
-                        else if (direction == "Y-")
-                            data = &massratebcYMinus_;
-                        else if (direction == "Y")
-                            data = &massratebcY_;
-                        else if (direction == "Z-")
-                            data = &massratebcZMinus_;
-                        else if (direction == "Z")
-                            data = &massratebcZ_;
-                        else
-                            throw std::logic_error("invalid direction for BC");
-
-                        const Evaluation rate = bc.getRecord(record).getItem("RATE").getSIDouble(0);
-                        for (int i = i1; i <= i2; ++i) {
-                            for (int j = j1; j <= j2; ++j) {
-                                for (int k = k1; k <= k2; ++k) {
-                                    std::array<int, 3> tmp = {i,j,k};
-                                    size_t elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
                                     (*data)[elemIdx][compIdx] = rate;
-                                }
                             }
                         }
-                    } else if (type == "FREE") {
-                        std::vector<bool>* data = 0;
-                        if (direction == "X-")
-                            data = &freebcXMinus_;
-                        else if (direction == "X")
-                            data = &freebcX_;
-                        else if (direction == "Y-")
-                            data = &freebcYMinus_;
-                        else if (direction == "Y")
-                            data = &freebcY_;
-                        else if (direction == "Z-")
-                            data = &freebcZMinus_;
-                        else if (direction == "Z")
-                            data = &freebcZ_;
-                        else
-                            throw std::logic_error("invalid direction for BC");
-
-                        for (int i = i1; i <= i2; ++i) {
-                            for (int j = j1; j <= j2; ++j) {
-                                for (int k = k1; k <= k2; ++k) {
-                                    std::array<int, 3> tmp = {i,j,k};
-                                    size_t elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
-                                    (*data)[elemIdx] = true;
-                                }
-                            }
-                        }
-                        // TODO: either the real initial solution needs to be computed or read from the restart file
-                        const auto& eclState = simulator.vanguard().eclState();
-                        const auto& initconfig = eclState.getInitConfig();
-                        if (initconfig.restartRequested()) {
-                            throw std::logic_error("restart is not compatible with using free boundary conditions");
-                        }
-                    } else {
-                        throw std::logic_error("invalid type for BC. Use FREE or RATE");
                     }
+                } else if (type == BCType::FREE) {
+                    std::vector<bool>* data = nullptr;
+                    switch (bcface.dir) {
+                    case FaceDir::XMinus:
+                        data = &freebcXMinus_;
+                        break;
+                    case FaceDir::XPlus:
+                        data = &freebcX_;
+                        break;
+                    case FaceDir::YMinus:
+                        data = &freebcYMinus_;
+                        break;
+                    case FaceDir::YPlus:
+                        data = &freebcY_;
+                        break;
+                    case FaceDir::ZMinus:
+                        data = &freebcZMinus_;
+                        break;
+                    case FaceDir::ZPlus:
+                        data = &freebcZ_;
+                        break;
+                    }
+
+                    for (int i = bcface.i1; i <= bcface.i2; ++i) {
+                        for (int j = bcface.j1; j <= bcface.j2; ++j) {
+                            for (int k = bcface.k1; k <= bcface.k2; ++k) {
+                                std::array<int, 3> tmp = {i,j,k};
+                                auto elemIdx = cartesianToCompressedElemIdx[vanguard.cartesianIndex(tmp)];
+                                if (elemIdx >= 0)
+                                    (*data)[elemIdx] = true;
+                            }
+                        }
+                    }
+
+                    // TODO: either the real initial solution needs to be computed or read from the restart file
+                    const auto& eclState = simulator.vanguard().eclState();
+                    const auto& initconfig = eclState.getInitConfig();
+                    if (initconfig.restartRequested()) {
+                        throw std::logic_error("restart is not compatible with using free boundary conditions");
+                    }
+                } else {
+                    throw std::logic_error("invalid type for BC. Use FREE or RATE");
                 }
             }
         }

opm/simulators/utils/ParallelRestart.cpp

@@ -64,6 +64,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTracerProperties.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WList.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WListManager.hpp>
+#include <opm/parser/eclipse/EclipseState/SimulationConfig/BCConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
 #include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
@@ -464,6 +465,7 @@ HANDLE_AS_POD(WellBrineProperties)
 HANDLE_AS_POD(Welldims)
 HANDLE_AS_POD(WellFoamProperties)
 HANDLE_AS_POD(WellSegmentDims)
+HANDLE_AS_POD(BCConfig::BCFace)
 
 std::size_t packSize(const data::Well& data, Dune::MPIHelper::MPICommunicator comm)
 {
@@ -513,6 +515,12 @@ std::size_t packSize(const ThresholdPressure& data, Dune::MPIHelper::MPICommunic
           packSize(data.pressureTable(), comm);
 }
 
+std::size_t packSize(const BCConfig& bc, Dune::MPIHelper::MPICommunicator comm)
+{
+    return packSize(bc.faces(), comm);
+}
+
+
 std::size_t packSize(const NNC& data, Dune::MPIHelper::MPICommunicator comm)
 {
    return packSize(data.data(), comm);
@@ -604,6 +612,7 @@ std::size_t packSize(const InitConfig& data, Dune::MPIHelper::MPICommunicator co
 std::size_t packSize(const SimulationConfig& data, Dune::MPIHelper::MPICommunicator comm)
 {
     return packSize(data.getThresholdPressure(), comm) +
+           packSize(data.bcconfig(), comm) +
            packSize(data.useCPR(), comm) +
            packSize(data.hasDISGAS(), comm) +
            packSize(data.hasVAPOIL(), comm) +
@@ -2324,6 +2333,13 @@ void pack(const ThresholdPressure& data, std::vector<char>& buffer, int& positio
     pack(data.pressureTable(), buffer, position, comm);
 }
 
+
+void pack(const BCConfig& bc, std::vector<char>& buffer, int& position,
+    Dune::MPIHelper::MPICommunicator comm)
+{
+    pack(bc.faces(), buffer, position, comm);
+}
+
 void pack(const NNC& data, std::vector<char>& buffer, int& position,
           Dune::MPIHelper::MPICommunicator comm)
 {
@@ -2431,6 +2447,7 @@ void pack(const SimulationConfig& data, std::vector<char>& buffer, int& position
           Dune::MPIHelper::MPICommunicator comm)
 {
     pack(data.getThresholdPressure(), buffer, position, comm);
+    pack(data.bcconfig(), buffer, position, comm);
     pack(data.useCPR(), buffer, position, comm);
     pack(data.hasDISGAS(), buffer, position, comm);
     pack(data.hasVAPOIL(), buffer, position, comm);
@@ -4320,6 +4337,16 @@ void unpack(ThresholdPressure& data, std::vector<char>& buffer, int& position,
     data = ThresholdPressure(active, restart, thpTable, pTable);
 }
 
+
+void unpack(BCConfig& bc, std::vector<char>& buffer, int& position,
+            Dune::MPIHelper::MPICommunicator comm)
+{
+    std::vector<BCConfig::BCFace> faces;
+
+    unpack(faces, buffer, position, comm);
+    bc = BCConfig(faces);
+}
+
 void unpack(NNC& data, std::vector<char>& buffer, int& position,
             Dune::MPIHelper::MPICommunicator comm)
 {
@@ -4465,13 +4492,15 @@ void unpack(SimulationConfig& data, std::vector<char>& buffer, int& position,
             Dune::MPIHelper::MPICommunicator comm)
 {
     ThresholdPressure thresholdPressure;
+    BCConfig bc;
     bool useCPR, DISGAS, VAPOIL, isThermal;
     unpack(thresholdPressure, buffer, position, comm);
+    unpack(bc, buffer, position, comm);
     unpack(useCPR, buffer, position, comm);
     unpack(DISGAS, buffer, position, comm);
     unpack(VAPOIL, buffer, position, comm);
     unpack(isThermal, buffer, position, comm);
-    data = SimulationConfig(thresholdPressure, useCPR, DISGAS, VAPOIL, isThermal);
+    data = SimulationConfig(thresholdPressure, bc, useCPR, DISGAS, VAPOIL, isThermal);
 }
 
 void unpack(TimeMap& data, std::vector<char>& buffer, int& position,

opm/simulators/utils/ParallelRestart.hpp

@@ -75,6 +75,8 @@ namespace Action {
 }
 
 class Aqudims;
+class BCConfig;
+class BCConfig::BCFace;
 class BrineDensityTable;
 class ColumnSchema;
 class Connection;
@@ -656,6 +658,8 @@ ADD_PACK_PROTOTYPES(Action::ASTNode)
 ADD_PACK_PROTOTYPES(Action::Condition)
 ADD_PACK_PROTOTYPES(Action::Quantity)
 ADD_PACK_PROTOTYPES(Aqudims)
+ADD_PACK_PROTOTYPES(BCConfig)
+ADD_PACK_PROTOTYPES(BCConfig::BCFace)
 ADD_PACK_PROTOTYPES(BrineDensityTable)
 ADD_PACK_PROTOTYPES(ColumnSchema)
 ADD_PACK_PROTOTYPES(Connection)

tests/test_ParallelRestart.cpp

@@ -73,6 +73,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WList.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WListManager.hpp>
+#include <opm/parser/eclipse/EclipseState/SimulationConfig/BCConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
 #include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
 #include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
@@ -183,6 +184,11 @@ Opm::ThresholdPressure getThresholdPressure()
 }
 
 
+Opm::BCConfig getBCConfig()
+{
+    return Opm::BCConfig({{10,11,12,13,14,15,Opm::BCType::RATE,Opm::FaceDir::XPlus, Opm::BCComponent::GAS, 100.0}});
+}
+
 Opm::TableSchema getTableSchema()
 {
     Opm::OrderedMap<std::string, Opm::ColumnSchema> data;
@@ -768,7 +774,18 @@ BOOST_AUTO_TEST_CASE(InitConfig)
 BOOST_AUTO_TEST_CASE(SimulationConfig)
 {
 #if HAVE_MPI
-    Opm::SimulationConfig val1(getThresholdPressure(), false, true, false, true);
+    Opm::SimulationConfig val1(getThresholdPressure(), getBCConfig(), false, true, false, true);
+    auto val2 = PackUnpack(val1);
+    BOOST_CHECK(std::get<1>(val2) == std::get<2>(val2));
+    BOOST_CHECK(val1 == std::get<0>(val2));
+#endif
+}
+
+
+BOOST_AUTO_TEST_CASE(BCConfig)
+{
+#if HAVE_MPI
+    Opm::BCConfig val1({{10,11,12,13,14,15,Opm::BCType::RATE, Opm::FaceDir::XPlus, Opm::BCComponent::GAS, 100}});
     auto val2 = PackUnpack(val1);
     BOOST_CHECK(std::get<1>(val2) == std::get<2>(val2));
     BOOST_CHECK(val1 == std::get<0>(val2));