lin_evol_opt_output_time_integrated.m

function sim_data = lin_evol_opt_output_time_integrated(model, sim_data, varargin)
%function sim_data = lin_evol_opt_output_time_integrated(model, sim_data, varargin)
%
% This function returns the time integrated output
%
% J(u) = 1/Tlen \int_0^T\int_{\Omega} f u dOmega dt
%
% where f is a suitable kernel defined by model.output_ptr ...oder so
% aehnlich
%

% Markus Dihlmann 01.08.2012


%decide if it's a p-part model
if isfield(model,'p_part_ranges')
    %yes_p_part model
    model = model.base_model;
end


%calculate the length of the time interval
if isfield(model,'save_time_indices')
    Tlen = (model.stime_stop-model.stime_start)/length(model.save_time_indices);
elseif isfield(model,'base_model')&&(isfield(model.base_model,'save_time_indices'))
    Tlen = (model.stime_stop-model.stime_start)/length(model.base_model.save_time_indices);
else
    Tlen = (model.stime_stop - model.stime_start)/model.nt;
end
    

if ~isfield(model, 'starting_time_step')
    model.starting_time_step = 0;
end
if ~isfield(model,'stopping_time_step')
    model.stopping_time_step = model.nt;
end
if isfield(model,'initial_value_output')
    initial_value_output = model.initial_value_output;
    if isfield(model, 'initial_value_output_Delta_s')
        initial_value_output_Delta_s = model.initial_value_output_Delta_s;
        initial_value_output_Delta_s_l2 = model.initial_value_output_Delta_s_l2;
    end
    if model.compute_derivative_info
        initial_value_output_gradient = model.initial_value_output_gradient;
        if isfield(model, 'initial_value_output_Delta_s_grad')
            initial_value_output_Delta_s_grad = model.initial_value_output_Delta_s_grad;
            initial_value_output_Delta_s_grad_l2 = model.initial_value_output_Delta_s_grad_l2;
        end
    end
else
    initial_value_output = 0;
    initial_value_output_Delta_s = 0;
    initial_value_output_Delta_s_l2 = 0;
    if model.compute_derivative_info
        if isfield(sim_data,'c')
            initial_value_output_gradient = zeros(length(sim_data.c),1);
            initial_value_output_Delta_s_grad = zeros(length(sim_data.c),1);
            initial_value_output_Delta_s_grad_l2 = zeros(length(sim_data.c),1);
        elseif isfield(sim_data,'U_der')
            initial_value_output_gradient = zeros(length(sim_data.U_der),1);
        end
        
    end
end


if isfield(varargin{1},'grid')
    % --> model_data
    model_data = varargin{1};
    
    U=sim_data.U;
    nt = size(U,2);
    
    %distance between two time steps
    fac = model.save_time_indices(end-1) - model.save_time_indices(end-2);
    if ~(nt == ceil(model.nt/fac))
        warning(['dt  and fac in get output time integration is wrong... maybe... error',...
            num2str(nt-ceil(model.nt/fac))])
    end
    
    %in case of t-partition model...
    t_ind = find(((model.save_time_indices)>=model.starting_time_step)&...
        ((model.save_time_indices)<=model.stopping_time_step));
    %t_ind = t_ind;
    
    t_s_ind = find(((model.save_time_indices)>=model.stime_start)&...
        ((model.save_time_indices)<=model.stime_stop));
    t_eval=find((t_ind>=t_s_ind(1))&(t_ind<=t_s_ind(end)));
    nt_eval = length(t_eval);
    
    %output field
    model.decomp_mode =1;
    v = model.operators_output(model, model_data);
    
    y = initial_value_output.*ones(1,nt);
    
    %time integration
    if ~isempty(t_eval)
        y(t_eval(1)) = fac.*model.dt.*(v{1}')*U(:,t_eval(1))+ initial_value_output;
        for t = 2:nt_eval
            y(t_eval(t)) = y(t_eval(t)-1) + fac/Tlen.*model.dt.*(v{1}')*U(:,t_eval(t));
        end
        y(t_eval(end)+1:end) =y(t_eval(end)); 
        sim_data.y = y;
        if model.compute_derivative_info

            U_der = sim_data.U_der;
            y_der = repmat(initial_value_output_gradient,1,nt);
            for i=1:length(U_der)
                y_der(i,t_eval(1)) = fac/Tlen.*model.dt.*(v{1}')*U_der{i}(:,t_eval(1))+ initial_value_output_gradient(i);
                for t =2:nt_eval
                    y_der(i,t_eval(t)) = y_der(i,t_eval(t)-1)+fac.*model.dt.*(v{1}')*U_der{i}(:,t_eval(t));
                end
                y_der(i,t_eval(end)+1:end) = y_der(i,t_eval(end));
            end
            sim_data.y_der = y_der;
        end
    else
        sim_data.y = y;
        if model.compute_derivative_info
            sim_data.y_der = repmat(initial_value_output_gradient,1,nt);
        end
    end
    
else
     
    %--> reduced_data
    reduced_data  = varargin{1};
    a = sim_data.a;
    nt = size(a,2);
    s = initial_value_output.*ones(1,nt);
     
    %fix start end end time of time integration
    %t1 = model.stime_start;
    %t2 = model.stime_stop;
   
    %in case of t-part model this interval is not always included in the
    %t-partition interval
    %so find out which time indices will be simulated
    t_ind = model.starting_time_step:model.stopping_time_step;
    %t_ind = t_ind+1;
    if length(t_ind)~=nt
        error('Hier stimmt was nicht')
    end
    t_eval=find((t_ind>=model.stime_start)&(t_ind<=model.stime_stop));
    nt_eval = length(t_eval);
    
    %time integration
    if ~isempty(t_eval)
        s(t_eval(1)) = model.dt.*reduced_data.s_RB'*a(:,t_eval(1))+initial_value_output;
        for i =2:nt_eval
            s(t_eval(i)) = s(t_eval(i)-1) + model.dt/Tlen.*reduced_data.s_RB'*a(:,t_eval(i));
        end
        s(t_eval(end)+1:end) = s(t_eval(end));
        
        
        
        if isfield(model,'s_l2norm')
            Delta_s = initial_value_output_Delta_s.*ones(1,nt);
            Delta_s(t_eval(1)) = model.dt.*model.s_l1norm * sim_data.Delta(t_eval(1))+initial_value_output_Delta_s;
            for i = 2:nt_eval
                Delta_s(t_eval(i)) = Delta_s(t_eval(i)-1) + model.dt.*model.s_l1norm * sim_data.Delta(t_eval(i));
            end
            Delta_s(t_eval(end):end) = Delta_s(t_eval(end));
            
            Delta_s_l2 = (initial_value_output_Delta_s_l2)^2.*ones(1,nt);
            Delta_s_l2(t_eval(1)) = model.dt.*model.s_l2norm * (sim_data.Delta(t_eval(1)))^2+(initial_value_output_Delta_s_l2)^2;
            for i = 2:nt_eval
                Delta_s_l2(t_eval(i)) = Delta_s_l2(t_eval(i)-1) + model.dt.*model.s_l2norm * (sim_data.Delta(t_eval(i)))^2;
            end
            Delta_s_l2(t_eval(end):end) = Delta_s_l2(t_eval(end));
            Delta_s_l2 = sqrt(Delta_s_l2);
            
        else
            if (model.verbose >=5)
                disp('warning: no L2_norm for output defined in lin_evol_rb_simulation...')
            end
        end
        
        
        if model.compute_derivative_info
            %derivative output and error estimation of gradient
            %sum=0;%for l2-norm
            c = sim_data.c;
            Delta_der = sim_data.Delta_der;
            Delta_s_grad = zeros(length(c),nt);
            Delta_s_grad_l2 = zeros(length(c),nt);
            s_der = cell(1,length(c));
            for i=1:length(c)
                s_der{i} = initial_value_output_gradient(i).*ones(1,nt);
                s_der{i}(t_eval(1)) =  model.dt/Tlen.*reduced_data.s_RB_der{i}'*c{i}(:,t_eval(1))+initial_value_output_gradient(i);
                for t = 2:nt_eval
                    s_der{i}(t_eval(t)) = s_der{i}(t_eval(t)-1) + model.dt.*reduced_data.s_RB_der{i}'*c{i}(:,t_eval(t));
                end
                s_der{i}(t_eval(end):end) = s_der{i}(t_eval(end));
                %sum = sum + max(Delta_der(i,:))^2;
                
                Delta_s_grad(i,:) = initial_value_output_Delta_s_grad(i).*ones(1,nt);
                Delta_s_grad(i,t_eval(1)) = model.dt.*model.s_l1norm.*Delta_der(i,t_eval(1))+initial_value_output_Delta_s_grad(i);
                for t=2:nt_eval
                    Delta_s_grad(i,t_eval(t)) = Delta_s_grad(i,t_eval(t-1))+ model.dt.*model.s_l1norm.*Delta_der(i,t_eval(t));
                end
                Delta_s_grad(i,t_eval(end):end) = Delta_s_grad(i,t_eval(end));
                
                Delta_s_grad_l2(i,:) = (initial_value_output_Delta_s_grad_l2(i))^2.*ones(1,nt);
                Delta_s_grad_l2(i,t_eval(1)) = model.dt.*model.s_l2norm.*(Delta_der(i,t_eval(1)))^2+(initial_value_output_Delta_s_grad_l2(i))^2;
                for t=2:nt_eval
                    Delta_s_grad_l2(i,t_eval(t)) = Delta_s_grad_l2(i,t_eval(t-1))+ model.dt.*model.s_l2norm.*(Delta_der(i,t_eval(t)))^2;
                end
                Delta_s_grad_l2(i,t_eval(end):end) = Delta_s_grad_l2(i,t_eval(end));
                Delta_s_grad_l2(i,:) = sqrt(Delta_s_grad_l2(i,:));
                
            end
            sim_data.s_der = s_der;
            sim_data.Delta_s_grad = Delta_s_grad;
            sim_data.Delta_s_grad_l2 = Delta_s_grad_l2;
            %simulation_data.Delta_s_grad = model.s_l2norm * Delta_der;
        end
    else
        %no time integration in this simulation time interval
        s = initial_value_output.*ones(1,nt);
        Delta_s = initial_value_output_Delta_s.*ones(1,nt);
        Delta_s_l2 = initial_value_output_Delta_s_l2.*ones(1,nt);
        if model.compute_derivative_info
            Delta_s_grad = repmat(initial_value_output_Delta_s_grad,1,nt);
            Delta_s_grad_l2 = repmat(initial_value_output_Delta_s_grad_l2,1,nt);
            for i=1:length(initial_value_output_gradient)
                s_der{i} = initial_value_output_gradient(i).*ones(1,nt);
            end
        end
    end
    
    
    sim_data.s = s;
    sim_data.Delta_s = Delta_s;
    sim_data.Delta_s_l2 = Delta_s_l2;
    if model.compute_derivative_info
        sim_data.s_der = s_der;
        sim_data.Delta_s_grad = Delta_s_grad;
        sim_data.Delta_s_grad_l2 = Delta_s_grad_l2;
    end
end