eop_fd_operators.m

Go to the documentation of this file.

function [L_I, L_E, b, L_I_adj, L_E_adj] = eop_fd_operators(model, model_data)
%[L_I, L_E, b, L_I_adj, L_E_adj] = eop_fd_operators(model, model_data)
%
% function which computes either L_I, L_E and b the components of the
% linear evolution scheme or L_I_adj and L_E_adj the components of the
% corresponding dual problem (controlled via model.want_dual). 
% Supports affine decomposition via model.decomp_mode
%
% decomp_mode: operation mode of the function
%     - 0='complete' (default): no parameter dependence or decomposition is 
%                performed.
%     - 1='components': For each output argument a cell array of output
%                 arguments is returned representing the q-th component
%                 independent of the parameters given in mu_names  
%     - 2='coefficients': For each output argument a cell array of output
%                 arguments is returned representing the q-th coefficient
%                 dependent of the parameters given in mu_names
%
% the construction of L_I follows a 9-point-star described in my diploma
% thesis (chapter 4.2)

% Dominik Garmatter 02.04 2012

L_I = [];
L_E = [];
b = [];
L_I_adj = [];
L_E_adj = [];

if ~isfield(model, 'decomp_mode')
    model.decomp_mode = 0;
end
if ~isfield(model, 'want_dual')
    model.want_dual = 0;
end

decomp_mode = model.decomp_mode;

if decomp_mode == 0 % no affine decomposition
k = model_data.grid.nvertices;  % k = number of grid points

% building the parametermatrix, which is used later for L_I and L_I_adj
XI = [model.sigma1*model.sigma1,...
    2*model.rho/(1+model.rho*model.rho)*model.sigma1*model.sigma2;...
    2*model.rho/(1+model.rho*model.rho)*model.sigma1*model.sigma2,...
    model.sigma2*model.sigma2];

% L_I always has the form L_I = Id + \bar{L_I}. This also contains the
% dirichlet boundary values.
L_I_basic = speye(k); 

% NOTE: this version using spalloc is still sufficiently fast!

% build the entrys of the lower boundary (i.e. S2 = 0)
L_I_S1bound = spalloc(k,k,model.N1*9);
    for m = 2 : model.N1-1
        L_I_S1bound(m,m-1) = - model.deltaT*XI(1,1)*model_data.grid.X(m)*model_data.grid.X(m) / (2*model.h1*model.h1) + model.deltaT*model_data.grid.X(m)*model.r / (2*model.h1); % W of the '3-point-star' at the S2 = 0 boundary
        L_I_S1bound(m,m) = model.deltaT*XI(1,1)*model_data.grid.X(m)*model_data.grid.X(m) / (model.h1*model.h1) + model.deltaT*model.r; % Z of the '3-point-star' at the S2 = 0 boundary
        L_I_S1bound(m,m+1) = - model.deltaT*XI(1,1)*model_data.grid.X(m)*model_data.grid.X(m) / (2*model.h1*model.h1) - model.deltaT*model_data.grid.X(m)*model.r / (2*model.h1); % O of the '3-point-star' at the S2 = 0 boundary
    end
    
% build the entrys of the left boundary (i.e. S1 = 0)
L_I_S2bound = spalloc(k,k,model.N2*9);
    for n = 1 : model.N2-2
        L_I_S2bound(n*model.N1+1,(n-1)*model.N1+1) = - model.deltaT*XI(2,2)*model_data.grid.Y(n*model.N1+1)*model_data.grid.Y(n*model.N1+1) / (2*model.h2*model.h2) + model.deltaT*model_data.grid.Y(n*model.N1+1)*model.r / (2*model.h2); % S of the '3-point-star' at the S1 = 0 boundary
        L_I_S2bound(n*model.N1+1, n*model.N1+1) = model.deltaT*XI(2,2)*model_data.grid.Y(n*model.N1+1)*model_data.grid.Y(n*model.N1+1) / (model.h2*model.h2) + model.deltaT*model.r; % Z of the '3-point-star' at the S1 = 0 boundary
        L_I_S2bound(n*model.N1+1,(n+1)*model.N1+1) = - model.deltaT*XI(2,2)*model_data.grid.Y(n*model.N1+1)*model_data.grid.Y(n*model.N1+1) / (2*model.h2*model.h2) - model.deltaT*model_data.grid.Y(n*model.N1+1)*model.r / (2*model.h2); % N of the '3-point-star' at the S1 = 0 boundary
    end
    
% build the matrix entrys at the inner points
L_I_Inner = spalloc(k,k,model.N1*model.N2*9);
    for j = 1 : model.N2-2
        for i = j*model.N1+2 : (j+1)*model.N1-1
            L_I_Inner(i,i-model.N1-1) = - model.deltaT*XI(1,2)*model_data.grid.X(i)*model_data.grid.Y(i) / (4*model.h1*model.h2); % SW of the 9-point-star
            L_I_Inner(i,i-model.N1) = - model.deltaT*XI(2,2)*model_data.grid.Y(i)*model_data.grid.Y(i) / (2*model.h2*model.h2) + model.deltaT*model_data.grid.Y(i)*model.r / (2*model.h2); % S of the 9-point-star
            L_I_Inner(i,i-model.N1+1) = model.deltaT*XI(1,2)*model_data.grid.X(i)*model_data.grid.Y(i) / (4*model.h1*model.h2); % SO of the 9-point-star
            L_I_Inner(i,i-1) = - model.deltaT*XI(1,1)*model_data.grid.X(i)*model_data.grid.X(i) / (2*model.h1*model.h1) + model.deltaT*model_data.grid.X(i)*model.r / (2*model.h1); % W of the 9-point-star
            L_I_Inner(i,i) = model.deltaT*XI(1,1)*model_data.grid.X(i)*model_data.grid.X(i) / (model.h1*model.h1) + model.deltaT*XI(2,2)*model_data.grid.Y(i)*model_data.grid.Y(i) / (model.h2*model.h2) + model.deltaT*model.r; % Z of the 9-point-star
            L_I_Inner(i,i+1) = - model.deltaT*XI(1,1)*model_data.grid.X(i)*model_data.grid.X(i) / (2*model.h1*model.h1) - model.deltaT*model_data.grid.X(i)*model.r / (2*model.h1); % O of the 9-point-star
            L_I_Inner(i,i+model.N1-1) = model.deltaT*XI(1,2)*model_data.grid.X(i)*model_data.grid.Y(i) / (4*model.h1*model.h2); % NW of the 9-point-star
            L_I_Inner(i,i+model.N1) = - model.deltaT*XI(2,2)*model_data.grid.Y(i)*model_data.grid.Y(i) / (2*model.h2*model.h2) - model.deltaT*model_data.grid.Y(i)*model.r / (2*model.h2); % N of the 9-point-star
            L_I_Inner(i,i+model.N1+1) = - model.deltaT*XI(1,2)*model_data.grid.X(i)*model_data.grid.Y(i) / (4*model.h1*model.h2); % NO of the 9-point-star   
        end
    end

L_I = L_I_basic + L_I_S1bound + L_I_S2bound + L_I_Inner; % build L_I
L_E = speye(k); % L_E = Id, because our discretization is completely implicit
b = zeros(k,1); % no inhomogeneity

% simply transpose the components for the dual case
if model.want_dual
    L_I_adj = L_I';
    L_E_adj = L_E';
end


elseif decomp_mode == 1 % the parameterindependant components are generated
    
grid = model_data.grid;
k = grid.nvertices;

% L_I always has the form L_I = Id + \bar{L_I}. This also contains the
% dirichlet boundary values.
L_I_basic = speye(k);

% build the matrix to the coefficient sigma1*sigma1. It is tridiagonal with 
% diagonal position -1/0/1 and contains parts of the W/Z/O entry of the
% 9-point-star

% approach:
% always build a vector containing the entrys of the corresponding diagonal
% and write those vector entrys on the diagonal. Therefore the vector
% contains the important information. The length of the vector always
% indicates the length of the corresponding diagonal.

% Z entrys
XI_11_center_build = zeros(model.N1,1);
XI_11_center       = zeros(k,1);
i = 1 : model.N1-1; % the index set specifies which vector entrys are non-zero
XI_11_center_build(i,1) = model.deltaT*grid.X(i).*grid.X(i)/(model.h1*model.h1);
XI_11_center(1:model.N1*(model.N2-1),1) = repmat(XI_11_center_build,model.N2-1,1); % build the vector
% W entrys
XI_11_west_build = zeros(model.N1,1);
XI_11_west       = zeros(k-1,1);
i = 1 : model.N1-2; % the index set specifies which vector entrys are non-zero
XI_11_west_build(i,1) = -model.deltaT*grid.X(i+1).*grid.X(i+1)/(2*model.h1*model.h1);
XI_11_west(1:model.N1*(model.N2-1),1) = repmat(XI_11_west_build,model.N2-1,1); % build the vector
% O entrys
XI_11_ost_build = zeros(model.N1,1);
XI_11_ost       = zeros(k-1,1);
i = 2 : model.N1-1; % the index set specifies which vector entrys are non-zero
XI_11_ost_build(i,1) = -model.deltaT*grid.X(i).*grid.X(i)/(2*model.h1*model.h1);
XI_11_ost(1:model.N1*(model.N2-1),1) = repmat(XI_11_ost_build,model.N2-1,1); % build the vector
%Matrix assembly
L_I_XI_11 = sparse(diag(XI_11_center,0)+diag(XI_11_west,-1)+diag(XI_11_ost,1));


% build the matrix to the coefficient sigma1*sigma1. It is tridiagonal with 
% diagonal position -model.N1/0/model.N1 and contains parts of the S/Z/N 
% entry of the 9-point-star

% approach:
% again build a vector containing the entrys of the corresponding diagonal
% and write those vector entrys on the diagonal. Therefore the vector
% contains the important information. The length of the vector always
% indicates the length of the corresponding diagonal. The construction of
% the vector is a little more complex this time.

% Z entrys
XI_22_center_build = zeros(1,model.N2-1);
XI_22_center       = zeros(k,1);
i = 1 : model.N2-1; % the index set specifies which vector entrys are non-zero
% building the vector
XI_22_center_build(1,i) = model.deltaT*grid.Y(model.N1*(i-1)+1).*grid.Y(model.N1*(i-1)+1)/(model.h2*model.h2);
XI_22_center_build_2    = repmat(XI_22_center_build,model.N1,1);
XI_22_center(1:numel(XI_22_center_build_2),1) = XI_22_center_build_2(:);
XI_22_center(model.N1:model.N1:length(XI_22_center)) = 0;
% S entrys
XI_22_sued_build = zeros(1,model.N2-2);
XI_22_sued       = zeros(k-model.N1,1);
i = 1 : model.N2-2; % the index set specifies which vector entrys are non-zero
% building the vector
XI_22_sued_build(1,i) = -model.deltaT*grid.Y(model.N1*i+1).*grid.Y(model.N1*i+1)/(2*model.h2*model.h2);
XI_22_sued_build_2    = repmat(XI_22_sued_build,model.N1,1);
XI_22_sued(1:numel(XI_22_sued_build_2),1) = XI_22_sued_build_2(:);
XI_22_sued(model.N1:model.N1:length(XI_22_sued)) = 0; 
% N entrys
XI_22_nord_build = zeros(1,model.N2-1);
XI_22_nord       = zeros(k-model.N1,1);
i = 2 : model.N2-1; % the index set specifies which vector entrys are non-zero
% building the vector
XI_22_nord_build(1,i) = -model.deltaT*grid.Y(model.N1*(i-1)+1).*grid.Y(model.N1*(i-1)+1)/(2*model.h2*model.h2);
XI_22_nord_build_2    = repmat(XI_22_nord_build,model.N1,1);
XI_22_nord(1:numel(XI_22_nord_build_2),1) = XI_22_nord_build_2(:);
XI_22_nord(model.N1:model.N1:length(XI_22_nord)) = 0;
%Matrix assembly
L_I_XI_22 = sparse(diag(XI_22_center,0)+diag(XI_22_sued,-model.N1)+diag(XI_22_nord,model.N1));


% The matrix for the coefficient sigma1*sigma2*2*model.rho/(1+model.rho^2):
%
%
% XI_12 = spalloc(k,k,4*k);   %Maximal 4 Einträge pro Zeile (NW/NO/SW/SO)
%     for j = 2 : model.N2-1  %Doppel-for-Schleife über alle inneren Punkte
%         for i = (j-1)*model.N1+2 : j*model.N1-1
%              XI_12(i,i+model.N1-1) = model.deltaT*grid.X(i+model.N1-1)*grid.Y(i+model.N1-1)/(4*model.h1*model.h2);   %NW-Eintrag
%              XI_12(i,i+model.N1+1) = -model.deltaT*grid.X(i+model.N1+1)*grid.Y(i+model.N1+1)/(4*model.h1*model.h2);  %NO-Eintrag
%              XI_12(i,i-model.N1-1) = -model.deltaT*grid.X(i-model.N1-1)*grid.Y(i-model.N1-1)/(4*model.h1*model.h2);  %SW-Eintrag
%              XI_12(i,i-model.N1+1) = model.deltaT*grid.X(i-model.N1+1)*grid.Y(i-model.N1+1)/(4*model.h1*model.h2);  %SO-Eintrag 
%         end
%      end
% new version:
i = zeros(model.N2-2,model.N1-2); % number of inner points
% the loop generates an indexmatrix with the indices of all inner points
for j = 1 : model.N2-2 
    i(j,:) = j*model.N1+2 : (j+1)*model.N1-1;
end
indices = reshape(i,1,size(i,1)*size(i,2)); % the inidces are reshaped as a rowvector
Zeilen_build = repmat(indices,4,1);         % this vector is reproduced 4 times, since per point there are 4 (NO,NW,SO,SW) entrys
% with this rowvector we create indexvectors for the columns and rows in
% which the corresponding entrys are writeen
Zeilen = Zeilen_build(:);
Spalten = Zeilen + repmat([model.N1-1;model.N1+1;-model.N1-1;-model.N1+1],size(i,1)*size(i,2),1);
% build the matrix
L_I_XI_12 = sparse(Zeilen, Spalten, grid.X(Zeilen).*grid.Y(Zeilen).*repmat([model.deltaT;-model.deltaT;-model.deltaT;model.deltaT],size(i,1)*size(i,2),1)/(4*model.h1*model.h2) ,k ,k);


%Die Matrix zu model.r -> 5 Diagonalen an den Posiotionen
%-model.N1/-1/0/1/model.N1 mit den Resten der Süd/West/Z/Ost/Nord -
%Einträgen.

% build the matrix to the coefficient model.r. It has 5 diagonals with 
% diagonal position -model.N1/-1/0/1/model.N1 and contains parts of the 
% S/W/Z/O/N entry of the 9-point-star
% the procedure is similar to the ones above

% Z entrys
r_center_build = zeros(model.N1,1);
r_center       = zeros(k,1);
r_center_build(1:model.N1-1,1) = model.deltaT;
r_center(1:model.N1*(model.N2-1)) = repmat(r_center_build,model.N2-1,1);
r_center(1,1) = 0;
% S entrys
r_sued_build = zeros(1,model.N2-2);
r_sued       = zeros(k-model.N1,1);
i = 1 : model.N2-2;
r_sued_build(1,i) = model.deltaT*grid.Y(model.N1*i+1)/(2*model.h2);
r_sued_build_2    = repmat(r_sued_build,model.N1,1);
r_sued(1:numel(r_sued_build_2),1) = r_sued_build_2(:);
r_sued(model.N1:model.N1:length(r_sued)) = 0;
% N entrys
r_nord_build = zeros(1,model.N2-1);
r_nord       = zeros(k-model.N1,1);
i = 2 : model.N2-1;
r_nord_build(1,i) = -model.deltaT*grid.Y(model.N1*(i-1)+1)/(2*model.h2);
r_nord_build_2    = repmat(r_nord_build,model.N1,1);
r_nord(1:numel(r_nord_build_2),1) = r_nord_build_2(:);
r_nord(model.N1:model.N1:length(r_nord)) = 0;
% W entrys
r_west_build = zeros(model.N1,1);
r_west       = zeros(k-1,1);
i = 1 : model.N1-2;
r_west_build(i,1) = model.deltaT*grid.X(i+1)/(2*model.h1);
r_west(1:model.N1*(model.N2-1),1) = repmat(r_west_build,model.N2-1,1);
% O entrys
r_ost_build = zeros(model.N1,1);
r_ost       = zeros(k-1,1);
i = 2 : model.N1-1;
r_ost_build(i,1) = -model.deltaT*grid.X(i)/(2*model.h1);
r_ost(1:model.N1*(model.N2-1),1) = repmat(r_ost_build,model.N2-1,1);
%Matrix assembly
L_I_r     = sparse(diag(r_center,0)+diag(r_sued,-model.N1)+diag(r_nord,model.N1)+diag(r_west,-1)+diag(r_ost,1));

% the dual components are simply the transposed
if model.want_dual
    L_I_r_adj = L_I_r';
    L_I_XI_11_adj = L_I_XI_11';
    L_I_XI_22_adj = L_I_XI_22';
    L_I_XI_12_adj = L_I_XI_12';
    L_I_basic_adj = L_I_basic';
    L_I_adj = {L_I_r_adj; L_I_XI_11_adj; L_I_XI_22_adj; L_I_XI_12_adj; L_I_basic_adj};
    L_E_adj = {speye(k)'};
else
    L_I = {L_I_r; L_I_XI_11; L_I_XI_22; L_I_XI_12; L_I_basic};
    L_E = {speye(k)};
    b   = {zeros(k,1)};
end


elseif decomp_mode == 2 % parameterdependant coefficients are computed
    if model.want_dual
        L_I_adj = [model.r; model.sigma1*model.sigma1; model.sigma2*model.sigma2; 2*model.rho/(1+model.rho*model.rho)*model.sigma1*model.sigma2; 1];
        L_E_adj = 1;
    else
        L_I = [model.r; model.sigma1*model.sigma1; model.sigma2*model.sigma2; 2*model.rho/(1+model.rho*model.rho)*model.sigma1*model.sigma2; 1];
        L_E = 1;
        b   = 1;
    end

end