diff --git a/flow/flow_blackoil_alugrid.cpp b/flow/flow_blackoil_alugrid.cpp index cf2cfb7ee..6432a2efc 100644 --- a/flow/flow_blackoil_alugrid.cpp +++ b/flow/flow_blackoil_alugrid.cpp @@ -17,6 +17,9 @@ along with OPM. If not, see . */ #include +// The default in dune-ALUGrird is "DISABLE_ALUGRID_SFC_ORDERING is false" +#undef DISABLE_ALUGRID_SFC_ORDERING +#undef USE_ALUGRID_SFC_ORDERING #include #include diff --git a/opm/simulators/flow/AluGridVanguard.hpp b/opm/simulators/flow/AluGridVanguard.hpp index d65fd8bdf..807455b8d 100644 --- a/opm/simulators/flow/AluGridVanguard.hpp +++ b/opm/simulators/flow/AluGridVanguard.hpp @@ -325,7 +325,8 @@ protected: factory_ = std::make_unique(); grid_ = factory_->convert(*equilGrid_, cartesianCellId_, ordering_); - OpmLog::warning("Space Filling Curve Ordering is not yet supported: DISABLE_ALUGRID_SFC_ORDERING is enabled"); + OpmLog::warning("Space Filling Curve Ordering (SFCO) is used"); + OpmLog::warning("Add '#define DISABLE_ALUGRID_SFC_ORDERING 1' to disable SFC reordering"); equilGridToGrid_.resize(ordering_.size()); for (std::size_t index = 0; index < ordering_.size(); ++index) { equilGridToGrid_[ordering_[index]] = index; diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index a9d359bb7..62eadfb31 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -451,10 +451,11 @@ public: const auto& residual = this->model().linearizer().residual(); for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); globalDofIdx ++) { - this->drift_[globalDofIdx] = residual[globalDofIdx] * simulator.timeStepSize(); + int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(globalDofIdx); + this->drift_[sfcdofIdx] = residual[sfcdofIdx] * simulator.timeStepSize(); if constexpr (getPropValue()) { - this->drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx); + this->drift_[sfcdofIdx] *= this->model().dofTotalVolume(sfcdofIdx); } } } @@ -1375,14 +1376,15 @@ protected: const auto& fp = eclState.fieldProps(); const std::vector porvData = this -> fieldPropDoubleOnLeafAssigner_()(fp, "PORV"); for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { + int sfcdofIdx = simulator.vanguard().gridEquilIdxToGridIdx(dofIdx); Scalar poreVolume = porvData[dofIdx]; // we define the porosity as the accumulated pore volume divided by the // geometric volume of the element. Note that -- in pathetic cases -- it can // be larger than 1.0! - Scalar dofVolume = simulator.model().dofTotalVolume(dofIdx); + Scalar dofVolume = simulator.model().dofTotalVolume(sfcdofIdx); assert(dofVolume > 0.0); - this->referencePorosity_[/*timeIdx=*/0][dofIdx] = poreVolume/dofVolume; + this->referencePorosity_[/*timeIdx=*/0][sfcdofIdx] = poreVolume/dofVolume; } } diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index af1ecda84..c3c87fb8b 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -134,7 +134,7 @@ namespace Opm { const auto& fs = intQuants.fluidState(); // use pore volume weighted averages. const Scalar pv_cell = - simulator.model().dofTotalVolume(cellIdx) + simulator.model().dofTotalVolume(simulator.vanguard().gridEquilIdxToGridIdx(cellIdx)) * intQuants.porosity().value(); // only count oil and gas filled parts of the domain