detailed_simulation_galerkin.m

function sim_data = detailed_simulation_galerkin(dmodel, detailed_data)
% function sim_data = detailed_simulation(dmodel, detailed_data)
% executes a detailed simulation for a given parameter
%
% required fields of model:
% T             : final time
% nt            : number of time-steps
%                 `0 = t_0 \leq t_1 \leq \cdots \leq t^K+1 = T`, i.e. `K+1`.
% init_values_algorithm : function pointer for computing the initvalues-DOF
%                 with arguments '(model, model_data)'
%                  - example: fv_init_values()
%                  .
% L_E_local_ptr :  function pointer to the local space-discretization
%                  operator evaluation with syntax
%                  @code
%                   INC_local = L_local(U_local_ext, ind_local, grid_local_ext)
%                  @endcode
%                  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
%                  @code
%                      NU_local = U_local_ext(ind_local) - dt * INC_local
%                  @endcode
%                  - example: fv_conv_explicit_space()
%                  .
%
% optional fields of model:
% 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.
%
% Return values:
%   sim_data: structure holding the `H`-dimensional simulation data.
%   fl:       debug output holding the flux evaluations returned by the
%             'L_E_local_ptr' operator.
%
% Generated fields of sim_data:
%  U:        matrix of size `H \times K+1` holding for each time step the solution
%            dof vector.
%  Unewton:  matrix of size `H \times K+1+\sum_{k=0}^K n_{\text{Newton}}(k)`
%            holding for each time step and each Newton iteration the
%            intermediate solution dof vectors.
%


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

dd = detailed_data;
model_data = dd.model_data;

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

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

if isempty(model_data)
  model_data = gen_model_data(dmodel);
end;


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

model = dmodel.descr;

% initial values by midpoint evaluation
S0 = model.init_values.evaluate(model,model_data);

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');
sat_N = size(dd_sat.RB, 2);
vel_N = size(dd_vel.RB, 2);
prs_N = size(dd_prs.RB, 2);

S = zeros(length(S0(:)),model.nt+1);
U = zeros(model_data.gn_edges, model.nt+1);
P = zeros(length(S0(:)),model.nt+1);
sat_a   = zeros(sat_N, model.nt+1);   % matrix of saturation RB-coefficients
vel_a   = zeros(vel_N, model.nt+1);   % matrix of velocity   RB-coefficients
prs_a   = zeros(prs_N, model.nt+1);   % matrix of pressure   RB-coefficients
exflags = zeros(1, model.nt+1);
outtext = cell(1, model.nt+1);
Slength   = length(S0(:));
Ulength   = model_data.gn_edges;

%totlength = 2*Slength+length(grid.VI);
S(:,1) = S0(:);
W = model_data.W;
dW = model_data.diamondWinv;
dd_satl2 = dd_sat.RB';
dd_vell2 = dd_vel.RB';
dd_prsl2 = dd_prs.RB';

sat_a(:,1) = dd_satl2 * W * S(:,1);
%sim_data = evalin('base', 'sim_data');
%prs_a(:,1) = dd_prsl2 * W * sim_data.P(:,2);
%vel_a(:,1) = dd_vell2 * dW * sim_data.U(:,2);


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

%  sat_a(:,t+1) = dd_satl2 * W * sim_data.S(:,t+1);
%  vel_a(:,t+1) = dd_vell2 * dW * sim_data.U(:,t+1);
%  prs_a(:,t+1) = dd_prsl2 * W * sim_data.P(:,t+1);
  %  V = zeros(totlength, model.newton_steps+1);
  V_last = [sat_a(:,t); vel_a(:,t); prs_a(:,t)];

  model.t = model.t + model.dt;

%    Srhs   = model.S_local_ptr(model, model_data, V(1:Slength+Ulength,n), []);
%    Urhs   = model.U_local_ptr(model, model_data, V(Slength+1:end,n), []);
%    Prhs   = model.P_local_ptr(model, model_data, V(Slength+Ulength+1:end,n), []);
%    gradSM = model.S_op_matrix(model, model_data, V(Slength+1:Slength+Ulength,n), []);
%    gradUM = model.U_op_matrix(model, model_data, V(1:Slength,n), []);
%    gradPM = model.P_op_matrix(model, model_data, v(Slength+Ulength+1:end,n), []);

  fsolve_fun = @(X) sat_fun(model, model_data, sat_a(:,t), X, dd_sat, dd_vel, dd_prs, dd_satl2, dd_vell2, dd_prsl2);
  %[(X(1:Slength)-S(:,t))./model.dt;zeros(Ulength+Slength+1,1)] ...
  %  + fv_two_phase_space(model, model_data, X, []);

  %LB = [zeros(Slength,1);-Inf(Ulength,1);-Inf(Slength,1)];
  %UB = [ones(Slength,1)-100*eps;Inf(Ulength,1);Inf(Slength,1)];
  %[V_new, resnorm, FV, exitflag, output] = lsqnonlin(fsolve_fun, V_last, LB, UB,...
  %  optimset('Display', 'iter', ...
  %           'Algorithm', 'trust-region-reflective',...
  %           'Jacobian', 'on', ...
  %           'FinDiffType', 'central', ...
  %           'Diagnostics', 'on', ...
  %           'DerivativeCheck', 'off'));
  %  [V_new, FV, exitflag, output, jacobian] = fsolve(fsolve_fun, V_last, ...
  %     optimset('LargeScale', 'on', ... %'Algorithm', 'levenberg-marquardt', ...
  %              'Display', 'iter', ...
  %              'Jacobian', 'on', ...
  %              'Diagnostics', 'off', ...
  %              'DerivativeCheck', 'off'...
  %             ));

  gplmmopts.Px = @(X) box_projection(X, Slength);
%  gplmmopts.HessFun = @(params, jac, dummy, arg) jac' * (jac * arg);
  gplmmopts.Slength = Slength;
  [V_new, resnorm, FV, exitflag, output] = gplmm(fsolve_fun, V_last, gplmmopts);
%  [V_new, resnorm, FV, exitflag, output] = newton_raphson(fsolve_fun, V_last, gplmmopts);

%  sat_a(:,t) = dd_sat.RB' * model_data.W * sim_data.S(:,t);
%  vel_a(:,t) = dd_vel.RB' * model_data.diamondWinv * sim_data.U(:,t);
%  prs_a(:,t) = dd_prs.RB' * model_data.W * sim_data.P(:,t);
  V_new = real(V_new);
  disp(['t= ', num2str(model.t)]);

  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);
  S(:,t+1) = dd_sat.RB * sat_a(:,t+1);
  %  Stemp(:,2*t+1:2*t+2) = tempsols(1:Slength,:);
  U(:,t+1) = dd_vel.RB * vel_a(:,t+1);
  %  Utemp(:,2*t+1:2*t+2) = tempsols(Slength+1:Slength+Ulength,:);
  P(:,t+1) = dd_prs.RB * prs_a(:,t+1);
  %  Ptemp(:,2*t+1:2*t+2) = tempsols(Slength+Ulength+1:end,:);

end

sim_data.sat_a = sat_a;
sim_data.vel_a = vel_a;
sim_data.prs_a = prs_a;
sim_data.S = S;
%sim_data.Stemp = Stemp;
sim_data.U = U;
%sim_data.Utemp = Utemp;
sim_data.P = P;
%sim_data.Ptemp = Ptemp;
sim_data.exit_flags = exflags;
sim_data.outputs    = outtext;

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

function [ret, J] = sat_fun(descr, model_data, Sold, X, dd_sat, dd_vel, dd_prs, dd_satl2, dd_vell2, dd_prsl2)

  model_data = dd_sat.model_data;
  W = model_data.W;
  dW = model_data.diamondWinv;

  sat_N = size(dd_sat.RB, 2);
  vel_N = size(dd_vel.RB, 2);
  prs_N = size(dd_prs.RB, 2);
  Uoff = size(dd_sat.RB, 1);
  Poff = Uoff + size(dd_vel.RB, 1);

  %  Stemp(Stemp < 0) = 0+eps;
  %  Stemp(Stemp > 1) = 1-eps;
  arglength = Poff + size(dd_prs.RB, 1);
  L_sat_arg.S = dd_sat.RB * X(1:sat_N);
  L_sat_arg.U = dd_vel.RB * X(sat_N+1:sat_N+vel_N);
%  L_sat_arg.S = X(1:Uoff);
%  L_sat_arg.U = X(Uoff+1:Poff);
  L_sat_arg.Uoff = Uoff;
  L_sat_arg.m    = arglength;
  [sat_ret, dummy, sat_Jret] = descr.L_saturation.apply(descr, model_data, L_sat_arg);
  sat_ret = (X(1:sat_N) - Sold)./descr.dt + ( dd_satl2 * W * sat_ret );
%  sat_ret = ddsatl2 * W * ( (X(1:Uoff) - Sold)./descr.dt + sat_ret );

  L_vel_arg.S = dd_sat.RB * X(1:sat_N);
  L_vel_arg.P = dd_prs.RB * X(sat_N+vel_N+1:end);
%  L_vel_arg.S = X(1:Uoff);
%  L_vel_arg.P = X(Poff+1:end);
  L_vel_arg.Poff = Poff;
  L_vel_arg.n    = Poff-Uoff;
  L_vel_arg.m    = arglength;
  [vel_ret, dummy, vel_Jret] = descr.L_velocity.apply(descr, model_data, L_vel_arg);

  L_div_arg.U = dd_vel.RB * X(sat_N+1:sat_N+vel_N);
%  L_div_arg.U = X(Uoff+1:Poff);
  [div_Jret, div_bb] = descr.L_divergence.full_matrix(descr, model_data);

  % real_NBI_ind = model_data.grid.NBI > 0 & model_data.grid.NBI < repmat((1:model_data.grid.nelements)', 1, 4);
  % prs_mat = repmat(dd_prs.RB, 1, 4);
  % prs_on_edge1 = zeros(1, model_data.gn_edges);
  % prs_on_edge2 = zeros(1, model_data.gn_edges);
  % prs_on_edge1(model_data.gEI(real_NBI_ind)) = prs_mat(real_NBI_ind);
  % prs_on_edge2(model_data.gEI(real_NBI_ind)) = dd_prs.RB(model_data.grid.NBI(real_NB_ind));
  div_ret = div_Jret *  L_div_arg.U + div_bb;
  %  div_ret = div_Jret *  L_div_arg.U + div_bb;

  [prs_ret, dummy, prs_Jret] = descr.L_pressure.apply(descr, model_data, L_vel_arg);

  %  ret = [sat_ret; dd_vel.RB' * dW * (vel_ret) - X(sat_N+1:sat_N+vel_N); div_ret];
  ret = [sat_ret; dd_vell2 * dW * (vel_ret) - X(sat_N+1:sat_N+vel_N); div_ret; prs_ret];
  %  ret = [sat_ret; dd_vell2 * dW * (vel_ret - X(Uoff+1:Poff)); div_ret; prs_ret];

  assert(all(~isnan(ret)));
  assert(all(~isinf(ret)));

%  [ret2, Jret2] = fv_two_phase_space(descr, model_data, X);
%  retlength = length(ret2);

  if nargout > 1

    n = vel_N;
    m = sat_N + vel_N + prs_N;
    sp_minus_U = sparse(1:n, sat_N+1:sat_N+vel_N, -1, n, m);
    sat_Jret = dd_satl2 * W * ([sat_Jret(:,1:Uoff) * dd_sat.RB, sat_Jret(:,Uoff+1:Poff) * dd_vel.RB, zeros(Uoff, prs_N)]);
    vel_Jret = dd_vell2 * dW * ([vel_Jret(:,1:Uoff) * dd_sat.RB, zeros(Poff-Uoff, vel_N), vel_Jret(:,Poff+1:end) * dd_prs.RB]);
    div_Jret = ([zeros(Uoff, sat_N), div_Jret * dd_vel.RB, zeros(Uoff, prs_N)]);
    J = [ sparse(1:sat_N, 1:sat_N, 1/descr.dt * ones(1,sat_N), sat_N, m) ...
          + sat_Jret; ...
          vel_Jret + sp_minus_U; ...
          div_Jret; ...
          zeros(1, sat_N+vel_N), prs_Jret(:,Poff+1:end) * dd_prs.RB ];
%   n = Poff-Uoff;
%    m = arglength;
%    sp_minus_U = [zeros(vel_N, Uoff), dd_vell2 * dW * sparse(1:n, 1:n, -1, n, n), zeros(vel_N, arglength - Poff)];
%    sat_Jret = dd_satl2 * W * ([sat_Jret(:,1:Uoff), sat_Jret(:,Uoff+1:Poff), zeros(Uoff, arglength-Poff)]);
%    vel_Jret = dd_vell2 * dW * ([vel_Jret(:,1:Uoff), zeros(Poff-Uoff, vel_N), vel_Jret(:,Poff+1:end)]);
%    div_Jret = ([zeros(Uoff, Uoff), div_Jret, zeros(Uoff, Uoff)]);
%    J = [ sparse(1:sat_N, 1:sat_N, 1/descr.dt * ones(1,sat_N), sat_N, m) ...
%          + sat_Jret; ...
%          vel_Jret + sp_minus_U; ...
%          div_Jret; ...
%          zeros(1, sat_N+vel_N), prs_Jret(:,Poff+1:end) * dd_prs.RB ];

%    J2 = sparse(1:Uoff, 1:Uoff, ...
%               1/descr.dt * ones(1,Uoff), ...
%               retlength, arglength) ...
%          + Jret2;
%    keyboard;
  end

function Xproj = box_projection(X, Slength)

  Xproj = X;
%  SX = X(1:Slength);
%  SX(SX<0) = eps;
%  upper_bound = 1.0;%0.6359878341047967;
%  SX(SX>upper_bound) = upper_bound-eps;
%  Xproj = [SX; X(Slength+1:end)];

%  dd = evalin('base', 'dd');
%  sat_dd = get_by_description(dd.datatree.rb, 'saturation');
%  vel_dd = get_by_description(dd.datatree.rb, 'velocity');
%  prs_dd = get_by_description(dd.datatree.rb, 'pressure');
%  aaa = sat_dd.RB' * dd.model_data.W * SX;
%  bbb = vel_dd.RB' * dd.model_data.diamondW * X(Slength+1:Slength+size(vel_dd.RB,1));
%  ccc = prs_dd.RB' * dd.model_data.W * X(Slength+size(vel_dd.RB,1)+1:end);
%  Xproj = [sat_dd.RB * aaa; vel_dd.RB * bbb; prs_dd.RB * ccc];