add changes to fix the order

enabling reordering as default
OPM · Sep 25, 2024 · 6c2e957 · 6c2e957
1 parent 02555ac
commit 6c2e957
Show file tree

Hide file tree

Showing 4 changed files with 12 additions and 6 deletions.
diff --git a/flow/flow_blackoil_alugrid.cpp b/flow/flow_blackoil_alugrid.cpp
@@ -17,6 +17,9 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 #include <config.h>
+// The default in dune-ALUGrird is "DISABLE_ALUGRID_SFC_ORDERING is false"
+#undef DISABLE_ALUGRID_SFC_ORDERING
+#undef USE_ALUGRID_SFC_ORDERING
 
 #include <dune/alugrid/grid.hh>
 #include <opm/simulators/flow/Main.hpp>

diff --git a/opm/simulators/flow/AluGridVanguard.hpp b/opm/simulators/flow/AluGridVanguard.hpp
@@ -325,7 +325,8 @@ class AluGridVanguard : public FlowBaseVanguard<TypeTag>
 
         factory_ = std::make_unique<Factory>();
         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
@@ -451,10 +451,11 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
             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<TypeTag, Properties::UseVolumetricResidual>()) {
-                    this->drift_[globalDofIdx] *= this->model().dofTotalVolume(globalDofIdx);
+                    this->drift_[sfcdofIdx] *= this->model().dofTotalVolume(sfcdofIdx);
                 }
             }
         }
@@ -1375,14 +1376,15 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
         const auto& fp = eclState.fieldProps();
         const std::vector<double> 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
@@ -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