/DataAssimilation/common/data_assimilation_routines.h
C Header | 1398 lines | 915 code | 248 blank | 235 comment | 87 complexity | a5c27bede34902184ef0e3c81debb85c MD5 | raw 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 MatAssemblyEnd (subfilter->S_x,MAT_FINAL_ASSEMBLY );
1028
1029 //diagnosis
1030 if(app_context->diagnosisMode & DIAGNOSIS_DISP_INITIAL_CONF){
1031 app_context_report_info(INFO_PETSC_MAT, 1, &subfilter->S_x,"subfilter->S_x_0", app_context);
1032 app_context_report_info(INFO_PETSC_MAT, 1, &subfilter->S_theta,"subfilter->S_theta_0", app_context);
1033 }
1034
1035 break;
1036 default: //others
1037 app_context_report_error(1,"filter->filter_type not implemented",app_context,__LINE__,__FUNCT__,__FILE__);
1038 }
1039
1040 return subfilter;
1041}
1042
1043void
1044app_context_create_and_initialize_dependently(APP_CONTEXT_TYPE* app_context, FILTER_TYPE* filter, MODEL_TYPE* model, MEAUSREMENT_TYPE* measurement)
1045{
1046 //S_bar
1047 app_context->S_bar = petsc_create_matrix(1, measurement->n_y, filter->subfilter->r,app_context->mpi_comm_global);
1048 //S_bar2
1049 MatDuplicate(app_context->S_bar, MAT_COPY_VALUES, &app_context->S_bar2);
1050 //tmp_r
1051 app_context->tmp_r = petsc_create_matrix(1, filter->subfilter->r, filter->subfilter->r,app_context->mpi_comm_global);
1052 //tmp_nz_by_r
1053 app_context->tmp_nz_by_r = petsc_create_matrix(1, app_context->n_x+app_context->n_theta, filter->subfilter->r,app_context->mpi_comm_global);
1054 //tmp_r_by_ny
1055 app_context->tmp_r_by_ny = petsc_create_matrix(1, filter->subfilter->r, measurement->n_y, app_context->mpi_comm_global);
1056 //tmp_y
1057 app_context->tmp_y = petsc_create_vector(1, measurement->n_y, app_context->mpi_comm_global);
1058
1059 //Kalman gain matrix: K
1060 filter->K = petsc_create_matrix(1, filter->n_x+filter->n_theta,measurement->n_y,app_context->mpi_comm_global);
1061
1062}
1063
1064// generate sigma points from filter->covariance and model->x_k_hat
1065MODEL_TYPE*
1066filter_generate_sigma_points(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context)
1067{
1068 int i;
1069 PetscScalar localScale; //localization parameter
1070
1071 app_context_assert(model->n_ensem>=1, "", app_context);
1072 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1073 app_context_report_info(INFO_PETSC_INT, 1, &model->k,"filter_generate_sigma_points, model->k", app_context);
1074 app_context_report_info(INFO_PETSC_REAL, 1, &model->t,"filter_generate_sigma_points, model->t", app_context);
1075 app_context_report_info(INFO_PETSC_INT, 1, &i,"i", app_context);
1076 app_context_report_info(INFO_PETSC_VEC, 1, &model->x_k_hat,"filter_generate_sigma_points, model->x_k_hat", app_context);
1077 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k_hat,"filter_generate_sigma_points, model->theta_k_hat", app_context);
1078 app_context_report_info(INFO_PETSC_MAT, 1, &filter->subfilter->S_x,"/*filter_generate_sigma_points, filter->subfilter->S_x*/", app_context);
1079 app_context_report_info(INFO_PETSC_MAT, 1, &filter->subfilter->S_theta,"filter_generate_sigma_points, filter->subfilter->S_theta", app_context);
1080 }
1081
1082 if(0){ //temp: assume filter->covariance is 0
1083 for(i=0; i<model->n_ensem; i++){ //model->n_ensem>=1
1084 VecCopy(model->x_k_hat, model->x_k[i]);
1085 VecCopy(model->theta_k_hat, model->theta_k[i]);
1086 }
1087 }
1088
1089 if(1){ //generate sigma points from filter->covariance and model->x_k_hat
1090 if(filter->subfilter->filter_type!=FILTER_TYPE_RUKF){
1091 //calculate the square-root matrix of filter->covariance
1092 //1. allocate memory
1093 //2. do svd, storge them in filter->subfilter->S_x, filter->subfilter->S_theta
1094 }
1095
1096 //3. piling up first r eig vectors weighted by eig values
1097 localScale = 1; //hard coded 1, localization parameter
1098 for(i=0; i<filter->subfilter->r; i++){
1099 MatGetColumnVector(filter->subfilter->S_x,app_context->tmp_x,i);
1100 VecWAXPY(model->x_k[i], 1*localScale, app_context->tmp_x, model->x_k_hat);
1101 VecWAXPY(model->x_k[i+filter->subfilter->r], -1*localScale, app_context->tmp_x, model->x_k_hat);
1102
1103 //TODO!!Arguments are incompatible!
1104 //[0]PETSC ERROR: Incompatible vector global lengths!
1105 //solution: generalize MatGetColumnVector-->get the global vector instead of only the local portion!
1106
1107 MatGetColumnVector(filter->subfilter->S_theta,app_context->tmp_theta,i);
1108 VecWAXPY(model->theta_k[i], 1*localScale, app_context->tmp_theta, model->theta_k_hat);
1109 VecWAXPY(model->theta_k[i+filter->subfilter->r], -1*localScale, app_context->tmp_theta, model->theta_k_hat);
1110
1111 //diagnosis
1112 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1113 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k[i],"filter_generate_sigma_points, model->theta_k[i]", app_context);
1114 app_context_report_info(INFO_PETSC_VEC,1, &model->x_k[i],"filter_generate_sigma_points, model->x_k[i]", app_context);
1115 }
1116 }
1117 }
1118
1119 return model;
1120}
1121
1122FILTER_TYPE*
1123filter_time_update(FILTER_TYPE* filter, MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context)
1124{
1125 int i=0;
1126 PetscErrorCode ierr;
1127 PetscInt nEigs; //number of eig vectors to be retrieved
1128
1129 //sigma points
1130 filter_generate_sigma_points(filter, model, app_context);
1131
1132 //model evaluation
1133 model->f(model, app_context);
1134
1135 //compute mean
1136 ierr=VecAXPY(model->x_k_hat,-1, model->x_k_hat); //nullize
1137 for(i=0; i<model->n_ensem; i++){
1138 ierr=VecAXPY(model->x_k_hat,1, model->x_k[i]);
1139 }
1140 ierr=VecScale(model->x_k_hat, 1.0/model->n_ensem );
1141
1142 ierr=VecAXPY(model->theta_k_hat,-1, model->theta_k_hat);//nullize
1143 for(i=0; i<model->n_ensem; i++){
1144 ierr=VecAXPY(model->theta_k_hat,1, model->theta_k[i]);
1145 }
1146 ierr=VecScale(model->theta_k_hat, 1.0/model->n_ensem );
1147
1148
1149 //1. compute covariance (filter->subfilter->S_x and filter->subfilter->S_theta)
1150 //todo: add inflation parameter later
1151 //1.1. join them together (use MatSetColumnVec)
1152 for(i=0; i<model->n_ensem; i++){
1153 ierr=VecWAXPY(app_context->tmp_x, -1.0, model->x_k_hat, model->x_k[i]); app_context_report_error(ierr,"VecWAXPY",app_context,__LINE__,__FUNCT__,__FILE__);
1154 ierr=MatSetColumnVec(&app_context->S_ensemble, &app_context->tmp_x, app_context->tmp_x_indices, i, 0, 0, model->n_x);
1155 ierr=VecWAXPY(app_context->tmp_theta, -1.0, model->theta_k_hat, model->theta_k[i]);
1156 app_context_report_error(ierr,"VecWAXPY",app_context,__LINE__,__FUNCT__,__FILE__);
1157 ierr=MatSetColumnVec(&app_context->S_ensemble, &app_context->tmp_theta, app_context->tmp_theta_indices, i, model->n_x, 0, model->n_theta);
1158 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1159 app_context_report_info(INFO_PETSC_VEC,1, &model->x_k[i],"compute covariance, filter_time_update, model->x_k[i]", app_context);
1160 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k[i],"compute covariance, filter_time_update, model->theta_k[i]", app_context);
1161 }
1162 }
1163 MatAssemblyBegin(app_context->S_ensemble,MAT_FINAL_ASSEMBLY );
1164 MatAssemblyEnd (app_context->S_ensemble,MAT_FINAL_ASSEMBLY );
1165 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1166 app_context_report_info(INFO_PETSC_MAT, 1, &app_context->S_ensemble,"compute covariance, filter_time_update, app_context->S_ensemble", app_context);
1167 }
1168
1169 //1.2. P_z=S_ensemble*S_ensemble'
1170 MatCopy(app_context->S_ensemble,app_context->S_ensemble2,SAME_NONZERO_PATTERN);
1171 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1172 app_context_report_info(INFO_PETSC_MAT, 1, &app_context->S_ensemble,"compute covariance, filter_time_update, app_context->S_ensemble", app_context);
1173 app_context_report_info(INFO_PETSC_MAT, 1, &app_context->S_ensemble2,"compute covariance, filter_time_update, app_context->S_ensemble2", app_context);
1174 }
1175 MatMultilXYT(app_context->S_ensemble, app_context->S_ensemble2, &app_context->P_z, app_context->mpi_comm_global);
1176 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1177 app_context_report_info(INFO_PETSC_MAT, 1, &app_context->P_z,"compute covariance, filter_time_update, app_context->P_z", app_context);
1178 }
1179//
1180 //2. do svd
1181 nEigs = filter->subfilter->r;
1182 MatSVD(&app_context->P_z, app_context->mpi_comm_global, &nEigs, NULL, &app_context->S_z, NULL,1);
1183 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1184 app_context_report_info(INFO_PETSC_INT, 1, &nEigs,"number of computed nEigs vectors, filter_time_update, nEigs", app_context);
1185 app_context_report_info(INFO_PETSC_MAT, 1, &app_context->S_z,"result of SVD, filter_time_update, app_context->S_z", app_context);
1186 }
1187 if(nEigs!=filter->subfilter->r){
1188 PetscPrintf(app_context->mpi_comm_global, "app_context->P_z only has %d eig vectors (filter->subfilter->r=%d)!!\n", nEigs, filter->subfilter->r);
1189 assert(0); //static rank
1190 }
1191
1192 //3. separate them (use MatSetColumnVec or MatSetMatBlock)
1193 MatSetMatBlock(app_context->mpi_comm_global, &filter->subfilter->S_x, app_context->S_z, 0, 0, 0, 0, model->n_x, filter->subfilter->r);
1194 MatSetMatBlock(app_context->mpi_comm_global, &filter->subfilter->S_theta, app_context->S_z, 0, 0, model->n_x, 0, model->n_theta,filter->subfilter->r);
1195 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1196 app_context_report_info(INFO_PETSC_MAT, 1, &filter->subfilter->S_x,"S_x, after SVD, filter_time_update, filter->subfilter->S_x", app_context);
1197 app_context_report_info(INFO_PETSC_MAT, 1, &filter->subfilter->S_theta,"S_theta, after SVD, filter_time_update, filter->subfilter->S_theta", app_context);
1198 app_context_report_info(INFO_PETSC_INT, 1, &nEigs,"==exit of filter_time_update==", app_context);
1199 }
1200
1201 app_context_report_error(ierr,"exit of filter_time_update",app_context,__LINE__,__FUNCT__,__FILE__);
1202 return filter;
1203}
1204
1205FILTER_TYPE*
1206filter_measurement_update(FILTER_TYPE* filter, MODEL_TYPE* model, MEAUSREMENT_TYPE* measurement, APP_CONTEXT_TYPE* app_context)
1207{
1208 PetscErrorCode ierr=0;
1209 app_context_report_error(ierr,"entry of filter_measurement_update",app_context,__LINE__,__FUNCT__,__FILE__);
1210
1211 if(0){ //generate synthetic data
1212 VecDuplicate(model->x_k_hat,&measurement->y);
1213 measurement_petsc_vector_disassemble(measurement,app_context);
1214 sprintf(app_context->measurement_data_file_name, "%s%04d.dat", app_context->measurement_dir_name, model->k); //assume model->k==model->t
1215 measurement_io_write_single_measurement_dat(app_context->measurement_data_file_name, measurement, app_context);
1216 return filter;
1217 }
1218
1219 //2. covariance update and gain matrix caculation----------
1220 // S_bar = H*S_minus;
1221 // inv_R_d = inv(R_d);
1222 // tmp = eye(npdof,npdof) + S_bar'*inv_R_d*S_bar;
1223 // inv_tmp = inv(tmp);
1224 // try
1225 // S_plus = S_minus * chol(inv_tmp,'lower');
1226 // catch ME
1227 // disp('exception:Matrix must be positive definite with real diagonal.');
1228 // keyboard;
1229 // end
1230 // P_plus = S_plus * S_plus';
1231 // K_x = S_minus * inv_tmp * S_bar'*inv_R_d;
1232
1233// S_minus (app_context->S_z) = [filter->subfilter->S_x; filter->subfilter->S_theta];
1234 for(int i=0; i<filter->subfilter->r; i++){
1235 ierr=MatGetColumnVector(filter->subfilter->S_x, app_context->tmp_x, i);
1236 app_context_report_error(ierr,"",app_context,__LINE__,__FUNCT__,__FILE__);
1237 ierr=MatSetColumnVec(&app_context->S_z, &app_context->tmp_x, app_context->tmp_x_indices, i, 0, 0, model->n_x);
1238 app_context_report_error(ierr,"",app_context,__LINE__,__FUNCT__,__FILE__);
1239 ierr=MatGetColumnVector(filter->subfilter->S_theta, app_context->tmp_theta, i);
1240 app_context_report_error(ierr,"",app_context,__LINE__,__FUNCT__,__FILE__);
1241 ierr=MatSetColumnVec(&app_context->S_ensemble, &app_context->tmp_theta, app_context->tmp_theta_indices, i, model->n_x, 0, model->n_theta);
1242 app_context_report_error(ierr,"",app_context,__LINE__,__FUNCT__,__FILE__);
1243 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1244 app_context_report_info(INFO_PETSC_VEC,1, &app_context->tmp_x,"filter_measurement_update, filter->subfilter->S_x(i,:)", app_context);
1245 app_context_report_info(INFO_PETSC_VEC, 1, &app_context->tmp_theta,"filter_measurement_update, filter->subfilter->S_theta(i,:)", app_context);
1246 }
1247 }
1248 MatAssemblyBegin(app_context->S_z,MAT_FINAL_ASSEMBLY );
1249 MatAssemblyEnd (app_context->S_z,MAT_FINAL_ASSEMBLY );
1250 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1251 app_context_report_info(INFO_PETSC_MAT, 1, &app_context->S_z,"filter_measurement_update, app_context->S_z", app_context);
1252 }
1253
1254// S_bar(app_context->S_bar) = H*S_minus;
1255 MatMultilXY(measurement->H, app_context->S_z, &app_context->S_bar, app_context->mpi_comm_global);
1256
1257// inv_R_d = inv(R_d);
1258 MatInverse(measurement->R_d, &measurement->inv_R_d, app_context->mpi_comm_global);
1259// tmp = eye(npdof,npdof) + S_bar'*inv_R_d*S_bar;
1260 MatMultilXY(measurement->inv_R_d, app_context->S_bar, &app_context->S_bar2, app_context->mpi_comm_global);
1261 MatMultilXTY(app_context->S_bar, app_context->S_bar2, &app_context->tmp_r, app_context->mpi_comm_global);
1262 MatShift(app_context->tmp_r,1);
1263// inv_tmp = inv(tmp);
1264 Mat tmp = app_context->tmp_r;
1265 MatInverse(app_context->tmp_r, &app_context->tmp_r, app_context->mpi_comm_global);
1266 MatDestroy(tmp);
1267// K_x = S_minus * inv_tmp * S_bar'*inv_R_d;
1268 MatMultilXY(app_context->S_z, app_context->tmp_r, &app_context->tmp_nz_by_r, app_context->mpi_comm_global);
1269 MatMultilXTY(app_context->S_bar, measurement->inv_R_d, &app_context->tmp_r_by_ny, app_context->mpi_comm_global);
1270 MatMultilXY(app_context->tmp_nz_by_r,app_context->tmp_r_by_ny, &filter->K, app_context->mpi_comm_global);
1271// S_plus(app_context->S_z) = S_minus * chol(inv_tmp,'lower');
1272 MatFactorize(app_context->tmp_r, app_context->mpi_comm_global, 1);
1273 MatMultilXY (app_context->S_z, app_context->tmp_r, &app_context->tmp_nz_by_r, app_context->mpi_comm_global);
1274 MatCopy(app_context->tmp_nz_by_r, app_context->S_z, SAME_NONZERO_PATTERN);
1275// [filter->subfilter->S_x; filter->subfilter->S_theta]=S_minus (app_context->S_z)
1276 MatSetMatBlock(app_context->mpi_comm_global, &filter->subfilter->S_x , app_context->S_z,0,0,0, 0, model->n_x, filter->subfilter->r);
1277 MatSetMatBlock(app_context->mpi_comm_global, &filter->subfilter->S_theta, app_context->S_z,0,0,model->n_x,0, model->n_theta, filter->subfilter->r);
1278
1279//(ignore) P_plus = S_plus * S_plus';
1280
1281// clean up
1282 MatDestroy(measurement->inv_R_d);
1283
1284
1285 //3. update state---------------------------------------
1286 // y = model.u_t(:); //synthetic measurement
1287 // correction =measurementAmplifyFactor * K_x*(H*z-y);
1288 // tmp = model.parameters(:) + correction(ndof+1:nadof);
1289
1290// y = model.u_t(:); //synthetic measurement
1291 //get the measurement
1292 int ret = measurement_get(model->k, measurement, app_context);
1293 if(ret==0) measurement_petsc_vector_assemble(measurement,app_context);
1294
1295// z(app_context->tmp_z) = [model->x_k_hat;model->theta_k_hat];
1296 VecSetVecBlock(app_context->mpi_comm_global, &app_context->tmp_z, model->x_k_hat, 0, 0, model->n_x);
1297 VecSetVecBlock(app_context->mpi_comm_global, &app_context->tmp_z, model->theta_k_hat, model->n_x, 0, model->n_theta);
1298// H*z
1299 MatMultilXY(measurement->H, app_context->S_z, &app_context->S_bar, app_context->mpi_comm_global);
1300
1301// correction(app_context->tmp_z) =measurementAmplifyFactor * K_x*(H*z-y);
1302 VecScale(measurement->y, -1);
1303 MatMultAdd(measurement->H, app_context->tmp_z,measurement->y,app_context->tmp_y);
1304 MatMult (filter->K, app_context->tmp_y, app_context->tmp_z);
1305// [model->x_k_hat;model->theta_k_hat] = z(app_context->tmp_z);
1306 VecSetVecBlock(app_context->mpi_comm_global, &model->x_k_hat, app_context->tmp_z, 0, 0, model->n_x);
1307 VecSetVecBlock(app_context->mpi_comm_global, &model->theta_k_hat, app_context->tmp_z, 0, model->n_x, model->n_theta);
1308
1309 //4. filter monitor---------------------------------------
1310 if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
1311 app_context_report_info(INFO_PETSC_VEC, 1, &model->x_k_hat,"x_k_hat, after filter_measurement_update, model->x_k_hat", app_context);
1312 app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k_hat,"theta_k_hat, after filter_measurement_update, model->theta_k_hat", app_context);
1313 app_context_report_info(INFO_PETSC_MAT, 1, &filter->subfilter->S_x,"S_x, after filter_measurement_update, filter->subfilter->S_x", app_context);
1314 app_context_report_info(INFO_PETSC_MAT, 1, &filter->subfilter->S_theta,"S_theta, after filter_measurement_update, filter->subfilter->S_theta", app_context);
1315 app_context_report_info(INFO_PETSC_MAT, 1, &filter->K,"K, after filter_measurement_update, filter->K", app_context);
1316 }
1317
1318 app_context_report_error(ierr,"exit of filter_measurement_update",app_context,__LINE__,__FUNCT__,__FILE__);
1319 return filter;
1320}
1321
1322//-------------------------------------------------------------------------------------------------
1323//---optimizer routines (not implemented yet. TODO!!)---
1324
1325OPTIMIZER_TYPE*
1326optimizer_create_and_initialize(APP_CONTEXT_TYPE* app_context)
1327{
1328 OPTIMIZER_TYPE* optimizer = NULL;
1329 app_context->ierr = PetscMalloc(sizeof(OPTIMIZER_TYPE), (void**) &optimizer);
1330 app_context_report_error(app_context->ierr,"filter_create_and_initialize: PetscMalloc",app_context,__LINE__,__FUNCT__,__FILE__);
1331
1332 optimizer_set_intial_condition(optimizer,app_context);
1333
1334 return optimizer;
1335}
1336
1337OPTIMIZER_TYPE*
1338optimizer_set_intial_condition(OPTIMIZER_TYPE* optimizer, APP_CONTEXT_TYPE* app_context)
1339{
1340 return optimizer;
1341}
1342
1343
1344//---- Pesudo-code (matlab) of reduced-order UKF ----------------------------
1345//--- global variable declaration and initialization
1346
1347//configurationfile='';
1348//estimationFromRealData = 0; %defualt value
1349//estimator = 'seik';
1350
1351
1352//---configuration
1353// load(configurationfile);
1354// save configuration state
1355
1356
1357// set up simulation folders and default variables
1358
1359
1360// -------------- filter variables---------
1361//nadof = ndof+npdof;
1362//r=npdof;
1363//u_0 = zeros(ndof,1);%only the size matter, no the values
1364//H = [eye(ndof,ndof) zeros(ndof,npdof)];
1365//R_d = eye(ndof,ndof)*R_scale;
1366//model.u_t = u_0;
1367
1368//parameters_truth = normalize('normalize', parameters_truth,model.parametersRange);
1369
1370//if(exist('pertubulation','var')==1)
1371// model.parameters = parameters_truth + (rand(size(parameters_truth))-0.5)*2 .* pertubulation;
1372//else
1373// model.parameters = normalize('normalize', model.parameters,model.parametersRange);
1374//end
1375
1376//S_plus_theta = diag(min(model.parameters, 1-model.parameters)) * S_plus_theta_scale;
1377//S_plus = [zeros(ndof,npdof); S_plus_theta];
1378
1379
1380//io_save_variable([objective '/-----configuration-----'],[]); %mark this simulation
1381
1382//io_dynamical_simulation_snapshot('initialize',t:Delta_t/nIter:N+Delta_t*(nIter-1)/nIter);
1383
1384
1385// looping
1386//while 1
1387//---time update
1388
1389 // model simulation
1390 // filter equations
1391//---meausrment update
1392 //load measurement
1393 //covariance update and gain matrix caculation
1394 //update state
1395
1396 //post-processing (disp,save)
1397
1398//save parameter history