diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt index ccd1a0468..beee85537 100644 --- a/lib/CMakeLists.txt +++ b/lib/CMakeLists.txt @@ -87,7 +87,6 @@ add_library( ecl/ecl_file_kw.cpp ecl/ecl_file_view.cpp ecl/ecl_grav.cpp - ecl/ecl_grav_calc.cpp ecl/ecl_smspec.cpp ecl/ecl_unsmry_loader.cpp ecl/ecl_sum_data.cpp diff --git a/lib/ecl/ecl_grav_calc.cpp b/lib/ecl/ecl_grav_calc.cpp deleted file mode 100644 index 1a3c50c39..000000000 --- a/lib/ecl/ecl_grav_calc.cpp +++ /dev/null @@ -1,136 +0,0 @@ -#include -#include -#include - -#include - -#include -#include -#include -#include - -/** - This file contains one function, ecl_grav_phase_deltag() which - calculates the change in local gravitational strength (in units of - micro Gal) in a point, based on base and monitor values for pore - volume, saturation and density. Should typically be called several - times, one time for each phase. -*/ - -/** - Check that the rporv values are in the right ballpark. For - ECLIPSE version 2008.2 they are way off. Check PORV - versus RPORV for ten 'random' locations in the grid. -*/ - -static void ecl_grav_check_rporv(const ecl_grid_type *ecl_grid, - const ecl_kw_type *rporv1_kw, - const ecl_kw_type *rporv2_kw, - const ecl_kw_type *init_porv_kw) { - int active_index; - int active_delta; - int active_size; - - ecl_grid_get_dims(ecl_grid, NULL, NULL, NULL, &active_size); - active_delta = active_size / 12; - for (active_index = active_delta; active_index < active_size; - active_index += active_delta) { - int global_index = ecl_grid_get_global_index1A(ecl_grid, active_index); - double init_porv = ecl_kw_iget_as_double( - init_porv_kw, global_index); /* NB - this uses global indexing. */ - double rporv1 = ecl_kw_iget_as_double(rporv1_kw, active_index); - double rporv2 = ecl_kw_iget_as_double(rporv2_kw, active_index); - double rporv12 = 0.5 * (rporv1 + rporv2); - double fraction = util_double_min(init_porv, rporv12) / - util_double_max(init_porv, rporv12); - - if (fraction < 0.50) { - fprintf(stderr, "--------------------------------------------------" - "---------------\n"); - fprintf(stderr, "INIT PORV: %g \n", init_porv); - fprintf(stderr, "RPORV1 : %g \n", rporv1); - fprintf(stderr, "RPORV2 : %g \n", rporv2); - fprintf(stderr, "Hmmm - the RPORV values extracted from the " - "restart file seem to be \n"); - fprintf(stderr, "veeery different from the initial rporv value. " - "This might indicate\n"); - fprintf(stderr, "an ECLIPSE bug. Version 2007.2 is known to be ok " - "in this respect, \n"); - fprintf(stderr, - "whereas version 2008.2 is known to have a bug. \n"); - fprintf(stderr, "--------------------------------------------------" - "---------------\n"); - exit(1); - } - } -} - -double -ecl_grav_phase_deltag(double utm_x, double utm_y, double tvd, - const ecl_grid_type *grid, const ecl_file_type *init_file, - const ecl_kw_type *sat1_kw, const ecl_kw_type *rho1_kw, - const ecl_kw_type *porv1_kw, const ecl_kw_type *sat2_kw, - const ecl_kw_type *rho2_kw, const ecl_kw_type *porv2_kw) { - - double deltag = 0; - const int *aquifern = NULL; - - const float *rho1 = ecl_kw_get_float_ptr(rho1_kw); - const float *rho2 = ecl_kw_get_float_ptr(rho2_kw); - const float *sat1 = ecl_kw_get_float_ptr(sat1_kw); - const float *sat2 = ecl_kw_get_float_ptr(sat2_kw); - const float *porv1 = ecl_kw_get_float_ptr(porv1_kw); - const float *porv2 = ecl_kw_get_float_ptr(porv2_kw); - - if (ecl_file_has_kw(init_file, "AQUIFERN")) { - const ecl_kw_type *aquifern_kw = - ecl_file_iget_named_kw(init_file, "AQUIFERN", 0); - aquifern = ecl_kw_get_int_ptr(aquifern_kw); - } - - { - const ecl_kw_type *init_porv_kw = - ecl_file_iget_named_kw(init_file, "PORV", 0); - ecl_grav_check_rporv(grid, porv1_kw, porv2_kw, init_porv_kw); - } - - { - int active_index; - for (active_index = 0; active_index < ecl_grid_get_active_size(grid); - active_index++) { - if (aquifern != NULL && aquifern[active_index] != 0) - continue; /* This is a numerical aquifer cell - skip it. */ - else { - double mas1, mas2; - double xpos, ypos, zpos; - - mas1 = rho1[active_index] * porv1[active_index] * - sat1[active_index]; - mas2 = rho2[active_index] * porv2[active_index] * - sat2[active_index]; - - ecl_grid_get_xyz1A(grid, active_index, &xpos, &ypos, &zpos); - { - double dist_x = xpos - utm_x; - double dist_y = ypos - utm_y; - double dist_z = zpos - tvd; - double dist_sq = - dist_x * dist_x + dist_y * dist_y + dist_z * dist_z; - - if (dist_sq == 0) - exit(1); - - /** - The Gravitational constant is 6.67E-11 N (m/kg)^2, we - return the result in microGal, i.e. we scale with 10^2 * - 10^6 => 6.67E-3. - */ - deltag += - 6.67E-3 * (mas2 - mas1) * dist_z / pow(dist_sq, 1.5); - } - } - } - } - - return deltag; -} diff --git a/lib/include/ert/ecl/ecl_grav_calc.h b/lib/include/ert/ecl/ecl_grav_calc.h deleted file mode 100644 index 04a35749b..000000000 --- a/lib/include/ert/ecl/ecl_grav_calc.h +++ /dev/null @@ -1,7 +0,0 @@ -/* - Warning: The libecl code has changed to be compiled as a C++ project. This - header file is retained for a period for compatibility, but you are encouraged - to switch to include the new hpp header directly in your code. -*/ - -#include diff --git a/lib/include/ert/ecl/ecl_grav_calc.hpp b/lib/include/ert/ecl/ecl_grav_calc.hpp deleted file mode 100644 index 95ec01843..000000000 --- a/lib/include/ert/ecl/ecl_grav_calc.hpp +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef ERT_ECL_GRAV_CALC_H -#define ERT_ECL_GRAV_CALC_H -#ifdef __cplusplus -extern "C" { -#endif - -#include -#include -#include - -double -ecl_grav_phase_deltag(double utm_x, double utm_y, double tvd, - const ecl_grid_type *grid, const ecl_file_type *init_file, - const ecl_kw_type *sat_kw1, const ecl_kw_type *rho_kw1, - const ecl_kw_type *porv_kw1, const ecl_kw_type *sat_kw2, - const ecl_kw_type *rho_kw2, const ecl_kw_type *porv_kw2); - -#ifdef __cplusplus -} -#endif -#endif diff --git a/python/ecl/gravimetry/__init__.py b/python/ecl/gravimetry/__init__.py index fbacda4c2..e7532f14d 100644 --- a/python/ecl/gravimetry/__init__.py +++ b/python/ecl/gravimetry/__init__.py @@ -7,6 +7,5 @@ functionality. """ -from .ecl_subsidence import EclSubsidence -from .ecl_grav_calc import phase_deltag, deltag from .ecl_grav import EclGrav +from .ecl_subsidence import EclSubsidence diff --git a/python/ecl/gravimetry/ecl_grav_calc.py b/python/ecl/gravimetry/ecl_grav_calc.py deleted file mode 100644 index 418a9e2c4..000000000 --- a/python/ecl/gravimetry/ecl_grav_calc.py +++ /dev/null @@ -1,79 +0,0 @@ -from ecl import EclPrototype - -__arglist = "double, double, double, " -__arglist += "ecl_grid, ecl_file, " -__arglist += "ecl_kw, ecl_kw, ecl_kw, ecl_kw, ecl_kw, ecl_kw" -_phase_deltag = EclPrototype("double ecl_grav_phase_deltag(%s)" % __arglist) - - -def phase_deltag(xyz, grid, init, sat1, rho1, porv1, sat2, rho2, porv2): - return _phase_deltag( - xyz[0], - xyz[1], - xyz[2], - grid.c_ptr, - init.c_ptr, - sat1.c_ptr, - rho1.c_ptr, - porv1.c_ptr, - sat2.c_ptr, - rho2.c_ptr, - porv2.c_ptr, - ) - - -def deltag(xyz, grid, init_file, restart_file1, restart_file2): - """ - 1. All restart files should have water, i.e. the SWAT keyword. - 2. All phases present in the restart file should also be present as densities, - in addition the model must contain one additional phase - which should have a density. - 3. The restart files can never contain oil saturation. - """ - - swat1 = restart_file1.iget_named_kw("SWAT", 0) - swat2 = restart_file2.iget_named_kw("SWAT", 0) - - phase_list = [(swat1, swat2)] - - if restart_file1.has_kw("SGAS"): - # This is a three phase Water / Gas / Oil system - sgas1 = restart_file1.iget_named_kw("SGAS", 0) - sgas2 = restart_file2.iget_named_kw("SGAS", 0) - - soil1 = 1 - (sgas1 + swat1) - soil2 = 1 - (sgas2 + swat2) - soil1.name = "SOIL" - soil2.name = "SOIL" - phase_list += [(sgas1, sgas2), (soil1, soil2)] - else: - # This is a two phase Water / xxx System. Will look for - # OIL_DEN and GAS_DEN keywords to determine whether it is a - # Water / Oil or Water / Gas system. - - if restart_file1.has_kw("OIL_DEN"): - # Oil / Water system - soil1 = 1 - swat1 - soil2 = 1 - swat2 - soil1.name = "SOIL" - soil2.name = "SOIL" - phase_list += [(soil1, soil2)] - else: - # Gas / Water system - sgas1 = 1 - swat1 - sgas2 = 1 - swat2 - sgas1.name = "SGAS" - sgas2.name = "SGAS" - phase_list += [(sgas1, sgas2)] - - porv1 = restart_file1.iget_named_kw("RPORV", 0) - porv2 = restart_file2.iget_named_kw("RPORV", 0) - - deltag = 0 - for sat1, sat2 in phase_list: - rho_name = "%s_DEN" % sat1.name[1:] - rho1 = restart_file1.iget_named_kw(rho_name, 0) - rho2 = restart_file2.iget_named_kw(rho_name, 0) - deltag += phase_deltag( - xyz, grid, init_file, sat1, rho1, porv1, sat2, rho2, porv2 - ) - return deltag