/DataAssimilation/DataAssimilationConstantSystem/src/data_assimilation_test01.cpp
C++ | 82 lines | 50 code | 18 blank | 14 comment | 4 complexity | c8fd2fae57a95678c4e4047a72315e49 MD5 | raw file
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}