Merge pull request #5630 from ElyesAhmed/fix_sfc_order_alugrid

add changes to fix the order

flow/flow_blackoil_alugrid.cpp

@@ -17,6 +17,15 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 #include <config.h>
+// By default, dune-ALUGrid uses Space-Filling Curve (SFC) ordering
+// for cells to optimize data access patterns for adaptive mesh refinements/coarsening.
+// However, if you want to use Cartesian ordering, as is used in OPM, you
+// can switch to it by defining the macro #define DISABLE_SFC_ORDERING 1.
+// This will change the default cell order to Cartesian.
+// Note that this option is not available for pre-built or installed versions of dune-ALUGrid.
+// To enable changig to Cartesian ordering, you will need to rebuild dune-ALUGrid from source, ensuring 
+// the build configuration allows disabling SFC ordering from OPM.
+// For more details, refer to the files gridfactory.hh and aluinline.hh located in the dune-alugrid/3d/
 
 #include <dune/alugrid/grid.hh>
 #include <opm/simulators/flow/Main.hpp>

opm/simulators/flow/AluGridVanguard.hpp

@@ -325,7 +325,7 @@ protected:
 
         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 (SFC) ordering is enabled: see flow_blackoil_alugrid for more informations on disabling/enabling SFC reordering");
         equilGridToGrid_.resize(ordering_.size());
         for (std::size_t index = 0; index < ordering_.size(); ++index) {
             equilGridToGrid_[ordering_[index]] = index;

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<TypeTag, Properties::UseVolumetricResidual>()) {
-                    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<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;
         }
     }
 

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

regressionTests.cmake

@@ -1496,19 +1496,13 @@ if(BUILD_FLOW_POLY_GRID)
                            DIR spe1)
 endif()
 
-if(dune-alugrid_FOUND AND BUILD_FLOW_ALU_GRID)
+if(dune-alugrid_FOUND AND BUILD_FLOW_ALU_GRID AND MPI_FOUND)
   add_test_compareECLFiles(CASENAME spe12_alugrid
                            FILENAME SPE1CASE2
                            SIMULATOR flow_blackoil_alugrid
                            ABS_TOL ${abs_tol}
                            REL_TOL ${coarse_rel_tol}
-                           DIR spe1
-                           TEST_ARGS --linear-solver=ilu0
-                                     --tolerance-mb=1e-6
-                                     --min-strict-cnv-iter=0
-                                     --local-well-solve-control-switching=false
-                                     --use-implicit-ipr=false
-                                     --enable-drift-compensation=true)
+                           DIR spe1)
 endif()
 
 if(BUILD_FLOW_FLOAT_VARIANTS)