%% 设置同化程序全局续算开关
flag_da_xusuan=0; % 设置是否进行全同化程序续算
                  % 0 - 不开启续算；
                  % 1 - 开启中断续算，可导入任意中断后保存.mat文件后自由调整同化参数，开始进行续算；
                  % 2 - 直接定位至指定同化轮数 p 的保存文件，开始续算；
%% 初始化操作
% 若不开启续算，则清理所有工作区变量
if (flag_da_xusuan==0)
    clear all
    clc
    clf
    close all
    flag_da_xusuan=0;
    p=0; % 初始化同化轮数 p
elseif (flag_da_xusuan==2)
    clear all
    flag_da_xusuan=2;
    p=1; % 定位至指定同化轮数 p，继续同化
end

% ------------  设置主要文件交互与并行相关参数 START ------------ %

% 设置数据交互路径
savename='Results_20210604_Case_S809_AoA_12.2_Nmax_120_P_199_DA_58_Beta';
ddir_recal=('Program_NWind2D_Parallel_Recal');  % 设置重计算路径
ddir_post=('Program_NWind2D_Parallel_Post');    % 设置结果后处理路径
new_output_folder=['.\',savename,'\'];

% 若不开启续算开关，则直接建立新文件夹，保存计算结果
if (flag_da_xusuan==0)
    mkdir(new_output_folder);
    mkdir(['.\Output_Figures\',savename]);
end
% 若开启指定续算开关：直接定位至指定同化轮数 p 的保存文件，开始续算（注意：指定同化轮数必须已经计算完毕）；
if (flag_da_xusuan==2)
    load([new_output_folder,savename,'_S_Va_pre_restart_',num2str(p,'%03d'),'.mat']);
end
% 重启/续算模块设置完毕

targetPath=new_output_folder;
initPath='Init_Ensemble_Member';

% 打开记录程序信息文件
fid_info=fopen([new_output_folder,'DA_CompInfo_Brief.dat'],'w+');
fid_all=fopen([new_output_folder,'DA_CompInfo_Detailed.dat'],'w+');

% 设置并行线程/分组数（可根据工作机线程设置）
PPNS=4;

% -------------  设置主要文件交互与并行相关参数 END ------------- %

% -------------  设置主要文件交互与并行相关参数 END ------------- %

% -------------  设置生成初始集合成员相关参数 START ------------- %

% 定义集合成员最大数量（注意：实际的 N_new ≤ N）
N_max=120;

% ----------  指定初始扰动变量、目标同化变量、同化观测变量模块 START ---------- %

% ① 指定初始扰动变量模块（目前默认为单变量扰动）
% 算例标准状态设置
% 定义初始扰动变量名称、文件路径、修改行数、修改格式（均用数组表示，一一对应）
% 初始扰动变量名称
S_Va_turb_Name=["Cb1","Sigma","Cv1","CKS"];
% S_Va_turb_Name=["Beta"];
S_Va_turb_CharName=char(S_Va_turb_Name);
% 初始扰动变量类型：0 - 常数变量及其组合；1 - 空间分布变量
S_Va_turb_Type=0;
% 初始扰动变量个数
S_Va_turb_Num=length(S_Va_turb_Name);
% 初始扰动变量文件路径（省略运行程序根目录）
% S_Va_turb_FilePath_raw=char('\Input\InputNS.dat','\Input\InputNS.dat','\Input\InputTurb.dat');
% S_Va_turb_FilePath=strip(S_Va_turb_FilePath_raw,'right');
S_Va_turb_FilePath=["\Input\InputSAconst.dat","\Input\InputSAconst.dat","\Input\InputSAconst.dat","\Input\InputSAconst.dat"];
% 初始扰动变量修改行数（对应文件内）
S_Va_turb_ModiLine=[3,5,6,7];
% 初始扰动变量修改格式（对应文件内）
S_Va_turb_ModiFormat=["%-16.5f","%-16.5f","%-16.5f","%-16.5f"];
% 初始扰动变量扰动区间
% 设置扰动方式：0 - 绝对值上下扰动；1 - 比例上下扰动（例如：设置50代表扰动区间为标准值的[50%,150%]）
S_Va_turb_way=[1,1,1,1];
S_Va_turb_Std=[0.1355,(2/3),7.1,0.41];  % 初始扰动变量标准值
S_Va_turb_Det_base_lower=[10,10,10,10];  % 初始扰动变量下区间
S_Va_turb_Det_base_upper=[10,10,10,10];  % 初始扰动变量上区间

S_Va_turb_Det_real_lower=S_Va_turb_Det_base_lower;
S_Va_turb_Det_real_upper=S_Va_turb_Det_base_upper;
for i=1:length(S_Va_turb_Std)
    % 若采用比例上下扰动方式
    if (S_Va_turb_way(i)==1)
        S_Va_turb_Det_real_lower(i)=S_Va_turb_Std(i)*(S_Va_turb_Det_base_lower(i)/100);
        S_Va_turb_Det_real_upper(i)=S_Va_turb_Std(i)*(S_Va_turb_Det_base_upper(i)/100);
    end
end
S_Va_turb_min=S_Va_turb_Std-S_Va_turb_Det_real_lower;
S_Va_turb_max=S_Va_turb_Std+S_Va_turb_Det_real_upper;
% 利用拉丁超立方进行初始扰动变量采样
Va_min=S_Va_turb_min;
Va_max=S_Va_turb_max;
% 初始同化变量个数
PN=length(Va_min);
S_Va_coef=lhsdesign(N_max,PN);
% 生成N*P列扰动变量向量
S_Va_init=(ones(N_max,1)*Va_min)+S_Va_coef.*(ones(N_max,1)*(Va_max-Va_min));
S_Va_init_ave=mean(S_Va_init);

% ② 指定目标同化变量模块
% 全部目标同化变量名称
S_Va_target_Name=["TMu"];
% S_Va_target_Name=["Beta"];
S_Va_target_CharName=char(S_Va_target_Name);
% 全部目标同化变量数量
S_Va_target_Num=length(S_Va_target_Name);
% 判断目标同化变量类型①（TypeTurb）：k - 是第k个初始扰动变量或分布（k~=0）；0 - 非初始扰动变量
% 注：若目标同化变量非初始扰动变量，则需要在更新步将值改为原始 STD 值
S_Va_target_TypeTurb=[0];
% 判断目标同化变量类型②（TypeDA）：1 - 单（值）同化变量；2 - 流场（多值）同化变量
S_Va_target_TypeDA=[2];
% 目标单值同化（且为初始扰动）的变量个数为 PT；流场（多值）同化变量个数 PS，即
% 例如在本算例中，如欲同化马赫数 Ma、迎角 AoA、全（部分）流场涡粘 TMu，则 PT = 2；
%                   欲仅同化迎角 AoA、全（部分）流场涡粘 TMu，则 PT = 1；
% PT=sum((S_Va_target_TypeTurb>0) & (S_Va_target_TypeDA==1));
PT=sum((S_Va_target_TypeDA==1));
PS=sum((S_Va_target_TypeDA==2));
% 标记初始扰动变量中同时满足[为目标同化]和[单值]的变量（e.g. "Ma" "AoA"），并保存在S_Va_target_NonTurb数组
Lia=ismember(S_Va_turb_Name,S_Va_target_Name);
S_Va_turb_target=double(Lia);

% ③ 指定同化观测变量模块
% 全部同化观测变量名称
S_Va_ref_Name=["Cp"];
% S_Va_ref_Name=["BL"];
S_Va_ref_CharName=char(S_Va_ref_Name);
% 全部同化观测变量数量
S_Va_ref_Num=length(S_Va_ref_Name);

% 若目标观测变量为Cp,Cf,yplus    
    % 设置参与同化过程的翼型数据点编号
    % num_da_min=1;    % 取样最小值
    % num_da_max=97;  % 取样最大值
    % num_da_gap=2;    % 取样间隔
    % 建立同化数据点序号数组
    % num_da_arr=(num_da_min:num_da_gap:num_da_max);
    % 可根据重点区域进行标志
    % num_da_arr=[(10:2:42),(48:64),(65:73),(74:2:92)];
    % num_da_arr=[(10:2:42),(48:92)];
    % num_da_arr=[(10:2:42),(48:73),(74:2:92)]; % RAE2822
    % num_da_arr=[(18:3:63),(90:2:168)]; % 56 Points Selected for S809 （P = 199）
    num_da_arr=[(20:2:60),(94:2:166)]; % 58 Points Selected for S809 （P = 199）
    % 记录同化数据点个数
    k=length(num_da_arr);

% 若目标观测变量为边界层速度型，则需要指定边界层站位数量、导入各站位观测信息（y/c,u/u_ref）、各站位处边界层同化位置
    % N_BL=4;
    % BL_Type=[1,1,1,1];  % 设置各站位边界层取点模式：0 - 指定起始、终止编号以及间隔；1 - 指定同化点编号
    % BL_S=[]; % 设置各站位计算值起始编号
    % BL_T=[]; % 设置各站位计算值终止编号
    % BL_G=[]; % 设置各站位计算值编号间隔
    % global BL_setting N_BL L_BL BL_num BL_lab BL_dot
    BL_setting=textread('BL_DA_Setting.dat'); % 导入边界层同化设置文件
    N_BL=size(BL_setting,1);
    L_BL=size(BL_setting,2);
    % BL_dot=(-1)*ones(N_BL,100); % 保存边界层同化位置信息
    BL_num=zeros(1,N_BL); % 保存各站位边界层同化点数
    BL_lab=zeros(1,N_BL); % 保存参与同化的边界层编号
    for i=1:N_BL
        BL_lab(i)=BL_setting(i,1); % 保存边界层站位编号
        BL_type=BL_setting(i,2); % 设置各站位边界层取点模式：0 - 指定起始、终止编号以及间隔；1 - 指定同化点编号
        if (BL_type==0)
            BL_S=BL_setting(i,3);
            BL_T=BL_setting(i,4);
            BL_G=BL_setting(i,5);
            cnt_dot=0;
            for j=BL_S:BL_G:BL_T
                cnt_dot=cnt_dot+1;
                BL_dot(i,cnt_dot)=j;
            end
            BL_num(i)=cnt_dot;
        elseif (BL_type==1)
            cnt_dot=0;
            for kkk=3:L_BL
                if (BL_setting(i,kkk)~=0)
                    cnt_dot=cnt_dot+1;
                    BL_dot(i,cnt_dot)=BL_setting(i,kkk);
                end
            end
            BL_num(i)=cnt_dot;
        end
    end
    
    % Pos_Obs=(-1)*ones(N_BL,100); % 保存边界层同化位置信息
    % 导入指定边界层的速度型信息（实验数据 + 位置信息）
    kb=sum(BL_num); % 统计参与同化的边界层速度型点数量
    Exp_Obs_Arr=zeros(1,kb);
    BL_obs_num=0;
    for i=1:N_BL
        Vel_Exp_Data=load(['.\Input_Data\Input_Exp_Hump_BL_Vel_',num2str(BL_lab(i),'%03d'),'.dat']);
        Vel_Shape_Data=load(['.\Input_Data\Input_Hump_BL_Profile_',num2str(BL_lab(i),'%03d'),'.dat']);
        
        % 整理指定编号同化点的位置信息（y/c）
        for j=1:BL_num(i)
            Pos_Obs(i,j)=Vel_Shape_Data(BL_dot(i,j));
            % 插值得到边界层指定位置的速度实验数据
            Exp_Obs(i,j)=interp1(Vel_Exp_Data(:,1),Vel_Exp_Data(:,2),Pos_Obs(i,j),'linear','extrap');
            BL_obs_num=BL_obs_num+1; % 同化观测点 + 1
            Exp_Obs_Arr(BL_obs_num)=Exp_Obs(i,j); % 保存所有插值出的实验观测数据，数量 = k
        end
        
    end
    

% -----------  指定初始扰动变量、目标同化变量、同化观测变量模块 END ----------- %

% --------------  设置生成初始集合成员相关参数 END -------------- %

% -----------------  设置主要求解迭代参数 START ----------------- %

MStepNS_cov=50000;   % 定常NS求解收敛最大迭代步数（集合生成过程）
Max_inner_iter=10;  % 设置内迭代操作最大次数
MStepNS_pre=50000;   % 设置内迭代预测步数
MStepNS_recal_ensemble=60000;  % 设置每代成员重计算步数（集合重计算过程）
CFL_Setting=[5,5,5]; % 分别设置集合成员生成、内迭代以及重计算过程的 CFL 数
Resi_target_log=[-9.0,-9.0,-9.0]; % 分别定义集合成员生成、内迭代、重计算过程的收敛精度（以 10 为底）
Resi_target=10.^(Resi_target_log);
Resi_target_lev=[-7.5,-6.0,-7.0]; % 规定各阶段（集合生成、内迭代、重计算）的截断残差水平，以确定集合成员捕捉范围
% Resi_target_gen_level=Resi_target_log(1)+1.5; % 规定集合生成截断残差水平为期望残差值 + f ，以扩大集合成员捕捉范围
JX_NS=[0,0,0]; % 更新（关闭）JX_NS 续算操纵开关，即设置流场是否续算：0 - 否（从初始流场开始计算）；1 - 是 （从前一轮计算流场开始计算）

JX_Turb=[3,3,3]; % 更新（关闭）JX_Turb 续算操纵开关，即设置涡粘场是否续算：0 - 否（不固定涡粘，直接采用湍流模型计算）；1 - 是 （固定涡粘进行计算）；3 - 引入 Beta 项进行计算
STD_Comp=[1,1,1]; % 是否进行当前算例与标准算例计算结果的比较（流场各变量 RU、U、V、P、TMu 的绝对误差场，以及总MSE、RMSE）
                  % 需要在Output文件夹中导入 FlowField_STD.PLT 文件，输出 FlowField_Err.PLT 和 FlowField_Err_Info.txt 文件
L2NORM_Setting=[0,0,0]; % 设置每轮 L2NORM 收敛要求水平，达到即退出当前计算（注：5 - 代表收敛截断水平为 5%；当且仅当湍流续算开关 JX_Turb = 1 时生效）


% 设置算例类型（0 - 翼型绕流：需导入上下翼型外形数据以及实验数据，共4个输入文件，并进行插值排序）
%               1 - 单壁内流，如后台阶（Backstep）、驼峰（Hump）等：需导入壁面外形参数以及实验数据，共2个输入文件，并进行插值排序
flag_flow_type=0;

flag_std=0;  % 设置是否进行初始标准算例：0 - 否；1 - 是
flag_da_part=0;  % 设置是否进行部分区域数据同化：0 - 否；1 - 是
flag_da_method=2; % 设置数据同化核心算法：1 - 集合卡尔曼滤波（EnKF）；2 - 集合变换卡尔曼滤波（ETKF）
flag_da_relative=1; % 设置数据同化观测值类型：0 - 相对误差，同化目标变为使相对误差变为零；1 - 绝对误差（传统方式）
flag_beta_turb_type=0; % 设置网格取样点规则（对于 Beta修正项 同化情形）：0 - 每次扰动点固定；1 - 每次扰动点随机；

% 如果初始扰动与同化目标中出现 "Beta" 项，则进行扰动设置
if ((strcmp(char(S_Va_target_Name(1)),'Beta')==1) || (strcmp(char(S_Va_target_Name(1)),'Beta')==1))
    flag_da_part=1; % 开启区域数据同化
    flag_da_beta=1; % 开启扰动项Beta(x)同化
    flag_Va_init=1; % 开启集合成员导入流场多维工作变量标志
    flag_sel_dot=2; % 设置扰动Beta(x)的方式
                    % 拟定规则 1 - 直接选取相关编号间隔，即在选取区域内每隔一定编号进行单点采样，选取散点并进行平面插值；
                    % 拟定规则 2 - 精细化采样均匀采样；
                    % 拟定规则 3 - 采样二维光滑函数进行采样（待开发......）；
    % 设置 Beta 扰动规则：方差 sigma_2，满足分布 Beta ~ N(std_beta,sigma_2)
    % 规则①取点间隔信息
    DA_NG_gap=200;
    % 规则①~③取点距离信息
    std_beta=1.0;
    sigma_beta=3.0;
    dis_ratio=5.0;
else
    % 若不进行 Beta扰动项 同化
    flag_da_beta=0; % 关闭扰动项Beta(x)同化
    flag_Va_init=0; % 关闭集合成员导入流场多维工作变量，使用普通单变量组合扰动
end

% 区域数据同化设置
NNDA=4; % NNDA -- 同化区域节点个数，可为任意多边形（须满足>=3），下面输入 NNDA 个坐标（XDA,YDA）
% XDA=[1.05,1.2,1.2,1.05]; % XDA(NNDA) -- 同化区域节点横坐标
% YDA=[0.1,0.1,0.0,0.0]; % YDA(NNDA) -- 同化区域节点纵坐标
% XDA=[0.8,1.5,1.5,0.8]; % XDA(NNDA) -- 同化区域节点横坐标
% YDA=[0.1,0.1,0.02,0.02]; % YDA(NNDA) -- 同化区域节点纵坐标
XDA=[0.5499846528026755,1.3977168074273609,1.3977168074273609,0.5499846528026755]; % XDA(NNDA) -- 同化区域节点横坐标
YDA=[0.11283900866252453,0.16714063248086963,0.005,0.005]; % YDA(NNDA) -- 同化区域节点纵坐标

if (flag_flow_type==0)
    % 导入上下翼面的外形数据（自动插值）
    shape_upper=load('.\Input_Data\Input_Data_Shape_S809_upper.dat');
    shape_lower=load('.\Input_Data\Input_Data_Shape_S809_lower.dat');
    % 导入上下翼面的实验数据
    exp_upper_data=load('.\Input_Data\s809_200w_aoa12.2_exp_down.dat');
    exp_lower_data=load('.\Input_Data\s809_200w_aoa12.2_exp_up.dat');
elseif (flag_flow_type==1)
    % 导入外形数据（自动插值）
    shape_outline=load('.\Input_Data\Input_Data_Shape_Hump.dat');
    % 导入实验数据
    exp_outline=load('.\Input_Data\Input_Data_Exp_Hump.dat');
end
    
% ------------------  设置主要求解迭代参数 END ------------------ %

% -----------  设置数据同化实验值与误差估计参数 START ----------- %

% 设置观测值协方差（方差），并假定这些测量误差不相关
% !!!!!!! 方差具体设置值需要再考虑（大一些？自适应方差调整模型？） !!!!!!!
% 2021/04/27 更新模块
% 设置是否开启 自适应观测值标准差调节模型 - Adaptive observation standard deviation adjustment model 
%               abbr. AOSDAM（否 - 0；固定放缩模式 - 1；区间平均模式 - 2；区间适应模式 - 3）
% 算法说明：若关闭自适应模型，则设定所有观测值标准差为一固定值 sigma_obs；
%           若开启自适应模型（1 - 固定放缩模式）：则每轮同化前所有观测值标准差设置为上一轮标准差乘以放缩因子；
%           若开启自适应模型（2 - 区间平均模式），则每轮同化前根据上一轮所有观测值对应的计算值包络区间平均宽度，将所有观测值的标准差调整为平均期望标准差水平（乘以放缩因子）；
%           若开启自适应模型（3 - 区间适应模式），则每轮同化前根据上一轮各观测值对应的计算值包络区间宽度，将各观测值的标准差调整为相应期望标准差水平（乘以放缩因子）；
% 变量说明：flag_sigma_obs_adapt - 模式开关
%           flag_sigma_obs_coeff - 放缩因子
%           sigma_obs_init       - 设定初始观测值标准差水平
%           sigma_obs_min        - 设定观测值标准差下限，即调整前后的标准差均不可低于下限
%           sigma_obs_max        - 设定观测值标准差上限，即调整前后的标准差均不可高于上限
flag_sigma_obs_adapt=0;
flag_sigma_obs_coeff=0.80;
sigma_obs_ini=0.15; % 设置初始观测标准差水平
sigma_obs_min=0.02; % 指定观测标准差下限为0.05，对应方差为0.0025
sigma_obs_max=0.80; % 指定观测标准差上限为0.80，对应方差为0.6400

% 计算数据同化初始集合分析值误差协方差
% C_pp_e_a=cov(S_Va_ana');

% ----------  设置数据同化实验值与误差估计参数 END ---------- %

num_da_arr=[(12:7:75),(76:5:186),(208:5:288),(290:3:320),(321:5:376),(377:4:405)]; % 81 Points Selected for RAE2822 (P = 420)