gfr_parallel.m

%% GRACE gravity field recovery (GFR)
%
% GFR_PARALLEL is the main m-file preforming the gravity field estimation. It
% is optimised for parallel computing using local and global parameters for
% different arcs. All necessary parameters and data are specified in the
% following. You have to set a number of iterations you want to perform. With
% the m-file continue_gfr_par_from_iteration.m you can continue a GFR from an
% iteration. Therefore you need to save all necessary variables, e.g. use the
% file save_vars2continue_itr.m to save them in one .mat file. In the end some
% example plots are given you may change or extend.
%
% Written by Florian Wöske, ZARM Uni Bremen, and
% Neda Darbeheshti, AEI Hannover, 2018-07.

%%
format longg;

%% variable inputs

% add folder with all functions
Release1_Path = pwd;
addpath(genpath(Release1_Path));

% EBA_Path = fullfile(Release1_Path,'/../function_eba/');
% addpath(genpath(EBA_Path));
GFR_Path = fullfile(Release1_Path,'/../function_gfr/');
addpath(genpath(GFR_Path));

% folder to save intermediate results
mkdir('results/OutputArcParall');
mkdir('results/OutputNormalEqu');

% folder where observation data are stored
Input_Observation_Path = fullfile(Release1_Path,'/../input_data','MockData_do10_nod4');
addpath(genpath(Input_Observation_Path));

% Define the number of iterations for batch processor,
ItrNo_max = 5;

% order of spherical harmonic coefficients to be estimated
lmaxcs = 10;

% parameters of background gravity field
lmaxf = 10;

% number of days
nods = 4;

% number of observations in each arc (one day, 17280 is the standard)
mKBR = 17280;

%% fixed inputs

% n: size of gravity fields in vec format, estimated and background
n_est = ((lmaxcs+1)^2 + lmaxcs+1)/2;
n_back = ((lmaxf+1)^2 + lmaxf+1)/2;

% GM and Earth radius of gravity field
GM = 0.3986004415E+15;
ae = 0.6378136300E+07;

% reference for comparison
ggm05s = readGFC('../input_data/GGM05S.gfc');
ggm05s = ggm05s(1:n_back,:); % cut to desired length
ggm05s_cs = vec2cs([ggm05s(:,1) ggm05s(:,2)]',ggm05s(:,3),ggm05s(:,4));
% d_vec_i = ggm05s(:,3:4)-gravityf(:,3:4);
% d_vec_i = [ggm05s(:,1:2) d_vec_i];

% Load a-priori/ reference gravity field
egm96 = readGFC('../input_data/EGM96.gfc');
egm96 = egm96(1:n_back,:); % cut to desired length
egm96_cs = vec2cs([egm96(:,1) egm96(:,2)]',egm96(:,3),egm96(:,4));

% reformatting spherical harmonics coefficients
field_cs = egm96_cs;
field_sc = cs2sc(egm96_cs(1:lmaxf+1,1:lmaxf+1),0);
field = field_sc(:)';
[outC1,outS1,nm1] = cs2vec(egm96_cs,false);
field_vec0 = [nm1',outC1',outS1'];

err_vec = egm96;
err_vec(:,3:4) = ggm05s(:,3:4) - egm96(:,3:4);

% initialzes Legendre polynomial calculation
data_plm = initplm(lmaxf,2);

%% load observations data and initial states for all days from files

t0 = [2005 05 01 0 0 2];
date0 = time2str(t0);

KBRL1B = zeros(mKBR,4,nods);
state0 = zeros(12,nods);

% State deviation:
x0x = zeros(12,nods);
t0i = t0;

for Mday = 1:nods

    if Mday >1
        t0i(3)=t0i(3)+1;
    end

    date = time2str(t0i);

    % load KBR data
    data_KBR = ['KBR1B_', date, '_X_02.asc'];
    data_file = fullfile(Input_Observation_Path, date, data_KBR);
    KBRL1Bi = readKBR(data_file);
    KBRL1B(:,:,Mday) = KBRL1Bi(1:mKBR,1:4);

    % load GNV data GRACE A
    data_GNV = ['GNV1B_', date, '_A_02.asc'];
    data_file = fullfile(Input_Observation_Path, date, data_GNV);
    GNVi = readGNV(data_file);
    state0(1:6,Mday) = GNVi(1,2:7);

    % load GNV data GRACE B
    data_GNV = ['GNV1B_', date, '_B_02.asc'];
    data_file = fullfile(Input_Observation_Path, date, data_GNV);
    GNVi = readGNV(data_file);
    state0(7:12,Mday) = GNVi(1,2:7);
end

% save initial state vector for comparisons
state00 = state0;

%% iteration

% number of coefficients
% Minimum degree and order of coefficients to be estimated
lmincs=2;
noc=lmaxcs^2-lmincs^2+2*lmaxcs+1; %total number of coefficients to be estimated

xhat_save = zeros (12, nods, ItrNo_max+1);
chat_save = zeros (noc, ItrNo_max+1);
eps_save = zeros(mKBR,nods,ItrNo_max+1);
chat_dv_save = zeros (ItrNo_max+1,lmaxcs-1);
dv_save = zeros (ItrNo_max+1,lmaxcs-1);

for ItrNo = 0:ItrNo_max

    % batch over all days
    parfor Mday = 1:nods
        tic
        batch_processor_partitioned(Mday,Release1_Path,lmaxcs,mKBR,field,data_plm,GM,ae,lmaxf,state0(:,Mday),KBRL1B(:,1,Mday),KBRL1B(:,3,Mday),x0x(:,Mday))
        toc
    end

    %% solve normal equations

    % read data from arc parallel
    invMxx = zeros(12,12,nods);
    iMxc = zeros(12,noc,nods);
    iMcc = zeros(noc,noc,nods);
    iNx = zeros(12,1,nods);
    iNc = zeros(noc,1,nods);
    iL = zeros(noc,noc,nods);
    iN = zeros(noc,1,nods);
    iHx = zeros(mKBR,12,nods);
    iHc = zeros(mKBR,noc,nods);
    irhodot_deviation = zeros(mKBR,nods);
    for Mday = 1:nods
        FileName = ['arcmatrix',num2str(Mday,'%.2d'),'.mat'];
        File = fullfile(Release1_Path,'results/OutputArcParall', FileName);
        arcstru=load(File);
        invMxx(:,:,Mday) = arcstru.invMxx;
        iMxc(:,:,Mday) = arcstru.Mxc;
        iMcc(:,:,Mday) = arcstru.Mcc;
        iNx(:,:,Mday) = arcstru.Nx;
        iNc(:,:,Mday) = arcstru.Nc;
        McxTinvMxx = iMxc(:,:,Mday)'*invMxx(:,:,Mday);
        iL(:,:,Mday) = McxTinvMxx*iMxc(:,:,Mday);
        iN(:,:,Mday) = McxTinvMxx*iNx(:,:,Mday);
        iHx(:,:,Mday) = arcstru.Hx_save;
        iHc(:,:,Mday) = arcstru.Hc_save;
        irhodot_deviation(:,Mday) = arcstru.rhodot_deviation;
    end

    % sum over the number of days
    sumiMcc = sum(iMcc,3);
    sumiNc = sum(iNc,3);
    sumiL = sum(iL,3);
    sumiN = sum(iN,3);
    L = sumiMcc-sumiL;
    N = sumiNc-sumiN;

    % estimate global parameters
    chat = L \ N;

    % estimate local parameters
    xhat = zeros(12,nods);
    yhat = zeros(mKBR,nods);
    for Mday = 1:nods
        xhat(:,Mday) = invMxx(:,:,Mday)*iNx(:,:,Mday)-invMxx(:,:,Mday)*iMxc(:,:,Mday)*chat;
        % estimate range rate residuals
        yhat(:,Mday)=iHx(:,:,Mday)*xhat(:,Mday)+iHc(:,:,Mday)*chat;
        eps = irhodot_deviation(:,Mday)-yhat(:,Mday);
        eps_save(:,Mday,ItrNo+1) = eps;
%         figure;plot(eps)
%         xlabel('Observation Number');
%         ylabel('residuals [m/s]');
%         title(['Range-rate residuals, Itr',num2str(ItrNo,'%.2d'),',Day',num2str(Mday,'%.2d')])
    end

    % save xhat and chat of iterations
    xhat_save(:,:,ItrNo+1) = xhat;
    chat_save(:,ItrNo+1) = chat;

    % put coefficients in right order for plotting and output
    ko=1;
    for i=1:lmaxcs+1
      for j=1:lmaxcs+1
        if i<=j
          ordering1(ko,1)=j-1;
          ordering1(ko,2)=i-1;
          ko=ko+1;
        end
      end
    end
    nocC=(noc+lmaxcs-1)/2;
    ordering1(3:end,3)=[chat(1:lmaxcs-1)', 0 ,chat(lmaxcs:nocC)']; % put C coeffs.
    ordering1(3+lmaxcs:end,4)=chat(nocC+1:end); % put S coeffs.
    ordering1_cs = vec2cs([ordering1(:,1) ordering1(:,2)]',ordering1(:,3),ordering1(:,4));

    chat_cs = ordering1_cs;
    % chat_vec = ordering1; % orderwise ordering
    [outC1,outS1,nm1] = cs2vec(chat_cs,false); %degreewise
    chat_vec = [nm1',outC1',outS1'];

    chat_dv_save(ItrNo+1,:) = dv_geoidn_no_plot(chat_vec,lmaxcs);
    % dv_geoidn(chat_vec,lmaxcs);

    % % save final monthly gravity solution in a file
    % FileName= ['gfr_NODs',num2str(nods,'%.2d'),'_DO',num2str(lmaxcs,'%.2d'),'_ItrNo',num2str(ItrNo,'%.2d'),'.gfc'];
    % File = fullfile(Release1_Path,'results/OutputNormalEqu', FileName);
    % fileID = fopen(File,'w');
    % for i=1:length(ordering2)
    %     fprintf(fileID,'%4i %4i %17.16e %17.16e\n',ordering2(i,:));
    % end
    % fclose(fileID);

    % set up values for next iteration
    state0 = state0 + xhat;
    x0x = x0x - xhat;

    % Add CS coeffs in cs format to reference field (just up to the order
    % that was estimated)
    field_cs(1:lmaxcs+1,1:lmaxcs+1) = field_cs(1:lmaxcs+1,1:lmaxcs+1) + chat_cs;
    field_sc = cs2sc(field_cs(1:lmaxf+1,1:lmaxf+1),0);
    field = field_sc(:)';

    [outC1,outS1,nm1] = cs2vec(field_cs,false);
    field_vec = [nm1',outC1',outS1'];

    % save dv of true gravity field minus estimated from each iteration
    d_vec = ggm05s(1:n_back,3:4)-field_vec(:,3:4);
    d_vec = [field_vec(:,1:2) d_vec];
    dv_save(ItrNo+1,:) = dv_geoidn_no_plot(d_vec,lmaxcs);

    fprintf('done with iteration %d \n', ItrNo)
end

% plot chat from all iterations:
figure;
for i = 1:ItrNo_max+1
    semilogy(2:lmaxcs,chat_dv_save(i,:),'.-','LineWidth',2,'MarkerSize',20);
    hold on
end

fs = 12;
set(gcf,'PaperPositionMode','auto')
set(gca,'FontSize',fs);
xlabel('Spherical harmonic degree n')
ylabel('Square root of degree variances')
title('chat from each iteration')
grid on

figure;
dv_geoidn(ggm05s,lmaxcs);
hold on
% dv_geoidn(d_vec_i,lmaxcs);

[outC1,outS1,nm1] = cs2vec(field_cs,false);
final_field_vec = [nm1',outC1',outS1'];
dv_geoidn(final_field_vec,lmaxcs);

d_vec = ggm05s(1:n_back,3:4)-final_field_vec(:,3:4);
d_vec = [final_field_vec(:,1:2) d_vec];
dv_geoidn(d_vec,lmaxcs);
legend('ggm05s', 'estimated field', 'ggm05 - estimated field')

% plot gravity field error from each iteration additional to initial fields and
% initial error

figure;
dv_geoidn(ggm05s,lmaxcs);
hold on
dv_geoidn(field_vec0,lmaxcs);
dv_geoidn(err_vec,lmaxcs);
legend('ggm05s','ggm05s + error field', 'error field')
for i = 1:ItrNo_max+1
    str = ['error iter ',num2str(i-1)];
    semilogy(2:lmaxcs,dv_save(i,:),'.-','LineWidth',2,'MarkerSize',16,'DisplayName', str);
end