rb_simulation_impes_impl.m

function rb_sim_data = rb_simulation_impes_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);

if nargin == 4
  dd_vel = get_by_description(dd.datatree.rb, 'velocity');
  dd_prs = get_by_description(dd.datatree.rb, 'pressure');
  vel_a(:,1) = dd_vel.RB' * dd.model_data.diamondW * sim_data.U(:,2);
  prs_a(:,1) = dd_prs.RB' * dd.model_data.W * sim_data.P(:,2);
end

% 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;
  fprintf('.');

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

  fsolve_fun_imp = @(X) sat_fun_imp(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.HessFun = @(params, jac, dummy, arg) jac' * (jac * arg);
  [UP_new, resnorm, FV, exitflag, output] = gplmm(fsolve_fun_imp, UP_last, gplmmopts);

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

  vel_a(:, t+1) = UP_new(1:vel_N);
  prs_a(:, t+1) = UP_new(vel_N+1:end);

  sat_a(:, t+1) = sat_a(:,t) - descr.dt * sat_fun_es(descr, rd_sat, vel_a(:, t+1), sat_a(:,t));

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 = sat_fun_es(descr, rd_sat, vel_a, sat_a)
  Xsat.S = rd_sat.RB_local_ext.S * sat_a;
  Xsat.U = rd_sat.RB_local_ext.U * vel_a;

  msatsat = size(rd_sat.RB_local_ext.S,1);
  Xsat.Uoff = msatsat;
  sat_ret = descr.L_saturation.apply(descr, rd_sat, Xsat, rd_sat.TM_local);
  ret = rd_sat.CE * sat_ret;
end

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

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


  if descr.linear_velocity
    vel_Jret = lincomb_sequence(rd_vel, descr);
    vel_ret  = vel_Jret * X(vel_N+1:end) - X(1:vel_N);
  else
    Xvel.S = rd_vel.RB_local_ext.S * Sold;
    %  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(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);

    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 * Sold;
    dd_Xvel.P = dd_prs.RB * X(vel_N+1:end);
    dd_ret = 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(1:vel_N);

%    vel_ret = rd_vel.CE * vel_ret - X(1:vel_N);
  end


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

  ret = [vel_ret; div_ret; prs_ret];

  if nargout > 1

    if descr.linear_velocity
      vel_Jret = [zeros(size(vel_Jret,1), vel_N), vel_Jret];
    else
      vel_Jret = rd_vel.CE * ...
                [ zeros(size(vel_Jret,1), vel_N), ...
                  vel_Jret(:,mvelsat+1:end) * rd_vel.RB_local_ext.P ];
    end

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

function Xproj = box_projection(X, Slength)

  Xproj = X;

end