-- fixed out of bound nnc

-- still observer some race condition on strange files with temperature

Author: hnil

opm/grid/MinpvProcessor.cpp

@@ -123,7 +123,7 @@ MinpvProcessor::process(const std::vector<double>& thickness,
                 bool c_active = actnum.empty() || actnum[c];
                 bool c_thin = (thickness[c] <= z_tolerance);
                 bool c_thin_inactive = !c_active && c_thin;
-                bool c_low_pv_active = pv[c] < minpvv[c] && c_active;
+                bool c_low_pv_active = (pv[c] < minpvv[c] && c_active) || (c_thin && c_active);
 
                 if (c_low_pv_active || c_thin_inactive) {
                     std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
@@ -166,7 +166,7 @@ MinpvProcessor::process(const std::vector<double>& thickness,
                     bool active = actnum.empty() || actnum[c_below];
                     bool thin = (thickness[c_below] <= z_tolerance);
                     bool thin_inactive = !active && thin;
-                    bool low_pv_active = pv[c_below] < minpvv[c_below] && active;
+                    bool low_pv_active = ((pv[c_below] < minpvv[c_below]) && active) || (thin && active);
 
 
                     while ( (thin_inactive || low_pv_active) && kk_iter < dims_[2] )
@@ -223,13 +223,26 @@ MinpvProcessor::process(const std::vector<double>& thickness,
                         active = actnum.empty() || actnum[c_below];
                         thin = (thickness[c_below] <= z_tolerance);
                         thin_inactive = (!actnum.empty() && !actnum[c_below]) && thin;
-                        low_pv_active = pv[c_below] < minpvv[c_below] && active;
+                        low_pv_active = ((pv[c_below] < minpvv[c_below]) && active) || (thin && active);
                     }
 
                     // create nnc if false or merge the cells if true
-                    if (mergeMinPVCells && c_low_pv_active) {
+                    if (mergeMinPVCells){// && c_low_pv_active) {
+                        // try to make a topological connected grid
+                        // in Flow this is currently called only for edge_conformal grids
+                        // however zcorn inbetween is not modified to make zcorn sorted
                         // Set lower k coordinates of cell below to upper cells's coordinates.
                         // i.e fill the void using the cell below
+                        if (kk==0  || kk_iter == dims_[2]) {
+                            kk = kk_iter;
+                            continue;
+                        }
+                        //bottom cell not active, hence no nnc is created
+                        if (!actnum.empty() && !actnum[c_below]) {
+                            kk = kk_iter;
+                            continue;
+                        }
+
                         std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk_iter, zcorn);
                         for (int count = 0; count < 4; ++count) {
                             cz_below[count] = cz[count];
@@ -258,15 +271,15 @@ MinpvProcessor::process(const std::vector<double>& thickness,
                         auto above_active = actnum.empty() || actnum[c_above];
                         auto above_inactive = !actnum.empty() && !actnum[c_above];
                         auto above_thin = thickness[c_above] < z_tolerance;
-                        auto above_small_pv = pv[c_above] < minpvv[c_above];
+                        auto above_small_pv = (pv[c_above] < minpvv[c_above]) ||  above_thin;
 
                         if ((above_inactive && above_thin) || (above_active && above_small_pv
                                                                && (!pinchNOGAP || above_thin) ) ) {
                             for (k_above = kk - 2; k_above > 0; --k_above) {
                                 c_above = ii + dims_[0] * (jj + dims_[1] * (k_above));
                                 above_active = actnum.empty() || actnum[c_above];
                                 above_inactive = !actnum.empty() && !actnum[c_above];
-                                auto above_significant_pv = pv[c_above] > minpvv[c_above];
+                                auto above_significant_pv = !((pv[c_above] < minpvv[c_above]) || (thickness[c_above] < z_tolerance));
                                 auto above_broad = thickness[c_above] > z_tolerance;
 
                                 // \todo if condition seems wrong and should be the negation of above?
@@ -291,25 +304,30 @@ MinpvProcessor::process(const std::vector<double>& thickness,
                         option4ALLZero = option4ALLZero || (!permz.empty() && permz[c_above] == 0.0) || multz(c_above) == 0.0;
                         nnc_allowed = nnc_allowed && (computeGap(cz_above, cz_below) < max_gap) && (!pinchOption4ALL || !option4ALLZero) ;
 
+                        //bool
+                        above_small_pv = (pv[c_above] < minpvv[c_above]) || (thickness[c_above] < z_tolerance);
+                        bool below_small_pv = (pv[c_below] < minpvv[c_below]) || (thickness[c_below] < z_tolerance);
                         if ( nnc_allowed &&
                              (actnum.empty() || (actnum[c_above] && actnum[c_below])) &&
-                             pv[c_above] > minpvv[c_above] && pv[c_below] > minpvv[c_below]) {
+                             !(above_small_pv) && !(below_small_pv) ){
                             result.add_nnc(c_above, c_below);
                         }
                         kk = kk_iter;
                     }
                 }
                 else
                 {
-                    if (kk < dims_[2] - 1 && (actnum.empty() || actnum[c]) && pv[c] > minpvv[c] &&
+                    if (kk < dims_[2] - 1 && (actnum.empty() || actnum[c]) && !((pv[c] < minpvv[c]) || (thickness[c] < z_tolerance)) &&
                         multz(c) != 0.0)
                     {
                         // Check whether there is a gap to the neighbor below whose thickness is less
                         // than MAX_GAP. In that case we need to create an NNC if there is a gap between the two cells.
                         int kk_below = kk + 1;
                         int c_below = ii + dims_[0] * (jj + dims_[1] * kk_below);
 
-                        if ((actnum.empty() || actnum[c_below]) && pv[c_below] > minpvv[c_below])
+                        if ( (actnum.empty() || actnum[c_below])
+                             && 
+                             !((pv[c_below] < minpvv[c_below])  || (thickness[c_below] < z_tolerance) ) ) 
                         {
                             // Check MAX_GAP threshold
                             std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);

opm/grid/common/GridPartitioning.cpp

@@ -576,7 +576,7 @@ int addOverlapLayer([[maybe_unused]] const CpGrid& grid,
         int index = ix.index(*it);
         auto owner = cell_part[index];
         exportProcs.insert(std::make_pair(owner, 0));
-        if ( trans ) {
+        if ( trans && false) {
             addOverlapLayerNoZeroTrans(grid, index, *it, owner, cell_part, exportList, addCornerCells, layers-1, trans, level);
         }
         else {

opm/grid/cpgrid/processEclipseFormat.cpp

@@ -204,8 +204,10 @@ namespace cpgrid
                     return transMult.getMultiplier(cartindex, ::Opm::FaceDir::ZPlus) *
                         transMult.getMultiplier(cartindex, ::Opm::FaceDir::ZMinus);
                 };
+                bool mergeMinPVCells = edge_conformal;// try to make geometrically connected grids
+
                 minpv_result = mp.process(thickness, z_tolerance, ecl_grid.getPinchMaxEmptyGap(),
-                                          poreVolume, ecl_grid.getMinpvVector(), actnumData, false,
+                                          poreVolume, ecl_grid.getMinpvVector(), actnumData, edge_conformal,
                                           zcornData.data(), nogap, pinchOptionALL,
                                           permZ, multZ, tolerance_unique_points);
                 if (!minpv_result.nnc.empty()) {
@@ -223,6 +225,9 @@ namespace cpgrid
 
             // Add PINCH NNCs.
             std::vector<Opm::NNCdata> pinchedNNCs;
+            if(edge_conformal){
+                assert(minpv_result.nnc.size() == 0);
+            }
 
             for (const auto& [cell1, cell2] : minpv_result.nnc) {
                 nnc_cells[PinchNNC].insert({cell1, cell2});
@@ -416,12 +421,19 @@ namespace cpgrid
         }
         else {
             // Make the grid.
+            auto pinchActive_copy = pinchActive;
+            if(edge_conformal){
+                // In edge conformal grids we need to set the pinchActive flag to true
+                pinchActive_copy = true;
+                assert(nnc_cells[PinchNNC].empty());
+                assert(nnc_cells[ExplicitNNC].empty());
+            }
             this->processEclipseFormat(g,
                                        ecl_state,
                                        nnc_cells,
                                        false,
                                        turn_normals,
-                                       pinchActive,
+                                       pinchActive_copy,
                                        tolerance_unique_points,
                                        edge_conformal);
         }