lin_stat_rb_simulation.m

function rb_sim_data = lin_stat_rb_simulation(model,reduced_data)
%function rb_sim_data = lin_stat_rb_simulation(model,reduced_data)
% 
% function performing a reduced simulation

% B. Haasdonk 22.2.2011

rb_sim_data =[];

old_mode = model.decomp_mode;
model.decomp_mode = 2; % coefficients
[A_coeff,f_coeff] = model.operators(model,[]);

N = reduced_data.N;
% for cell arrays:
%AN = lincomb_sequence(reduced_data.AN_comp, A_coeff);
%fN = lincomb_sequence(reduced_data.fN_comp, f_coeff);
AN = reshape(reduced_data.AN_comp * A_coeff,N,N);
fN = reduced_data.fN_comp * f_coeff;

% solution variable:
%uh = femdiscfunc([],model_data.df_info);
%uh.dofs = A\r;    
uN = AN\fN;

% return results:
rb_sim_data.uN = uN;

% plus error estimator
% res_norm = ... % residual norm

if isfield(model, 'RB_error_indicator') ...
        && strcmp(model.RB_error_indicator, 'estimator')
    
    % for elliptic compliant case, X-norm (=H10-norm) error estimator:
    Q_r = size(reduced_data.G,1);
    neg_auN_coeff = -uN * A_coeff'; 
    %neg_auN_coeff = -A_coeff * uN'; 
    res_coeff = [f_coeff; neg_auN_coeff(:)];
    res_norm_sqr = res_coeff' * reduced_data.G * res_coeff;
    
    % direct computation (expensive):
    if 0
        % for debugging:
        neg_auN_coeff = neg_auN_coeff(:);
        Q_f = length(f_coeff);
        res_norm_sqr_ff = f_coeff' * reduced_data.G(1:Q_f,1:Q_f) * f_coeff;
        res_norm_sqr_fAu = f_coeff' * reduced_data.G(1:Q_f,(Q_f+1):end) * neg_auN_coeff;
        res_norm_sqr_Auf = neg_auN_coeff' * reduced_data.G((Q_f+1):end,1:Q_f) * f_coeff;
        res_norm_sqr_AuAu = neg_auN_coeff' * reduced_data.G((Q_f+1):end,(Q_f+1):end) * neg_auN_coeff;
        
        model_data = gen_model_data(model);
        detailed_data = gen_detailed_data(model,model_data);
        model.decomp_mode = 0;
        [A,f] = model.operators(model,model_data);
        model.decomp_mode = 2;
        resAu = A * (detailed_data.RB(:,1:model.N) * uN);
        resf = -f;
        res = resAu + resf;
        % residuum functional is res * v 
        % riesz representant: v_r' K v = (v_r,v) = res*v
        % so res = K * v_r;
        K = model.get_inner_product_matrix(detailed_data);
        v_r = K \ res;
        v_rAu = K \ resAu;
        v_rf= K \ resf;
        % res_norm_sqr = (v_r,v_r) = v_r' K v_r = v_r' * res;
        res_norm_sqr2 = v_r' * res;
        res_norm_sqr2_ff = v_rf' * resf;
        res_norm_sqr2_fAu = v_rf' * resAu;
        res_norm_sqr2_Auf = v_rAu' * resf;
        res_norm_sqr2_AuAu = v_rAu' * resAu;
        keyboard;
    end;
    
    % prevent possibly negative numerical disturbances:
    res_norm_sqr = max(res_norm_sqr,0);    
    res_norm = sqrt(res_norm_sqr);
    % if the SCM is being used perform an online-phase and use the resulting
    % lower bound. Otherwise use the old code.
    if model.use_scm == 0
        rb_sim_data.Delta = ...
            res_norm/model.coercivity_alpha(model); 
        rb_sim_data.Delta_s = ...
            res_norm_sqr/model.coercivity_alpha(model); 
    elseif model.use_scm == 1
        if ~isfield(reduced_data, 'scm_offline_data')
            error('There are no scm_offline_data to work with. They must be included in the reduced_data!')
        end
        constant_LB = scm_lower_bound(model, reduced_data);
        rb_sim_data.Delta = ...
            res_norm/constant_LB; 
        rb_sim_data.Delta_s = ...
            res_norm_sqr/constant_LB; 
    end
    
end % error estimation

if model.compute_output_functional
  % get operators
 l_coeff = model.operators_output(model,reduced_data);
%  lN = lincomb_sequence(reduced_data.lN_comp,l_coeff);
  lN = reduced_data.lN_comp * l_coeff;
  rb_sim_data.s = lN(:)' * rb_sim_data.uN;
end;

model.decomp_mode = old_mode;