ReducedData.m

classdef ReducedData < Greedy.User.IReducedDataNode
  % Reduced data implementation for linear evolution problems with finite
  % volume discretizations
  %
  % See @ref HO08a for details on the implementation of the reduced matrices
  % and vectors.


  properties
    % Dof vectors of projections `\left\{{ \cal P }_{\text{red}}[u_0^q]
    % \right\}_{q=1}^{Q_{u_0}}`.
    %
    % @sa LinEvol.ReducedModel.rb_init_values() for details.
    a0;

    % `Q_I` reduced matrix for implicit operator `{\cal L}_{h,I}`
    %
    % The matrix entries are:
    % ``({ \bf L}_I^q)_{ij} = \int \varphi_i {\cal L}^q_{h,I}[\varphi_j]``
    % for `i,j=1,\ldots,N` and `q=1,\ldots,Q_I`.
    LL_I;

    % `Q_E`-sequence of reduced matrices for explicit operator `{\cal L}_{h,E}`
    %
    % The matrix entries are:
    % ``({ \bf L}_E^q)_{ij} = \int \varphi_i {\cal L}^q_{h,E}[\varphi_j]``
    % for `i,j=1,\ldots,N` and `q=1,\ldots,Q_E`.
    LL_E;

    % `Q_b`-sequence of reduced vectors for constant contributions `b_h`
    %
    % The matrix entries are:
    % ``({ \bf b})_{i} = \int \varphi_i b^q_h``
    % for `i=1,\ldots,N` and `q=1,\ldots,Q_b`.
    bb;

    % `Q_I \cdot Q_I`-sequence of gram matrices
    %
    % The matrix entries are:
    % ``({ \bf K}_{II}^{q_1,q_2})_{ij} = \int {\cal L}^{q_1}_{h,I}[\varphi_i] {\cal L}^{q_2}_{h,I}[\varphi_j]``
    % for `i,j=1,\ldots,N` and `q_1,q_2=1,\ldots,Q_I`.
    K_II;

    % `Q_I \cdot Q_E`-sequence of gram matrices
    %
    % The matrix entries are:
    % ``({ \bf K}_{IE}^{q_1,q_2})_{ij} = \int {\cal L}^{q_1}_{h,I}[\varphi_i] {\cal L}^{q_2}_{h,E}[\varphi_j]``
    % for `i,j=1,\ldots,N` and `q_1=1,\ldots,Q_I, q_2=1,\ldots,Q_E`.
    K_IE;

    % `Q_E \cdot Q_E`-sequence of gram matrices
    %
    % The matrix entries are:
    % ``({ \bf K}_{EE}^{q_1,q_2})_{ij} = \int {\cal L}^{q_1}_{h,E}[\varphi_i] {\cal L}^{q_2}_{h,E}[\varphi_j]``
    % for `i,j=1,\ldots,N` and `q_1,q_2=1,\ldots,Q_E`.
    K_EE;

    % `Q_I \cdot Q_b`-sequence of vectors
    %
    % The vector entries are:
    % ``({ \bf m}_I^{q_1,q_2})_i = \int {\cal L}^{q_1}_{h,I}[\varphi_i] b_h^{q_2}``
    % for `i=1,\ldots,N` and `q_1=1,\ldots,Q_I, q_2=1,\ldots,Q_b`.
    m_I;

    % `Q_E \cdot Q_b`-sequence of vectors
    %
    % The vector entries are:
    % ``({ \bf m}_E^{q_1,q_2})_i = \int {\cal L}^{q_1}_{h,E}[\varphi_i] b_h^{q_2}``
    % for `i=1,\ldots,N` and `q_1=1,\ldots,Q_E, q_2=1,\ldots,Q_b`.
    m_E;

    % `Q_b \cdot Q_b`-sequence of scalars
    %
    % The scalars are:
    % ``{ \bf m}^{q_1,q_2} = \int b_h^{q_1} b_h^{q_2}``
    %  and `q_1,q_2=1,\ldots,Q_b`.
    m;

    % in case of output functional this is a `Q_s`-sequence of vectors.
    %
    % The vector entries are:
    % ``({ \bf s}^q)_i = s(\varphi_i)``
    % for `i=1,\ldots,N` and `q_1=1,\ldots,Q_s`
    s_RB;

    % in case of output functional this is a `Q_s`-sequence of scalars.
    %
    % The scalars are:
    % ``{ \bf s}_{l^2}^q = \|s(\varphi_1), \ldots, s(\varphi_N)\|_{l^2}``
    % for `q_1=1,\ldots,Q_s`
    s_l2norm;
  end

  properties (SetAccess = private, Dependent)
    % number of reduced basis vectors stored in this data node.
    N;
  end

  properties (SetAccess = private)
    % number of collateral reduced basis vectors stored in this data node.
    M = 0;
    % number of collateral reduced basis vectors used for error estimation
    Mstrich = 0;
  end

  methods

    function rd = ReducedData(rmodel, detailed_data)
      % Constructor for the generation of the reduced matrices and vectors.
      %
      % Parameters:
      %  rmodel: of type LinEvol.ReducedModel
      %  detailed_data: of type IDetailedData
      %
      %
      if nargin == 0
        error('LinEvol.ReducedData constructor needs an argument');
      elseif nargin == 2 && isa(rmodel, 'IReducedModel') ...
          && isa(detailed_data, 'LinEvol.ReducedData')
        copy = detailed_data;
        copy_extract(rd, copy, rmodel);
      elseif nargin == 2 && isa(rmodel, 'IReducedModel') && isa(detailed_data, 'IDetailedData')
        fill_fields(rd, rmodel, detailed_data);
      else
        error('Did not find constructor for your arguments');
      end
    end

    function N = get.N(this)
      if ~isempty(this.a0)
        N = size(this.a0{1},2);
      else
        N = 0;
      end
    end

    function yesno = needs_subset_copy(this, rmodel)
      % function yesno = needs_subset_copy(this, rmodel)
      % @copybrief .Greedy.User.IReducedDataNode.needs_subset_copy()
      %
      % @copydetails .Greedy.User.IReducedDataNode.needs_subset_copy()
      %
      % Parameters:
      %   rmodel: of type LinEvol.ReducedModel

      yesno = this.N ~= rmodel.N;
    end


    function conds = get_conds(this)
      disp('get_conds needs to be implemented for LinEvol.ReducedData');
      conds = [];
    end

  end

  methods (Access=private)
    function fill_fields(this, rmodel, detailed_data)
      % function fill_fields(this, rmodel, detailed_data)

      this.a0 = rmodel.descr.mexptr('rb_init_values', 1);
      [this.LL_I, this.LL_E, this.bb, ...
       this.K_II, this.K_IE, this.K_EE, ...
       this.m_I, this.m_E, this.m ] ...
        = LinEvolDune.ReducedData.rb_operators(rmodel, detailed_data, 1);
    end

    function copy_extract(this, copy, rmodel)
      % function copy_extract(this, copy, rmodel)
      % @copybrief IReducedModelextract_reduced_data_subset()
      %
      % Parameters:
      %  rmodel: reduced model of type LinEvol.ReducedModel indicating the
      %          number of reduced basis vectors to be used in the
      %          'reduced_data_subset'

      % extract correct N-sized submatrices and subvectors from reduced_data
      if rmodel.N > copy.N
        error('N too large for current size of reduced basis!');
      end

      N = rmodel.N;
      this.a0   = subblock_sequence(copy.a0,1:N);
      this.LL_I = subblock_sequence(copy.LL_I,1:N,1:N);
      this.LL_E = subblock_sequence(copy.LL_E,1:N,1:N);
      this.bb   = subblock_sequence(copy.bb,1:N);
      this.K_EE = subblock_sequence(copy.K_EE,1:N,1:N);
      this.K_IE = subblock_sequence(copy.K_IE,1:N,1:N);
      this.K_II = subblock_sequence(copy.K_II,1:N,1:N);
      this.m_I  = subblock_sequence(copy.m_I,1:N);
      this.m_E  = subblock_sequence(copy.m_E,1:N);
      this.m    = copy.m;
      %        if ~isempty(copy.T)
      %          this.T    = copy.T(1:N,1:N);
      %        end

      model = rmodel.detailed_model;
      if isfield(model,'name_output_functional')
        this.s_RB = copy.s_RB(1:N);
        this.s_l2norm = copy.s_l2norm;
      end

    end
  end

  methods(Static)
    %
    %
    % Parameters:
    %   rmodel        : of type LinEvol.ReducedModel
    %   detailed_data : of type Greedy.DataTree.Detailed.RBLeafNode
    function [LL_I, LL_E, bb, K_II, K_IE, K_EE, m_I, m_E, m] = rb_operators(rmodel, detailed_data, decomp_mode)

      descr = rmodel.descr;

      if decomp_mode == 2
        mu = get_mu(rmodel);
        dLL_I = descr.coeff_ops.LL_I_ptr(mu);
        dLL_E = descr.coeff_ops.LL_E_ptr(mu);
        dbb   = descr.coeff_ops.bb_ptr(mu);
        %    model.mexptr('set_mu', model.mu);
        %    foo   = model.mexptr('rb_operators', decomp_mode);
        %    disp('start');
        %    dLL_I == foo{1}
        %    dLL_E == foo{2}
        %    dbb   == foo{3}
        %    kron(dLL_I,dLL_E) == foo{6}
        %    kron(dLL_E,dbb) == foo{8}
        %    kron(dbb,dbb) == foo{9}
        LL_I = [1, dLL_I];
        LL_E = [1, dLL_E];
        bb   = dbb;
        K_II = [1, dLL_I, dLL_I, kron(dLL_I,dLL_I) ];
        K_EE = [1, dLL_E, dLL_E, kron(dLL_E,dLL_E) ];
        K_IE = [1, dLL_I, dLL_E, kron(dLL_I,dLL_E) ];
        m_I  = [dbb, kron(dLL_I,dbb) ];
        m_E  = [dbb, kron(dLL_E,dbb) ];
        m    = kron(dbb,dbb);
      else
        foo   = descr.mexptr('rb_operators', decomp_mode);
        dLL_I = foo{1};
        dLL_E = foo{2};
        dbb   = foo{3};
        dK_II = foo{4};
        dK_EE = foo{5};
        dK_IE = foo{6};
        dm_I  = foo{7};
        dm_E  = foo{8};
        dm    = foo{9};

        if decomp_mode == 1
          func2  = @(X) descr.dt * X;
          func2t = @(X) descr.dt * X';
          func3  = @(X) descr.dt.^2 * X;

          identity = speye(size(dLL_I{1}));

          LL_I = [{identity}; ...
            cellfun(func2 , dLL_I, 'UniformOutput', false)];

          LL_E = [{identity}; ...
            cellfun(func2 , dLL_E, 'UniformOutput', false)];

          bb   =  cellfun(func2 , dbb  , 'UniformOutput', false);

          K_II = [{identity}; ...
            cellfun(func2t, dLL_I, 'UniformOutput', false); ...
            cellfun(func2 , dLL_I, 'UniformOutput', false); ...
            cellfun(func3 , dK_II, 'UniformOutput', false)];

          K_EE = [{identity}; ...
            cellfun(func2t, dLL_E, 'UniformOutput', false); ...
            cellfun(func2 , dLL_E, 'UniformOutput', false); ...
            cellfun(func3 , dK_EE, 'UniformOutput', false)];

          K_IE = [{identity}; ...
            cellfun(func2t, dLL_I, 'UniformOutput', false); ...
            cellfun(func2 , dLL_E, 'UniformOutput', false); ...
            cellfun(func3 , dK_IE, 'UniformOutput', false)];

          m_I  = [bb; ...
            cellfun(func3 , dm_I , 'UniformOutput', false)];

          m_E  = [bb; ...
            cellfun(func3 , dm_E , 'UniformOutput', false)];

          m    =  cellfun(func3 , dm   , 'UniformOutput', false);

        else % decomp_mode == 0
          assert(all(get_mu(rmodel) == get_mu(rmodel.detailed_model)));
          LL_I = 1 + model.dt * dLL_I;
          LL_E = 1 + model.dt * dLL_E;
          bb   =     model.dt * dbb;
          K_II = 1 + model.dt * dLL_I' + model.dt * dLL_I + model.dt^2 * dK_II;
          K_EE = 1 + model.dt * dLL_E' + model.dt * dLL_E + model.dt^2 * dK_EE;
          K_IE = 1 + model.dt * dLL_I' + model.dt * dLL_E + model.dt^2 * dK_IE;
          m_I  =     model.dt * dbb                       + model.dt^2 * dm_I;
          m_E  =     model.dt * dbb                       + model.dt^2 * dm_E;
          m    =                                            model.dt^2 * dm;
        end
      end
    end
  end
end