PageRenderTime 120ms CodeModel.GetById 23ms app.highlight 85ms RepoModel.GetById 1ms app.codeStats 1ms

/DataAssimilation/common/data_assimilation_routines.h

http://github.com/xyan075/examples
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