rb_simulation_impl.m

function rb_sim_data = rb_simulation_impl(rmodel, reduced_data)
%function rb_sim_data = rb_simulation_impl(rmodel, reduced_data)
% function, which performs a reduced basis online simulation for
% the parameter vector `\mu \in {\cal P} \subset \mathbb{R}^p`, which is
% assumed to be set by IDetailedModel.set_mu()
%
% - allowed dependency of data: Nmax, Mmax, N, M, mu
% - not allowed dependency of data: H
% - allowed dependency of computation: Nmax, Mmax, N, M, mu
% - not allowed dependency of computation: H
% - Unknown at this stage: ---
%
% The behaviour of this simulation is controlled by the ModelDescr structure
% 'rmodel.descr'.
%
% Parameters:
%   reduced_data: of type NonlinEvol.EiRbReducedDataNode
%
% Required fields of descr:
%  mu_names       : the cell array of names of mu-components and
%                   for each of these stringt, a corresponding
%                   field in model is expected.
%  T              : end-time of simulation
%  nt             : number of timesteps to compute
%  L_E_local_name : name of the local space-discretization operator 
%                   evaluation with syntax
%                   @code
%                   INC_local = L_local(U_local_ext, ind_local,
%                                 grid_local_ext, model) @endcode
%                   where '*_local' indicates results on a set of
%                   elements, '*_local_ext' includes information including
%                   their neighbours, i.e. the grid must be known also for the
%                   neighbours, the values of the previous timstep must be
%                   known on the neighbours. The subset of elements, which are
%                   to be computed are given in ind_local and the values
%                   produced are 'INC_local'.
%                   Subsequent time-step operator would be defined
%                   as 'NU_local = U(ind_local) - dt * INC_local'
%  rb_operators   : function pointer to a method assembling the decomposed
%                   discretization operators
%
% optional fields of model:
%  data_const_in_time : if this flag is set, then only operators
%                       for first time instance are computed
%
% generated fields of rb_sim_data:
%  a     : 'a(:,k)' is the coefficient vector of the reduced simulation at time
%          `t^{k-1}`
%

% Bernard Haasdonk 15.5.2007

if ~(nargin == 2)
  error('wrong number of arguments: should be 2');
end

if rmodel.compute_conditions
  LIconds = zeros(1, rmodel.descr.nt);
  gradLIconds = [];
end

%% compute a posteriori error estimator (c.f. diploma thesis Martin Drohmann
%   p.  26ff.)
if ~rmodel.enable_error_estimator
  rmodel.Mstrich = 0;
end

if rmodel.verbose >= 10
  disp(['performing simulation for mu = ',mat2str(get_mu(rmodel))]);
end;

if ~isfield(rmodel, 'filecache_velocity_matrixfile_extract')
  filecache_velocity_matrixfile_extract = 0;
else
  filecache_velocity_matrixfile_extract = rmodel.filecache_velocity_matrixfile_extract;
end

% set init data
rmodel.decomp_mode = 2;

sa0 = rb_init_values(rmodel,[], 2);

a0  = lincomb_sequence(reduced_data.a0, sa0);

model = rmodel.descr;


% perform simulation loop in time
model.dt = model.T/model.nt;
a        = zeros(rmodel.N, model.nt+1);   % matrix of RB-coefficients
a(:,1)   = a0;

% in case of velocity-file access, the correct filename must be set here
% the following is only relevant in case of use of a 
% model.use_velocitymatrix_file and filecaching mode 2
if filecache_velocity_matrixfile_extract == 2;
  % change global velocity file to MMax velocity file generated in offline
  rmodel.descr.velocity_matrixfile = ...
      ['M',num2str(model.M(1)),'_','gridpart_Mmax',num2str(model.Mmax),'_',...
       model.velocity_matrixfile];
end

%disp('temporary true mass matrix is considered');
%M = reduced_data.M;

if rmodel.enable_error_estimator
  % Lipschitz constants
  if ~isfield(model, 'C_E')
    C_E = 1;
  else
    C_E = model.C_E;
  end
  if ~isfield(model, 'C_I')
    C_I = 1;
  else
    C_I = model.C_I;
  end
  % initialize error estimator
  Delta = 0;
  % initialize error estimator log (strictly increasing)
  Deltalog = zeros(1, model.nt-1);
  % initialize log of residuum for each time step
  reslog   = zeros(1, model.nt-1);
  % initialize part for ei estimation for each time step
  eilog    = zeros(1, model.nt-1);
  % initialize part for newton error for each time step
  newtonlog    = zeros(1, model.nt);
end

if isa(rmodel.M, 'DataTree.INode')
  datanode_mode = true;
else
  datanode_mode = false;
end

Mstrich = rmodel.Mstrich;
M       = get_current_M(rmodel);
% M plus Mstrich
MM      = M + Mstrich;

impl_tstop = 0;
expl_tstop = 0;

% model for explicit part (has different affine decomposition)
emodel             = model;
emodel.decomp_mode = 0;

% loop over time steps: computation of a(t+1)==a^t
for t = 1:model.nt
  model.t      = (t-1)*model.dt;
  emodel.t     = model.t;
  emodel.tstep = t;
  if rmodel.verbose >= 10
    disp(['entered time-loop step ',num2str(t)]);
  end;

  % update reduced data for operators if necessary.
  if t > expl_tstop
    rmodel.t     = emodel.t;
    rmodel.tstep = emodel.tstep;
    [expl_reduced_data, expl_tstop, Mexpl, Mred_expl]...
      = get_reduced_data_for_operator('explicit', rmodel, reduced_data, datanode_mode);
    rmodel.t     = model.t;
    rmodel.tstep = model.tstep;
  end
  if t >= impl_tstop
    imodel = model;
    imodel.tstep = t+1;
    rmodel.t     = imodel.t;
    rmodel.tstep = imodel.tstep;
    [impl_reduced_data, impl_tstop, Mimpl, Mred_impl]...
      = get_reduced_data_for_operator('implicit', rmodel, reduced_data, datanode_mode);
    rmodel.t     = model.t;
    rmodel.tstep = model.tstep;
  end

  % asssemble affine parameter compositions of implicit parts
  % only once, if data is const in time
  if ~model.newton_solver && ((t==1) || (~model.data_const_in_time))
    if model.implicit_nonlinear
      tmpmodel = model;
      tmpmodel.decomp_mode = 0;
      clear_all_caches;
      [L_I_red, b_I_red] = ...
        model.implicit_operators_algorithm(tmpmodel, impl_reduced_data,...
                                           impl_reduced_data.TM_local); % tii
      clear_all_caches;

      LL_I_red = L_I_red * impl_reduced_data.RB_local_ext; % tii
      LL_I     = impl_reduced_data.CE(:,1:Mred_impl) * ... % tii
                    LL_I_red;

      bb_I     = impl_reduced_data.CE(:,1:Mred_impl) * ...
                   b_I_red;  %tii

    else
      [sLL_I, sbb_I] = ...
        NonlinEvol.ReducedData.rb_operators(rmodel,[], 2); % detailed_data omitted
      if ~isempty(reduced_data.LL_I)
        LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
      else
        LL_I = sparse([],[],[],rmodel.N,rmodel.N);
      end;
      if ~isempty(reduced_data.bb_I)
        bb_I = lincomb_sequence(reduced_data.bb_I, sbb_I);
      else
        bb_I = sparse([],[],[],rmodel.N,1);
      end;
    end;
  end;

  % setup linear system     A * a(:,t+1) = rhs
  % with A = Id + dt * LL_I
  % rhs = a(:,t) - dt * BB_I - dt * CE * ll_E(a(:,t))
  %                  where ll_E is the local operator evaluation

  if ~model.newton_solver
    A = LL_I * model.dt + speye(size(LL_I));
    if rmodel.compute_conditions
      LIconds(t) = condest(A);
    end
  else
    A    = 0;
    bb_I = 0;
  end

  % determine local values by reconstructing U values from reduced
  % basis parts
  U_local_ext = expl_reduced_data.RB_local_ext * a(:,t);  % ti

  % the following performs
  %  U_local = L_local(U_local_ext, TM_local, [], ...
  %                    reduced_data.grid_local_ext);
  %
  % note the use of emodel

  le_local_inc = emodel.L_E_local_ptr(emodel, expl_reduced_data, ...
                                      U_local_ext, ...
                                      expl_reduced_data.TM_local);                    % ti
%  le_local_inc(U_local_ext(reduced_data.TM_local{crbInd_exp,tii}) + le_local_inc > 1) = 0;

  if model.newton_solver
    nmax = -1;
    updelta_a = zeros(rmodel.N, model.newton_steps+1);
    updelta_a(:,1) = a(:,t);
    n=0;

    for n=1:model.newton_steps
      if rmodel.verbose > 10
        disp(['Computing time step no. ', num2str(t), ...
              ', Newton step no. ', num2str(n), ...
              '... Residual norm: ']);
      end
      model.decomp_mode = 0;
      updelta_rec = impl_reduced_data.RB_local_ext * updelta_a(:,n);                  % tii
      if rmodel.debug && any(updelta_rec < -10*eps)
        disp(num2str(min(updelta_rec)));
        warning('RBMatlab:rb_simulation:runtime', ...
                'reconstructed solutions at interpolation points is negative!');
      end
      model.t = model.t + model.dt;
      % updelta_rec(updelta_rec < 0) = 0.001;
      li_local_inc = model.L_I_local_ptr(model, impl_reduced_data, ...
                                         updelta_rec, ...
                                         impl_reduced_data.TM_local);                % tii

      [gradL_I, gL_I_offset] = ...
        model.implicit_gradient_operators_algorithm(model, impl_reduced_data, ...
                                          updelta_rec,...
                           impl_reduced_data.TM_local(1:Mred_impl));                 % tii
      model.t = model.t - model.dt;
      if any(gL_I_offset ~= 0)
        warning('offset has non-zero entries!');
      end

      gradL_I  = impl_reduced_data.CE_red * gradL_I ...
                   * impl_reduced_data.RB_local_ext;

      gradL_I  = gradL_I * model.dt + speye(size(gradL_I));

      opts.LT = true;

      if impl_reduced_data ~= expl_reduced_data
        theta_update_expl = linsolve(expl_reduced_data.BM(1:MM,1:MM), ...
                                    le_local_inc(1:MM), ...
                                    opts);
        theta_update_impl = linsolve(impl_reduced_data.BM(1:MM,1:MM), ...
                                    li_local_inc(1:MM), ...
                                    opts);

        arhs = expl_reduced_data.DE(:,1:Mred_expl) ...
                          * theta_update_expl(1:Mred_expl) ...
               + impl_reduced_data.DE(:,1:Mred_impl) ...
                          * theta_update_impl(1:Mred_impl);
      else
        theta_update = linsolve(expl_reduced_data.BM(1:MM,1:MM), ...
                                le_local_inc(1:MM) + li_local_inc(1:MM), ...
                                opts);

        arhs = expl_reduced_data.DE(:,1:Mred_expl) ...
                          * theta_update(1:Mred_expl);
      end

      full_rhs = - updelta_a(:,n) + a(:,t) - model.dt * arhs;

      delta_a = gradL_I \ full_rhs;

      if rmodel.compute_conditions
        gradLIconds = [gradLIconds, condest(gradL_I)];
      end

      updelta_a(:,n+1) = updelta_a(:,n) + delta_a;
      nres_a = updelta_a(:,n+1) -a(:,t) + model.dt * arhs;
      nres = sqrt(nres_a' * nres_a);
      if model.verbose > 10
        disp(num2str(nres));
      end
      if nres < model.newton_epsilon
        nmax = max(nmax, n);
        break;
      end
    end
    a(:,t+1) = updelta_a(:,n+1);
  else % no newton solver...
    % compute coefficients of update in CRB Space `\W_{M+M'}`.
    opts.LT = true;
    theta_update_expl = linsolve(expl_reduced_data.BM(1:MM,1:MM), ...
                                 le_local_inc, opts);
    theta_update = theta_update_expl;

    rhs = a(:,t) - model.dt * bb_I ...
                 - model.dt * expl_reduced_data.DE(:,1:Mred_expl) ...
                            * theta_update_expl(1:Mred_expl);

    % check if an LES is to be solved, or problem is an explicit problem
    if isequal(A, speye(size(A)))
      a(:,t+1) = rhs;
    else
%        disp('check symmetry and choose solver accordingly!');
%       keyboard;
       % nonsymmetric solvers:
       [a(:,t+1), flag, relres, iter] = bicgstab(A,rhs,1e-10,1000);
       % [a(:,t+1), flag] = cgs(Li,rhs,[],1000); 
       % symmetric solver, non pd:
       % see bug_symmlq for a very strange bug: cannot solve identity system!
       % [a(:,t+1), flag] = symmlq(Li,rhs,[],1000);
       %  [a(:,t+1), flag] = minres(Li,rhs,[],1000); 
       % symmetric solver, pd:
       %[a(:,t+1), flag] = bicgstab(A,rhs,[],1000); 
       if(t==1 && model.debug)
         disp(cond(A));
         disp(all(A'==A));
         disp(eig(full(A)));
       end
       if(t==1) % && model.verbose > 10)
         disp(['relres=', num2str(relres), ' iter=', num2str(iter)]);
       end
       % exact solver:
       % a(:,t+1)=A \ rhs;
       if flag>0
         disp(['error in system solution, solver return flag = ', ...
               num2str(flag)]);
         keyboard;
       end
    end
  end

  if rmodel.enable_error_estimator
    if ~isinf(impl_tstop)
      error('error estimator does not work with varying crbs for different time slices yet');
    end
    ainc           = a(:,t) - a(:,t+1);

    if model.newton_solver
      lipschitzConst = C_I^(model.nt - t + 2)*C_E^(model.nt - t + 1);
      newtonlog(t)   = nres + model.newton_epsilon;
    else
      lipschitzConst = C_E^(model.nt - t + 1);
      newtonlog(t)   = 0;
    end
    if impl_reduced_data == expl_reduced_data
      % note that in case of one basis for both operators all values are the
      % same for both operators, i.e. Mexpl == Mimpl, impl_tstop == expl_tstop,
      % expl_reduced_data.Mmass === impl_reduced_data.Mmass
      %
      % compute the residuum `\delta t norm(R^k)`
      % TODO revise the following
      residuum  = sqrt( abs( ainc' * ainc ...
                      - 2*model.dt * ainc' ...
                                   * ( expl_reduced_data.DE(:,1:Mexpl) ...
                                       * theta_update(1:Mexpl) ) ...
                      + model.dt^2 * theta_update(1:Mexpl)' ...
                                   * expl_reduced_data.Mmass(1:Mexpl,1:Mexpl) ...
                                   * theta_update(1:Mexpl)) );
    else
%    theta_update   = theta_update_expl + theta_update_impl;
      error('residuum can only be computed for the case when _one_ reduced basis space for both operators has been computed.');
    end
    reslog(t) = residuum;

    % compute the empirical interpolation error estimation
    eilog(t) = model.dt*lipschitzConst ...
                     * sqrt(abs( ...
                     theta_update(M+1:MM)' ...
                        * impl_reduced_data.Mmass(M+1:MM,M+1:MM) ...
                     * theta_update(M+1:MM) ));

    Delta    = Delta + eilog(t) ...
                     + lipschitzConst*newtonlog(t) ...
                     + lipschitzConst*residuum;
    if ~isreal(Delta)
      warning('RBMatlab:rb_simulation:runtime',...
              'error estimator returned complex result');
      Delta = real(Delta);
    end
    Deltalog(t) = Delta;
  end
end; % timeloop

if model.newton_solver && rmodel.verbose > 11
  disp(['computation ready: We needed at most ', num2str(nmax), ' Newton steps']);
end

if rmodel.enable_error_estimator
  rb_sim_data.Delta  = [0, Deltalog];
  rb_sim_data.reslog = [0, reslog];
  rb_sim_data.eilog  = [0, eilog];
end

% prepare return parameters
rb_sim_data.a = a;

if rmodel.compute_conditions
  rb_sim_data.conds.LI     = LIconds;
  rb_sim_data.conds.gradLI = gradLIconds;
end

end

function [reduced_data_op, tstop, Mop, Mred] = get_reduced_data_for_operator(id, rmodel, reduced_data, datanode_mode)

  if datanode_mode
    % Mind = get_index(rmodel.M, id, model.t);
    Mind = get_index(rmodel.M, id, get_mu(rmodel), rmodel.t);
    M    = get(rmodel.M, [Mind,1]);
  else
    M = rmodel.M;
  end
  MM = M + rmodel.Mstrich;

  %crbInd = get_index(reduced_data, id, rmodel);
  crbInd = get_index(reduced_data, id, get_mu(rmodel), rmodel.t);
  reduced_data_op = get(reduced_data, crbInd);
  tstop = index_valid_till(reduced_data, crbInd);
  Mop = MM;
  assert(Mop <= length(reduced_data_op.TM_local));

  Mred = min(Mop, M);

  if isempty(reduced_data_op.CE) ...
       || any(size(reduced_data_op.CE) ~= size(reduced_data_op.DE))
    reduced_data_op.CE     = reduced_data_op.DE / reduced_data_op.BM;
    reduced_data_op.CE_red = reduced_data_op.DE(:,1:Mred) ...
      / reduced_data_op.BM(1:Mred, 1:Mred);
  end
end