lin_evol_opt_rb_simulation.m

function simulation_data = lin_evol_opt_rb_simulation(model, reduced_data)
%function simulation_data = lin_evol_opt_rb_simulation(model, reduced_data)
%
% function, which performs a reduced basis online simulation for
% the parameter vector mu, which is assumed to be set in the model
% class, i.e. a previous model = model.set_mu(model,[...]) is assumed
% to have taken place.
%
% allowed dependency of data: Nmax, N, M, mu
% not allowed dependency of data: H
% allowed dependency of computation: Nmax, N, M, mu
% not allowed dependency of computation: H
% Unknown at this stage: ---
%
% Required fields of model as required by the numerical operators and
%  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
%  error_norm     : 'l2' or 'energy'
%  L_I_inv_norm_bound: upper bound on implicit operator inverse norm
%                  constant in time, a single scalar
%  L_E_norm_bound: upper bounds on explicit operator
%                  constant in time, a single scalar
%  energy_norm_gamma: gamma >= 0 defining the weight in the energy norm
%  coercivity_bound_name : name of the function, which can
%                  computes a lower bound on the coercivity, e.g.
%                  fv_coercivity_bound
%
% plus fields required by coercivity_bound_name.m
%
% optional fields of model:
%  data_const_in_time : if this flag is set, then only operators
%                       for first time instance are computed
%  name_output_functional: if this field is existent, then an
%        output estimation is performed and error.estimtations
%
% generated fields of simulation_data:
%   a: time sequence of solution coefficients, columns are
%      'a(:,k)' = `a^(k-1)`
%   Delta:time sequence of `L^2`-posteriori error estimates
%            'Delta(k)' = `\Delta^{k-1}_N`
%         or energy-norm-posterior error estimates
%            'Delta_energy(k)' = `\bar \Delta^{k-1}_N`
%         depending on the field "error_norm" in model.
%
% if model.name_output_functional is set, then additionally, a
% sequence of output estimates 's(U(:,))' and error bound Delta_s is returned
%
% see the ***_gen_reduced_data() routine for specifications of the
% fields of reduced_data

% Markus Dihlmann 11.06.2010

if model.verbose >= 10
 printf('performing simulation for mu = %s', mat2str(get_mu(model)));
end;

if (~isequal(model.error_norm,'l2')) && ...
      (~isequal(model.error_norm,'energy'))
  error('error_norm unknown!!');
end;

%cgrid = grid_geometry(params);
nRB = length(reduced_data.a0{1});
a = zeros(nRB,model.nt+1);   % matrix of RB-coefficients
Delta = zeros(1,model.nt+1); % vector of error estimators
model.dt = model.T/model.nt;

% initial data projection and linear combination
model.decomp_mode = 2;
sa0 = rb_init_values(model,[]);
a0 = lincomb_sequence(reduced_data.a0,sa0);
a(:,1) = a0;

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

  % linear combination of all quantities
  % only once, if data is const in time

  if (t==1) || (~model.data_const_in_time)
    [sLL_I, sLL_E, sbb, sM_E, sM_b, sM_EE, sM_Eb, sM_bb, sM_I, sM_II, sM_IE, sM_Ib, ...
        sM_EdEd, sM_IdId, sM_bdbd, sM_Ed, sM_Id, sM_bd, sM_IEd, sM_IId, sM_Ibd, sM_EEd, ...
        sM_EId, sM_Ebd, sM_EdId, sM_Edbd, sM_Idbd] = rb_operators(model,[]);
    LL_I = lincomb_sequence(reduced_data.LL_I,sLL_I);
    LL_E = lincomb_sequence(reduced_data.LL_E,sLL_E);
    bb   = lincomb_sequence(reduced_data.bb,sbb);
    M_E  = lincomb_sequence(reduced_data.M_E,sM_E);
    M_b  = lincomb_sequence(reduced_data.M_b,sM_b);
    M_EE = lincomb_sequence(reduced_data.M_EE,sM_EE);
    M_Eb = lincomb_sequence(reduced_data.M_Eb,sM_Eb);
    M_bb = lincomb_sequence(reduced_data.M_bb,sM_bb);
    M_I  = lincomb_sequence(reduced_data.M_I,sM_I);
    M_II = lincomb_sequence(reduced_data.M_II,sM_II);
    M_Ib = lincomb_sequence(reduced_data.M_Ib,sM_Ib);
    M_IE = lincomb_sequence(reduced_data.M_IE,sM_IE);
    if model.data_const_in_time
      % inverse of implicit matrix! L_I' = Id, L_E' = L_I^(-1) *
      % L_E, bb = L_I^(-1) bb
%      ILL_I = inv(LL_I);
      LL_E = LL_I \ LL_E;
      bb = LL_I \ bb;
      LL_I = speye(size(LL_I));
      LL_I_is_speye = 1;
    elseif isequal(LL_I, speye(size(LL_I)))
      LL_I_is_speye = 1;
    end;
  end;

  % solve LL_I* a(:,t+1) = LL_E * a(:,t) + bb
  rhs = LL_E * a(:,t) + bb;
%  warning('this does not work for pure space operators (independent of delta t). In that case use the commented line below.');
  %rhs = a(:,t) + model.dt * (LL_E * a(:,t) + bb);
  % check whether pure explicit problem:
  if (LL_I_is_speye == 1)
    a(:,t+1) = a(:,t)+model.dt*rhs;
  elseif (LL_I_is_speye == -1) && (max(max(LL_I-speye(size(LL_I))))<1.0e-3)
    a(:,t+1) = a(:,t)+model.dt*rhs;
  else % solve linear system
       %    disp('check symmetry and choose solver accordingly!');
       %    keyboard;
       % nonsymmetric solvers:
       %  [U(:,t+1), flag] = bicgstab(Li,rhs,[],1000);
       %  [U(:,t+1), flag] = cgs(Li,rhs,[],1000);
       % symmetric solver, non pd:
       % see bug_symmlq for a very strange bug: cannot solve identity system!
       % [U(:,t+1), flag] = symmlq(Li,rhs,[],1000);
       %  [U(:,t+1), flag] = minres(Li,rhs,[],1000);
       % symmetric solver, pd:
       % [a(:,t+1), flag] = pcg(LL_I,rhs,[],1000);
       [a(:,t+1), flag] = bicgstab(LL_I,rhs,[],1000);
       if flag>0
         disp(['error in system solution, solver return flag = ', ...
               num2str(flag)]);
         keyboard;
       end;
  end;

  % compute error estimate recursively
  % 

  
  res_norm_times_dt_sqr = a(:,t+1)'*a(:,t+1) - 2*a(:,t+1)' * a(:,t) + a(:,t)'*a(:,t) ...
                            +2*model.dt*( a(:,t)'*M_E*a(:,t) - a(:,t)'*M_E*a(:,t+1) ...
                                        +M_b*a(:,t) - M_b*a(:,t+1) ...
                                        +a(:,t+1)'*M_I*a(:,t) - a(:,t+1)'*M_I*a(:,t+1)) ...
                            +((model.dt)^2)* ( a(:,t)'*M_EE*a(:,t) + a(:,t+1)'*M_II*a(:,t+1) + M_bb ...
                                               +2*( a(:,t+1)'*M_IE*a(:,t) +a(:,t)'*M_Eb + a(:,t+1)'*M_Ib));
           
  % ensure real estimators (could be complex due to numerics)
  if res_norm_times_dt_sqr>=0
    res_norm_times_dt = sqrt(res_norm_times_dt_sqr);
  else
    res_norm_times_dt = 0;
  end;

  if isequal(model.error_norm,'l2')
    Delta(t+1) =    res_norm_times_dt + model.L_E_norm_bound * Delta(t);
  else
    % only sum up the squares of the residuals, multiplication with coeff
    % and squareroots are taken at the end.
    disp('warning: error calculation not defined!!! See lin_evol_opt_rb_simulation line 161');
    %Rkplus1sqr = max(res_norm_sqr/(params.dt.^2), 0);
    %Delta(t+1) = Delta(t)+Rkplus1sqr;
  end;
end;

if isequal(model.error_norm,'energy')
  disp('warning: error calculation for energynorm  not defined!!! See lin_evol_opt_rb_simulation line 168');  
  % now perform scaling and squareroot of energy-error-estimate:
  %alpha = model.coercivity_bound_ptr(model);
  %if max(model.L_E_norm_bound>1) || (alpha == 0)
    % no reasonable estimation possible
   % Delta(:) = nan;
  %else
  %  C = (sqrt(1-model.L_E_norm_bound.^2)+1) * 0.5;
  %  Coeff = model.dt /(4*alpha*C*(1-model.energy_norm_gamma * C));
  %  Delta = sqrt(Delta * Coeff);
  %end;
end;

% return simulation result
simulation_data       = [];
simulation_data.a     = a;
simulation_data.Delta = Delta;
simulation_data.LL_I  = LL_I;
simulation_data.LL_E  = LL_E;

if isfield(model,'name_output_functional')
    if isequal(model.error_norm,'l2')
    s = a' * reduced_data.s_RB;
    if isfield(reduced_data,'s_l2norm')
        Delta_s = reduced_data.s_l2norm * Delta;
        simulation_data.Delta_s = Delta_s(:);
    else
        if(model.verbose>=10)
            disp('warning: no L2_norm for output defined in lin_evol_opt_rb_simulation')
        end
    end
    simulation_data.s = s(:);
    
  else
    error(['output estimation not implemented concerning energy-norm.' ...
           ' choose l2 as error_norm!']);
  end;
end;

%| \docupdate