/DataAssimilation/DataAssimilationConstantSystem/src/data_assimilation_test01.cpp

 1#define MEASUREMENT_DIR  "./" //where to find the measurement files, relative to the path of running the exe file
 2
 3#include "data_assimilation_routines.h"
 4
 5//model operator 01: simplest
 6MODEL_TYPE* 
 7model_opeartor_01(MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context);
 8
 9//model operator 01: constant system
10//evalute model with (model->n_ensem) sigma points, 
11//forwad the system state (model->theta_k[i],model->x_k[i]) and previous states (model->x_k_minus_1[i],model->theta_k_minus_1[i])  
12//by 1 (model->k) time step or 0.01s (model->t)
13MODEL_TYPE* 
14model_opeartor_01(MODEL_TYPE* model, APP_CONTEXT_TYPE* app_context)
15{
16  PetscInt i=0;
17
18  for(i=0; i<model->n_ensem; i++){
19    //1. update the system state vector
20    VecCopy(model->x_k[i],model->x_k_minus_1[i]);
21    
22    //2. update the system parameter vector
23    VecCopy(model->theta_k[i],model->theta_k_minus_1[i]);
24    //VecCopy(model->theta_k[i],model->x_k[i]);
25    
26    //debug diagnosis
27    if(app_context->diagnosisMode & DIAGNOSIS_DISP_INITIAL_CONF){
28      app_context_report_info(INFO_PETSC_INT, 1, &model->k,"===simulation==== \nmodel_opeartor_01, model->k", app_context);
29      app_context_report_info(INFO_PETSC_REAL, 1, &model->t,"model_opeartor_01, model->t", app_context);
30
31      app_context_report_info(INFO_PETSC_INT, 1, &i,"i",  app_context);
32      
33      app_context_report_info(INFO_PETSC_INT, 1, &model->n_x,"model_opeartor_01, model->n_x", app_context);
34      app_context_report_info(INFO_PETSC_INT, 1, &model->n_theta,"model_opeartor_01, model->n_theta", app_context);
35      
36      app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k[i],"model_opeartor_01, model->theta_k_hat", app_context);    
37      app_context_report_info(INFO_PETSC_VEC, 1, &model->theta_k_minus_1[i],"model_opeartor_01, model->theta_k_minus_1_hat", app_context);    
38    }
39    if(app_context->diagnosisMode & DIAGNOSIS_DISP_EVERYTHING){
40      app_context_report_info(INFO_PETSC_VEC,1,  &model->x_k[i],"model_opeartor_01, model->x_k_hat", app_context);        
41      app_context_report_info(INFO_PETSC_VEC, 1, &model->x_k_minus_1[i],"model_opeartor_01, model->x_k_minus_1_hat", app_context);    
42    }
43  }
44  
45  //update the system time stamp
46  model->k=model->k+1;
47  model->t=model->t+0.01;
48  
49  return model;
50}
51
52
53//-------------------------------------------------------------------------------------------------
54//---test routines---
55int
56test01(int argc, char **args)
57{
58  APP_CONTEXT_TYPE* app_context = app_context_create_and_initialize(argc, args);
59  MEAUSREMENT_TYPE* measurement = measurement_create_and_initialize(app_context);
60  MODEL_TYPE*       model       = model_create_and_initialize      (&model_opeartor_01, app_context);
61  FILTER_TYPE*      filter      = filter_create_and_initialize     (app_context, model);
62  
63  app_context_create_and_initialize_dependently(app_context, filter, model, measurement);
64  filter_set_intial_condition(filter, model, app_context);
65  
66  for(int k=0; k<3; k++){
67    filter_time_update       (filter, model,              app_context);
68    filter_measurement_update(filter, model, measurement, app_context);
69  }
70  
71  app_context_destroy_and_finalize(app_context);
72
73  return 0;
74}
75
76//-------------------------------------------------------------------------------------------------
77//---driver routines---
78int 
79main(int argc,char **args)
80{
81  return test01(argc,args);
82}