OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Updated ERT to 946e512ac17c1e2469892072a925428c0a012fa1

Browse Source
This commit is contained in:
Magne Sjaastad
2014-10-09 20:13:04 +02:00
parent a8ce90e404
commit b93053b089
1188 changed files with 49680 additions and 11622 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
ThirdParty/Ert/devel/libanalysis/CMakeLists.txt vendored
View File

@@ -4,11 +4,11 @@ add_subdirectory( script )
add_subdirectory( src )
add_subdirectory( modules )
if (BUILD_TESTS)
add_subdirectory( tests )
endif()
if (BUILD_APPLICATIONS)
add_subdirectory( applications )
endif()
if (BUILD_TESTS)
add_subdirectory( tests )
endif()

14
ThirdParty/Ert/devel/libanalysis/applications/CMakeLists.txt vendored
View File

@@ -5,4 +5,16 @@ if (USE_RUNPATH)
add_runpath( ert_module_test )
endif()
install(TARGETS ert_module_test DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
set (destination ${CMAKE_INSTALL_PREFIX}/bin)
install(TARGETS ert_module_test DESTINATION ${destination})
if (INSTALL_GROUP)
install(CODE "EXECUTE_PROCESS(COMMAND chgrp ${INSTALL_GROUP} ${destination}/ert_module_test)")
install(CODE "EXECUTE_PROCESS(COMMAND chmod g+w ${destination}/ert_module_test)")
endif()
if (BUILD_TESTS)
ert_module_name( VAR_RML rml_enkf ${LIBRARY_OUTPUT_PATH} )
add_test( analysis_module_test_RML ${EXECUTABLE_OUTPUT_PATH}/ert_module_test ${VAR_RML})
endif()

14
ThirdParty/Ert/devel/libanalysis/include/ert/analysis/analysis_module.h vendored
View File

@@ -24,6 +24,7 @@ extern "C" {
#include <ert/util/type_macros.h>
#include <ert/util/matrix.h>
#include <ert/util/bool_vector.h>
/*
@@ -94,13 +95,15 @@ typedef enum {
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D );
void analysis_module_init_update( analysis_module_type * module ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D );
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
const matrix_type * D );
const char * analysis_module_get_lib_name( const analysis_module_type * module);
@@ -114,6 +117,7 @@ typedef enum {
bool analysis_module_has_var( const analysis_module_type * module , const char * var );
double analysis_module_get_double( const analysis_module_type * module , const char * var);
int analysis_module_get_int( const analysis_module_type * module , const char * var);
bool analysis_module_get_bool( const analysis_module_type * module , const char * var);
void * analysis_module_get_ptr( const analysis_module_type * module , const char * var);
const char * analysis_module_flag_enum_iget( int index, int * value);

5
ThirdParty/Ert/devel/libanalysis/include/ert/analysis/analysis_table.h vendored
View File

@@ -8,7 +8,7 @@ extern "C" {
#include <ert/util/matrix.h>
#include <ert/util/rng.h>
#include <ert/util/bool_vector.h>
typedef void (analysis_updateA_ftype) (void * module_data ,
matrix_type * A ,
@@ -38,6 +38,7 @@ extern "C" {
typedef void (analysis_init_update_ftype) (void * module_data,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
@@ -51,6 +52,7 @@ extern "C" {
typedef bool (analysis_has_var_ftype) (const void * module_data , const char * var_name);
typedef int (analysis_get_int_ftype) (const void * module_data , const char * var_name );
typedef double (analysis_get_double_ftype) (const void * module_data , const char * var_name );
typedef bool (analysis_get_bool_ftype) (const void * module_data , const char * var_name );
typedef void * (analysis_get_ptr_ftype) (const void * module_data , const char * var_name );
/*****************************************************************/
@@ -77,6 +79,7 @@ typedef struct {
analysis_has_var_ftype * has_var;
analysis_get_int_ftype * get_int;
analysis_get_double_ftype * get_double;
analysis_get_bool_ftype * get_bool;
analysis_get_ptr_ftype * get_ptr;
} analysis_table_type;

2
ThirdParty/Ert/devel/libanalysis/include/ert/analysis/cv_enkf.h vendored
View File

@@ -18,6 +18,7 @@
#include <ert/util/rng.h>
#include <ert/util/matrix.h>
#include <ert/util/bool_vector.h>
typedef struct cv_enkf_data_struct cv_enkf_data_type;
@@ -25,6 +26,7 @@ void * cv_enkf_data_alloc( rng_type * rng );
void cv_enkf_data_free( void * arg );
void cv_enkf_init_update( void * arg ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,

11
ThirdParty/Ert/devel/libanalysis/include/ert/analysis/enkf_linalg.h vendored
View File

@@ -3,6 +3,7 @@
#include <ert/util/matrix_lapack.h>
#include <ert/util/matrix.h>
#include <ert/util/double_vector.h>
void enkf_linalg_get_PC( const matrix_type * S0,
@@ -10,7 +11,8 @@ void enkf_linalg_get_PC( const matrix_type * S0,
double truncation,
int ncomp,
matrix_type * PC,
matrix_type * PC_obs );
matrix_type * PC_obs ,
double_vector_type * singular_values);
void enkf_linalg_init_stdX( matrix_type * X ,
@@ -96,16 +98,11 @@ void enkf_linalg_checkX(const matrix_type * X , bool bootstrap);
void enkf_linalg_rml_enkfX1(matrix_type *X1, matrix_type * Udr ,matrix_type * S ,matrix_type *R);
void enkf_linalg_rml_enkfX2(matrix_type *X2, double *Wdr, matrix_type * X1 ,double a , int nsign);
void enkf_linalg_rml_enkfX3(matrix_type *X3, matrix_type *VdTr, double *Wdr,matrix_type *X2, int nsign);
void enkf_linalg_rml_enkfdA(matrix_type *dA1,matrix_type *Dm,matrix_type *X3);
double enkf_linalg_data_mismatch(matrix_type *D , matrix_type *R , matrix_type *Sk);
void enkf_linalg_Covariance(matrix_type *Cd, const matrix_type *E, double nsc ,int nrobs);
void enkf_linalg_rml_enkfAm(matrix_type * Um, double * Wm,int nsign1);
void enkf_linalg_rml_enkfAm(matrix_type * Um, const double * Wm,int nsign1);
void enkf_linalg_rml_enkfX4 (matrix_type *X4,matrix_type *Am, matrix_type *A);
void enkf_linalg_rml_enkfX5(matrix_type * X5,matrix_type * Am, matrix_type * X4);
void enkf_linalg_rml_enkfX6(matrix_type * X6, matrix_type * Dk, matrix_type * X5);
void enkf_linalg_rml_enkfX7(matrix_type * X7, matrix_type * VdT, double * Wdr, double a,matrix_type * X6);
void enkf_linalg_rml_enkfXdA2(matrix_type * dA2,matrix_type * Dk,matrix_type * X7);
#endif

6
ThirdParty/Ert/devel/libanalysis/include/ert/analysis/std_enkf.h vendored
View File

@@ -13,8 +13,6 @@ extern "C" {
#define DEFAULT_ENKF_TRUNCATION_ 0.98
#define ENKF_TRUNCATION_KEY_ "ENKF_TRUNCATION"
#define ENKF_NCOMP_KEY_ "ENKF_NCOMP"
#define ENKF_LAMBDA0_KEY_ "LAMBDA0"
#define ENKF_ITER_KEY_ "ITER"
typedef struct std_enkf_data_struct std_enkf_data_type;
@@ -27,7 +25,9 @@ extern "C" {
void std_enkf_set_truncation( std_enkf_data_type * data , double truncation );
void std_enkf_set_subspace_dimension( std_enkf_data_type * data , int subspace_dimension);
void std_enkf_set_lambda0( std_enkf_data_type * data , double lambda0 );
bool std_enkf_has_var( const void * arg, const char * var_name);
int std_enkf_get_int( const void * arg, const char * var_name);
double std_enkf_get_double( const void * arg, const char * var_name);
double std_enkf_get_truncation( std_enkf_data_type * data );
void * std_enkf_data_alloc( rng_type * rng);

47
ThirdParty/Ert/devel/libanalysis/modules/CMakeLists.txt vendored
View File

@@ -1,12 +1,43 @@
set( args "--silent --exclude-ert -I${PROJECT_SOURCE_DIR}/libanalysis/include -I${PROJECT_SOURCE_DIR}/libert_util/include -I${CMAKE_CURRENT_SOURCE_DIR} -I${PROJECT_BINARY_DIR}/libert_util/include")
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
set( RML_SOURCE_FILES
rml_enkf.c
set( RML_SOURCE_FILES
rml_enkf.c
rml_enkf_common.c )
set( RMLI_SOURCE_FILES
rml_enkf_imodel.c
rml_enkf_common.c )
set( header_files analysis_module.h enkf_linalg.h analysis_table.h std_enkf.h rml_enkf_common.h)
add_library( rml_enkf SHARED ${RML_SOURCE_FILES} )
set_target_properties( rml_enkf PROPERTIES VERSION 1.0 SOVERSION 1.0 PREFIX "")
target_link_libraries( rml_enkf analysis )
target_link_libraries( rml_enkf dl )
if (USE_RUNPATH)
add_runpath( rml_enkf )
endif()
if (BUILD_TESTS)
add_subdirectory( tests )
endif()
if (INSTALL_ERT)
install(TARGETS rml_enkf DESTINATION ${CMAKE_INSTALL_LIBDIR})
endif()
#-----------------------------------------------------------------
# Alternative script based build:
#if (BUILD_TESTS)
# if (BUILD_APPLICATIONS)
#set( args "--silent --exclude-ert -I${PROJECT_SOURCE_DIR}/libanalysis/include -I${PROJECT_SOURCE_DIR}/libert_util/include -I${CMAKE_CURRENT_SOURCE_DIR} -I${PROJECT_BINARY_DIR}/libert_util/include")
#set( RML_SOURCE_FILES
# rml_enkf.c
# rml_enkf_common.c )
#ert_module( ${LIBRARY_OUTPUT_PATH}/rml_enkf ${args} "${RML_SOURCE_FILES}")
ert_module( ${LIBRARY_OUTPUT_PATH}/rml_enkf ${args} "${RML_SOURCE_FILES}")
ert_module( ${LIBRARY_OUTPUT_PATH}/rmli_enkf ${args} "${RMLI_SOURCE_FILES}")

427
ThirdParty/Ert/devel/libanalysis/modules/deprecated/rml_enkf.c vendored Normal file
View File

@@ -0,0 +1,427 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
The file 'rml_enkf.c' is part of ERT - Ensemble based Reservoir Tool.
ERT is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ERT is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
for more details.
*/
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <ert/util/type_macros.h>
#include <ert/util/util.h>
#include <ert/util/rng.h>
#include <ert/util/matrix.h>
#include <ert/util/matrix_blas.h>
#include <ert/util/bool_vector.h>
#include <ert/analysis/analysis_module.h>
#include <ert/analysis/analysis_table.h>
#include <ert/analysis/enkf_linalg.h>
#include <ert/analysis/std_enkf.h>
#include <rml_enkf_common.h>
/*
A random 'magic' integer id which is used for run-time type checking
of the input data.
*/
#define RML_ENKF_TYPE_ID 261123
/*
Observe that only one of the settings subspace_dimension and
truncation can be valid at a time; otherwise the svd routine will
fail. This implies that the set_truncation() and
set_subspace_dimension() routines will set one variable, AND
INVALIDATE THE OTHER. For most situations this will be OK, but if
you have repeated calls to both of these functions the end result
might be a surprise.
*/
#define INVALID_SUBSPACE_DIMENSION -1
#define INVALID_TRUNCATION -1
#define DEFAULT_SUBSPACE_DIMENSION INVALID_SUBSPACE_DIMENSION
/*
The configuration data used by the rml_enkf module is contained in a
rml_enkf_data_struct instance. The data type used for the rml_enkf
module is quite simple; with only a few scalar variables, but there
are essentially no limits to what you can pack into such a datatype.
All the functions in the module have a void pointer as the first
argument, this will immediately be casted to a rml_enkf_data_type
instance, to get some type safety the UTIL_TYPE_ID system should be
used (see documentation in util.h)
The data structure holding the data for your analysis module should
be created and initialized by a constructor, which should be
registered with the '.alloc' element of the analysis table; in the
same manner the desctruction of this data should be handled by a
destructor or free() function registered with the .freef field of
the analysis table.
*/
typedef struct rml_enkf_data_struct rml_enkf_data_type;
struct rml_enkf_data_struct {
UTIL_TYPE_ID_DECLARATION;
double truncation; // Controlled by config key: ENKF_TRUNCATION_KEY
int subspace_dimension; // Controlled by config key: ENKF_NCOMP_KEY (-1: use Truncation instead)
long option_flags;
int iteration_nr; // Keep track of the outer iteration loop
double lambda; // parameter to control the search direction in Marquardt levenberg optimization
double lambda0; // Initial lambda value
double Sk; // Objective function value
double Std; // Standard Deviation of the Objective function
matrix_type *state;
bool_vector_type * ens_mask;
};
/*
This is a macro which will expand to generate a function:
rml_enkf_data_type * rml_enkf_data_safe_cast( void * arg ) {}
which is used for runtime type checking of all the functions which
accept a void pointer as first argument.
*/
static UTIL_SAFE_CAST_FUNCTION( rml_enkf_data , RML_ENKF_TYPE_ID )
static UTIL_SAFE_CAST_FUNCTION_CONST( rml_enkf_data , RML_ENKF_TYPE_ID )
double rml_enkf_get_truncation( rml_enkf_data_type * data ) {
return data->truncation;
}
int rml_enkf_get_subspace_dimension( rml_enkf_data_type * data ) {
return data->subspace_dimension;
}
void rml_enkf_set_truncation( rml_enkf_data_type * data , double truncation ) {
data->truncation = truncation;
if (truncation > 0.0)
data->subspace_dimension = INVALID_SUBSPACE_DIMENSION;
}
void rml_enkf_set_lambda0(rml_enkf_data_type * data , double lambda0 ) {
data->lambda0 = lambda0;
}
void rml_enkf_set_subspace_dimension( rml_enkf_data_type * data , int subspace_dimension) {
data->subspace_dimension = subspace_dimension;
if (subspace_dimension > 0)
data->truncation = INVALID_TRUNCATION;
}
void rml_enkf_set_iteration_number( rml_enkf_data_type *data , int iteration_number ) {
data->iteration_nr = iteration_number;
}
void * rml_enkf_data_alloc( rng_type * rng) {
rml_enkf_data_type * data = util_malloc( sizeof * data );
UTIL_TYPE_ID_INIT( data , RML_ENKF_TYPE_ID );
rml_enkf_set_truncation( data , DEFAULT_ENKF_TRUNCATION_ );
rml_enkf_set_subspace_dimension( data , DEFAULT_SUBSPACE_DIMENSION );
data->option_flags = ANALYSIS_NEED_ED + ANALYSIS_UPDATE_A + ANALYSIS_ITERABLE + ANALYSIS_SCALE_DATA;
data->iteration_nr = 0;
data->Std = 0;
data->state = matrix_alloc(1,1); // This will be resized under use; but we need a valid instance
data->lambda0 = -1.0;
data->ens_mask = bool_vector_alloc(0,false);
return data;
}
void rml_enkf_data_free( void * module_data ) {
rml_enkf_data_type * data = rml_enkf_data_safe_cast( module_data );
matrix_free( data->state );
bool_vector_free(data->ens_mask);
free( data );
}
/*
About the matrix Cd: The matrix Cd is calculated based on the content
of the E input matrix. In the original implementation this matrix was
only calculated in the first iteration, and then reused between subsequent
iterations.
Due to deactivating outliers the number of active observations can change
from one iteration to the next, if the matrix Cd is then reused between
iterations we will get a matrix size mismatch in the linear algebra. In the
current implementation the Cd matrix is recalculated based on the E input
for each iteration.
*/
void rml_enkf_updateA(void * module_data ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D) {
rml_enkf_data_type * data = rml_enkf_data_safe_cast( module_data );
double truncation = data->truncation;
double Sk_new;
double Std_new;
int ens_size = matrix_get_columns( S );
int nrobs = matrix_get_rows( S );
matrix_type * Cd = matrix_alloc( nrobs , nrobs);
double nsc = 1/sqrt(ens_size-1);
matrix_type * Skm = matrix_alloc(matrix_get_columns(D),matrix_get_columns(D));
FILE *fp = util_fopen("rml_enkf_output","a");
int nrmin = util_int_min( ens_size , nrobs);
matrix_type * Ud = matrix_alloc( nrobs , nrmin ); /* Left singular vectors. */
matrix_type * VdT = matrix_alloc( nrmin , ens_size ); /* Right singular vectors. */
double * Wd = util_calloc( nrmin , sizeof * Wd );
Cd = matrix_alloc( nrobs, nrobs );
enkf_linalg_Covariance(Cd ,E ,nsc, nrobs);
matrix_inv(Cd);
if (data->iteration_nr == 0) {
Sk_new = enkf_linalg_data_mismatch(D,Cd,Skm); //Calculate the intitial data mismatch term
Std_new = matrix_diag_std(Skm,Sk_new);
rml_enkf_common_store_state( data->state , A , data->ens_mask );
if (data->lambda0 < 0)
data->lambda = pow(10,floor(log10(Sk_new/(2*nrobs))));
else
data->lambda = data->lambda0;
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lambda,Ud,Wd,VdT);
data->Sk = Sk_new;
data->Std = Std_new;
printf("Prior Objective function value is %5.3f \n", data->Sk);
fprintf(fp,"Iteration number\t Lamda Value \t Current Mean (OB FN) \t Old Mean\t Current Stddev\n");
fprintf(fp, "\n\n");
fprintf(fp,"%d \t\t NA \t %5.5f \t \t %5.5f \n",data->iteration_nr, Sk_new, Std_new);
} else {
Sk_new = enkf_linalg_data_mismatch(D , Cd , Skm); //Calculate the intitial data mismatch term
Std_new= matrix_diag_std(Skm,Sk_new);
printf(" Current Objective function value is %5.3f \n\n",Sk_new);
printf("The old Objective function value is %5.3f \n", data->Sk);
if ((Sk_new< (data->Sk)) && (Std_new< (data->Std)))
{
if ( (1- (Sk_new/data->Sk)) < .0001) // check convergence ** model change norm has to be added in this!!
data-> iteration_nr = 16;
fprintf(fp,"%d \t\t %5.5f \t %5.5f \t %5.5f \t %5.5f \n",data->iteration_nr,data->lambda, Sk_new,data->Sk, Std_new);
data->lambda = data->lambda / 10 ;
data->Std = Std_new;
rml_enkf_common_store_state( data->state , A , data->ens_mask );
data->Sk = Sk_new;
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lambda,Ud,Wd,VdT);
}
else if((Sk_new< (data->Sk)) && (Std_new > (data->Std)))
{
if ( (1- (Sk_new/data->Sk)) < .0001) // check convergence ** model change norm has to be added in this!!
data-> iteration_nr = 16;
fprintf(fp,"%d \t\t %5.5f \t %5.5f \t %5.5f \t %5.5f \n",data->iteration_nr,data->lambda, Sk_new,data->Sk, Std_new);
data->Std=Std_new;
rml_enkf_common_store_state( data->state , A , data->ens_mask );
data->Sk = Sk_new;
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lambda,Ud,Wd,VdT);
}
else {
fprintf(fp,"%d \t\t %5.5f \t %5.5f \t %5.5f \t %5.5f \n",data->iteration_nr,data->lambda, Sk_new,data->Sk, Std_new);
printf("The previous step is rejected !!\n");
data->lambda = data ->lambda * 4;
rml_enkf_common_recover_state( data->state , A , data->ens_mask );
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lambda,Ud,Wd,VdT);
data->iteration_nr--;
}
}
data->iteration_nr++;
// setting the lower bound for lambda
if (data->lambda <.01)
data->lambda= .01;
printf ("The current iteration number is %d \n ", data->iteration_nr);
matrix_free(Cd);
matrix_free(Ud);
matrix_free(VdT);
matrix_free(Skm);
free(Wd);
fclose(fp);
}
void rml_enkf_init_update(void * arg ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
const matrix_type * D ) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
bool_vector_memcpy( module_data->ens_mask , ens_mask );
}
bool rml_enkf_set_double( void * arg , const char * var_name , double value) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , ENKF_TRUNCATION_KEY_) == 0)
rml_enkf_set_truncation( module_data , value );
else if (strcmp( var_name , ENKF_LAMBDA0_KEY_) == 0)
rml_enkf_set_lambda0( module_data , value );
else
name_recognized = false;
return name_recognized;
}
}
bool rml_enkf_set_int( void * arg , const char * var_name , int value) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , ENKF_NCOMP_KEY_) == 0)
rml_enkf_set_subspace_dimension( module_data , value );
else if(strcmp( var_name , ENKF_ITER_KEY_) == 0)
rml_enkf_set_iteration_number( module_data , value );
else
name_recognized = false;
return name_recognized;
}
}
long rml_enkf_get_options( void * arg , long flag ) {
rml_enkf_data_type * module_data = rml_enkf_data_safe_cast( arg );
{
return module_data->option_flags;
}
}
bool rml_enkf_has_var( const void * arg, const char * var_name) {
bool ret = false;
if ((strcmp(var_name , ENKF_ITER_KEY_) == 0) ||
(strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0) ||
(strcmp(var_name , ENKF_LAMBDA0_KEY_) == 0)) {
ret = true;
}
return ret;
}
int rml_enkf_get_int( const void * arg, const char * var_name) {
const rml_enkf_data_type * module_data = rml_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_ITER_KEY_) == 0)
return module_data->iteration_nr;
else
return -1;
}
}
double rml_enkf_get_double( const void * arg, const char * var_name) {
const rml_enkf_data_type * module_data = rml_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0)
return module_data->truncation;
else if (strcmp(var_name , ENKF_LAMBDA0_KEY_) == 0)
return module_data->lambda0;
else
return -1.0;
}
}
/**
gcc -fpic -c <object_file> -I?? <src_file>
gcc -shared -o <lib_file> <object_files>
*/
#ifdef INTERNAL_LINK
#define SYMBOL_TABLE rml_enkf_symbol_table
#else
#define SYMBOL_TABLE EXTERNAL_MODULE_SYMBOL
#endif
analysis_table_type SYMBOL_TABLE = {
.alloc = rml_enkf_data_alloc,
.freef = rml_enkf_data_free,
.set_int = rml_enkf_set_int ,
.set_double = rml_enkf_set_double ,
.set_bool = NULL ,
.set_string = NULL ,
.get_options = rml_enkf_get_options ,
.initX = NULL,
.updateA = rml_enkf_updateA ,
.init_update = rml_enkf_init_update ,
.complete_update = NULL,
.has_var = rml_enkf_has_var,
.get_int = rml_enkf_get_int,
.get_double = rml_enkf_get_double,
.get_ptr = NULL,
};

913
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf.c vendored
View File

File diff suppressed because it is too large Load Diff

115
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_common.c vendored
View File

@@ -1,5 +1,3 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
@@ -34,69 +32,60 @@
#include <rml_enkf_common.h>
/* This program contains common functions to both rml_enkf & rml_enkf_imodel*/
void rml_enkf_common_initA__( matrix_type * A ,
matrix_type * S ,
matrix_type * Cd ,
matrix_type * E ,
matrix_type * D ,
double truncation,
double lamda,
matrix_type * Udr,
double * Wdr,
matrix_type * VdTr) {
int nrobs = matrix_get_rows( S );
int ens_size = matrix_get_columns( S );
double a = lamda + 1;
matrix_type *tmp = matrix_alloc (nrobs, ens_size);
double nsc = 1/sqrt(ens_size-1);
printf("The lamda Value is %5.5f\n",lamda);
printf("The Value of Truncation is %4.2f \n",truncation);
matrix_subtract_row_mean( S ); /* Shift away the mean in the ensemble predictions*/
matrix_inplace_diag_sqrt(Cd);
matrix_dgemm(tmp, Cd, S,false, false, 1.0, 0.0);
matrix_scale(tmp, nsc);
printf("The Scaling of data matrix completed !\n ");
// Explanation
// zzz_enkf_common_store_state( state , A ,ens_mask) assigns A to state. RESIZES state to rows(A)-by-LEN(ens_mask)
// zzz_enkf_common_recover_state(state , A ,ens_mask) assigns state to A. RESIZES A to rows(state)-by-SUM(ens_mask)
// SVD(S) = Ud * Wd * Vd(T)
int nsign = enkf_linalg_svd_truncation(tmp , truncation , -1 , DGESVD_MIN_RETURN , Wdr , Udr , VdTr);
/* After this we only work with the reduced dimension matrices */
printf("The number of siginificant ensembles are %d \n ",nsign);
matrix_type * X1 = matrix_alloc( nsign, ens_size);
matrix_type * X2 = matrix_alloc (nsign, ens_size );
matrix_type * X3 = matrix_alloc (ens_size, ens_size );
// Compute the matrices X1,X2,X3 and dA
enkf_linalg_rml_enkfX1(X1, Udr ,D ,Cd ); //X1 = Ud(T)*Cd(-1/2)*D -- D= -(dk-d0)
enkf_linalg_rml_enkfX2(X2, Wdr ,X1 ,a, nsign); //X2 = ((a*Ipd)+Wd^2)^-1 * X1
matrix_free(X1);
enkf_linalg_rml_enkfX3(X3, VdTr ,Wdr,X2, nsign); //X3 = Vd *Wd*X2
printf("The X3 matrix is computed !\n ");
matrix_type *dA1= matrix_alloc (matrix_get_rows(A), ens_size);
matrix_type * Dm = matrix_alloc_copy( A );
matrix_subtract_row_mean( Dm ); /* Remove the mean from the ensemble of model parameters*/
matrix_scale(Dm, nsc);
enkf_linalg_rml_enkfdA(dA1, Dm, X3); //dA = Dm * X3
matrix_inplace_add(A,dA1); //dA
matrix_free(X3);
matrix_free(Dm);
matrix_free(dA1);
void rml_enkf_common_store_state( matrix_type * state , const matrix_type * A , const bool_vector_type * ens_mask ) {
matrix_resize( state , matrix_get_rows( A ) , bool_vector_size( ens_mask ) , false);
{
const int ens_size = bool_vector_size( ens_mask );
int active_index = 0;
for (int iens = 0; iens < ens_size; iens++) {
if (bool_vector_iget( ens_mask , iens )) {
matrix_copy_column( state , A , iens , active_index );
active_index++;
} else
matrix_set_const_column( state , iens , 0);
}
}
}
void rml_enkf_common_recover_state( const matrix_type * state , matrix_type * A , const bool_vector_type * ens_mask ) {
const int ens_size = bool_vector_size( ens_mask );
const int active_size = bool_vector_count_equal( ens_mask , true );
const int rows = matrix_get_rows( state );
matrix_resize( A , rows , active_size , false );
{
int active_index = 0;
for (int iens = 0; iens < ens_size; iens++) {
if (bool_vector_iget( ens_mask , iens )) {
matrix_copy_column( A , state , active_index , iens );
active_index++;
}
}
}
}
// Scale rows by the entries in the vector Csc
void rml_enkf_common_scaleA(matrix_type *A , const double * Csc, bool invert ){
int nrows = matrix_get_rows(A);
if (invert) {
for (int i=0; i< nrows ; i++) {
double sc= 1/Csc[i];
matrix_scale_row(A, i, sc);
}
} else {
for (int i=0; i< nrows ; i++) {
double sc= Csc[i];
matrix_scale_row(A, i, sc);
}
}
}

16
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_common.h vendored
View File

@@ -1,4 +1,3 @@
#ifndef __RML_ENKF_COMMON_H__
#define __RML_ENKF_COMMON_H__
@@ -6,18 +5,11 @@
#include <ert/util/matrix.h>
#include <ert/util/rng.h>
#include <ert/util/bool_vector.h>
void rml_enkf_common_initA__( matrix_type * A ,
matrix_type * S ,
matrix_type * Cd ,
matrix_type * E ,
matrix_type * D ,
double truncation,
double lamda,
matrix_type * Ud,
double * Wd,
matrix_type * VdT);
void rml_enkf_common_store_state( matrix_type * state , const matrix_type * A , const bool_vector_type * ens_mask );
void rml_enkf_common_recover_state( const matrix_type * state , matrix_type * A , const bool_vector_type * ens_mask );
void rml_enkf_common_scaleA(matrix_type *A , const double * Csc, bool invert );
#endif

532
ThirdParty/Ert/devel/libanalysis/modules/rml_enkf_imodel.c vendored
View File

@@ -1,532 +0,0 @@
/*
Copyright (C) 2011 Statoil ASA, Norway.
The file 'rml_enkf_imodel.c' is part of ERT - Ensemble based Reservoir Tool.
ERT is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ERT is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
for more details.
*/
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <ert/util/util.h>
#include <ert/util/type_macros.h>
#include <ert/util/rng.h>
#include <ert/util/matrix.h>
#include <ert/util/matrix_blas.h>
#include <ert/analysis/analysis_module.h>
#include <ert/analysis/analysis_table.h>
#include <ert/analysis/enkf_linalg.h>
#include <ert/analysis/std_enkf.h>
#include <rml_enkf_common.h>
/*
A random 'magic' integer id which is used for run-time type checking
of the input data.
*/
#define RML_ENKF_IMODEL_TYPE_ID 261123
typedef struct rml_enkf_imodel_data_struct rml_enkf_imodel_data_type;
/*
Observe that only one of the settings subspace_dimension and
truncation can be valid at a time; otherwise the svd routine will
fail. This implies that the set_truncation() and
set_subspace_dimension() routines will set one variable, AND
INVALIDATE THE OTHER. For most situations this will be OK, but if
you have repeated calls to both of these functions the end result
might be a surprise.
*/
#define INVALID_SUBSPACE_DIMENSION -1
#define INVALID_TRUNCATION -1
#define DEFAULT_SUBSPACE_DIMENSION INVALID_SUBSPACE_DIMENSION
/*
The configuration data used by the rml_enkf_imodel module is contained in a
rml_enkf_imodel_data_struct instance. The data type used for the rml_enkf_imodel
module is quite simple; with only a few scalar variables, but there
are essentially no limits to what you can pack into such a datatype.
All the functions in the module have a void pointer as the first
argument, this will immediately be casted to a rml_enkf_imodel_data_type
instance, to get some type safety the UTIL_TYPE_ID system should be
used (see documentation in util.h)
The data structure holding the data for your analysis module should
be created and initialized by a constructor, which should be
registered with the '.alloc' element of the analysis table; in the
same manner the desctruction of this data should be handled by a
destructor or free() function registered with the .freef field of
the analysis table.
*/
struct rml_enkf_imodel_data_struct {
UTIL_TYPE_ID_DECLARATION;
double truncation; // Controlled by config key: ENKF_TRUNCATION_KEY
int subspace_dimension; // Controlled by config key: ENKF_NCOMP_KEY (-1: use Truncation instead)
long option_flags;
int iteration_nr; // Keep track of the outer iteration loop
double lamda; // parameter to control the search direction in Marquardt levenberg optimization
double Sk; // Objective function value
double Std; // Standard Deviation of the Objective function
double * Csc;
matrix_type * Am;
matrix_type *prior;
matrix_type *state;
};
/*
This is a macro which will expand to generate a function:
rml_enkf_imodel_data_type * rml_enkf_imodel_data_safe_cast( void * arg ) {}
which is used for runtime type checking of all the functions which
accept a void pointer as first argument.
*/
static UTIL_SAFE_CAST_FUNCTION( rml_enkf_imodel_data , RML_ENKF_IMODEL_TYPE_ID )
static UTIL_SAFE_CAST_FUNCTION_CONST( rml_enkf_imodel_data , RML_ENKF_IMODEL_TYPE_ID )
double rml_enkf_imodel_get_truncation( rml_enkf_imodel_data_type * data ) {
return data->truncation;
}
int rml_enkf_imodel_get_subspace_dimension( rml_enkf_imodel_data_type * data ) {
return data->subspace_dimension;
}
void rml_enkf_imodel_set_truncation( rml_enkf_imodel_data_type * data , double truncation ) {
data->truncation = truncation;
if (truncation > 0.0)
data->subspace_dimension = INVALID_SUBSPACE_DIMENSION;
}
void rml_enkf_imodel_set_subspace_dimension( rml_enkf_imodel_data_type * data , int subspace_dimension) {
data->subspace_dimension = subspace_dimension;
if (subspace_dimension > 0)
data->truncation = INVALID_TRUNCATION;
}
void * rml_enkf_imodel_data_alloc( rng_type * rng) {
rml_enkf_imodel_data_type * data = util_malloc( sizeof * data);
UTIL_TYPE_ID_INIT( data , RML_ENKF_IMODEL_TYPE_ID );
rml_enkf_imodel_set_truncation( data , DEFAULT_ENKF_TRUNCATION_ );
rml_enkf_imodel_set_subspace_dimension( data , DEFAULT_SUBSPACE_DIMENSION );
data->option_flags = ANALYSIS_NEED_ED + ANALYSIS_UPDATE_A + ANALYSIS_ITERABLE + ANALYSIS_SCALE_DATA;
data->iteration_nr = 0;
data->Std = 0;
return data;
}
void rml_enkf_imodel_data_free( void * data ) {
free( data );
}
void rml_enkf_imodel_init1__( matrix_type * A,
rml_enkf_imodel_data_type * data,
double truncation,
double nsc) {
int nstate = matrix_get_rows( A );
int ens_size = matrix_get_columns( A );
int nrmin = util_int_min( ens_size , nstate);
matrix_type * Dm = matrix_alloc_copy( A );
matrix_type * Um = matrix_alloc( nstate , nrmin ); /* Left singular vectors. */
matrix_type * VmT = matrix_alloc( nrmin , ens_size ); /* Right singular vectors. */
double * Wm = util_calloc( nrmin , sizeof * Wm );
matrix_subtract_row_mean(Dm);
//This routine only computes the SVD of Ensemble State matrix
for (int i=0; i<nstate; i++){
double sc = nsc/ (data->Csc[i]);
matrix_scale_row (Dm,i, sc);
}
int nsign1 = enkf_linalg_svd_truncation(Dm , truncation , -1 , DGESVD_MIN_RETURN , Wm , Um , VmT);
printf("The significant Eigen values are %d\n",nsign1);
enkf_linalg_rml_enkfAm(Um, Wm, nsign1);
data->Am=matrix_alloc_copy(Um);
printf("\n Init1 completed\n");
matrix_free(Um);
matrix_free(VmT);
matrix_free(Dm);
free(Wm);
}
void rml_enkf_imodel_Create_Csc(rml_enkf_imodel_data_type * data){
// Create the scaling matrix based on the state vector
int nstate = matrix_get_rows( data->prior );
int ens_size = matrix_get_columns( data->prior );
for (int i=0; i < nstate; i++) {
double sumrow = matrix_get_row_sum(data->prior , i);
double tmp = sumrow / ens_size;
if (abs(tmp)< 1)
data->Csc[i]=0.05;
else
data->Csc[i]= 1;
}
}
void rml_enkf_imodel_scalingA(matrix_type *A, double * Csc, bool invert ){
int nrows = matrix_get_rows(A);
if (invert)
for (int i=0; i< nrows ; i++)
{
double sc= 1/Csc[i];
matrix_scale_row(A, i, sc);
}
else
for (int i=0; i< nrows ; i++)
{
double sc= Csc[i];
matrix_scale_row(A, i, sc);
}
}
void rml_enkf_imodel_init2__( rml_enkf_imodel_data_type * data,
matrix_type *A,
matrix_type *Acopy,
double * Wdr,
double nsc,
matrix_type * VdTr) {
int nstate = matrix_get_rows( Acopy );
int ens_size = matrix_get_columns( Acopy );
matrix_type * Dk = matrix_alloc_copy( Acopy );
double a = data->lamda + 1;
matrix_type *Am= matrix_alloc_copy(data->Am);
matrix_type *Apr= matrix_alloc_copy(data->prior);
double *Csc = util_calloc(nstate , sizeof * Csc );
for (int i=0; i< nstate ; i++)
{
Csc[i]= data->Csc[i];
}
int nsign1= matrix_get_columns(data->Am);
matrix_type * X4 = matrix_alloc(nsign1,ens_size);
matrix_type * X5 = matrix_alloc(nstate,ens_size);
matrix_type * X6 = matrix_alloc(ens_size,ens_size);
matrix_type * X7 = matrix_alloc(ens_size,ens_size);
matrix_type * dA2 = matrix_alloc(nstate, ens_size);
//Compute dA2
printf("\n Starting init 2 \n");
matrix_inplace_sub(Dk, Apr);
rml_enkf_imodel_scalingA(Dk,Csc,true);
enkf_linalg_rml_enkfX4(X4, Am, Dk);
matrix_free(Dk);
enkf_linalg_rml_enkfX5(X5, Am, X4);
printf("\nMatrix X5 computed\n");
matrix_type * Dk1 = matrix_alloc_copy( Acopy );
matrix_subtract_row_mean(Dk1);
rml_enkf_imodel_scalingA(Dk1,Csc,true);
matrix_scale(Dk1,nsc);
enkf_linalg_rml_enkfX6(X6, Dk1,X5);
printf("Matrix X6 computed!\n");
enkf_linalg_rml_enkfX7(X7, VdTr ,Wdr, a, X6);
printf("Matrix X7 computed!\n");
rml_enkf_imodel_scalingA(Dk1,Csc,false);
printf("Matrix Dk1 Scaling done!\n");
enkf_linalg_rml_enkfXdA2(dA2,Dk1,X7);
printf("Matrix dA2 computed!\n");
matrix_inplace_sub(A, dA2);
free(Csc);
matrix_free(Am);
matrix_free(Apr);
matrix_free(X4);
matrix_free(X5);
matrix_free(X6);
matrix_free(X7);
matrix_free(dA2);
matrix_free(Dk1);
}
void rml_enkf_imodel_updateA(void * module_data ,
matrix_type * A ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D) {
rml_enkf_imodel_data_type * data = rml_enkf_imodel_data_safe_cast( module_data );
double truncation = data->truncation;
double Sk_new;
double Std_new;
matrix_type * Skm = matrix_alloc(matrix_get_columns(D),matrix_get_columns(D));
FILE *fp = util_fopen("rml_enkf_imodel_output","a");
int nstate = matrix_get_rows(A);
int ens_size = matrix_get_columns( S );
int nrobs = matrix_get_rows( S );
matrix_type * Cd = matrix_alloc( nrobs, nrobs );
double nsc = 1/sqrt(ens_size-1);
int nrmin = util_int_min( ens_size , nrobs);
matrix_type * Ud = matrix_alloc( nrobs , nrmin ); /* Left singular vectors. */
matrix_type * VdT = matrix_alloc( nrmin , ens_size ); /* Right singular vectors. */
double * Wd = util_calloc( nrmin , sizeof * Wd );
enkf_linalg_Covariance(Cd ,E ,nsc, nrobs);
matrix_inv(Cd);
if (data->iteration_nr == 0) {
Sk_new = enkf_linalg_data_mismatch(D,Cd,Skm); //Calculate the intitial data mismatch term
Std_new= matrix_diag_std(Skm,Sk_new);
data->lamda =pow(10,floor(log10(Sk_new/(2*nrobs))));
data->prior = matrix_alloc_copy(A);
data->state = matrix_alloc_copy(A);
data->Csc = util_calloc(nstate , sizeof * data->Csc);
rml_enkf_imodel_Create_Csc(data);
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lamda,Ud,Wd,VdT);
printf("\n model scaling matrix computed\n");
rml_enkf_imodel_init1__(data->prior, data, truncation, nsc);
data->Sk = Sk_new;
data->Std= Std_new;
printf("Prior Objective function value is %5.3f \n", data->Sk);
fprintf(fp,"Iteration number\t Lamda Value \t Current Mean (OB FN) \t Old Mean\t Current Stddev\n");
fprintf(fp, "\n\n");
fprintf(fp,"%d \t\t NA \t %5.5f \t \t %5.5f \n",data->iteration_nr, Sk_new, Std_new);
}
else
{
matrix_type *Acopy = matrix_alloc_copy (A);
Sk_new = enkf_linalg_data_mismatch(D,Cd,Skm); //Calculate the intitial data mismatch term
Std_new= matrix_diag_std(Skm,Sk_new);
printf(" Current Objective function value is %5.3f \n\n",Sk_new);
printf("The old Objective function value is %5.3f \n", data->Sk);
if ((Sk_new< (data->Sk)) && (Std_new<= (data->Std)))
{
if ( (1- (Sk_new/data->Sk)) < .0001) // check convergence ** model change norm has to be added in this!!
data-> iteration_nr = 16;
fprintf(fp,"%d \t\t %5.5f \t %5.5f \t %5.5f \t %5.5f \n",data->iteration_nr,data->lamda, Sk_new,data->Sk, Std_new);
data->lamda = data->lamda / 10 ;
data->Std = Std_new;
data->state = matrix_alloc_copy(A);
data->Sk = Sk_new;
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lamda,Ud,Wd,VdT);
rml_enkf_imodel_init2__(data,A,Acopy,Wd,nsc,VdT);
}
else if((Sk_new< (data->Sk)) && (Std_new > (data->Std)))
{
if ( (1- (Sk_new/data->Sk)) < .0001) // check convergence ** model change norm has to be added in this!!
data-> iteration_nr = 16;
fprintf(fp,"%d \t\t %5.5f \t %5.5f \t %5.5f \t %5.5f \n",data->iteration_nr,data->lamda, Sk_new,data->Sk, Std_new);
data->lamda = data->lamda;
data->Std=Std_new;
data->state = matrix_alloc_copy(A);
data->Sk = Sk_new;
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lamda,Ud,Wd,VdT);
rml_enkf_imodel_init2__(data,A,Acopy,Wd,nsc,VdT);
}
else {
fprintf(fp,"%d \t\t %5.5f \t %5.5f \t %5.5f \t %5.5f \n",data->iteration_nr,data->lamda, Sk_new,data->Sk, Std_new);
printf("The previous step is rejected !!\n");
data->lamda = data ->lamda * 4;
matrix_assign( A , data->state );
printf("matrix copied \n");
rml_enkf_common_initA__(A,S,Cd,E,D,truncation,data->lamda,Ud,Wd,VdT);
rml_enkf_imodel_init2__(data,A,Acopy,Wd,nsc,VdT);
data->iteration_nr--;
}
matrix_free(Acopy);
}
data->iteration_nr++;
//setting the lower bound for lamda
if (data->lamda <.01)
data->lamda= .01;
printf ("The current iteration number is %d \n ", data->iteration_nr);
// free(data->Csc);
matrix_free(Ud);
matrix_free(VdT);
free(Wd);
matrix_free(Skm);
matrix_free(Cd);
fclose(fp);
}
bool rml_enkf_imodel_set_double( void * arg , const char * var_name , double value) {
rml_enkf_imodel_data_type * module_data = rml_enkf_imodel_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , ENKF_TRUNCATION_KEY_) == 0)
rml_enkf_imodel_set_truncation( module_data , value );
else
name_recognized = false;
return name_recognized;
}
}
bool rml_enkf_imodel_set_int( void * arg , const char * var_name , int value) {
rml_enkf_imodel_data_type * module_data = rml_enkf_imodel_data_safe_cast( arg );
{
bool name_recognized = true;
if (strcmp( var_name , ENKF_NCOMP_KEY_) == 0)
rml_enkf_imodel_set_subspace_dimension( module_data , value );
else
name_recognized = false;
return name_recognized;
}
}
long rml_enkf_imodel_get_options( void * arg , long flag ) {
rml_enkf_imodel_data_type * module_data = rml_enkf_imodel_data_safe_cast( arg );
{
return module_data->option_flags;
}
}
bool rml_enkf_imodel_has_var( const void * arg, const char * var_name) {
{
if (strcmp(var_name , ENKF_ITER_KEY_) == 0)
return true;
else
return false;
}
}
int rml_enkf_imodel_get_int( const void * arg, const char * var_name) {
const rml_enkf_imodel_data_type * module_data = rml_enkf_imodel_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_ITER_KEY_) == 0)
return module_data->iteration_nr;
else
return -1;
}
}
/**
gcc -fpic -c <object_file> -I?? <src_file>
gcc -shared -o <lib_file> <object_files>
*/
#ifdef INTERNAL_LINK
#define SYMBOL_TABLE rml_enkf_imodel_symbol_table
#else
#define SYMBOL_TABLE EXTERNAL_MODULE_SYMBOL
#endif
analysis_table_type SYMBOL_TABLE = {
.alloc = rml_enkf_imodel_data_alloc,
.freef = rml_enkf_imodel_data_free,
.set_int = rml_enkf_imodel_set_int ,
.set_double = rml_enkf_imodel_set_double ,
.set_bool = NULL ,
.set_string = NULL ,
.get_options = rml_enkf_imodel_get_options ,
.initX = NULL,
.updateA = rml_enkf_imodel_updateA ,
.init_update = NULL,
.complete_update = NULL,
.has_var = rml_enkf_imodel_has_var,
.get_int = rml_enkf_imodel_get_int,
.get_double = NULL,
.get_ptr = NULL,
};

5
ThirdParty/Ert/devel/libanalysis/modules/tests/CMakeLists.txt vendored Normal file
View File

@@ -0,0 +1,5 @@
add_executable(analysis_rml_enkf_common analysis_rml_enkf_common.c ../rml_enkf_common.c)
target_link_libraries( analysis_rml_enkf_common analysis util test_util )
add_test( analysis_rml_enkf_common ${EXECUTABLE_OUTPUT_PATH}/analysis_rml_enkf_common )

133
ThirdParty/Ert/devel/libanalysis/modules/tests/analysis_rml_enkf_common.c vendored Normal file
View File

@@ -0,0 +1,133 @@
/*
Copyright (C) 2014 Statoil ASA, Norway.
The file 'analysis_rml_common.c' is part of ERT - Ensemble based Reservoir Tool.
ERT is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
ERT is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
for more details.
*/
#include <ert/util/test_util.h>
#include <ert/util/rng.h>
#include <ert/util/mzran.h>
#include <ert/util/matrix.h>
#include <ert/util/bool_vector.h>
#include <rml_enkf_common.h>
void test_store_recover_state() {
rng_type * rng = rng_alloc( MZRAN , INIT_DEFAULT );
int ens_size = 10;
int active_size = 8;
int rows = 100;
matrix_type * state = matrix_alloc(1,1);
bool_vector_type * ens_mask = bool_vector_alloc(ens_size , false);
matrix_type * A = matrix_alloc( rows , active_size);
matrix_type * A2 = matrix_alloc( rows, active_size );
matrix_type * A3 = matrix_alloc( 1,1 );
for (int i=0; i < active_size; i++)
bool_vector_iset( ens_mask , i + 1 , true );
matrix_random_init(A , rng);
rml_enkf_common_store_state( state , A , ens_mask );
test_assert_int_equal( matrix_get_rows( state ) , rows );
test_assert_int_equal( matrix_get_columns( state ) , ens_size );
{
int g;
int a = 0;
for (g=0; g < ens_size; g++) {
if (bool_vector_iget( ens_mask , g )) {
test_assert_true( matrix_columns_equal( state , g , A , a ));
a++;
}
}
}
rml_enkf_common_recover_state( state , A2 , ens_mask);
rml_enkf_common_recover_state( state , A3 , ens_mask);
test_assert_true( matrix_equal( A , A2 ));
test_assert_true( matrix_equal( A , A3 ));
bool_vector_free( ens_mask );
matrix_free( state );
matrix_free( A );
}
void test_scaleA() {
const int N = 10;
matrix_type * m1 = matrix_alloc(N , N);
matrix_type * m2 = matrix_alloc(N , N);
double * csc = util_calloc( N , sizeof * csc );
rng_type * rng = rng_alloc( MZRAN , INIT_DEFAULT );
matrix_random_init( m1 , rng );
matrix_assign(m2 , m1);
test_assert_true( matrix_equal(m1 , m2));
{
for (int i=0; i < N; i++)
csc[i] = (i + 2);
}
rml_enkf_common_scaleA( m1 , csc , false );
{
int row,col;
for (row = 0; row < N; row++) {
for (col=0; col < N; col++) {
double v1 = matrix_iget(m1 , row , col);
double v2 = matrix_iget(m2 , row , col);
test_assert_double_equal( v1 , v2 * csc[row] );
}
}
}
rml_enkf_common_scaleA( m2 , csc , false );
test_assert_true( matrix_equal(m1 , m2));
rml_enkf_common_scaleA( m2 , csc , true );
{
int row,col;
for (row = 0; row < N; row++) {
for (col=0; col < N; col++) {
double v1 = matrix_iget(m1 , row , col);
double v2 = matrix_iget(m2 , row , col);
test_assert_double_equal( v1 , v2 * csc[row] );
}
}
}
rml_enkf_common_scaleA( m1 , csc , true );
test_assert_true( matrix_equal(m1 , m2));
rng_free( rng );
matrix_free(m1);
matrix_free(m2);
free( csc );
}
int main(int argc , char ** argv) {
test_store_recover_state();
test_scaleA();
exit(0);
}

9
ThirdParty/Ert/devel/libanalysis/script/CMakeLists.txt vendored
View File

@@ -1 +1,8 @@
install(PROGRAMS ert_module DESTINATION ${CMAKE_INSTALL_PREFIX}/bin)
set (destination ${CMAKE_INSTALL_PREFIX}/bin)
install(PROGRAMS ert_module DESTINATION ${destination})
if (INSTALL_GROUP)
install(CODE "EXECUTE_PROCESS(COMMAND chgrp ${INSTALL_GROUP} ${destination}/ert_module)")
install(CODE "EXECUTE_PROCESS(COMMAND chmod g+w ${destination}/ert_module)")
endif()

41
ThirdParty/Ert/devel/libanalysis/script/ert_module vendored
View File

@@ -9,7 +9,7 @@ ert_root = os.path.realpath( os.path.join(os.path.dirname( os.path.realpath( os.
#-----------------------------------------------------------------
default_lib_list = ["analysis" , "ert_util"]
default_define_list = ["HAVE_PTHREAD"]
default_define_list = ["HAVE_PTHREAD" , "COMPILER_GCC"]
CFLAGS = "-std=gnu99 -O2 -Wall -fpic -g"
@@ -169,10 +169,10 @@ Link: %s %s %s
parser = OptionParser( usage )
parser.add_option("--ert-root" , dest="ert_root" , action="store")
parser.add_option("-I" , dest = "include_path_list", action = "append")
parser.add_option("-D" , dest = "define_list" , action = "append")
parser.add_option("-L" , dest = "lib_path_list" , action = "append")
parser.add_option("-l" , dest = "lib_list" , action = "append")
parser.add_option("-I" , dest = "include_path_list", action = "append" , default = [])
parser.add_option("-D" , dest = "define_list" , action = "append" , default = [])
parser.add_option("-L" , dest = "lib_path_list" , action = "append" , default = [])
parser.add_option("-l" , dest = "lib_list" , action = "append" , default = [])
parser.add_option("--exclude-ert" , dest = "exclude_ert" , action="store_true" , default = False)
parser.add_option("--use-rpath" , dest="use_rpath" , action="store_true" , default = False)
parser.add_option("--silent" , dest="silent" , action="store_true" , default = False)
@@ -185,28 +185,21 @@ if options.ert_root:
ert_root = options.ert_root
if options.exclude_ert:
include_path_list = ["./"]
lib_path_list = []
define_list = []
lib_list = []
default_include_path_list = ["./"]
default_lib_path_list = []
default_lib_list = []
else:
include_path_list = ["./" , "%s/include" % ert_root]
lib_path_list = ["%s/lib64" % ert_root]
define_list = default_define_list
lib_list = default_lib_list
default_include_path_list = ["./" , "%s/include" % ert_root]
default_lib_path_list = ["%s/lib64" % ert_root]
default_lib_list = default_lib_list
if options.include_path_list:
include_path_list += options.include_path_list
if options.define_list:
define_list += options.define_list
if options.lib_list:
lib_list += options.lib_list
if options.lib_path_list:
lib_path_list += options.lib_path_list
# What is supplied as commandline options should take presedence, this
# implies that it should be first OR last on the commandline.
include_path_list = options.include_path_list + default_include_path_list
define_list = default_define_list + options.define_list
lib_list = options.lib_list + default_lib_list
lib_path_list = options.lib_path_list + default_lib_path_list
verbose = not options.silent

2
ThirdParty/Ert/devel/libanalysis/src/CMakeLists.txt vendored
View File

@@ -1,6 +1,6 @@
# Common libanalysis library
set( source_files analysis_module.c enkf_linalg.c std_enkf.c sqrt_enkf.c cv_enkf.c bootstrap_enkf.c null_enkf.c fwd_step_enkf.c )
set( header_files analysis_module.h enkf_linalg.h analysis_table.h std_enkf.h)
set( header_files analysis_module.h enkf_linalg.h analysis_table.h std_enkf.h fwd_step_enkf.h)
add_library( analysis SHARED ${source_files} )
set_target_properties( analysis PROPERTIES COMPILE_DEFINITIONS INTERNAL_LINK)
set_target_properties( analysis PROPERTIES VERSION 1.0 SOVERSION 1.0 )

29
ThirdParty/Ert/devel/libanalysis/src/analysis_module.c vendored
View File

@@ -55,6 +55,7 @@ struct analysis_module_struct {
analysis_has_var_ftype * has_var;
analysis_get_int_ftype * get_int;
analysis_get_double_ftype * get_double;
analysis_get_bool_ftype * get_bool;
analysis_get_ptr_ftype * get_ptr;
bool internal;
@@ -81,6 +82,7 @@ static analysis_module_type * analysis_module_alloc_empty( const char * user_nam
module->has_var = NULL;
module->get_int = NULL;
module->get_double = NULL;
module->get_bool = NULL;
module->get_ptr = NULL;
module->alloc = NULL;
@@ -121,6 +123,7 @@ static analysis_module_type * analysis_module_alloc__( rng_type * rng ,
module->has_var = table->has_var;
module->get_int = table->get_int;
module->get_double = table->get_double;
module->get_bool = table->get_bool;
module->get_ptr = table->get_ptr;
if (module->alloc)
@@ -270,13 +273,14 @@ void analysis_module_updateA(analysis_module_type * module ,
void analysis_module_init_update( analysis_module_type * module ,
matrix_type * S ,
matrix_type * R ,
matrix_type * dObs ,
matrix_type * E ,
matrix_type * D ) {
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
const matrix_type * E ,
const matrix_type * D ) {
if (module->init_update != NULL)
module->init_update( module->module_data , S , R , dObs , E , D);
module->init_update( module->module_data , ens_mask , S , R , dObs , E , D);
}
@@ -407,6 +411,19 @@ int analysis_module_get_int( const analysis_module_type * module , const char *
}
bool analysis_module_get_bool( const analysis_module_type * module , const char * var) {
if (analysis_module_has_var( module , var )) {
if (module->get_bool != NULL)
return module->get_bool( module->module_data , var );
else
util_exit("%s: Tried to get bool variable:%s from module:%s - get_int() method not implemented for this module\n" , __func__ , var , module->user_name);
} else
util_exit("%s: Tried to get bool variable:%s from module:%s - module does not support this variable \n" , __func__ , var , module->user_name);
return false;
}
double analysis_module_get_double( const analysis_module_type * module , const char * var) {
if (analysis_module_has_var( module , var )) {
if (module->get_double != NULL)

307
ThirdParty/Ert/devel/libanalysis/src/bootstrap_enkf.c vendored
View File

@@ -56,281 +56,7 @@ typedef struct {
static UTIL_SAFE_CAST_FUNCTION( bootstrap_enkf_data , BOOTSTRAP_ENKF_TYPE_ID )
/*****************************************************************/
//CV: static void lowrankCinv_pre_cv(const matrix_type * S ,
//CV: const matrix_type * R ,
//CV: matrix_type * V0T ,
//CV: matrix_type * Z ,
//CV: double * eig ,
//CV: matrix_type * U0,
//CV: double truncation ,
//CV: int ncomp) {
//CV:
//CV: const int nrobs = matrix_get_rows( S );
//CV: const int nrens = matrix_get_columns( S );
//CV: const int nrmin = util_int_min( nrobs , nrens );
//CV:
//CV: double * inv_sig0 = util_malloc( nrmin * sizeof * inv_sig0 , __func__);
//CV:
//CV: enkf_linalg_svdS( S , truncation , ncomp , DGESVD_MIN_RETURN , inv_sig0 , U0 , V0T);
//CV:
//CV: {
//CV: matrix_type * B = matrix_alloc( nrmin , nrmin );
//CV: enkf_linalg_Cee( B , nrens , R , U0 , inv_sig0); /* B = Xo = (N-1) * Sigma0^(+) * U0'* Cee * U0 * Sigma0^(+') (14.26)*/
//CV: /*USE SVD INSTEAD*/
//CV: matrix_dgesvd(DGESVD_MIN_RETURN , DGESVD_NONE, B , eig, Z , NULL);
//CV: matrix_free( B );
//CV: }
//CV:
//CV: {
//CV: int i,j;
//CV: /* Lambda1 = (I + Lambda)^(-1) */
//CV:
//CV: for (i=0; i < nrmin; i++)
//CV: eig[i] = 1.0 / (1 + eig[i]);
//CV:
//CV: for (j=0; j < nrmin; j++)
//CV: for (i=0; i < nrmin; i++)
//CV: matrix_imul(Z , i , j , inv_sig0[i]); /* Z2 = Sigma0^(+) * Z; */
//CV: }
//CV: }
//CV:
//CV:
//CV:
//CV: void enkf_analysis_invertS_pre_cv(double truncation,
//CV: int ncomp ,
//CV: const matrix_type * S ,
//CV: const matrix_type * R ,
//CV: matrix_type * V0T ,
//CV: matrix_type * Z ,
//CV: double * eig ,
//CV: matrix_type * U0 ) {
//CV:
//CV: lowrankCinv_pre_cv( S , R , V0T , Z , eig , U0 , truncation , ncomp );
//CV: }
//CV:
//CV:
//CV:
//CV: static void getW_pre_cv(matrix_type * W ,
//CV: const matrix_type * V0T,
//CV: const matrix_type * Z ,
//CV: double * eig ,
//CV: const matrix_type * U0 ,
//CV: int nfolds_CV,
//CV: const matrix_type * A,
//CV: int unique_bootstrap_components ,
//CV: rng_type * rng,
//CV: bool pen_press) {
//CV:
//CV: const int nrobs = matrix_get_rows( U0 );
//CV: const int nrens = matrix_get_columns( V0T );
//CV: const int nrmin = util_int_min( nrobs , nrens );
//CV:
//CV: int i,j;
//CV:
//CV:
//CV: /* Vector with random permutations of the itegers 1,...,nrens */
//CV: int * randperms = util_malloc( sizeof * randperms * nrens, __func__);
//CV: int * indexTest = util_malloc( sizeof * indexTest * nrens, __func__);
//CV: int * indexTrain = util_malloc( sizeof * indexTrain * nrens, __func__);
//CV:
//CV: if(nrens != unique_bootstrap_components)
//CV: nfolds_CV = util_int_min( nfolds_CV , unique_bootstrap_components-1);
//CV:
//CV:
//CV: matrix_type * cvError = matrix_alloc( nrmin,nfolds_CV );
//CV:
//CV: /*Copy Z */
//CV: matrix_type * workZ = matrix_alloc_copy( Z );
//CV:
//CV: int optP;
//CV:
//CV:
//CV: /* start cross-validation: */
//CV:
//CV: const int maxp = matrix_get_rows(V0T);
//CV:
//CV: /* draw random permutations of the integers 0,...,nrens-1 */
//CV: for (i=0; i < nrens; i++)
//CV: randperms[i] = i;
//CV: rng_shuffle_int( rng , randperms , nrens );
//CV:
//CV:
//CV:
//CV: /*need to init cvError to all zeros (?) */
//CV: for (i = 0; i < nrmin; i++){
//CV: for( j = 0; j> nfolds_CV; j++){
//CV: matrix_iset( cvError , i , j , 0.0 );
//CV: }
//CV: }
//CV:
//CV: int ntest, ntrain, k;
//CV: printf("\nStarting cross-validation\n");
//CV: for (i = 0; i < nfolds_CV; i++) {
//CV: printf(".");
//CV:
//CV: ntest = 0;
//CV: ntrain = 0;
//CV: k = i;
//CV: /*extract members for the training and test ensembles */
//CV: for (j = 0; j < nrens; j++) {
//CV: if (j == k) {
//CV: indexTest[ntest] = randperms[j];
//CV: k += nfolds_CV;
//CV: ntest++;
//CV: } else {
//CV: indexTrain[ntrain] = randperms[j];
//CV: ntrain++;
//CV: }
//CV: }
//CV:
//CV: //This one instaed: ?? enkf_analysis_get_cv_error_prin_comp( cv_data , cvError , A , indexTest , indexTrain, ntest, ntrain , i , maxP);
//CV: enkf_analysis_get_cv_error( cvError , A , V0T , workZ , eig , indexTest , indexTrain, ntest, ntrain , i );
//CV: }
//CV: printf("\n");
//CV: /* find optimal truncation value for the cv-scheme */
//CV: optP = enkf_analysis_get_optimal_numb_comp( cvError , maxp, nfolds_CV, pen_press);
//CV:
//CV: printf("Optimal number of components found: %d \n",optP);
//CV: printf("\n");
//CV: FILE * compSel_log = util_fopen("compSel_log_local_cv" , "a");
//CV: fprintf( compSel_log , " %d ",optP);
//CV: fclose( compSel_log);
//CV:
//CV:
//CV: /*free cvError vector and randperm */
//CV: matrix_free( cvError );
//CV: free( randperms );
//CV: free( indexTest );
//CV: free( indexTrain );
//CV:
//CV: /* need to update matrices so that we only use components 1,...,optP */
//CV: /* remove non-zero entries of the z matrix (we do not want to recompute sigma0^(+') * z */
//CV: /* this can surely be done much more efficiently, but for now we want to minimize the
//CV: number of potential bugs in the code for now */
//CV: for (i = optP; i < nrmin; i++) {
//CV: for (j = 0; j < nrmin; j++) {
//CV: matrix_iset(workZ , i , j, 0.0);
//CV: }
//CV: }
//CV:
//CV:
//CV: /*fix the eig vector as well: */
//CV: {
//CV: int i;
//CV: /* lambda1 = (i + lambda)^(-1) */
//CV: for (i=optP; i < nrmin; i++)
//CV: eig[i] = 1.0;
//CV: }
//CV:
//CV: matrix_matmul(W , U0 , workZ); /* x1 = w = u0 * z2 = u0 * sigma0^(+') * z */
//CV:
//CV:
//CV:
//CV:
//CV: /*end cross-validation */
//CV: }
// CV: void enkf_analysis_initX_pre_cv(matrix_type * X ,
// CV: int nfolds_CV,
// CV: enkf_mode_type enkf_mode,
// CV: bool bootstrap ,
// CV: matrix_type * S ,
// CV: matrix_type * R ,
// CV: matrix_type * E ,
// CV: matrix_type * D ,
// CV: rng_type * rng ,
// CV: meas_data_type * meas_data ,
// CV: obs_data_type * obs_data ,
// CV: const matrix_type * randrot ,
// CV: const matrix_type * A ,
// CV: const matrix_type * V0T ,
// CV: const matrix_type * Z ,
// CV: const double * eig ,
// CV: const matrix_type * U0 ,
// CV: meas_data_type * fasit ,
// CV: int unique_bootstrap_components) {
// CV:
// CV: int ens_size = meas_data_get_ens_size( meas_data );
// CV: {
// CV: int nrobs = obs_data_get_active_size(obs_data);
// CV: int nrmin = util_int_min( ens_size , nrobs);
// CV:
// CV:
// CV: /*
// CV: 1: Allocating all matrices
// CV: */
// CV: /*Need a copy of A, because we need it later */
// CV: matrix_type * workA = matrix_alloc_copy( A ); /* <- This is a massive memory requirement. */
// CV: matrix_type * innov = enkf_linalg_alloc_innov( dObs , S );
// CV:
// CV: matrix_type * W = matrix_alloc(nrobs , nrmin);
// CV: double * workeig = util_malloc( sizeof * workeig * nrmin , __func__);
// CV:
// CV:
// CV: // Really unclear whether the row mean should have been shifted from S????
// CV:
// CV: /*copy entries in eig:*/
// CV: {
// CV: int i;
// CV: for (i = 0 ; i < nrmin ; i++)
// CV: workeig[i] = eig[i];
// CV: }
// CV:
// CV: /* Subtracting the ensemble mean of the state vector ensemble */
// CV: matrix_subtract_row_mean( workA );
// CV:
// CV: /*
// CV: 2: Diagonalize the S matrix; singular vectors are stored in W
// CV: and singular values (after some massage) are stored in eig.
// CV: W = X1, eig = inv(I+Lambda1),(Eq.14.30, and 14.29, Evensen, 2007, respectively)
// CV: */
// CV:
// CV: getW_pre_cv(W , V0T , Z , workeig , U0 , nfolds_CV , workA , unique_bootstrap_components , rng);
// CV:
// CV: /*
// CV: 3: actually calculating the X matrix.
// CV: */
// CV: switch (enkf_mode) {
// CV: case(ENKF_STANDARD):
// CV: enkf_linalg_init_stdX( X , S , D , W , workeig , bootstrap);
// CV: break;
// CV: case(ENKF_SQRT):
// CV: enkf_linalg_init_sqrtX(X , S , randrot , dObs , W , workeig , bootstrap );
// CV: break;
// CV: default:
// CV: util_abort("%s: INTERNAL ERROR \n",__func__);
// CV: }
// CV:
// CV: matrix_free( W );
// CV: matrix_free( R );
// CV: matrix_free( S );
// CV: matrix_free( workA );
// CV: matrix_free( innov );
// CV: free( workeig );
// CV:
// CV: if (enkf_mode == ENKF_STANDARD) {
// CV: matrix_free( E );
// CV: matrix_free( D );
// CV: }
// CV:
// CV: enkf_analysis_checkX(X , bootstrap);
// CV: }
// CV: }
/*****************************************************************/
static UTIL_SAFE_CAST_FUNCTION_CONST( bootstrap_enkf_data , BOOTSTRAP_ENKF_TYPE_ID )
void bootstrap_enkf_set_doCV( bootstrap_enkf_data_type * data , bool doCV) {
@@ -376,6 +102,7 @@ void bootstrap_enkf_data_free( void * arg ) {
std_enkf_data_free( boot_data->std_enkf_data );
cv_enkf_data_free( boot_data->cv_enkf_data );
}
free( boot_data );
}
@@ -443,7 +170,8 @@ void bootstrap_enkf_updateA(void * module_data ,
}
if (bootstrap_data->doCV) {
cv_enkf_init_update( bootstrap_data->cv_enkf_data , S_resampled , R , dObs , E , D);
const bool_vector_type * ens_mask = NULL;
cv_enkf_init_update( bootstrap_data->cv_enkf_data , ens_mask , S_resampled , R , dObs , E , D);
cv_enkf_initX( bootstrap_data->cv_enkf_data , X , A_resampled , S_resampled , R , dObs , E , D);
} else
std_enkf_initX(bootstrap_data->std_enkf_data , X , NULL , S_resampled,R, dObs, E,D );
@@ -518,6 +246,26 @@ bool bootstrap_enkf_set_bool( void * arg , const char * var_name , bool value) {
}
bool bootstrap_enkf_has_var( const void * arg, const char * var_name) {
const bootstrap_enkf_data_type * module_data = bootstrap_enkf_data_safe_cast_const( arg );
{
return std_enkf_has_var(module_data->std_enkf_data, var_name);
}
}
double bootstrap_enkf_get_double( const void * arg, const char * var_name) {
const bootstrap_enkf_data_type * module_data = bootstrap_enkf_data_safe_cast_const( arg );
{
return std_enkf_get_double( module_data->std_enkf_data , var_name);
}
}
int bootstrap_enkf_get_int( const void * arg, const char * var_name) {
const bootstrap_enkf_data_type * module_data = bootstrap_enkf_data_safe_cast_const( arg );
{
return std_enkf_get_int( module_data->std_enkf_data , var_name);
}
}
@@ -543,8 +291,9 @@ analysis_table_type SYMBOL_TABLE = {
.updateA = bootstrap_enkf_updateA,
.init_update = NULL,
.complete_update = NULL,
.has_var = NULL,
.get_int = NULL,
.get_double = NULL,
.has_var = bootstrap_enkf_has_var,
.get_int = bootstrap_enkf_get_int,
.get_double = bootstrap_enkf_get_double,
.get_bool = NULL,
.get_ptr = NULL,
};

57
ThirdParty/Ert/devel/libanalysis/src/cv_enkf.c vendored
View File

@@ -63,6 +63,7 @@ struct cv_enkf_data_struct {
static UTIL_SAFE_CAST_FUNCTION( cv_enkf_data , CV_ENKF_TYPE_ID )
static UTIL_SAFE_CAST_FUNCTION_CONST( cv_enkf_data , CV_ENKF_TYPE_ID )
void cv_enkf_set_truncation( cv_enkf_data_type * data , double truncation ) {
@@ -114,6 +115,7 @@ void cv_enkf_data_free( void * arg ) {
matrix_safe_free( cv_data->Rp );
matrix_safe_free( cv_data->Dp );
}
free( cv_data );
}
@@ -122,6 +124,7 @@ void cv_enkf_data_free( void * arg ) {
void cv_enkf_init_update( void * arg ,
const bool_vector_type * ens_mask ,
const matrix_type * S ,
const matrix_type * R ,
const matrix_type * dObs ,
@@ -643,6 +646,53 @@ long cv_enkf_get_options( void * arg , long flag) {
}
}
bool cv_enkf_has_var( const void * arg, const char * var_name) {
{
if (strcmp(var_name , ENKF_NCOMP_KEY_) == 0)
return true;
else if (strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0)
return true;
else if (strcmp(var_name , NFOLDS_KEY) == 0)
return true;
else if (strcmp(var_name , CV_PEN_PRESS_KEY) == 0)
return true;
else
return false;
}
}
double cv_enkf_get_double( const void * arg, const char * var_name) {
const cv_enkf_data_type * module_data = cv_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0)
return module_data->truncation;
else
return -1;
}
}
int cv_enkf_get_int( const void * arg, const char * var_name) {
const cv_enkf_data_type * module_data = cv_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_NCOMP_KEY_) == 0)
return module_data->subspace_dimension;
else if (strcmp(var_name , NFOLDS_KEY) == 0)
return module_data->nfolds;
else
return -1;
}
}
bool cv_enkf_get_bool( const void * arg, const char * var_name) {
const cv_enkf_data_type * module_data = cv_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , CV_PEN_PRESS_KEY) == 0)
return module_data->penalised_press;
else
return false;
}
}
@@ -665,8 +715,9 @@ analysis_table_type SYMBOL_TABLE = {
.updateA = NULL,
.init_update = cv_enkf_init_update ,
.complete_update = cv_enkf_complete_update ,
.has_var = NULL,
.get_int = NULL,
.get_double = NULL,
.has_var = cv_enkf_has_var,
.get_int = cv_enkf_get_int,
.get_double = cv_enkf_get_double,
.get_bool = cv_enkf_get_bool,
.get_ptr = NULL,
};

196
ThirdParty/Ert/devel/libanalysis/src/enkf_linalg.c vendored
View File

@@ -73,19 +73,19 @@ void enkf_linalg_genX2(matrix_type * X2 , const matrix_type * S , const matrix_t
/*This function is similar to enkf_linalg_svdS but it returns the eigen values without its inverse and also give the matrices truncated U VT and Sig0*/
// Trunc.SVD(S) = U0 * Sig0 * V0T
int enkf_linalg_svd_truncation(const matrix_type * S ,
double truncation ,
int ncomp ,
dgesvd_vector_enum store_V0T ,
double * sig0,
matrix_type * U0 ,
matrix_type * V0T) {
double truncation ,
int ncomp ,
dgesvd_vector_enum store_V0T ,
double * sig0,
matrix_type * U0 ,
matrix_type * V0T) {
int num_significant = -1;
int nrows = matrix_get_rows(S);
int ncolumns= matrix_get_columns(S);
if (((truncation > 0) && (ncomp < 0)) ||
((truncation < 0) && (ncomp > 0))) {
@@ -95,7 +95,6 @@ int enkf_linalg_svd_truncation(const matrix_type * S ,
matrix_dgesvd(DGESVD_MIN_RETURN , store_V0T , workS , sig0 , U0 , V0T);
matrix_free( workS );
}
int i;
if (ncomp > 0)
@@ -125,7 +124,7 @@ int enkf_linalg_svd_truncation(const matrix_type * S ,
matrix_resize(U0 , nrows , num_significant , true);
matrix_resize(V0T , num_significant , ncolumns , true);
}
else
else
util_abort("%s: truncation:%g ncomp:%d - invalid ambigous input.\n",__func__ , truncation , ncomp );
return num_significant;
}
@@ -590,7 +589,8 @@ void enkf_linalg_get_PC( const matrix_type * S0,
double truncation,
int ncomp,
matrix_type * PC,
matrix_type * PC_obs ) {
matrix_type * PC_obs,
double_vector_type * singular_values) {
const int nrobs = matrix_get_rows( S0 );
const int nrens = matrix_get_columns( S0 );
@@ -598,10 +598,12 @@ void enkf_linalg_get_PC( const matrix_type * S0,
matrix_type * U0 = matrix_alloc( nrobs , nrens );
matrix_type * S = matrix_alloc_copy( S0 );
double * inv_sig0 = util_calloc( nrmin , sizeof * inv_sig0 );
double * inv_sig0;
double_vector_iset( singular_values , nrmin - 1 , 0 );
matrix_subtract_row_mean( S );
ncomp = util_int_min( ncomp , nrmin );
inv_sig0 = double_vector_get_ptr( singular_values );
{
matrix_type * S_mean = matrix_alloc( nrobs , 1 );
int num_PC = enkf_linalg_svdS(S , truncation , ncomp, DGESVD_NONE , inv_sig0 , U0 , NULL);
@@ -629,10 +631,11 @@ void enkf_linalg_get_PC( const matrix_type * S0,
matrix_dgemm( PC_obs , U0 , S_mean , true , false , 1.0 , 0.0 );
}
for (int i=0; i < double_vector_size( singular_values ); i++)
inv_sig0[i] = 1.0 / inv_sig0[i];
matrix_free( S_mean );
}
free( inv_sig0 );
matrix_free( S );
matrix_free( U0 );
}
@@ -641,147 +644,114 @@ void enkf_linalg_get_PC( const matrix_type * S0,
void enkf_linalg_rml_enkfX1(matrix_type *X1, matrix_type *Udr, matrix_type *D, matrix_type *R)
{
/*This routine computes X1 for RML_EnKF module as X1 = Ud(T)*Cd(-1/2)*D -- D= (dk-do)
here the negative sign cancels with one needed in X3 matrix computation*/
/*
This routine computes X1 for RML_EnKF module as X1 = Ud(T)*Cd(-1/2)*D -- D= (dk-do)
here the negative sign cancels with one needed in X3 matrix computation
*/
matrix_type * tmp = matrix_alloc(matrix_get_columns(Udr),matrix_get_rows(R));
matrix_dgemm( tmp, Udr, R, true, false, 1.0, 0.0);
printf("tmp matrix is computed");
matrix_dgemm( X1 , tmp, D,false,false, 1.0, 0.0);
printf("X1 is computed");
matrix_free(tmp);
matrix_matmul_with_transpose( tmp, Udr, R, true, false);
matrix_matmul( X1 , tmp, D);
matrix_free(tmp);
}
void enkf_linalg_rml_enkfX2(matrix_type *X2, double *Wdr, matrix_type * X1 ,double a, int nsign)
{
/*This routine computes X2 for RML_EnKF module as X2 = ((a*Ipd)+Wd^2)^-1 * X1 */
/* Since a+Ipd & Wd are diagonal in nature the computation is reduced to array operations*/
double * tmp = util_calloc(nsign , sizeof * tmp );
for (int i=0; i< nsign ; i++)
{
tmp[i] = 1/ (a+ (Wdr[i]*Wdr[i]));
matrix_scale_row(X1,i,tmp[i]);
}
matrix_assign(X2,X1);
void enkf_linalg_rml_enkfX2(matrix_type *X2 , double *Wdr , matrix_type * X1 , double a , int nsign)
{
/*
This routine computes X2 for RML_EnKF module as X2 = ((a*Ipd)+Wd^2)^-1 * X1
Since a+Ipd & Wd are diagonal in nature the computation is reduced to array operations
*/
for (int i=0; i< nsign ; i++) {
double scale_factor = 1 / (a + (Wdr[i]*Wdr[i]));
matrix_scale_row(X1 , i , scale_factor);
}
matrix_assign(X2,X1);
}
void enkf_linalg_rml_enkfX3(matrix_type *X3, matrix_type *VdTr, double *Wdr, matrix_type *X2, int nsign)
{
/*This routine computes X3 for RML_EnKF module as X3 = Vd *Wd*X2 */
matrix_type *tmp = matrix_alloc_copy(VdTr);
matrix_matlab_dump(tmp, "matrixVdTrcopy.dat");
printf("matrix X3 started");
for (int i=0; i< nsign ; i++)
{
printf("The value of Wd(%d) is = %5.2f \n",i, Wdr[i]);
matrix_scale_row(tmp, i, Wdr[i]);
}
/*
This routine computes X3 for RML_EnKF module as X3 = Vd *Wd*X2
*/
matrix_dgemm( X3 , tmp , X2 , true, false, 1.0, 0.0);
printf("\nWd: ");
matrix_type *tmp = matrix_alloc_copy(VdTr);
for (int i=0; i< nsign ; i++) {
printf("%5.2f ", Wdr[i]);
matrix_scale_row(tmp, i, Wdr[i]);
}
printf("\n\n");
matrix_matmul_with_transpose( X3 , tmp , X2 , true, false);
matrix_free(tmp);
}
void enkf_linalg_rml_enkfdA(matrix_type * dA1, matrix_type * Dm, matrix_type * X3) //dA = Dm * X3
{
matrix_dgemm (dA1, Dm, X3,false, false, 1.0, 0.0 );
}
double enkf_linalg_data_mismatch(matrix_type *D , matrix_type *R , matrix_type *Sk)
{
matrix_type * tmp = matrix_alloc (matrix_get_columns(D), matrix_get_columns(R));
double mean;
printf("-----------------------------------------------------------------\n");
printf("%s:%d Calling matrix_dgemm() \n",__func__ , __LINE__);
// C A B
matrix_dgemm(tmp, D, R,true, false, 1.0, 0.0);
printf("-----------------------------------------------------------------\n");
printf("%s:%d Calling matrix_dgemm() \n",__func__ , __LINE__);
matrix_dgemm(Sk, tmp, D, false, false, 1.0, 0.0);
printf("-----------------------------------------------------------------\n");
printf("The data mismatch computed");
double mismatch;
matrix_matmul_with_transpose(tmp, D, R,true, false); // tmp = D' * R, i.e. N-by-p
matrix_matmul(Sk, tmp, D); // Sk = D' * R * D
// Calculate the mismatch
mean = matrix_trace(Sk)/(matrix_get_columns(D));
return mean;
mismatch = matrix_trace(Sk)/(matrix_get_columns(D));
return mismatch;
}
// Cd = SampCov(E) (including (N-1) normalization)
void enkf_linalg_Covariance(matrix_type *Cd, const matrix_type *E, double nsc ,int nrobs)
{
matrix_matlab_dump( E, "matrixE.dat");
printf("Starting Dgemm for EE(T)\n");
matrix_dgemm(Cd, E, E,false,true, 1.0, 0.0);
for (int i=0; i< nrobs; i++)
{
for (int j=0;j< nrobs; j++)
{
if (i!=j)
matrix_iset(Cd, i, j, 0.0);
}
matrix_matmul_with_transpose(Cd, E, E,false,true);
for (int i=0; i< nrobs; i++) {
for (int j=0;j< nrobs; j++) {
if (i!=j)
matrix_iset(Cd, i, j, 0.0);
}
nsc= nsc*nsc;
printf("nsc = %5.3f\n",nsc);
}
nsc = nsc*nsc;
matrix_scale(Cd,nsc);
}
void enkf_linalg_rml_enkfAm(matrix_type * Um, double * Wm,int nsign1){
for (int i=0; i< nsign1 ; i++)
{
double sc = 1/ Wm[i];
matrix_scale_column(Um, i, sc);
}
// Scale columns (not rows!) of Um by entries in diagonal Wm
void enkf_linalg_rml_enkfAm(matrix_type * Um, const double * Wm,int nsign1){
for (int i=0; i< nsign1 ; i++) {
double sc = 1 / Wm[i];
matrix_scale_column(Um, i, sc);
}
}
void enkf_linalg_rml_enkfX4 (matrix_type *X4,matrix_type *Am, matrix_type *A){
matrix_dgemm(X4, Am, A, true, false, 1.0, 0.0);
}
void enkf_linalg_rml_enkfX5(matrix_type * X5,matrix_type * Am, matrix_type * X4){
matrix_dgemm (X5, Am, X4, false,false,1.0,0.0);
}
void enkf_linalg_rml_enkfX6(matrix_type * X6, matrix_type * Dk, matrix_type * X5){
matrix_dgemm(X6, Dk, X5, true, false, 1.0, 0.0);
}
void enkf_linalg_rml_enkfX7(matrix_type * X7, matrix_type * VdT, double * Wdr, double a, matrix_type * X6){
int nsign = matrix_get_rows(VdT);
int ens_size= matrix_get_columns(VdT);
matrix_type *tmp1= matrix_alloc_copy(VdT);
matrix_type *tmp2= matrix_alloc(ens_size,ens_size);
int ens_size = matrix_get_columns(VdT);
matrix_type *tmp1 = matrix_alloc_copy(VdT);
matrix_type *tmp2 = matrix_alloc(ens_size,ens_size);
for (int i=0; i< nsign ; i++)
{
double tmp = 1/ ( a + (Wdr[i]*Wdr[i]));
matrix_scale_row(tmp1, i ,tmp);
}
matrix_dgemm(tmp2, tmp1, VdT, true, false, 1.0, 0.0);
for (int i=0; i < nsign ; i++) {
double scale_factor = 1 / ( a + (Wdr[i]*Wdr[i]));
matrix_scale_row( tmp1 , i , scale_factor);
}
matrix_matmul_with_transpose(tmp2, tmp1, VdT, true, false);
matrix_matmul(X7, tmp2, X6);
matrix_free(tmp1);
matrix_dgemm(X7, tmp2, X6, false,false, 1.0, 0.0);
matrix_free(tmp2);
}
void enkf_linalg_rml_enkfXdA2(matrix_type * dA2, matrix_type * Dk,matrix_type * X7){
matrix_dgemm(dA2, Dk, X7, false, false, 1.0, 0.0);
}

44
ThirdParty/Ert/devel/libanalysis/src/fwd_step_enkf.c vendored
View File

@@ -20,7 +20,8 @@
#include <string.h>
#include <math.h>
#include <stdio.h>
#include <ert/util/type_macros.h>
#include <ert/util/util.h>
#include <ert/util/rng.h>
#include <ert/util/matrix.h>
@@ -50,7 +51,7 @@ struct fwd_step_enkf_data_struct {
};
static UTIL_SAFE_CAST_FUNCTION_CONST( fwd_step_enkf_data , FWD_STEP_ENKF_TYPE_ID )
static UTIL_SAFE_CAST_FUNCTION( fwd_step_enkf_data , FWD_STEP_ENKF_TYPE_ID )
@@ -177,6 +178,7 @@ void fwd_step_enkf_data_free( void * arg ) {
}
}
}
free( fwd_step_data );
}
@@ -217,6 +219,37 @@ long fwd_step_enkf_get_options( void * arg , long flag) {
}
}
bool fwd_step_enkf_has_var( const void * arg, const char * var_name) {
{
if (strcmp(var_name , NFOLDS_KEY) == 0)
return true;
else if (strcmp(var_name , R2_LIMIT_KEY ) == 0)
return true;
else
return false;
}
}
double fwd_step_enkf_get_double( const void * arg, const char * var_name) {
const fwd_step_enkf_data_type * module_data = fwd_step_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , R2_LIMIT_KEY ) == 0)
return module_data->r2_limit;
else
return -1;
}
}
int fwd_step_enkf_get_int( const void * arg, const char * var_name) {
const fwd_step_enkf_data_type * module_data = fwd_step_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , NFOLDS_KEY) == 0)
return module_data->nfolds;
else
return -1;
}
}
@@ -239,9 +272,10 @@ analysis_table_type SYMBOL_TABLE = {
.updateA = fwd_step_enkf_updateA,
.init_update = NULL ,
.complete_update = NULL ,
.has_var = NULL ,
.get_int = NULL ,
.get_double = NULL ,
.has_var = fwd_step_enkf_has_var,
.get_int = fwd_step_enkf_get_int ,
.get_double = fwd_step_enkf_get_double ,
.get_bool = NULL ,
.get_ptr = NULL
};

29
ThirdParty/Ert/devel/libanalysis/src/sqrt_enkf.c vendored
View File

@@ -48,6 +48,7 @@ typedef struct {
static UTIL_SAFE_CAST_FUNCTION( sqrt_enkf_data , SQRT_ENKF_TYPE_ID )
static UTIL_SAFE_CAST_FUNCTION_CONST( sqrt_enkf_data , SQRT_ENKF_TYPE_ID )
void * sqrt_enkf_data_alloc( rng_type * rng ) {
@@ -164,6 +165,27 @@ void sqrt_enkf_complete_update( void * arg ) {
}
}
bool sqrt_enkf_has_var( const void * arg, const char * var_name) {
const sqrt_enkf_data_type * module_data = sqrt_enkf_data_safe_cast_const( arg );
{
return std_enkf_has_var(module_data->std_data, var_name);
}
}
double sqrt_enkf_get_double( const void * arg, const char * var_name) {
const sqrt_enkf_data_type * module_data = sqrt_enkf_data_safe_cast_const( arg );
{
return std_enkf_get_double( module_data->std_data , var_name);
}
}
int sqrt_enkf_get_int( const void * arg, const char * var_name) {
const sqrt_enkf_data_type * module_data = sqrt_enkf_data_safe_cast_const( arg );
{
return std_enkf_get_int( module_data->std_data , var_name);
}
}
/*****************************************************************/
@@ -189,9 +211,10 @@ analysis_table_type SYMBOL_TABLE = {
.init_update = sqrt_enkf_init_update,
.complete_update = sqrt_enkf_complete_update,
.get_options = sqrt_enkf_get_options,
.has_var = NULL,
.get_int = NULL,
.get_double = NULL,
.has_var = sqrt_enkf_has_var,
.get_int = sqrt_enkf_get_int,
.get_double = sqrt_enkf_get_double,
.get_bool = NULL,
.get_ptr = NULL
};

39
ThirdParty/Ert/devel/libanalysis/src/std_enkf.c vendored
View File

@@ -86,6 +86,8 @@ struct std_enkf_data_struct {
long option_flags;
};
static UTIL_SAFE_CAST_FUNCTION_CONST( std_enkf_data , STD_ENKF_TYPE_ID )
/*
This is a macro which will expand to generate a function:
@@ -229,6 +231,36 @@ long std_enkf_get_options( void * arg , long flag ) {
}
}
bool std_enkf_has_var( const void * arg, const char * var_name) {
{
if (strcmp(var_name , ENKF_NCOMP_KEY_) == 0)
return true;
else if (strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0)
return true;
else
return false;
}
}
double std_enkf_get_double( const void * arg, const char * var_name) {
const std_enkf_data_type * module_data = std_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_TRUNCATION_KEY_) == 0)
return module_data->truncation;
else
return -1;
}
}
int std_enkf_get_int( const void * arg, const char * var_name) {
const std_enkf_data_type * module_data = std_enkf_data_safe_cast_const( arg );
{
if (strcmp(var_name , ENKF_NCOMP_KEY_) == 0)
return module_data->subspace_dimension;
else
return -1;
}
}
/**
@@ -257,9 +289,10 @@ analysis_table_type SYMBOL_TABLE = {
.updateA = NULL,
.init_update = NULL,
.complete_update = NULL,
.has_var = NULL,
.get_int = NULL,
.get_double = NULL,
.has_var = std_enkf_has_var,
.get_int = std_enkf_get_int,
.get_double = std_enkf_get_double,
.get_bool = NULL,
.get_ptr = NULL,
};

21
ThirdParty/Ert/devel/libanalysis/tests/CMakeLists.txt vendored
View File

@@ -1,10 +1,19 @@
if (BUILD_APPLICATIONS)
add_test( analysis_module_test_RML ${EXECUTABLE_OUTPUT_PATH}/ert_module_test ${LIBRARY_OUTPUT_PATH}/rml_enkf.so )
add_test( analysis_module_test_RMLI ${EXECUTABLE_OUTPUT_PATH}/ert_module_test ${LIBRARY_OUTPUT_PATH}/rmli_enkf.so )
endif()
ert_module_name( VAR_RML rml_enkf ${LIBRARY_OUTPUT_PATH} )
add_executable(analysis_test_external_module analysis_test_external_module.c )
target_link_libraries( analysis_test_external_module analysis util test_util )
add_test( analysis_module_rml ${EXECUTABLE_OUTPUT_PATH}/analysis_test_external_module "RML_ENKF" ${LIBRARY_OUTPUT_PATH}/rml_enkf.so 41 ITER LAMBDA0)
add_test( analysis_module_rmli ${EXECUTABLE_OUTPUT_PATH}/analysis_test_external_module "RMLI_ENKF" ${LIBRARY_OUTPUT_PATH}/rmli_enkf.so 41 ITER)
add_test( analysis_module_rml ${EXECUTABLE_OUTPUT_PATH}/analysis_test_external_module "RML_ENKF" ${VAR_RML} 41
ITER:45
USE_PRIOR:False
LAMBDA_REDUCE:0.10
LAMBDA_INCREASE:2.5
ENKF_TRUNCATION:0.77
LAMBDA0:0.25
LAMBDA_MIN:0.01
LOG_FILE:LogFile.txt
CLEAR_LOG:True
LAMBDA_RECALCULATE:True )

49
ThirdParty/Ert/devel/libanalysis/tests/analysis_test_external_module.c vendored
View File

@@ -26,6 +26,49 @@
void test_set_get(analysis_module_type * module , const char * var_value) {
char * var , *string_value;
util_binary_split_string( var_value , ":" , false , &var , &string_value);
if (var && string_value) {
printf("Testing variable:%s \n",var);
while (true) {
int int_value;
double double_value;
bool bool_value;
test_assert_true(analysis_module_has_var( module , var ));
if (util_sscanf_int( string_value , &int_value)) {
test_assert_true(analysis_module_set_var( module , var , string_value ));
test_assert_int_equal( int_value , analysis_module_get_int( module , var ));
break;
}
if (util_sscanf_double( string_value , &double_value)) {
test_assert_true(analysis_module_set_var( module , var , string_value ));
test_assert_double_equal( double_value , analysis_module_get_double( module , var ));
break;
}
if (util_sscanf_bool( string_value , &bool_value)) {
test_assert_true(analysis_module_set_var( module , var , string_value ));
test_assert_bool_equal( bool_value , analysis_module_get_bool( module , var ));
break;
}
test_assert_true(analysis_module_set_var( module , var , string_value ));
test_assert_string_equal( string_value , (const char *) analysis_module_get_ptr( module , var ));
break;
}
} else {
fprintf(stderr,"Invalid test input data: %s -> could not split in var:value\n" , var_value);
exit(1);
}
}
void load_module( rng_type * rng , const char * user_name , const char * lib_name, const char * options_str , int nvar , const char ** var_list) {
long flags = strtol(options_str , NULL , 10);
analysis_module_type * analysis_module = analysis_module_alloc_external(rng , user_name , lib_name);
@@ -36,9 +79,9 @@ void load_module( rng_type * rng , const char * user_name , const char * lib_nam
test_assert_string_equal( lib_name , analysis_module_get_lib_name(analysis_module));
test_assert_true( analysis_module_is_instance( analysis_module));
for (int i=0; i < nvar; i++) {
printf("has_var(%s) \n" , var_list[i]);
test_assert_true( analysis_module_has_var(analysis_module , var_list[i]));
{
for (int i=0; i < nvar; i++)
test_set_get( analysis_module , var_list[i] );
}
test_assert_false( analysis_module_has_var(analysis_module , "DoesNotHaveThisVariable"));