/DataAssimilation/common/data_assimilation_routines.h
C Header | 1398 lines | 915 code | 248 blank | 235 comment | 87 complexity | a5c27bede34902184ef0e3c81debb85c MD5 | raw file
Large files files are truncated, but you can click here to view the full file
1 2//command: mpiexec -n 1 ./01 -measurement-noise-scale 0.1 -global-noise-scale 0.1 -initialization-error 100 -measurement-noise-amount 100 3static char help[] = "rUKF approach.\n\ 4 ... \n\n"; 5 6#undef __FUNCT__ 7#define __FUNCT__ "main" 8#define PETSC_OPERATIONS_DEBUG 1 9 10#include "slepcsvd.h" 11#include "petscksp.h" 12#include "petscmat.h" 13#include <stdlib.h> 14#include <assert.h> 15#include <gsl/gsl_errno.h> 16#include <gsl/gsl_spline.h> 17#include "data_assimilation_utility_routines.h" 18 19#define ERROR_TYPE_COMMON 1 20#define INFO_PETSC_INT 1 21#define INFO_PETSC_REAL 2 22#define INFO_PETSC_VEC 3 23#define INFO_PETSC_MAT 4 24#define INFO_PETSC_INTARRAY 5 25#define INFO_PETSC_REALARRAY 6 26 27#define FILTER_TYPE_RUKF 7 28 29#define MEAUSREMENT_TYPE_FROM_HEADER_AND_DAT_FILES 1 30 31#define DIAGNOSIS_DISP_EVERYTHING 0x8000 32#define DIAGNOSIS_DISP_MEDIUM 0x4000 33#define DIAGNOSIS_DISP_MINIMUM 0x2000 34#define DIAGNOSIS_DISP_INITIAL_CONF 0x0001 35 36#define app_context_assert(x,msg,app_context) {if(!(x)) app_context_report_error(ERROR_TYPE_COMMON, msg, app_context, __LINE__,__FUNCT__,__FILE__);} 37 38#ifndef MEASUREMENT_DIR 39 #define MEASUREMENT_DIR "./" //where to find the measurement files 40#endif 41 42#define MAX_SIZE_PATH_NAME 256 43#define MAX_SIZE_FILE_NAME 128 44 45//------------------------------------------------------------------------------------------------- 46//---structures definitions--- 47struct _APP_CONTEXT_TYPE{ 48 MPI_Comm mpi_comm_global;//process in the work group 49 MPI_Comm mpi_comm_local; //this process 50 PetscViewer pv_global; 51 PetscViewer pv_local; 52 PetscErrorCode ierr; 53 PetscInt diagnosisMode; 54 55 PetscReal measurement_noise_scale; // for R 56 PetscReal global_noise_scale; // for both R and P_0 57// PetscReal initialization_error; // for theta_0 58 PetscReal measurement_noise_amount; // for y_k 59 60 PetscInt n_x;//ndof 61 PetscInt n_theta; //npdof 62 PetscReal* x_0; 63 PetscReal x_scale_constant; 64 PetscReal* x_scale; //currently initialized from x_scale_constant; 65 PetscReal* theta_truth; 66 PetscReal* theta_0; 67 PetscReal* theta_scale; 68 69 //temporary storage for computation 70 //-SVD related 71 Mat S_ensemble; //ensemble matrix, only for storage purpose (for doing SVD), app_context->n_x+app_context->n_theta by model->n_ensem 72 Mat S_ensemble2; //ensemble matrix, only for storage purpose (for S_ensemble*S_ensemble2^T) 73 Mat P_z; //BIG! error covariance for the augmented state vecotr, used by rUKF for doing SVD,app_context->n_x+app_context->n_theta by app_context->n_x+app_context->n_theta 74 Mat S_z; //SVD result of P_z, app_context->n_x+app_context->n_theta by filter->subfilter->r 75 Mat S_bar; //used in the measurement update by rUKF, measurement->n_y by filter->subfilter->r 76 Mat S_bar2; //temperory storage for a matrix of the size measurement->n_y by filter->subfilter->r 77 //-general temperory usage 78 Mat tmp_r; // r by r matrix 79 Mat tmp_nz_by_r; 80 Mat tmp_r_by_ny; 81 Vec tmp_z; //n_x+n_theta by 1 vector 82 Vec tmp_x; //n_x by 1 vector 83 Vec tmp_theta; //n_theta by 1 vector 84 Vec tmp_y; //measurement->n_y by 1 vector 85 PetscInt* tmp_x_indices; 86 PetscInt* tmp_theta_indices; 87 88 //measurement info 89 char* measurement_dir_name; 90 char* measurement_header_file_name; 91 char* measurement_data_file_name; 92 char* measurement_matrix_file_name; 93 char* measurement_error_covariance_matrix_file_name; 94 int measurementType; 95}; 96typedef struct _APP_CONTEXT_TYPE APP_CONTEXT_TYPE; 97 98 99struct _MODEL_TYPE{ 100 struct _MODEL_TYPE* (*f)(struct _MODEL_TYPE*, APP_CONTEXT_TYPE*); //f: model forward operator : n (number of ensembles) states at k --> n states at k+1 (n>=1) 101 102 PetscInt n_x;//ndof 103 PetscInt n_theta; 104 PetscInt n_z; 105 106 Vec x_k_hat; 107 Vec theta_k_hat; 108 Vec x_k_minus_1_hat; 109 Vec theta_k_minus_1_hat; 110 111 PetscInt k; //time step 112 PetscReal t; //time 113 114 PetscInt n_ensem; //number of ensembles 115 Vec* x_k_minus_1; 116 Vec* theta_k_minus_1; 117 Vec* x_k; 118 Vec* theta_k; 119// measurementDisplacementScale=[]; 120//ipfibreNums=[]; 121//fixed_contractility_factor=[]; 122//phase=[]; 123//Cai_i_max=[]; 124//parameter_permutation_matrx =[]; 125//repeatScaleFacotrs =[]; 126//repeatNumber=[]; 127}; 128typedef struct _MODEL_TYPE MODEL_TYPE; 129 130struct _FILTER_SUBTYPE{ 131 PetscInt filter_type; 132 PetscInt r; //rank ( r>=(n_theat>0)?n_theat:n_x, r<= n_theta+n_x; 133 Mat P_x; //BIG! error covariance for the state vecotr, not used by rUKF 134 Mat P_theta; //error covariance for the parameter vector, not used by rUKF 135 136 Mat S_x; //square-root of P_x 137 Mat S_theta; //square-root of P_theta 138}; 139 140typedef struct _FILTER_SUBTYPE FILTER_SUBTYPE; 141 142struct _FILTER_TYPE{ 143 PetscInt n_x; //ndof, n_x>=0 144 PetscInt n_theta; //ndof, n_theta>=0 145 PetscInt n_z; 146 Mat K; //Kalman gain matrix, n_x + n_theta by measurement->n_y 147// Vec x; 148// Vec theta; 149// PetscInt filter_type; 150 FILTER_SUBTYPE* subfilter; //FILTER_SUBTYPE 151}; 152typedef struct _FILTER_TYPE FILTER_TYPE; 153 154struct _OPTIMIZER_TYPE{ 155 PetscInt n_theta; //ndof 156}; 157typedef struct _OPTIMIZER_TYPE OPTIMIZER_TYPE; 158 159 160struct _MEAUSREMENT_TYPE{ 161 //Vec k_y; 162 //Vec t_y; 163 Vec y; //measurement_petsc_vector, current measurement in use(maybe interpolated), n_y by 1 164 Mat H; //measurement matrix, n_y by n_z 165 Mat R_d; //measurement error covariance, n_y by n_y 166 Mat inv_R_d; //measurement error covariance inverse, n_y by n_y 167 168 int n_y;//length of single measurement vector 169 int n_z;//length of state vector 170 int K; //number of measurements 171 int nMeta; //default 2 172 173 double* map_index2time; //nMeta(2, index, time) by K 174 double* all_measurements; //n_y by K 175 double* single_measurement; //n_y by 1 176 double* H_values; //n_y by n_z; 177 double* R_d_values;//n_y by n_y 178// Vec theta_y; //ground-truth parameter 179//xiCoordsFileName=[]; 180//sampleType=[]; 181//nSampleOrder=[]; 182//elementsID=[]; 183//facesID=[]; 184//tempDirID =[]; 185//experimentID =[]; 186 187 188}; 189typedef struct _MEAUSREMENT_TYPE MEAUSREMENT_TYPE; 190 191 192//------------------------------------------------------------------------------------------------- 193//---all routines--- 194 195PetscErrorCode 196app_context_report_error(PetscErrorCode errorcode, const char* msg, APP_CONTEXT_TYPE* app_context, int lineNum, const char* funct, const char* file); 197 198PetscErrorCode 199app_context_report_info(PetscInt info_type, PetscInt info_size, void* info, const char* msg, APP_CONTEXT_TYPE* app_context); 200 201int 202app_context_destroy_and_finalize(APP_CONTEXT_TYPE* app_context); 203 204APP_CONTEXT_TYPE* 205app_context_create_and_initialize(int argc,char **args); 206 207MODEL_TYPE* 208model_create_and_initialize(MODEL_TYPE* (*model_opeartor)(MODEL_TYPE*, APP_CONTEXT_TYPE*), APP_CONTEXT_TYPE* app_context); 209 210MODEL_TYPE* 211model_set_opeartor_and_intial_condition(MODEL_TYPE* model, MODEL_TYPE* (*model_opeartor)(MODEL_TYPE*, APP_CONTEXT_TYPE*), APP_CONTEXT_TYPE* app_context); 212 213FILTER_TYPE* 214filter_create_and_initialize(APP_CONTEXT_TYPE* app_context,MODEL_TYPE* model); 215 216FILTER_TYPE* 217filter_set_intial_condition(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context); 218 219OPTIMIZER_TYPE* 220optimizer_create_and_initialize(APP_CONTEXT_TYPE* app_context); 221 222OPTIMIZER_TYPE* 223optimizer_set_intial_condition(OPTIMIZER_TYPE* optimizer, APP_CONTEXT_TYPE* app_context); 224 225//create and initialize the filter 226FILTER_SUBTYPE* 227filter_subtype_create_and_initialize(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context); 228 229// Initialize the filter subtype 230// PetscInt filter_type; 231// PetscInt r; //rank ( r>=(n_theat>0)?n_theat:1, r<= n_theta+n_x; 232// Mat S_x; 233// Mat S_theta; 234FILTER_SUBTYPE* 235filter_subtype_set_intial_condition(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context); 236 237//------------------------------------------------------------------------------------------------- 238//---app_context routines--- 239PetscErrorCode 240app_context_report_error(PetscErrorCode errorcode, const char* msg, APP_CONTEXT_TYPE* app_context, int lineNum, const char* funct, const char* file) 241{ 242 static int previous_lineNum=0; 243 if(errorcode){ 244 PetscPrintf(app_context->mpi_comm_global, "== fatal error: %s: at line %d, function %s(), file %s. prevous line number: %d \n", msg, lineNum, funct, file,previous_lineNum); 245 //CHKERRQ(app_context->ierr); //if errorcode is non-zero, 246 SETERRQ(errorcode,msg); 247 }else{ 248 previous_lineNum = lineNum; 249 return 0; 250 } 251} 252 253PetscErrorCode 254app_context_report_info(PetscInt info_type, PetscInt info_size, void* info, const char* msg, APP_CONTEXT_TYPE* app_context) 255{ 256 app_context->ierr = PetscPrintf(app_context->mpi_comm_global, "== info: %s: ", msg); 257 //app_context_report_error(app_context->ierr,"app_context_report_info: PetscPrintf",app_context,__LINE__,__FUNCT__,__FILE__); 258 259 switch(info_type) 260 { 261 case INFO_PETSC_INT: 262 case INFO_PETSC_INTARRAY: 263 app_context->ierr = PetscIntView(info_size, (PetscInt*)info, app_context->pv_global); 264 app_context_report_error(app_context->ierr,"app_context_report_info: INFO_PETSC_INT",app_context,__LINE__,__FUNCT__,__FILE__); 265 break; 266 case INFO_PETSC_REAL: 267 case INFO_PETSC_REALARRAY: 268 app_context->ierr = PetscRealView(info_size, (PetscReal*)info, app_context->pv_global); 269 app_context_report_error(app_context->ierr,"app_context_report_info: INFO_PETSC_REAL",app_context,__LINE__,__FUNCT__,__FILE__); 270 break; 271 case INFO_PETSC_VEC: 272 app_context->ierr = PetscPrintf(app_context->mpi_comm_global, "\n", msg); 273 app_context_report_error(app_context->ierr,"app_context_report_info: INFO_PETSC_VEC",app_context,__LINE__,__FUNCT__,__FILE__); 274 app_context->ierr = VecView( *((Vec*)info), app_context->pv_global); 275 app_context_report_error(app_context->ierr,"app_context_report_info: INFO_PETSC_VEC",app_context,__LINE__,__FUNCT__,__FILE__); 276 break; 277 case INFO_PETSC_MAT: 278 app_context->ierr = PetscPrintf(app_context->mpi_comm_global, "\n", msg); 279 app_context_report_error(app_context->ierr,"app_context_report_info: INFO_PETSC_MAT",app_context,__LINE__,__FUNCT__,__FILE__); 280 app_context->ierr = MatView( *((Mat*)info), app_context->pv_global); 281 app_context_report_error(app_context->ierr,"app_context_report_info: INFO_PETSC_MAT",app_context,__LINE__,__FUNCT__,__FILE__); 282 break; 283 default: 284 app_context->ierr = PetscPrintf(app_context->mpi_comm_global, "== wrong info_type: %d\n", info_type); 285 app_context_report_error(app_context->ierr,"app_context_report_info: default",app_context,__LINE__,__FUNCT__,__FILE__); 286 break; 287 } 288 app_context->ierr = PetscPrintf(app_context->mpi_comm_global, "\n"); 289 app_context_report_error(app_context->ierr,"app_context_report_info: PetscPrintf",app_context,__LINE__,__FUNCT__,__FILE__); 290 return 0; 291} 292 293APP_CONTEXT_TYPE* 294app_context_create_and_initialize(int argc,char **args) 295{ 296 APP_CONTEXT_TYPE* app_context = NULL; 297 PetscTruth flg = PETSC_FALSE, flg_1=PETSC_FALSE; 298 int i=0; 299// char options[2][PETSC_MAX_PATH_LEN]; /* input file name */ 300 301 //app_context->ierr = PetscMalloc(sizeof(APP_CONTEXT_TYPE), (void **) &app_context); 302 app_context = (APP_CONTEXT_TYPE*) malloc(sizeof(APP_CONTEXT_TYPE)); 303 if(app_context==NULL){ 304 PetscPrintf(MPI_COMM_WORLD, "== fatal error: app_context_create_and_initialize: app_context==NULL"); 305 PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,ERROR_TYPE_COMMON,1,"app_context_create_and_initialize: app_context==NULL"); 306 return NULL; 307 } 308 309 //app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscMalloc",app_context); 310 311 //program initialization 312 313// PetscInitialize(&argc,&args,(char *)0,help); 314 SlepcInitialize(&argc,&args,(char*)0,(char*) 0); 315 316 app_context->mpi_comm_local=MPI_COMM_SELF; 317 app_context->mpi_comm_global=MPI_COMM_WORLD; 318 app_context->pv_global=PETSC_VIEWER_STDOUT_WORLD; 319 app_context->pv_local=PETSC_VIEWER_STDOUT_SELF; 320 PetscViewerSetFormat(app_context->pv_local, PETSC_VIEWER_ASCII_MATLAB) ; //PETSC_VIEWER_ASCII_INDEX); 321 PetscViewerSetFormat(app_context->pv_global, PETSC_VIEWER_ASCII_MATLAB) ;//PETSC_VIEWER_ASCII_INDEX); 322 323 // --- program parameter initialization ---- 324 325 //default values 326 app_context->diagnosisMode = DIAGNOSIS_DISP_INITIAL_CONF | DIAGNOSIS_DISP_EVERYTHING; 327 app_context->measurement_noise_scale = 0.25; 328 app_context->global_noise_scale = 0.0001; 329// app_context->initialization_error = 0.5; 330 app_context->measurement_noise_amount = 0.0; 331 332 app_context->n_x = 2; 333 334 app_context->x_scale_constant=1.0; 335 336 PetscMalloc(sizeof(PetscReal)*app_context->n_x, (void**) &app_context->x_0); 337 app_context->x_0[0] = 1.0; 338 app_context->x_0[1] = 1.0; 339 340 app_context->n_theta = 2 ; 341 PetscMalloc(sizeof(PetscReal)*app_context->n_theta, (void**) &app_context->theta_truth); 342 app_context->theta_truth[0] = 1.0; 343 app_context->theta_truth[1] = 1.0; 344 345 PetscMalloc(sizeof(PetscReal)*app_context->n_theta, (void**) &app_context->theta_0); 346 app_context->theta_0[0] = 1.5; 347 app_context->theta_0[1] = 0.5; 348 349 PetscMalloc(sizeof(PetscReal)*app_context->n_theta, (void**) &app_context->theta_scale); 350 app_context->theta_scale[0] = 1.0; 351 app_context->theta_scale[1] = 1.0; 352 353 //user-defined values 354 app_context->ierr = PetscOptionsGetReal(PETSC_NULL,"-measurement-noise-scale",&app_context->measurement_noise_scale,&flg); 355 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetReal",app_context,__LINE__,__FUNCT__,__FILE__); 356 357 app_context->ierr = PetscOptionsGetReal(PETSC_NULL,"-global-noise-scale",&app_context->global_noise_scale,&flg); 358 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetReal",app_context,__LINE__,__FUNCT__,__FILE__); 359 360// app_context->ierr = PetscOptionsGetReal(PETSC_NULL,"-initialization-error",&app_context->initialization_error,&flg); 361// app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetReal",app_context,__LINE__,__FUNCT__,__FILE__); 362 363 app_context->ierr = PetscOptionsGetReal(PETSC_NULL,"-measurement-noise-amount",&app_context->measurement_noise_amount,&flg); 364 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetReal",app_context,__LINE__,__FUNCT__,__FILE__); 365 366 //app_context->x_0 367 app_context->ierr = PetscOptionsGetInt(PETSC_NULL,"-n-x",&app_context->n_x,&flg); 368 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetInt",app_context,__LINE__,__FUNCT__,__FILE__); 369 if(flg){ //resize the array 370 PetscFree(app_context->x_0); 371 app_context->ierr=PetscMalloc(sizeof(PetscReal)*app_context->n_x, (void **) &app_context->x_0); 372 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__); 373 } 374 flg_1 = flg; 375 PetscOptionsHasName(PETSC_NULL,"-x-0",&flg); 376 if( flg_1 | flg ){ 377 app_context_report_error(ERROR_TYPE_COMMON,"app_context_create_and_initialize: -n-x and -x_0 should be given together",app_context,__LINE__,__FUNCT__,__FILE__); 378 } 379 PetscOptionsHasName(PETSC_NULL,"-x-0",&flg); 380 if(flg){ 381 app_context->ierr = PetscOptionsGetRealArray(PETSC_NULL,"-x-0",app_context->x_0, &app_context->n_x,&flg); 382 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetInt",app_context,__LINE__,__FUNCT__,__FILE__); 383 } 384 385 //x_scale_constant 386 app_context->ierr = PetscOptionsGetReal(PETSC_NULL,"-x-scale",&app_context->x_scale_constant, &flg); 387 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetReal x_scale",app_context,__LINE__,__FUNCT__,__FILE__); 388 389 //x-scale 390 app_context->ierr = PetscMalloc(sizeof(PetscReal)*app_context->n_x, (void **) &app_context->x_scale); 391 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__); 392 for(i=0; i<app_context->n_x; i++) 393 app_context->x_scale[i]=app_context->x_scale_constant; 394 395 //app_context->theta_truth 396 app_context->ierr = PetscOptionsGetInt(PETSC_NULL,"-n-theta-truth",&app_context->n_theta,&flg); 397 app_context_assert(app_context->n_theta>=0, "", app_context); 398 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetInt",app_context,__LINE__,__FUNCT__,__FILE__); 399 if(flg){ //resize the array 400 PetscFree(app_context->x_0); 401 PetscMalloc(sizeof(PetscReal)*app_context->n_theta, (void **) &app_context->theta_truth); 402 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__); 403 } 404 flg_1 = flg; 405 PetscOptionsHasName(PETSC_NULL,"-theta-truth",&flg); 406 if( app_context->n_theta>0 && (flg_1 | flg) ){ 407 app_context_report_error(ERROR_TYPE_COMMON,"app_context_create_and_initialize: -n-theta-truth and -theta_truth should be given together",app_context,__LINE__,__FUNCT__,__FILE__); 408 } 409 if(app_context->n_theta>0 && flg){ 410 app_context->ierr = PetscOptionsGetRealArray(PETSC_NULL,"-theta-truth",app_context->theta_truth,&app_context->n_theta, &flg); 411 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetRealArray",app_context,__LINE__,__FUNCT__,__FILE__); 412 } 413 414 //app_context->theta_0 415 app_context->ierr = PetscOptionsGetInt(PETSC_NULL,"-n-theta-0",&app_context->n_theta,&flg); 416 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetInt",app_context,__LINE__,__FUNCT__,__FILE__); 417 if(flg){ //resize the array 418 if(app_context->n_theta != sizeof(app_context->theta_truth)/sizeof(PetscReal)) 419 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: -n-theta-0 is different from -n-theta-truth",app_context,__LINE__,__FUNCT__,__FILE__); 420 PetscFree(app_context->x_0); 421 PetscMalloc(sizeof(PetscReal)*app_context->n_theta, (void **) &app_context->theta_0); 422 } 423 flg_1 = flg; 424 PetscOptionsHasName(PETSC_NULL,"-theta-0",&flg); 425 if( flg_1 | flg ){ 426 app_context_report_error(ERROR_TYPE_COMMON,"app_context_create_and_initialize: -n-theta-0 and -theta_0 should be given together",app_context,__LINE__,__FUNCT__,__FILE__); 427 } 428 if(flg){ 429 app_context->ierr = PetscOptionsGetRealArray(PETSC_NULL,"-theta-0",app_context->theta_0,&app_context->n_theta, &flg); 430 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetIntArray",app_context,__LINE__,__FUNCT__,__FILE__); 431 } 432 433 //app_context->theta_scale 434 app_context->ierr = PetscOptionsGetInt(PETSC_NULL,"-n-theta-scale",&app_context->n_theta,&flg); 435 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetInt",app_context,__LINE__,__FUNCT__,__FILE__); 436 if(flg){ //resize the array 437 if(app_context->n_theta != sizeof(app_context->theta_truth)/sizeof(PetscReal)) 438 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: -n-theta-scale is different from -n-theta",app_context,__LINE__,__FUNCT__,__FILE__); 439 PetscFree(app_context->theta_scale); 440 PetscMalloc(sizeof(PetscReal)*app_context->n_theta, (void **) &app_context->theta_scale); 441 } 442 flg_1 = flg; 443 PetscOptionsHasName(PETSC_NULL,"-theta-scale",&flg); 444 if( flg_1 | flg ){ 445 app_context_report_error(ERROR_TYPE_COMMON,"app_context_create_and_initialize: n_theta_scale and -theta-scale should be given together",app_context,__LINE__,__FUNCT__,__FILE__); 446 } 447 if(flg){ 448 app_context->ierr = PetscOptionsGetRealArray(PETSC_NULL,"-theta-scale",app_context->theta_scale,&app_context->n_theta, &flg); 449 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscOptionsGetRealArray",app_context,__LINE__,__FUNCT__,__FILE__); 450 } 451 452 //echo configuration information 453 if(app_context->diagnosisMode & DIAGNOSIS_DISP_INITIAL_CONF){ 454 app_context->ierr = app_context_report_info(INFO_PETSC_REAL,1, &app_context->measurement_noise_scale,"measurement_noise_scale", app_context); 455 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 456 457 app_context->ierr = app_context_report_info(INFO_PETSC_REAL,1, &app_context->global_noise_scale,"global_noise_scale", app_context); 458 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 459 460 //app_context->ierr = app_context_report_info(INFO_PETSC_REAL,1, &app_context->initialization_error,"initialization_error", app_context); 461 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 462 463 app_context->ierr = app_context_report_info(INFO_PETSC_INT,1, &app_context->n_x,"n_x", app_context); 464 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 465 466 app_context->ierr = app_context_report_info(INFO_PETSC_REALARRAY, app_context->n_x, app_context->x_0,"x_0", app_context); 467 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 468 469 app_context->ierr = app_context_report_info(INFO_PETSC_REALARRAY, app_context->n_x, app_context->x_scale,"x_scale", app_context); 470 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 471 472 app_context->ierr = app_context_report_info(INFO_PETSC_INT,1, &app_context->n_theta,"n_theta", app_context); 473 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 474 475 app_context->ierr = app_context_report_info(INFO_PETSC_REALARRAY, app_context->n_theta, app_context->theta_truth,"theta_truth", app_context); 476 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 477 478 app_context->ierr = app_context_report_info(INFO_PETSC_REALARRAY, app_context->n_theta, app_context->theta_0,"theta_0", app_context); 479 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 480 481 app_context->ierr = app_context_report_info(INFO_PETSC_REALARRAY, app_context->n_theta, app_context->theta_scale,"theta_scale", app_context); 482 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: app_context_report_info",app_context,__LINE__,__FUNCT__,__FILE__); 483 } 484 485 //--allocate temp space 486 VecCreate(app_context->mpi_comm_global, &app_context->tmp_x); 487 VecCreate(app_context->mpi_comm_global, &app_context->tmp_theta); 488 VecCreate(app_context->mpi_comm_global, &app_context->tmp_z); 489 490 VecSetSizes(app_context->tmp_x, PETSC_DECIDE,app_context->n_x); 491 VecSetSizes(app_context->tmp_theta, PETSC_DECIDE,app_context->n_theta); 492 VecSetSizes(app_context->tmp_z, PETSC_DECIDE,app_context->n_x+app_context->n_theta); 493 app_context->ierr=PetscMalloc(sizeof(PetscInt)*app_context->n_x, &app_context->tmp_x_indices); 494 app_context_report_error(app_context->ierr,"PetscMalloc: app_context->tmp_x_indices",app_context,__LINE__,__FUNCT__,__FILE__); 495 app_context->ierr=PetscMalloc(sizeof(PetscInt)*app_context->n_theta, &app_context->tmp_theta_indices); 496 app_context_report_error(app_context->ierr,"PetscMalloc: app_context->tmp_theta_indices",app_context,__LINE__,__FUNCT__,__FILE__); 497 498 VecSetFromOptions(app_context->tmp_x); 499 VecSetFromOptions(app_context->tmp_theta); 500 VecSetFromOptions(app_context->tmp_z); 501 502 //measurement 503 app_context->measurement_dir_name = new char[MAX_SIZE_PATH_NAME+1]; 504 app_context->measurement_data_file_name = new char[MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME+1]; 505 app_context->measurement_header_file_name = new char[MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME+1]; 506 app_context->measurement_matrix_file_name = new char[MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME+1]; 507 app_context->measurement_error_covariance_matrix_file_name = new char[MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME+1]; 508 509 PetscOptionsGetString(PETSC_NULL,"-measurement_dir_name", app_context->measurement_dir_name,MAX_SIZE_PATH_NAME,&flg); 510 if(flg==PETSC_FALSE) sprintf(app_context->measurement_dir_name, "%smeasurement/", MEASUREMENT_DIR); 511 512 PetscOptionsGetString(PETSC_NULL,"-measurement_header_file_name", app_context->measurement_header_file_name,MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME,&flg); 513 if(flg==PETSC_FALSE) 514 sprintf(app_context->measurement_header_file_name, "%sheader", app_context->measurement_dir_name); 515 else if(!strchr(app_context->measurement_header_file_name,'/')) 516 sprintf(app_context->measurement_header_file_name, "%s%s", app_context->measurement_dir_name,app_context->measurement_header_file_name); 517 else{} 518 519 PetscOptionsGetString(PETSC_NULL,"-measurement_matrix_file_name", app_context->measurement_matrix_file_name,MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME,&flg); 520 if(flg==PETSC_FALSE) 521 sprintf(app_context->measurement_matrix_file_name, "%sH.dat", app_context->measurement_dir_name); 522 else if(!strchr(app_context->measurement_matrix_file_name,'/')) 523 sprintf(app_context->measurement_matrix_file_name, "%s%s", app_context->measurement_dir_name,app_context->measurement_matrix_file_name); 524 else{} 525 526 PetscOptionsGetString(PETSC_NULL,"-measurement_error_covariance_matrix_file_name", app_context->measurement_error_covariance_matrix_file_name,MAX_SIZE_FILE_NAME+MAX_SIZE_PATH_NAME,&flg); 527 if(flg==PETSC_FALSE) 528 sprintf(app_context->measurement_error_covariance_matrix_file_name, "%sR_d.dat", app_context->measurement_dir_name); 529 else if(!strchr(app_context->measurement_error_covariance_matrix_file_name,'/')) 530 sprintf(app_context->measurement_error_covariance_matrix_file_name, "%s%s", app_context->measurement_dir_name,app_context->measurement_error_covariance_matrix_file_name); 531 else{} 532 533 PetscOptionsGetInt(PETSC_NULL,"-measurementType",&app_context->measurementType,&flg); 534 if(flg==PETSC_FALSE) app_context->measurementType=MEAUSREMENT_TYPE_FROM_HEADER_AND_DAT_FILES; 535 536 return app_context; 537} 538 539int 540app_context_destroy_and_finalize(APP_CONTEXT_TYPE* app_context) 541{ 542 if(app_context==NULL){ 543 PetscPrintf(MPI_COMM_WORLD, "== fatal error: app_context_create_and_initialize: app_context==NULL"); 544 PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,ERROR_TYPE_COMMON,1,"app_context_finalize: app_context==NULL"); 545 return 1; 546 } 547 548// app_context->ierr = PetscFinalize(); CHKERRQ(app_context->ierr); 549 app_context->ierr = SlepcFinalize(); CHKERRQ(app_context->ierr); 550 551 return 0; 552} 553 554//------------------------------------------------------------------------------------------------- 555//---measurement routines--- 556 557//load single mesurement from .dat file with simple format: % dimension length(V) 1\n V \n 558// into measurement->single_measurement 559int 560measurement_io_read_single_measurement_dat(char* filename, MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 561{ 562 int measurement_size_1 = 0, measurement_size_2 = 0; 563 564 if(measurement->single_measurement) free(measurement->single_measurement); 565 measurement->single_measurement = load_matrix_2D_ascii(filename, measurement_size_1, measurement_size_2); 566 assert(measurement_size_2==1); assert(measurement_size_1 > 1); 567 measurement->n_y = measurement_size_1; 568 569 return 0; 570} 571 572//write .dat file with simple format: length(V) \n V \n 573int 574measurement_io_write_single_measurement_dat(char* filename, MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 575{ 576 int ret = 0; 577 578 if(measurement->single_measurement) ret = save_matrix_2D_ascii(filename, measurement->single_measurement, measurement->n_y, 1, NULL); 579 assert(ret==0); 580 581 return 0; 582} 583 584//load all meausrment into measurement data structure 585int 586measurement_io_read_all_measurements(MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 587{ 588 int meta_info_size_1 = 0, meta_info_size_2 = 0; 589 590 //header + dat files 591 if(app_context->measurementType==MEAUSREMENT_TYPE_FROM_HEADER_AND_DAT_FILES){ 592 593 //load header file 594 measurement->map_index2time = load_matrix_2D_ascii(app_context->measurement_header_file_name, meta_info_size_1, meta_info_size_2); 595 assert(meta_info_size_1 == 2); assert(meta_info_size_2 > 2); 596 measurement->K = meta_info_size_2; 597 598 //load each measurement 599 measurement->all_measurements = new double[measurement->n_y*measurement->K]; 600 for (int i=0; i<measurement->K; i++){ 601 sprintf(app_context->measurement_data_file_name, "%s%04d.dat", app_context->measurement_dir_name, (int)measurement->map_index2time[i*2+0]); 602 measurement_io_read_single_measurement_dat(app_context->measurement_data_file_name, measurement, app_context); 603 for (int j=0; j<measurement->n_y; j++) measurement->all_measurements[j*measurement->K+i]=measurement->single_measurement[j]; 604 } 605 } 606 607 return 0; 608} 609 610//create (Vec) measurement->y from (double*) measurement->single_measurement 611int 612measurement_petsc_vector_assemble(MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 613{ 614 static int* index = NULL; 615 if(index==NULL) { 616 index = new int[measurement->n_y]; 617 VecCreate(app_context->mpi_comm_global, &measurement->y); 618 VecSetSizes(measurement->y, PETSC_DECIDE, measurement->n_y); 619 VecSetFromOptions(measurement->y); 620 for (int i = 0; i<measurement->n_y; i++) index[i]=i; 621 } 622 app_context->ierr=VecSetValues(measurement->y,measurement->n_y, index, measurement->single_measurement,ADD_VALUES); 623 app_context_report_error(app_context->ierr,"measurement_petsc_vector_assemble: VecSetValues",app_context,__LINE__,__FUNCT__,__FILE__); 624 app_context->ierr=VecAssemblyBegin(measurement->y); 625 app_context_report_error(app_context->ierr,"measurement_petsc_vector_assemble: VecAssemblyBegin",app_context,__LINE__,__FUNCT__,__FILE__); 626 app_context->ierr=VecAssemblyEnd(measurement->y); 627 app_context_report_error(app_context->ierr,"measurement_petsc_vector_assemble: VecAssemblyEnd",app_context,__LINE__,__FUNCT__,__FILE__); 628 629 return 0; 630} 631 632//create (double*) measurement->single_measurement from (Vec) measurement->y 633int 634measurement_petsc_vector_disassemble(MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 635{ 636 static int* index = NULL; 637 if(index==NULL) { 638 index = new int[measurement->n_y]; 639 for (int i = 0; i<measurement->n_y; i++) index[i]=i; 640 } 641 app_context->ierr=VecGetValues(measurement->y,measurement->n_y, index, measurement->single_measurement); 642 app_context_report_error(app_context->ierr,"measurement_petsc_vector_disassemble: VecGetValues",app_context,__LINE__,__FUNCT__,__FILE__); 643 app_context->ierr=VecAssemblyBegin(measurement->y); 644 app_context_report_error(app_context->ierr,"measurement_petsc_vector_disassemble: VecAssemblyBegin",app_context,__LINE__,__FUNCT__,__FILE__); 645 app_context->ierr=VecAssemblyEnd(measurement->y); 646 app_context_report_error(app_context->ierr,"measurement_petsc_vector_disassemble: VecAssemblyEnd",app_context,__LINE__,__FUNCT__,__FILE__); 647 648 return 0; 649} 650 651//initialize the interpolation for all meausrements, from all_measurements and map_index2time 652int 653measurement_interpolation_initialize(MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 654{ 655 interp_v(0, measurement->n_y, measurement->K, &measurement->map_index2time[1*measurement->K+0], measurement->all_measurements, -1, NULL); 656 return 0; 657} 658 659//finalize the interpolation for all meausrements 660int 661measurement_interpolation_finalize(MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 662{ 663 interp_v(2, measurement->n_y, measurement->K, NULL, NULL, -1, NULL); 664 return 0; 665} 666 667//get the meausrement at any given time point, stored in single_measurement 668int 669measurement_get(double time, MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 670{ 671 if(1){ //interpolation 672 interp_v(1, measurement->n_y, measurement->K, NULL, NULL,time, measurement->single_measurement); 673 return 0; 674 } 675} 676 677//load measurement matrix from a data file into measurement->H_values; 678int 679measurement_io_load_measurement_matrix_dat(char* filename, MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 680{ 681 int measurement_matrix_size_1 = 0, measurement_matrix_size_2 = 0; 682 683 if(measurement->H_values){PetscPrintf(app_context->mpi_comm_global, "measurement matrix was already loaded (measurement->H_values=%x), so no need to load it from %s again\n", measurement->H_values,filename); assert(0);} 684 685 measurement->H_values = load_matrix_2D_ascii(filename, measurement_matrix_size_1, measurement_matrix_size_2); 686 assert(measurement_matrix_size_2==app_context->n_x+app_context->n_theta); 687 assert(measurement_matrix_size_1==measurement->n_y); 688 measurement->n_z = measurement_matrix_size_2; 689 690 return 0; 691} 692//load measurement error covariance matrix from a data file into measurement->R_d_values; 693int 694measurement_io_load_measurement_error_covariance_matrix_dat(char* filename, MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context) 695{ 696 int measurement_matrix_size_1 = 0, measurement_matrix_size_2 = 0; 697 698 if(measurement->R_d_values){PetscPrintf(app_context->mpi_comm_global, "measurement error covariance matrix was already loaded (measurement->R_d_values=%x), so no need to load it from %s again\n", measurement->R_d_values,filename); assert(0);} 699 700 measurement->R_d_values = load_matrix_2D_ascii(filename, measurement_matrix_size_1, measurement_matrix_size_2); 701 assert(measurement_matrix_size_2==measurement->n_y); 702 assert(measurement_matrix_size_1==measurement->n_y); 703 return 0; 704} 705//prepare measurement (IO, create petsc data structure etc.) 706MEAUSREMENT_TYPE* 707measurement_create_and_initialize(APP_CONTEXT_TYPE* app_context) 708{ 709 MEAUSREMENT_TYPE* measurement = new MEAUSREMENT_TYPE; 710 memset(measurement, 0, sizeof(MEAUSREMENT_TYPE)); 711 if(0){ 712 measurement->n_y = 2; 713 measurement->n_z = 4; 714 measurement->single_measurement=new double[measurement->n_y]; 715 return measurement; 716 } 717 718 //load all measurements 719 measurement_io_read_all_measurements(measurement, app_context); 720 721 //load measurement matrix 722 measurement_io_load_measurement_matrix_dat(app_context->measurement_matrix_file_name, measurement, app_context); 723 { //create measurement->H (Mat type) 724 int* index_row = new int[measurement->n_y]; 725 int* index_col = new int[measurement->n_z]; 726 MatCreate(app_context->mpi_comm_global, &measurement->H); 727 MatSetSizes(measurement->H, PETSC_DECIDE, PETSC_DECIDE, measurement->n_y, measurement->n_z); 728 MatSetFromOptions(measurement->H); 729 for (int i = 0; i<measurement->n_y; i++) index_row[i]=i; 730 for (int i = 0; i<measurement->n_z; i++) index_col[i]=i; 731 732 app_context->ierr=MatSetValues(measurement->H,measurement->n_y,index_row,measurement->n_z, index_col,measurement->H_values,ADD_VALUES); 733 MatAssemblyBegin(measurement->H,MAT_FINAL_ASSEMBLY); 734 MatAssemblyEnd(measurement->H,MAT_FINAL_ASSEMBLY); 735 app_context_report_error(app_context->ierr,"measurement_create_and_initialize: MatSetValues",app_context,__LINE__,__FUNCT__,__FILE__); 736 free(index_row); free(index_col); 737 } 738 739 740 //load measurement error covariance matrix 741 measurement_io_load_measurement_error_covariance_matrix_dat(app_context->measurement_error_covariance_matrix_file_name, measurement, app_context); 742 { //create measurement->R_d (Mat type) 743 int* index_row = new int[measurement->n_y]; 744 MatCreateMPIDense(app_context->mpi_comm_global, PETSC_DECIDE, PETSC_DECIDE, measurement->n_y, measurement->n_y, PETSC_NULL, &measurement->R_d); 745 //MatCreate(app_context->mpi_comm_global, &measurement->R_d); 746 //MatSetSizes(measurement->R_d, PETSC_DECIDE, PETSC_DECIDE, measurement->n_y, measurement->n_y); 747 MatSetFromOptions(measurement->R_d); 748 for (int i = 0; i<measurement->n_y; i++) index_row[i]=i; 749 750 app_context->ierr=MatSetValues(measurement->R_d,measurement->n_y,index_row,measurement->n_y, index_row,measurement->R_d_values,ADD_VALUES); 751 MatAssemblyBegin(measurement->R_d,MAT_FINAL_ASSEMBLY); 752 MatAssemblyEnd(measurement->R_d,MAT_FINAL_ASSEMBLY); 753 app_context_report_error(app_context->ierr,"measurement_create_and_initialize: MatSetValues",app_context,__LINE__,__FUNCT__,__FILE__); 754 free(index_row); 755 } 756 757 //Create measurement error covariance matrix inverse inv_R_d 758 MatDuplicate(measurement->R_d, MAT_COPY_VALUES, &measurement->inv_R_d); 759 760 //prepare the interpolation 761 measurement_interpolation_initialize(measurement, app_context); 762 763 //usage: 764 //ret = measurement_get(0, measurement, app_context); 765 //if(ret==0) measurement_petsc_vector_assemble(measurement,app_context); 766 767 //load the measurement error covariance R_d 768 769 770 return measurement; 771} 772 773 774//------------------------------------------------------------------------------------------------- 775//---models routines--- 776MODEL_TYPE* 777model_create_and_initialize(MODEL_TYPE* (*model_opeartor)(MODEL_TYPE*, APP_CONTEXT_TYPE* ), APP_CONTEXT_TYPE* app_context) 778{ 779 MODEL_TYPE* model = NULL; 780 app_context->ierr = PetscMalloc(sizeof(MODEL_TYPE), (void**) &model); 781 app_context_report_error(app_context->ierr,"app_context_create_and_initialize: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__); 782 783 model_set_opeartor_and_intial_condition(model, model_opeartor,app_context); 784 785 return model; 786} 787 788 789// Initialize 790//1) model operator (model->f), 791//2) dimension (model->n_x, model->n_theta), 792//3) time stamp (), 793//4) state (model->x_k_hat,model->theta_k_hat) , 794//5) previous state (model->x_k_minus_1_hat,model->theta_k_minus_1_hat) 795MODEL_TYPE* 796model_set_opeartor_and_intial_condition(MODEL_TYPE* model,MODEL_TYPE* (*model_opeartor)(MODEL_TYPE*, APP_CONTEXT_TYPE*), APP_CONTEXT_TYPE* app_context) 797{ 798 PetscInt i=0; 799 800 if(1){ //model 01: 2 parameter, 2 state, simplest model 801 802 model->f = *model_opeartor; 803 model->k = 0; 804 model->t = 0.0; 805 806 model->n_x=app_context->n_x; 807// app_context->ierr=VecCreateMPIWithArray(app_context->mpi_comm_global, model->n_x, PETSC_DECIDE, app_context->x_0, &model->x_k_hat); 808// app_context->ierr=VecCreateMPIWithArray(app_context->mpi_comm_global, model->n_x, model->n_x, app_context->x_0, &model->x_k_hat); 809 810 app_context->ierr=VecCreate(app_context->mpi_comm_global, &model->x_k_hat); 811 app_context->ierr=VecSetSizes(model->x_k_hat, PETSC_DECIDE,app_context->n_x); 812 app_context->ierr=VecSetFromOptions(model->x_k_hat); 813 for(i=0;i<model->n_x; i++) app_context->tmp_x_indices[i]=i; 814 VecSetValues(model->x_k_hat,model->n_x,app_context->tmp_x_indices,app_context->x_0,ADD_VALUES); 815 VecAssemblyBegin(model->x_k_hat); 816 VecAssemblyEnd(model->x_k_hat); 817 818 app_context_report_error(app_context->ierr,"model_set_opeartor_and_intial_condition: VecCreateMPIWithArray",app_context,__LINE__,__FUNCT__,__FILE__); 819 820 model->n_theta=app_context->n_theta; 821//app_context->ierr=VecCreateMPIWithArray(app_context->mpi_comm_global, model->n_theta, PETSC_DECIDE, app_context->theta_0, &model->theta_k_hat); 822// app_context->ierr=VecCreateMPIWithArray(app_context->mpi_comm_global, model->n_theta, model->n_theta, app_context->theta_0, &model->theta_k_hat); 823 app_context->ierr=VecCreate(app_context->mpi_comm_global, &model->theta_k_hat); 824 app_context->ierr=VecSetSizes(model->theta_k_hat, PETSC_DECIDE,app_context->n_theta); 825 app_context->ierr=VecSetFromOptions(model->theta_k_hat); 826 for(i=0;i<model->n_theta; i++) app_context->tmp_theta_indices[i]=i; 827 VecSetValues(model->theta_k_hat,model->n_theta,app_context->tmp_theta_indices,app_context->theta_0,ADD_VALUES); 828 VecAssemblyBegin(model->theta_k_hat); 829 VecAssemblyEnd(model->theta_k_hat); 830 app_context_report_error(app_context->ierr,"model_set_opeartor_and_intial_condition: VecCreateMPIWithArray",app_context,__LINE__,__FUNCT__,__FILE__); 831 832 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){ 833 app_context_report_info(INFO_PETSC_REALARRAY,app_context->n_x, app_context->x_0,"model_set_opeartor_and_intial_condition, app_context->x_0", app_context); 834 app_context_report_info(INFO_PETSC_REALARRAY, app_context->n_theta, app_context->theta_0,"model_set_opeartor_and_intial_condition, app_context->theta_0", app_context); 835 app_context_report_info(INFO_PETSC_VEC,1, &model->x_k_hat,"model_set_opeartor_and_intial_condition, model->x_0_hat", app_context); 836 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k_hat,"model_set_opeartor_and_intial_condition, model->theta_0_hat", app_context); 837 } 838 839 840 }else if(0){ //model 02: 841 842 app_context_report_error(app_context->ierr,"model_set_opeartor_and_intial_condition: not implemented yet",app_context,__LINE__,__FUNCT__,__FILE__); 843 }else{ 844 845 app_context_report_error(app_context->ierr,"model_set_opeartor_and_intial_condition: not implemented yet",app_context,__LINE__,__FUNCT__,__FILE__); 846 } 847 848 //other standard setups 849 VecDuplicate(model->x_k_hat, &model->x_k_minus_1_hat); 850 VecDuplicate(model->theta_k_hat, &model->theta_k_minus_1_hat); 851 852 //diagnosis 853 if(app_context->diagnosisMode & DIAGNOSIS_DISP_INITIAL_CONF){ 854 app_context_report_info(INFO_PETSC_REAL, 1, &model->t,"model_set_opeartor_and_intial_condition, model->t", app_context); 855 app_context_report_info(INFO_PETSC_INT, 1, &model->k,"model_set_opeartor_and_intial_condition, model->k", app_context); 856 857 app_context_report_info(INFO_PETSC_INT, 1, &model->n_x,"model_set_opeartor_and_intial_condition, model->n_x", app_context); 858 app_context_report_info(INFO_PETSC_INT, 1, &model->n_theta,"model_set_opeartor_and_intial_condition, model->n_theta", app_context); 859 860 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k_hat,"model_set_opeartor_and_intial_condition, model->theta_0_hat", app_context); 861 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k_minus_1_hat,"model_set_opeartor_and_intial_condition, model->theta_0_minus_1_hat", app_context); 862 } 863 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){ 864 app_context_report_info(INFO_PETSC_VEC,1, &model->x_k_hat,"model_set_opeartor_and_intial_condition, model->x_0_hat", app_context); 865 app_context_report_info(INFO_PETSC_VEC, 1, &model->x_k_minus_1_hat,"model_set_opeartor_and_intial_condition, model->x_0_minus_1_hat", app_context); 866 } 867 868 return model; 869} 870 871 872//------------------------------------------------------------------------------------------------- 873//---filter routines--- 874 875//create and initialize the filter 876FILTER_TYPE* 877filter_create_and_initialize(APP_CONTEXT_TYPE* app_context,MODEL_TYPE* model) 878{ 879 FILTER_TYPE* filter = NULL; 880 app_context->ierr = PetscMalloc(sizeof(FILTER_TYPE), (void**) &filter); 881 app_context_report_error(app_context->ierr,"filter_create_and_initialize: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__); 882 883 filter->n_x = app_context->n_x; 884 filter->n_theta = app_context->n_theta; 885 filter = filter_set_intial_condition(filter,model, app_context); 886 887 return filter; 888} 889 890 891// Initialize 892//1) covariance for state, measurement (portion_afilter->), 893//2) dimension (filter->n_x, filter->n_theta), 894//3) time stamp (filter->), 895//4) state (K) 896FILTER_TYPE* 897filter_set_intial_condition(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context) 898{ 899 int i=0; 900 901 //allocate space and initialize covariance matrix (filter subtype) 902 filter->subfilter = filter_subtype_create_and_initialize(filter,model, app_context); 903 904 //model->n_ensem = 2; 905 //set in filter_subtype_set_intial_condition() 906 907 //allocate space for sigma points (in model struture) 908 PetscMalloc(sizeof(Vec*)*model->n_ensem, &model->x_k_minus_1); 909 PetscMalloc(sizeof(Vec*)*model->n_ensem, &model->theta_k_minus_1); 910 PetscMalloc(sizeof(Vec*)*model->n_ensem, &model->x_k); 911 PetscMalloc(sizeof(Vec*)*model->n_ensem, &model->theta_k); 912 913 for(i=0; i<model->n_ensem; i++){ 914 VecCreate(app_context->mpi_comm_global, &model->x_k_minus_1[i]); 915 VecCreate(app_context->mpi_comm_global, &model->theta_k_minus_1[i]); 916 VecCreate(app_context->mpi_comm_global, &model->x_k[i]); 917 VecCreate(app_context->mpi_comm_global, &model->theta_k[i]); 918 919 VecSetSizes(model->x_k_minus_1[i], PETSC_DECIDE,app_context->n_x); 920 VecSetSizes(model->theta_k_minus_1[i], PETSC_DECIDE,app_context->n_theta); 921 VecSetSizes(model->x_k[i], PETSC_DECIDE,app_context->n_x); 922 VecSetSizes(model->theta_k[i], PETSC_DECIDE,app_context->n_theta); 923 924 VecSetFromOptions(model->x_k_minus_1[i]); 925 VecSetFromOptions(model->theta_k_minus_1[i]); 926 VecSetFromOptions(model->x_k[i]); 927 VecSetFromOptions(model->theta_k[i]); 928 } 929 930 return filter; 931} 932 933//create and initialize the filter 934FILTER_SUBTYPE* 935filter_subtype_create_and_initialize(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context) 936{ 937 filter->subfilter = NULL; 938 app_context->ierr = PetscMalloc(sizeof(FILTER_SUBTYPE), (void**) &filter->subfilter); 939 app_context_report_error(app_context->ierr,"filter_set_intial_condition: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__); 940 941 filter->subfilter = filter_subtype_set_intial_condition(filter,model, app_context); 942 return filter->subfilter; 943} 944 945// Initialize the filter subtype 946// PetscInt filter_type; 947// PetscInt r; //rank ( r>=(n_theat>0)?n_theat:1, r<= n_theta+n_x; 948// Mat S_x; 949// Mat S_theta; 950FILTER_SUBTYPE* 951filter_subtype_set_intial_condition(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context) 952{ 953 //local values 954 PetscScalar S_x_value; //all 1 default 955 PetscScalar* S_theta_values; //1,1,1... default 956 PetscInt i; //temperory usage 957 958 //check pointer 959 app_context_report_error(filter==NULL,"FILTER_TYPE* filter is null",app_context,__LINE__,__FUNCT__,__FILE__); 960 app_context_report_error(model==NULL,"MODEL_TYPE* model is null",app_context,__LINE__,__FUNCT__,__FILE__); 961 962 FILTER_SUBTYPE* subfilter=filter->subfilter; 963 app_context_report_error(subfilter==NULL,"filter->subfilter is null",app_context,__LINE__,__FUNCT__,__FILE__); 964 965 //initialize the subfilter 966 subfilter->filter_type = FILTER_TYPE_RUKF; //hard-code rUKF 967 968 // Initialize 969 // Mat S_x; 970 // Mat S_theta; 971 switch(subfilter->filter_type) 972 { 973 case FILTER_TYPE_RUKF: //reduced-order UKF 974 //choose rank 975 if(model->n_theta==0){ //state estimation only 976 subfilter->r = model->n_x - 1; 977 }//chanllenge it by setting it to be not equal to n_x 978 else{ //joint estimation 979 subfilter->r = model->n_theta + 0 ; //chanllenge it by setting it to be not equal to n_theta 980 } 981 model->n_ensem = subfilter->r * 2; 982 983 //allocate space for S_x, S_theta; app_context->S_ensemble, app_context->P_z, app_context->S_z 984 app_context->ierr=MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, model->n_x,subfilter->r, PETSC_NULL, &subfilter->S_x); 985 app_context_report_error(app_context->ierr,"subfilter->S_x MatCreateMPIDense",app_context,__LINE__,__FUNCT__,__FILE__); 986 app_context->ierr=MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, model->n_theta,subfilter->r, PETSC_NULL, &subfilter->S_theta); 987 app_context_report_error(app_context->ierr,"subfilter->S_theta MatCreateMPIDense",app_context,__LINE__,__FUNCT__,__FILE__); 988 app_context->ierr=MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, model->n_x+model->n_theta,model->n_ensem, PETSC_NULL, &app_context->S_ensemble); 989 app_context_report_error(app_context->ierr,"app_context->S_ensemble MatCreateMPIDense",app_context,__LINE__,__FUNCT__,__FILE__); 990 app_context->ierr=MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, model->n_x+model->n_theta,model->n_ensem, PETSC_NULL, &app_context->S_ensemble2); 991 app_context_report_error(app_context->ierr,"app_context->S_ensemble2 MatCreateMPIDense",app_context,__LINE__,__FUNCT__,__FILE__); 992 app_context->ierr=MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, model->n_x+model->n_theta,model->n_x+model->n_theta, PETSC_NULL, &app_context->P_z); 993 app_context_report_error(app_context->ierr,"app_context->P_z MatCreateMPIDense",app_context,__LINE__,__FUNCT__,__FILE__); 994 995 app_context->ierr=MatCreateMPIDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, model->n_x+model->n_theta,subfilter->r, PETSC_NULL, &app_context->S_z); 996 app_context_report_error(app_context->ierr,"app_context->S_z MatCreateMPIDense",app_context,__LINE__,__FUNCT__,__FILE__); 997 MatAssemblyBegin(app_context->S_z ,MAT_FINAL_ASSEMBLY ); //for later svd usage 998 MatAssemblyEnd (app_context->S_z ,MAT_FINAL_ASSEMBLY ); //for later svd usage 999 1000 MatSetFromOptions( subfilter->S_x ); 1001 MatSetFromOptions( subfilter->S_theta ); 1002 MatSetFromOptions( app_context->S_ensemble ); 1003 MatSetFromOptions( app_context->P_z ); 1004 MatSetFromOptions( app_context->S_z ); 1005 1006 1007 //initialization S_x and S_theta 1008 //values 1009 app_context->ierr=PetscMalloc(sizeof(PetscScalar)*model->n_theta, &S_theta_values); 1010 for(i=0;i<model->n_theta; i++){ 1011 S_theta_values[i] = 1.0; 1012 } 1013 S_x_value = 1.0; 1014 1015 //set (set diagonal element of S_theta first) 1016 for(i=0; i<model->n_theta; i++){ 1017 MatSetValue(subfilter->S_theta, i,i, S_theta_values[i], INSERT_VALUES); 1018 } 1019 MatAssemblyBegin(subfilter->S_theta,MAT_FINAL_ASSEMBLY ); 1020 MatAssemblyEnd (subfilter->S_theta,MAT_FINAL_ASSEMBLY ); 1021 1022 //(then set diagonal element of S_x ) 1023 for(i=0; i < subfilter->r-model->n_theta; i++){ 1024 MatSetValue(subfilter->S_x, i, i+model->n_theta, S_x_value, INSERT_VALUES); 1025 } 1026 MatAssemblyBegin(subfilter->S_x,MAT_FINAL_ASSEMBLY ); 1027 MatA…
Large files files are truncated, but you can click here to view the full file