lin_evol_rb_derivative_simulation_t_part.m

function simulation_data = lin_evol_rb_derivative_simulation_t_part(model, reduced_data)
%function simulation_data = lin_evol_rb_derivative_simulation(model, reduced_data)
%
%performs in parallel a reduced simulation and calculation of the
%derivative @f$\partial_{\mu_i}u@f$
%
% This function is an extension of the \\c lin_evol_rb_simulation function.
% Except that in this method the time step is not included in the
% operators, but the time integration is conducted here by an explicit Euler.
% In parallel to the reduced basis online simulation a calculation of the
% derivative @f$\partial_{\mu_i}u@f$.
%
% It is assumed that the right mu is already set in the model.
% (setmu(model,mu))
%
% The calculation of the derivative is made in parallel by solving in each
% time step the ODE @f$\partial_t \partial_{\mu_i} a + L_E \partial_{\mu_i}
% a + \sum \partial_{\mu_i}(\Theta^q(\mu))L^q a -\partial_{mu_i} b(x,t,\mu) =
% 0@f$ plus derivative of the initial conditions.
% All information calculating the derivative is already contained in
% reduced_data or is already calculated for the normal rb_simulation except
% @f$\partial_{\mu_i}\Theta(\mu)@f$ and @f$\partial_{\mu_i}\Theta_b @f$
% which are obtained calling the field \\c model.rb_derivative_operators.
% The field \\c model.rb_derivative_init_values points to a method
% providing derivative initial values.
%
% 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
%
% Required fields of model
%  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
%  L_E_der_norm_bound: upper bound on explicit derivative operator
%                  constant in time, a single scalar
%  L_I_der_norm_bound: upper bound on implicit derivative 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
%  compute_derivative_info: 0 if no derivative information is needed
%                           1 if derivative information is needed
%   starting_tim_step:  in t-partition case this is the starting time step
%                       for the actual t-partitioni time step
%   stopping_time_step: in t-partition this is the stopping time step for
%                       the actual t-partition
%
% 
% 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)`
%   c: matrix(3-dimensional) of time secquence representing the derivative of a for
%   @f$\mu_i@f$
%   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 01.06.2010, Oliver Zeeb 14.10.2010

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

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

%t_part_preparation
if ~isfield(model,'transition_model')
    model.transition_model = 0;
end

if isfield(model,'starting_time_step')
    t_ind_start = model.starting_time_step;
else
    t_ind_start = 0;
end

if isfield(model, 'stopping_time_step')
    t_ind_stop = model.stopping_time_step;
else
    t_ind_stop = model.nt;
end
if ~model.transition_model
    nRB = length(reduced_data.a0{1});
else
    nRB = size(reduced_data.LL_E{1},2);
end


%needed for derivative information
if ~isfield(model,'optimization')
    disp('Model has no field optimization.')
    model.optimization =[];
end
if isfield(model.optimization,'params_to_optimize')
    indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
else
    indx_mu_i = 1:length(get_mu(model));
end
nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized

if isfield(model,'compute_derivative_indices')
    indx_derivatives = find(model.compute_derivative_indices);
else
    indx_derivatives = 1:length(get_mu(model));
end

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

%debug
res_u     = zeros(1,t_ind_stop-t_ind_start+1);
res_u_der     = zeros(nr_indx_mu_i,t_ind_stop-t_ind_start+1);
%end debug

if model.compute_derivative_info
    Delta_der = zeros(nr_indx_mu_i, t_ind_stop-t_ind_start+1); % vector of derivative error estimators
    c         = zeros(nRB, t_ind_stop-t_ind_start+1,nr_indx_mu_i); %matrix of RB-derivative-coefficients
end

%initial erro in t-part simulations
if isfield(reduced_data,'Delta0')
    Delta(1) = reduced_data.Delta0;
end
if isfield(reduced_data,'Delta_der0')
    Delta_der(:,1) = reduced_data.Delta_der0;
end


% 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;



if model.compute_derivative_info
    sc_0 = model.rb_derivative_init_values_coefficients(model);
    if isfield(reduced_data,'c0')
        for i=1:size(sc_0,1)
            c_0(:,:,i) = lincomb_sequence(reduced_data.c0{i},sc_0(i,:)); 
        end
    else
        for i=1:size(sc_0,1)
            c_0(:,:,i) = lincomb_sequence(reduced_data.a0,sc_0(i,:)'); 
        end
    end
    %use only parameter derivatives for parameters which are to optimize
    c(:,1,:)=c_0(:,:,indx_mu_i);
    
end

% create dummy files for derivative matrices --> memory preallocation
if model.compute_derivative_info 
model.t=0;
[dummy_LL_I, dummy_LL_E, dummy_bb, dummy_M_E, dummy_M_b, dummy_M_EE, ...
    dummy_M_Eb, dummy_M_bb, dummy_M_I, dummy_M_II, dummy_M_IE, ...
    dummy_M_Ib, dummy_M_EdEd, dummy_M_IdId, dummy_M_bdbd, ...
    dummy_M_Ed, dummy_M_Id, dummy_M_bd, dummy_M_IEd, dummy_M_IId, ...
    dummy_M_Ibd, dummy_M_EEd, dummy_M_EId, dummy_M_Ebd, dummy_M_EdId,...
    ]= rb_operators(model,[]);%dummy_M_Edbd, dummy_M_Idbd] = rb_operators(model,[]); --> dummy_M_Edbd, dummy_M_Idbd not needed

    if isempty(dummy_M_EdEd)
        M_EdEd = zeros(1,1,nr_indx_mu_i);
    else
        M_EdEd = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_IdId)
        M_IdId = zeros(1,1,nr_indx_mu_i);
    else
        M_IdId = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_Ed)
        M_Ed = zeros(1,1,nr_indx_mu_i);
    else
        M_Ed = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_Id)
        M_Id = zeros(1,1,nr_indx_mu_i);
    else
        M_Id = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_IEd)
        M_IEd = zeros(1,1,nr_indx_mu_i);
    else
        M_IEd = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_IId)
        M_IId = zeros(1,1,nr_indx_mu_i);
    else
        M_Id = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_EEd)
        M_EEd = zeros(1,1,nr_indx_mu_i);
    else
        M_EEd = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_EId)
        M_EId = zeros(1,1,nr_indx_mu_i);
    else
        M_EId = zeros(nRB,nRB,nr_indx_mu_i);
    end

    if isempty(dummy_M_EdId)
        M_EdId = zeros(1,1,nr_indx_mu_i);
    else
        M_EdId = zeros(nRB,nRB,nr_indx_mu_i);
    end
    
    M_bdbd = zeros(1  ,1  ,nr_indx_mu_i);
    M_bd   = zeros(1  ,nRB,nr_indx_mu_i);
    M_Ibd  = zeros(nRB,1  ,nr_indx_mu_i);
    M_Ebd  = zeros(nRB,1  ,nr_indx_mu_i);
    M_Edbd = zeros(nRB,1  ,nr_indx_mu_i);
    M_Idbd = zeros(nRB,1  ,nr_indx_mu_i);
end %end "if compute_derivative_info"
% end of: create dummy files for derivative matrices --> memory
% preallocation


% loop over time steps: computation of a(t+1)==a^t
for t = (t_ind_start+1):t_ind_stop
  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==t_ind_start+1) || (~model.data_const_in_time)
% 10.11.10 original:      
%     [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,[]);
      % 10.11.10 NEU:
      if model.compute_derivative_info% NEU
        [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,[]);
      else
          [sLL_I, sLL_E, sbb, sM_E, sM_b, sM_EE, sM_Eb, sM_bb, sM_I, sM_II, sM_IE, sM_Ib] ...
           = rb_operators(model,[]);
      end
      % END 10.11.10 NEU:
    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.compute_derivative_info
        for i=1:nr_indx_mu_i
            M_EdEd(:,:,i) = lincomb_sequence(reduced_data.M_EdEd,sM_EdEd(:,i));
            M_IdId(:,:,i) = lincomb_sequence(reduced_data.M_IdId,sM_IdId(:,i));
            M_bdbd(:,:,i) = lincomb_sequence(reduced_data.M_bdbd,sM_bdbd(:,i));
            M_Ed(:,:,i)   = lincomb_sequence(reduced_data.M_Ed  ,sM_Ed(:,i));
            M_Id(:,:,i)   = lincomb_sequence(reduced_data.M_Id  ,sM_Id(:,i));
            M_bd(:,:,i)   = lincomb_sequence(reduced_data.M_bd  ,sM_bd(:,i));
            M_IEd(:,:,i)  = lincomb_sequence(reduced_data.M_IEd ,sM_IEd(:,i));
            M_IId(:,:,i)  = lincomb_sequence(reduced_data.M_IId ,sM_IId(:,i));
            M_Ibd(:,:,i)  = lincomb_sequence(reduced_data.M_Ibd ,sM_Ibd(:,i));
            M_EEd(:,:,i)  = lincomb_sequence(reduced_data.M_EEd ,sM_EEd(:,i));
            M_EId(:,:,i)  = lincomb_sequence(reduced_data.M_EId ,sM_EId(:,i));
            M_Ebd(:,:,i)  = lincomb_sequence(reduced_data.M_Ebd ,sM_Ebd(:,i));
            M_EdId(:,:,i) = lincomb_sequence(reduced_data.M_EdId,sM_EdId(:,i));
            M_Edbd(:,:,i) = lincomb_sequence(reduced_data.M_Edbd,sM_Edbd(:,i));
            M_Idbd(:,:,i) = lincomb_sequence(reduced_data.M_Idbd,sM_Idbd(:,i));
        end
    end


    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
      LL_E = (speye(size(LL_I)) - LL_I) \ LL_E;
      bb   = (speye(size(LL_I)) - LL_I) \ bb;
      LL_I = zeros(size(LL_I)); 
    end;
  end;


  % solve (I-dt*LL_I) a(:,t+1) = (I+dt*LL_E) * a(:,t) + bb
  rhs = (speye(size(LL_E)) + model.dt*LL_E) * a(:,t-t_ind_start) + model.dt*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 isequal(LL_I, zeros(size(LL_I))) 
    a(:,t+1-t_ind_start) = rhs;
    %debug 28.10.10: folgende zeile mit reingenommen, damit
    %res_norm_times_dt_sqr ein skalar wird
    M_Ib=zeros(nRB,1);
    % end debug 28.10.10
  else % solve linear system
    a(:,t+1-t_ind_start) = (speye(size(LL_I)) - model.dt*LL_I)\rhs;
  end;
  
    
  % compute error estimate recursively
  % Delta^k := |L_I^(k-1)|^-1 * (res_norm(k-1)+
  % |L_E^(k-1)|*Delta^(k-1))
  res_norm_times_dt_sqr =   a(:,t+1-t_ind_start)'*a(:,t+1-t_ind_start) - 2*a(:,t+1-t_ind_start)' * a(:,t-t_ind_start) + a(:,t-t_ind_start)'*a(:,t-t_ind_start) ...
                            +2*model.dt*( a(:,t-t_ind_start)'*M_E*a(:,t-t_ind_start) - a(:,t-t_ind_start)'*M_E*a(:,t+1-t_ind_start) ...
                                        +M_b*a(:,t-t_ind_start) - M_b*a(:,t+1-t_ind_start) ...
                                        +a(:,t+1-t_ind_start)'*M_I*a(:,t-t_ind_start) - a(:,t+1-t_ind_start)'*M_I*a(:,t+1-t_ind_start)) ...
                            +((model.dt)^2)* ( a(:,t-t_ind_start)'*M_EE*a(:,t-t_ind_start) + a(:,t+1-t_ind_start)'*M_II*a(:,t+1-t_ind_start) + M_bb ...
                                               +2*( a(:,t+1-t_ind_start)'*M_IE*a(:,t-t_ind_start) +a(:,t-t_ind_start)'*M_Eb + a(:,t+1-t_ind_start)'*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
%    keyboard;
    disp('residual < 0')
    res_norm_times_dt = 0;
  end;
    
  %debug
    res_u(t-t_ind_start) = res_norm_times_dt;
  %end debug
  
  if isequal(model.error_norm,'l2')
    Delta(t+1-t_ind_start) = res_norm_times_dt + model.L_E_norm_bound * Delta(t-t_ind_start);
  else
    % only sum up the squares of the residuals, multiplication with coeff
    % and squareroots are taken at the end.
    keyboard;
    Rkplus1sqr = max(res_norm_sqr/(params.dt.^2), 0);
    Delta(t+1-t_ind_start) = Delta(t-t_ind_start)+Rkplus1sqr;
  end;
  
   
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %Calculation of derivative \partial_{\mu}a
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if model.compute_derivative_info
    model.decomp_mode=2;
  
    if (t==t_ind_start) || (~model.data_const_in_time)
        [sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model);
        %sLe = -sLe;
        %debug 28.10.10
        if ~isempty(reduced_data.LL_I)
            LL_I_derived=zeros(size(reduced_data.LL_I{1},1),size(reduced_data.LL_I{1},2),nr_indx_mu_i);
        else
            for j=1:nr_indx_mu_i
                LL_I_derived(:,:,j);
            end
        end
        %end debug 28.10.10
        LL_E_derived=zeros(size(reduced_data.LL_E{1},1),size(reduced_data.LL_E{1},2),nr_indx_mu_i);
        bb_derived=zeros(length(reduced_data.bb{1}),nr_indx_mu_i);
    
        %LL_E stays the same
        for i=1:nr_indx_mu_i
            if ~isempty(sLi) %only do the linear combination if the coefficients sLi are not empty
                LL_I_derived(:,:,i) = lincomb_sequence(reduced_data.LL_I,sLi(indx_mu_i(i),:));
            end
            LL_E_derived(:,:,i) = lincomb_sequence(reduced_data.LL_E,sLe(indx_mu_i(i),:));
            bb_derived(:,i) = lincomb_sequence(reduced_data.bb,sb(indx_mu_i(i),:));
        end
    end  

    % solve d/dt c(:,:,i) = LL_E * c(:,:,i) + LL_E_derived a() +
    % bb_derived
    for i=1:nr_indx_mu_i
        rhs_der = (speye(size(LL_E)) + model.dt*LL_E) * c(:,t-t_ind_start,i) + model.dt*LL_E_derived(:,:,i) * a(:,t-t_ind_start) ...
                + model.dt*LL_I_derived(:,:,i) * a(:,t+1-t_ind_start) + model.dt*bb_derived(:,i);
                
         % check whether pure explicit problem:  
        if isequal(LL_I, zeros(size(LL_I))) 
            c(:,t+1-t_ind_start,i) = rhs_der;
        else %solve linear system
            c(:,t+1-t_ind_start,i) = (speye(size(LL_I)) - model.dt*LL_I)\rhs_der;
        end
      
      
        %compute error estrimator recuresivly
        res_norm_der_times_dt_sqr = c(:,t+1-t_ind_start,i)'*c(:,t+1-t_ind_start,i) ...
          + c(:,t-t_ind_start,i)'*c(:,t-t_ind_start,i) - 2*c(:,t+1-t_ind_start,i)'*c(:,t-t_ind_start,i) ...
          + 2*model.dt*(-c(:,t+1-t_ind_start,i)'*M_I'*c(:,t+1-t_ind_start,i) - c(:,t+1-t_ind_start,i)'*M_E'*c(:,t-t_ind_start,i) ...
                        -c(:,t+1-t_ind_start,i)'*M_Ed(:,:,i)'*a(:,t-t_ind_start) - c(:,t+1-t_ind_start,i)'*M_Id(:,:,i)'*a(:,t+1-t_ind_start)...
                        -c(:,t+1-t_ind_start,i)'*M_bd(:,:,i)' + c(:,t+1-t_ind_start,i)'*M_I*c(:,t-t_ind_start,i) ...
                        +c(:,t-t_ind_start,i)'*M_E*c(:,t-t_ind_start,i) + c(:,t-t_ind_start,i)'*M_Ed(:,:,i)'*a(:,t-t_ind_start) ...
                        +c(:,t-t_ind_start,i)'*M_Id(:,:,i)*a(:,t+1-t_ind_start) + c(:,t-t_ind_start,i)'*M_bd(:,:,i)' ) ... %13.10.10: bei M_bd aufpassen, nur ein Doppelpunkt?!
          + ((model.dt)^2) * (c(:,t+1-t_ind_start,i)'*M_II*c(:,t+1-t_ind_start,i) + c(:,t-t_ind_start,i)'*M_EE*c(:,t-t_ind_start,i) ...
                              + a(:,t-t_ind_start)'*M_EdEd(:,:,i)*a(:,t-t_ind_start) + a(:,t+1-t_ind_start)'*M_IdId(:,:,i)*a(:,t+1-t_ind_start)...
                              + M_bdbd(:,:,i) + 2*c(:,t+1-t_ind_start,i)'*M_IE*c(:,t-t_ind_start,i) ...
                              + 2*c(:,t+1-t_ind_start,i)'*M_IEd(:,:,i)*a(:,t-t_ind_start) + 2*c(:,t+1-t_ind_start,i)'*M_IId(:,:,i)*a(:,t+1-t_ind_start)...
                              + 2*c(:,t+1-t_ind_start,i)'*M_Ibd(:,:,i) + 2*c(:,t-t_ind_start,i)'*M_EEd(:,:,i)*a(:,t-t_ind_start) ...
                              + 2*c(:,t-t_ind_start,i)'*M_EId(:,:,i)*a(:,t+1-t_ind_start) + 2*c(:,t-t_ind_start,i)'*M_Ebd(:,:,i) ...
                              + 2*a(:,t-t_ind_start)'*M_EdId(:,:,i)*a(:,t+1-t_ind_start) + 2*a(:,t-t_ind_start)'*M_Edbd(:,:,i)...
                              + 2*a(:,t+1-t_ind_start)'*M_Idbd(:,:,i) );                          
   
    % ensure real estimators (could be complex due to numerics)
    if res_norm_der_times_dt_sqr>=0
        res_norm_der_times_dt = sqrt(res_norm_der_times_dt_sqr);
    else
        res_norm_der_times_dt = 0;
    end;
   
    %debug
        res_u_der(i,t-t_ind_start) = res_norm_der_times_dt;
    %end debug 
    
    if isequal(model.error_norm,'l2')
        if model.time_const
            t_old = model.t;
            model.t=0;
        end
       Delta_der(i,t+1-t_ind_start) = model.L_E_norm_bound*Delta_der(i,t-t_ind_start) + model.L_E_der_norm_bound(model,indx_derivatives(i))*Delta(t-t_ind_start) ...
           + model.L_I_der_norm_bound(model,indx_derivatives(i))*Delta(t+1-t_ind_start) + res_norm_der_times_dt ;
       if model.time_const
           model.t = t_old;
       end
    end
    
    
    end %for i=1:nr_indx_mu
  end %if model.compute_derivative_info
end;

if isequal(model.error_norm,'energy')
  % 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;
simulation_data.bb = bb;
simulation_data.nRB = nRB;

%debug
simulation_data.res_u = res_u;
simulation_data.res_u_der = res_u_der;
%end debug

if model.compute_derivative_info
    simulation_data.Delta_der = Delta_der;
    simulation_data.c = c;
end


%output calculation
if isfield(model,'name_output_functional')
  if isequal(model.error_norm,'l2')
        simulation_data = model.get_output(model, simulation_data,reduced_data);
%     s = reduced_data.s_RB'*a;
%     if isfield(model,'s_l2norm')
%         Delta_s = model.s_l2norm * Delta; %war ursprünglich in reduced_data.s_l2norm
%         simulation_data.Delta_s = Delta_s;
%     else
%         if (model.verbose >=10)
%             disp('warning: no L2_norm for output defined in lin_evol_rb_simulation...')
%         end
%     end
%     simulation_data.s = s;
%     
%     if model.compute_derivative_info
%         %derivative output and error estimation of gradient
%         %sum=0;%for l2-norm
%         for i=1:size(c,3)
%             s_der{i} = reduced_data.s_RB'*c(:,:,i);
%             %sum = sum + max(Delta_der(i,:))^2;
%         end
%         simulation_data.s_der = s_der;
%         simulation_data.Delta_s_grad = model.s_l2norm .*max(Delta_der,[],2); %Gradient error bound
%         %simulation_data.Delta_s_der = model.s_l2norm * Delta_der;
%     end
    
  else
    error(['output estimation not implemented concerning energy-norm.' ...
           ' choose l2 as error_norm!']);
  end;
end;

%| \docupdate