Merge pull request #445 from stefoss23/load_egrid_without_actnum

WIP: ecl_grid can be allocated with int * instead of actnum_kw

lib/CMakeLists.txt

@@ -352,6 +352,7 @@ foreach (name   ecl_alloc_cpgrid
                 ecl_grid_export
                 ecl_grid_init_fwrite
                 ecl_grid_reset_actnum
+                ecl_grid_ext_actnum
                 ecl_init_file
                 ecl_kw_cmp_string
                 ecl_kw_equal

lib/ecl/ecl_grid.cpp

@@ -21,6 +21,7 @@
 #include <stdio.h>
 #include <stdbool.h>
 #include <math.h>
+#include <vector>
 
 #include <ert/util/util.h>
 #include <ert/util/double_vector.hpp>
@@ -2689,9 +2690,9 @@ static ecl_grid_type * ecl_grid_alloc_GRDECL_kw__(ecl_grid_type * global_grid ,
                                                   const ecl_kw_type * gridhead_kw ,
                                                   const ecl_kw_type * zcorn_kw ,
                                                   const ecl_kw_type * coord_kw ,
-                                                  const ecl_kw_type * actnum_kw ,    /* Can be NULL */
                                                   const ecl_kw_type * mapaxes_kw ,   /* Can be NULL */
-                                                  const ecl_kw_type * corsnum_kw) {   /* Can be NULL */
+                                                  const ecl_kw_type * corsnum_kw,    /* Can be NULL */
+                                                  const int * actnum_data) {         /* Can be NULL */
    int gtype, nx,ny,nz, lgr_nr;
 
   gtype   = ecl_kw_iget_int(gridhead_kw , GRIDHEAD_TYPE_INDEX);
@@ -2714,15 +2715,11 @@ static ecl_grid_type * ecl_grid_alloc_GRDECL_kw__(ecl_grid_type * global_grid ,
 
   {
     const float * mapaxes_data = NULL;
-    const int   * actnum_data  = NULL;
     const int   * corsnum_data = NULL;
 
     if (mapaxes_kw != NULL)
       mapaxes_data = ecl_grid_get_mapaxes_from_kw__(mapaxes_kw);
 
-    if (actnum_kw != NULL)
-      actnum_data = ecl_kw_get_int_ptr(actnum_kw);
-
     if (corsnum_kw != NULL)
       corsnum_data = ecl_kw_get_int_ptr( corsnum_kw );
 
@@ -2753,17 +2750,22 @@ ecl_grid_type * ecl_grid_alloc_GRDECL_kw( int nx, int ny , int nz ,
                                           const ecl_kw_type * actnum_kw ,      /* Can be NULL */
                                           const ecl_kw_type * mapaxes_kw ) {   /* Can be NULL */
 
+  const int * actnum_data = NULL;
+  if (actnum_kw)
+    actnum_data = ecl_kw_get_int_ptr(actnum_kw);
+
   bool apply_mapaxes = true;
   ecl_kw_type * gridhead_kw = ecl_grid_alloc_gridhead_kw( nx, ny, nz, 0);
+
   ecl_grid_type * ecl_grid = ecl_grid_alloc_GRDECL_kw__(NULL,
                                                         FILEHEAD_SINGLE_POROSITY,
                                                         apply_mapaxes,
                                                         gridhead_kw,
                                                         zcorn_kw,
                                                         coord_kw,
-                                                        actnum_kw,
                                                         mapaxes_kw,
-                                                        NULL);
+                                                        NULL,
+                                                        actnum_data);
   ecl_kw_free( gridhead_kw );
   return ecl_grid;
 
@@ -2986,7 +2988,7 @@ static void ecl_grid_init_nnc_amalgamated(ecl_grid_type * main_grid, ecl_file_ty
 */
 
 
-static ecl_grid_type * ecl_grid_alloc_EGRID__( ecl_grid_type * main_grid , const ecl_file_type * ecl_file , int grid_nr, bool apply_mapaxes) {
+static ecl_grid_type * ecl_grid_alloc_EGRID__( ecl_grid_type * main_grid , const ecl_file_type * ecl_file , int grid_nr, bool apply_mapaxes, const int * ext_actnum) {
   ecl_kw_type * gridhead_kw  = ecl_file_iget_named_kw( ecl_file , GRIDHEAD_KW  , grid_nr);
   ecl_kw_type * zcorn_kw     = ecl_file_iget_named_kw( ecl_file , ZCORN_KW     , grid_nr);
   ecl_kw_type * coord_kw     = ecl_file_iget_named_kw( ecl_file , COORD_KW     , grid_nr);
@@ -3004,10 +3006,23 @@ static ecl_grid_type * ecl_grid_alloc_EGRID__( ecl_grid_type * main_grid , const
     eclipse_version = main_grid->eclipse_version;
   }
 
-
-  /** If ACTNUM is not present - that is is interpreted as - all active. */
-  if (ecl_file_get_num_named_kw(ecl_file , ACTNUM_KW) > grid_nr)
+  // If ACTNUM and ext_actnum are not present - that is is interpreted as all active. 
+  const int * actnum_data;
+  std::vector<int> actnum_product;
+  if (ecl_file_get_num_named_kw(ecl_file , ACTNUM_KW) > grid_nr) {
     actnum_kw = ecl_file_iget_named_kw( ecl_file , ACTNUM_KW    , grid_nr);
+    actnum_data = ecl_kw_get_int_ptr(actnum_kw);
+    if (ext_actnum) {
+      int size = ecl_kw_get_size(actnum_kw);
+      actnum_product.resize(size);
+      for (int i = 0; i < size; i++)
+        actnum_product[i] = actnum_data[i] * ext_actnum[i];
+      actnum_data = actnum_product.data();
+    }
+  }
+  else
+    actnum_data = ext_actnum;
+
 
   if (grid_nr == 0) {
     /* MAPAXES and COARSENING only apply to the global grid. */
@@ -3018,18 +3033,16 @@ static ecl_grid_type * ecl_grid_alloc_EGRID__( ecl_grid_type * main_grid , const
       corsnum_kw   = ecl_file_iget_named_kw( ecl_file , CORSNUM_KW , 0);
   }
 
-
-
   {
     ecl_grid_type * ecl_grid = ecl_grid_alloc_GRDECL_kw__( main_grid ,
                                                            dualp_flag ,
                                                            apply_mapaxes,
                                                            gridhead_kw ,
                                                            zcorn_kw ,
                                                            coord_kw ,
-                                                           actnum_kw ,
                                                            mapaxes_kw ,
-                                                           corsnum_kw );
+                                                           corsnum_kw, 
+                                                           actnum_data);
 
     if (ECL_GRID_MAINGRID_LGR_NR != grid_nr) ecl_grid_set_lgr_name_EGRID(ecl_grid , ecl_file , grid_nr);
     ecl_grid->eclipse_version = eclipse_version;
@@ -3038,8 +3051,13 @@ static ecl_grid_type * ecl_grid_alloc_EGRID__( ecl_grid_type * main_grid , const
 }
 
 
+static ecl_grid_type * ecl_grid_alloc_GRID_all_grids(const char * filename) {
+  util_abort("%s .GRID files - %s - not supported \n", __func__ , filename);
+  return NULL;
+}
 
-ecl_grid_type * ecl_grid_alloc_EGRID(const char * grid_file, bool apply_mapaxes) {
+
+static ecl_grid_type * ecl_grid_alloc_EGRID_all_grids(const char * grid_file, bool apply_mapaxes, const int * ext_actnum) {
   ecl_file_enum   file_type;
   file_type = ecl_util_get_file_type(grid_file , NULL , NULL);
   if (file_type != ECL_EGRID_FILE)
@@ -3048,11 +3066,13 @@ ecl_grid_type * ecl_grid_alloc_EGRID(const char * grid_file, bool apply_mapaxes)
     ecl_file_type * ecl_file   = ecl_file_open( grid_file , 0);
     if (ecl_file) {
       int num_grid               = ecl_file_get_num_named_kw( ecl_file , GRIDHEAD_KW );
-      ecl_grid_type * main_grid  = ecl_grid_alloc_EGRID__( NULL , ecl_file , 0 , apply_mapaxes);
+      ecl_grid_type * main_grid  = ecl_grid_alloc_EGRID__( NULL , ecl_file , 0 , apply_mapaxes, ext_actnum );
       int grid_nr;
 
       for ( grid_nr = 1; grid_nr < num_grid; grid_nr++) {
-        ecl_grid_type * lgr_grid = ecl_grid_alloc_EGRID__( main_grid , ecl_file , grid_nr , false);  /* The apply_mapaxes argument is ignored for LGR - it inherits from parent anyway. */
+        // The apply_mapaxes argument is ignored for LGR - 
+        //   it inherits from parent anyway.
+        ecl_grid_type * lgr_grid = ecl_grid_alloc_EGRID__( main_grid , ecl_file , grid_nr , false, NULL );
         ecl_grid_add_lgr( main_grid , lgr_grid );
         {
           ecl_grid_type * host_grid;
@@ -3077,6 +3097,11 @@ ecl_grid_type * ecl_grid_alloc_EGRID(const char * grid_file, bool apply_mapaxes)
 }
 
 
+ecl_grid_type * ecl_grid_alloc_EGRID(const char * grid_file, bool apply_mapaxes) {
+  return ecl_grid_alloc_EGRID_all_grids(grid_file, apply_mapaxes, NULL);
+}
+
+
 
 
 
@@ -3607,6 +3632,23 @@ ecl_grid_type * ecl_grid_alloc(const char * grid_file ) {
 }
 
 
+// This function is used to override use of the keyword ACTNUM from the EGRID file.
+// ext_actnum must have size equal to the number of cells in the main grid
+// if ext_actnum = NULL, actnum is taken from file, otherwise ext_actnums
+// determines which cells are active.
+ecl_grid_type * ecl_grid_alloc_ext_actnum(const char * grid_file, const int * ext_actnum) {
+  ecl_file_enum file_type = ecl_util_get_file_type(grid_file , NULL ,  NULL);
+  if (file_type == ECL_EGRID_FILE)
+    return ecl_grid_alloc_EGRID_all_grids(grid_file, true, ext_actnum);
+  else if (file_type == ECL_GRID_FILE)
+    ecl_grid_alloc_GRID_all_grids(grid_file);
+  else
+    util_abort("%s must have .EGRID file - %s not recognized \n", __func__ , grid_file);
+
+  return NULL;
+}
+
+
 static void ecl_grid_file_nactive_dims( fortio_type * data_fortio , int * dims) {
   if (data_fortio) {
     if (ecl_kw_fseek_kw( INTEHEAD_KW , false , false , data_fortio )) {

lib/ecl/tests/ecl_grid_ext_actnum.cpp

@@ -0,0 +1,64 @@
+#include <vector>
+
+#include <ert/util/test_work_area.hpp>
+#include <ert/util/test_util.hpp>
+
+#include <ert/ecl/ecl_endian_flip.hpp>
+#include <ert/ecl/ecl_kw.hpp>
+#include <ert/ecl/ecl_kw_magic.hpp>
+#include <ert/ecl/ecl_file.hpp>
+#include <ert/ecl/ecl_grid.hpp>
+
+
+void test_1() {
+
+  test_work_area_type * work_area = test_work_area_alloc("ext_actnum_main_grid");
+  {
+    const char * filename = "FILE.EGRID";
+
+    ecl_grid_type * grid_write = ecl_grid_alloc_rectangular(2,3,1,   1,1,1,NULL);
+    ecl_grid_fwrite_EGRID( grid_write , filename, true);
+    ecl_grid_free(grid_write);
+
+    ecl_file_type * ecl_file = ecl_file_open(filename, 0);
+    ecl_kw_type * filehead_kw  = ecl_file_iget_named_kw( ecl_file , FILEHEAD_KW  , 0);
+    int * filehead_data = ecl_kw_get_int_ptr(filehead_kw);
+    filehead_data[FILEHEAD_DUALP_INDEX] = FILEHEAD_DUAL_POROSITY;
+
+    ecl_kw_type * actnum_kw = ecl_file_iget_named_kw( ecl_file, ACTNUM_KW, 0);
+    int * actnum_data = ecl_kw_get_int_ptr(actnum_kw);
+    actnum_data[0] = 1;
+    actnum_data[1] = 2;
+    actnum_data[2] = 2;
+    actnum_data[3] = 0;
+    actnum_data[4] = 1;
+    actnum_data[5] = 1;
+    const char * filename1 = "FILE1.EGRID";
+    fortio_type * f = fortio_open_writer( filename1, false, ECL_ENDIAN_FLIP );
+    ecl_file_fwrite_fortio(ecl_file, f, 0);
+    fortio_fclose( f );
+    ecl_file_close(ecl_file);
+
+    std::vector<int> ext_actnum = {0, 1, 0, 1, 1, 1};
+    ecl_grid_type * grid = ecl_grid_alloc_ext_actnum(filename1, ext_actnum.data());
+    test_assert_int_equal( 2, ecl_grid_get_nactive(grid) );
+    test_assert_int_equal( 1, ecl_grid_get_nactive_fracture(grid) );
+    test_assert_true( !ecl_grid_cell_active1(grid, 0) );
+
+    test_assert_true( !ecl_grid_cell_active1(grid, 2) );
+    test_assert_true( !ecl_grid_cell_active1(grid, 3) );
+    test_assert_true(  ecl_grid_cell_active1(grid, 4) );
+    test_assert_true(  ecl_grid_cell_active1(grid, 5) );
+
+    ecl_grid_free( grid );
+
+  }
+  test_work_area_free( work_area );
+
+
+}
+
+
+int main() {
+  test_1();
+}

lib/include/ert/ecl/ecl_grid.hpp

@@ -112,6 +112,7 @@ extern "C" {
   ecl_grid_type * ecl_grid_alloc_GRDECL_data(int , int , int , const float *  , const float *  , const int * , bool apply_mapaxes , const float * mapaxes);
   ecl_grid_type * ecl_grid_alloc_GRID_data(int num_coords , int nx, int ny , int nz , int coords_size , int ** coords , float ** corners , bool apply_mapaxes, const float * mapaxes);
   ecl_grid_type * ecl_grid_alloc(const char * );
+  ecl_grid_type * ecl_grid_alloc_ext_actnum(const char * , const int * ext_actnum);
   ecl_grid_type * ecl_grid_load_case( const char * case_input );
   ecl_grid_type * ecl_grid_load_case__( const char * case_input , bool apply_mapaxes);
   ecl_grid_type * ecl_grid_alloc_rectangular( int nx , int ny , int nz , double dx , double dy , double dz , const int * actnum);