rb_simulation_impl.m

function rb_sim_data = rb_simulation_impl(rmodel, reduced_data, sim_data, dd)
%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 TwoPhaseFlow.EiRbReducedDataNode
%
% generated fields of rb_sim_data:
%  a     : 'a(:,k)' is the coefficient vector of the reduced simulation at time
%          `t^{k-1}`
%

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

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

% set init data
rmodel.decomp_mode = 2;

sat0Ind = get_index(reduced_data, 'sat_init', get_mu(rmodel), rmodel.t);
rd_sat0 = get(reduced_data, sat0Ind);

descr = rmodel.descr;

sat_a0  = lincomb_sequence(rd_sat0, descr);

% perform simulation loop in time
descr.dt = descr.T/descr.nt;

sat_N = get_N(rmodel, 'saturation');
vel_N = get_N(rmodel, 'velocity');
prs_N = get_N(rmodel, 'pressure');
descr.prs_N = prs_N;

sat_a   = zeros(sat_N, descr.nt+1);   % matrix of saturation RB-coefficients
vel_a   = zeros(vel_N, descr.nt+1);   % matrix of velocity   RB-coefficients
prs_a   = zeros(prs_N, descr.nt+1);   % matrix of pressure   RB-coefficients
exflags = zeros(1, descr.nt+1);
outtext = cell(1, descr.nt+1);

sat_a(:,1) = sat_a0(1:sat_N);

[rd_sat] = get_reduced_data_for_operator('L_saturation', rmodel, reduced_data);
[rd_vel] = get_reduced_data_for_operator('L_velocity', rmodel, reduced_data);
[rd_prs] = get_reduced_data_for_operator('L_pressure_mean', rmodel, reduced_data);
rd_divInd = get_index(reduced_data, 'L_divergence', get_mu(rmodel), rmodel.t);
rd_div = get(reduced_data, rd_divInd);

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

  V_last = [sat_a(:,t); vel_a(:,t); prs_a(:,t)];
  descr.t = descr.t + descr.dt;

  fsolve_fun = @(X) sat_fun(descr, sat_a(:,t), X, rd_sat, rd_vel, rd_div, rd_prs);

  gplmmopts.Px = @(X) box_projection(X, sat_N);
  gplmmopts.debug = false;
  gplmmopts.TolRes = 1e-7;
  gplmmopts.TolChange = 1e-10;
%  gplmmopts.HessFun = @(params, jac, dummy, arg) jac' * (jac * arg);
  [V_new, resnorm, FV, exitflag, output] = gplmm(fsolve_fun, V_last, gplmmopts);

  exflags(t) = exitflag;
  outtext{t} = output;

  sat_a(:, t+1) = V_new(1:sat_N);
  vel_a(:, t+1) = V_new(sat_N+1:sat_N+vel_N);
  prs_a(:, t+1) = V_new(sat_N+vel_N+1:end);

end

if rmodel.enable_error_estimator
  error('error estimator is not implemented yet.');
  rb_sim_data.Delta = 0;
end

% prepare return parameters
rb_sim_data.sat_a   = sat_a;
rb_sim_data.vel_a   = vel_a;
rb_sim_data.prs_a   = prs_a;
rb_sim_data.exit_flags = exflags;
rb_sim_data.outputs    = outtext;

end

function [reduced_data_op] = get_reduced_data_for_operator(id, rmodel, reduced_data)

  crbInd = get_index(reduced_data, id, get_mu(rmodel), rmodel.t);
  reduced_data_op = get(reduced_data, crbInd);


  if isa(reduced_data_op, 'TwoPhaseFlow.EiRbReducedDataNode') ...
       && ( isempty(reduced_data_op.CE) ...
            || any(size(reduced_data_op.CE) ~= size(reduced_data_op.DE)) )

    M = reduced_data_op.M;
    reduced_data_op.CE     = reduced_data_op.DE / reduced_data_op.BM;
    reduced_data_op.CE_red = reduced_data_op.DE(:,1:M) ...
      / reduced_data_op.BM(1:M, 1:M);
  end
end

function [ret, J] = sat_fun(descr, Sold, X, rd_sat, rd_vel, rd_div, rd_prs)

  sat_N = length(Sold);
  vel_N = size(rd_sat.RB_local_ext.U, 2);
  prs_N = descr.prs_N;

  Xsat.S = rd_sat.RB_local_ext.S * X(1:sat_N);
  Xsat.U = rd_sat.RB_local_ext.U * X(sat_N+1:sat_N+vel_N);
  msatsat = size(rd_sat.RB_local_ext.S,1);
  Xsat.Uoff = msatsat;
  [sat_ret, dummy, sat_Jret] = descr.L_saturation.apply(descr, rd_sat, Xsat, rd_sat.TM_local);
  sat_ret = rd_sat.CE * sat_ret;

  if descr.linear_velocity
    vel_Jret = lincomb_sequence(rd_vel, descr);
    vel_ret  = vel_Jret * X(sat_N+vel_N+1:end) - X(sat_N+1:sat_N+vel_N);
  else

    % hier koennte man die empirische interpolation auch lediglich fuer die
    % totale mobilitaet machen, statt fuer den ganzen velocity operator.
    Xvel.S = rd_vel.RB_local_ext.S * X(1:sat_N);
    %  Xvel.U = rd_vel.RB_local_ext.U * X(sat_N+1:sat_N+vel_N);
    Xvel.P = rd_vel.RB_local_ext.P * X(sat_N+vel_N+1:end);
    mvelsat = size(rd_vel.RB_local_ext.S,1);
    Xvel.Poff = mvelsat;
    Xvel.n    = vel_N;
    %  Xvel.Uoff = mvelsat;
    [vel_ret, dummy, vel_Jret] = descr.L_velocity.apply(descr, rd_vel, Xvel, rd_vel.TM_local);
    vel_ret = rd_vel.CE * vel_ret - X(sat_N+1:sat_N+vel_N);

%    dd = evalin('base', 'dd');
%    dd_sat = get_by_description(dd.datatree.rb, 'saturation');
%    dd_vel = get_by_description(dd.datatree.rb, 'velocity');
%    dd_prs = get_by_description(dd.datatree.rb, 'pressure');
%    dd_Lvel = get_by_description(dd.datatree.ei, 'L_velocity');
%    dd_Xvel.S = dd_sat.RB * X(1:sat_N);
%    dd_Xvel.P = dd_prs.RB * X(sat_N+vel_N+1:end);
%    [dd_ret, dummy, dd_Jret] = descr.L_velocity.apply(descr, dd.model_data, dd_Xvel);
%    vel_ret_recon_sigma = dd_Lvel.BM \ vel_ret;
%    vel_ret_recon = dd_Lvel.QM * vel_ret_recon_sigma;
%    debug_diff = vel_ret_recon - dd_ret;
%    disp(['linfty-norm: ', num2str(max(abs(debug_diff))), ...
%      '    l2-norm: ', num2str(sqrt(debug_diff' * dd.model_data.diamondW * debug_diff))]);
%%    disp(['TM-linfty-norm: ', num2str(max(abs(dd_ret(dd_Lvel.TM) - vel_ret)))]);
%    vel_ret = dd_vel.RB' * dd.model_data.diamondW * vel_ret_recon - X(sat_N+1:sat_N+vel_N);

  end


  [div_Jret, div_bb] = lincomb_sequence(rd_div, descr);
  prs_Jret           = lincomb_sequence(rd_prs, descr);
  div_ret  = div_Jret * (X(sat_N + 1:sat_N + vel_N)) + div_bb;
  prs_ret = prs_Jret * (X(sat_N+vel_N+1:end));

  sat_ret = (X(1:sat_N) - Sold)./descr.dt + sat_ret;
  ret = [sat_ret; vel_ret; div_ret; prs_ret];

  if nargout > 1

    sat_Jret = rd_sat.CE * ...
           [sat_Jret(:,1:msatsat) * rd_sat.RB_local_ext.S, ...
            sat_Jret(:,msatsat+1:end) * rd_sat.RB_local_ext.U, ...
            zeros(size(sat_Jret,1), prs_N)];
    ntime = size(sat_Jret, 1);
    Jtime = sparse(1:ntime, 1:ntime, 1/descr.dt * ones(1, ntime),...
                   ntime, size(sat_Jret,2));

    if descr.linear_velocity
      vel_Jret = [zeros(size(vel_Jret,1), sat_N+vel_N), vel_Jret];
    else
      vel_Jret = rd_vel.CE * ...
                [ vel_Jret(:,1:mvelsat) * rd_vel.RB_local_ext.S, ...
                  zeros(size(vel_Jret,1), vel_N), ...
                  vel_Jret(:,mvelsat+1:end) * rd_vel.RB_local_ext.P ];
%       vel_Jret = rd_vel.CE * ...
%           [ dd_Jret(dd_Lvel.TM,1:400) * dd_sat.RB, ...
%             zeros(size(vel_Jret, 1), vel_N), ...
%             dd_Jret(dd_Lvel.TM,401:end) * dd_prs.RB ];
    end

    sp_minus_U = sparse(1:vel_N, sat_N+1:sat_N+vel_N, -1, vel_N, size(vel_Jret,2));
    J = [ Jtime + sat_Jret;...
          vel_Jret + sp_minus_U;...
          zeros(size(div_Jret,1),sat_N), div_Jret, zeros(size(div_Jret,1),prs_N);
          zeros(1, sat_N+vel_N), prs_Jret ];
  end
end

function Xproj = box_projection(X, Slength)

  Xproj = X;

end