nonlin_evol_detailed_rb_proj_simulation.m

function [sim_data, fl] = nonlin_evol_detailed_rb_proj_simulation(model,detailed_data, refU)
%function sim_data = nonlin_evol_detailed_rb_proj_simulation(model,detailed_data)
%
% required fields of model:
% T             : final time
% nt            : number of time-intervals until T, i.e. nt+1
%                 solution slices are computed
% init_values_algorithm : name of function for computing the
%                 initvalues-DOF with arguments (model, model_data)
%                 example: fv_init_values
% data_const_in_time : if this optional field is 1, the time
%                 evolution is performed with constant operators,
%                 i.e. only the initial-time-operators are computed
%                 and used throughout the time simulation.
% L_E_local_ptr :  function pointer to the local space-discretization
%                  operator evaluation with syntax
%
%                   INC_local = L_local(U_local_ext, ind_local,
%                                      grid_local_ext)
%                   where *_local indicates results on a set of
%                   elements, *_local_ext includes information including
%                   their neighbours, i.e. the grid must be known
%                   also for the neighbours, the values of the
%                   previous timstep must be known on the
%                   neighbours. The subset of elements, which are
%                   to be computed are given in ind_local and the
%                   values produced are INC_local.
%                   A time-step can then be realized by
%                        NU_local = U_local_ext(ind_local) - dt * INC_local
%                   example: fv_conv_explicit_space

% Martin Drohmann 14.04.2007 based on detailed_nonlin_evol_simulation.m by
% Bernard Haasdonk

if nargin ~= 3
  error('wrong number of parameters!');
end;

if model.verbose >= 9
  disp('entered detailed simulation ');
end;

% for a detailed simulation we do not need the affine parameter
% decomposition
model.decomp_mode = 0;

% initial values by midpoint evaluation
U0 = model.init_values_algorithm(model,detailed_data);

digits(60);
%digits(60);
%U0 = vpa(U0);

%Nm = model.Nm;
%Np = model.Np;
%N  = model.Np + model.Nm;
%RBp = detailed_data.RBp;
%RBm = detailed_data.RBm;
N = model.N;
%RB = [detailed_data.RBp,detailed_data.RBm];
RB = detailed_data.RB;
%RB = RBp + RBm;
disp(size(RB));
W  = model.get_inner_product_matrix(detailed_data);

%U = vpa(zeros(length(U0(:)),model.nt+1));
U = zeros(length(U0(:)),model.nt+1);
a = zeros(N, model.nt+1);

%a(1:Np,1) = RBp' * W * U0(:);
%U0min = U0(:) - RBp * a(1:Np,1);
%a(Np+1:N,1) = RBm' * W * U0min;

a = RB' * W * U0(:);

%a = RBp' * W * U0(:) + RBm' * W * U0(:);

U(:,1) = RB * a(:,1);

%disp(norm(U(:,1)-U0(:)));

if ~isfield(model,'data_const_in_time')
  model.data_const_in_time = 0;
end,

fl = cell(model.nt, 1);
% loop over fv-steps
if model.fv_impl_diff_weight ~= 0
  % diffusive part is computed implicitly.
  clear_all_caches;
  [L_diff_I, const_diff_I] = model.implicit_operators_algorithm(model, detailed_data, []);
  clear_all_caches;
end

model.dt = model.T / model.nt;

for t = 1:model.nt
  model.t = (t-1)*model.dt;
  model.tstep = t;
  if model.verbose > 19
    disp(['entered time-loop step ',num2str(t), ' for non linear evolution equation']);
  elseif model.verbose > 9
    fprintf('.');
  end;

  [inc, fluxes] = model.L_E_local_ptr(model, detailed_data, U(:,t), []);
  [refinc, dummy] = model.L_E_local_ptr(model, detailed_data, refU(:,t), []);

  if model.debug == 1
    fl{t} = fluxes;
  end
  if all([model.fv_impl_conv_weight, model.fv_impl_diff_weight, model.fv_impl_react_weight] == [0, 0, 0])
    % purely explicit scheme
%    a(:,t+1) = a(:,t) - model.dt * ainc;
    a(:,t+1) = a(:,t) - model.dt * inc;
  else % scheme with implicit parts

    inc    = inc + const_diff_I;
%    inc = double(single(inc));
    refinc = refinc + const_diff_I;
%    refres(t) = (inc-refinc)' * W * (inc-refinc);

    ainc = RB' * W * inc;
%    ainc = RBp' * W * inc + RBm' * W * inc;
%    ainc = zeros(N,1);

%    ainc(1:Np) = RBp' * W * inc;
%    incminus = inc - RBp * ainc(1:Np);
%    ainc(Np+1:N) = RBm' * W * incminus;

    L_I = L_diff_I;
    L_I = L_I * model.dt + speye(size(L_I));
%    L_I = sparse(double(single(full(L_I))));
    aL_I = RB' * W * L_I * RB;
%    aL_Ip = RBp' * W * L_I * RBp + RBm' * W * L_I * RBm;
%    aL_Im = RBm' * W * L_I * RBp + RBp' * W * L_I * RBm;
%    aL_I = [aL_Ip, aL_Im; aL_Im, aL_Ip];
    %aL_I = RBp' * W * L_I * RBp + RBp' * W * L_I * RBm + RBm' * W * L_I * RBp + RBm' * W * L_I * RBm;
%    L_I = speye(size(L_diff_I));

    arhs = a(:,t) - model.dt * ainc; % + model.dt * L_diff_I * U(:,t);
%    rhs = U(:,t) - model.dt * RB * ainc; % + model.dt * L_diff_I * U(:,t);
%    rhsdiff = U(:,t) - model.dt * inc - refU(:,t) + model.dt * refinc;
%    rhsres(t) = rhsdiff' * W * rhsdiff;
    if size(RB,2) == 20
%      keyboard;
    end
%    U(:,t+1) = rhs;
%    U(:,t+1) = L_I \ rhs;
%    atemp = aL_I \ [arhs; zeros(size(arhs))];
%    a(:,t+1) = atemp(1:length(arhs)) + atemp(length(arhs)+1:end);
    a(:,t+1) = aL_I \ arhs;
%    [a(:,t+1), flag] = bicgstab(aL_I, arhs, 100*eps, 1000);
%    a(:,t+1) = aL_I \ arhs;
%    LIresi = RB * a(:,t+1) - (RB*RB')*W*U(:,t+1);
%     a(:,t+1) = aL_I \ arhs;
%    residuum = aL_I * a(:,t+1) - arhs;
%    residuum = residuum' * residuum;
%    disp(residuum);
%    U(:,t+1) = L_I \ rhs;
%    disp(condest(L_I));
    U(:,t+1) = RB * a(:,t+1);
%    residuum = L_I * U(:,t+1) - U(:,t) + model.dt * inc;
%    res(t) = sqrt(residuum' * W * residuum);
%    LIres(t) = sqrt(LIresi' * W * LIresi);
%    Urefres(t+1) = sqrt((U(:,t+1)-refU(:,t+1))' * W * (U(:,t+1)-refU(:,t+1)));
%    disp(['res: ', num2str(res(t))]);
%    disp(['inc: ', num2str(refres(t))]);
%    disp(['rhs: ', num2str(rhsres(t))]);
%    disp(['LI: ', num2str(LIres(t))]);
%    disp(['err: ', num2str(Urefres(t))]);
%    pause
%    U(U(:,t+1)<0,t+1);
%    a(:,t+1) = RB' * W * U(:,t+1);
%    U(:,t+1) = RB * a(:,t+1);
    if flag>0
      disp(['error in system solution, solver return flag = ', ...
        num2str(flag)]);
      keyboard;
    end;
  end
end

%    disp(['res: max=', num2str(max(res)), ' sum=', num2str(sum(res))]);
%    disp(['inc: max=', num2str(max(refres)), ' sum=', num2str(sum(refres))]);
%    disp(['rhs: max=', num2str(max(rhsres)), ' sum=', num2str(sum(rhsres))]);
%    disp(['LI: max=', num2str(max(LIres)), ' sum=', num2str(sum(LIres))]);

sim_data.U = U;
sim_data.a = a;

get_mu(model)

if model.verbose > 9
  fprintf('\n');
end;

%| \docupdate