detailed_simulation.m

function sim_data = detailed_simulation(dmodel, model_data)
% function sim_data = detailed_simulation(dmodel, model_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

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);

S = zeros(length(S0(:)),model.nt+1);
%Stemp = zeros(length(S0(:)),2*model.nt+2);
U = zeros(model_data.gn_edges, model.nt+1);
%Utemp = zeros(model_data.gn_edges, 2*model.nt+2);
P = zeros(length(S0(:)),model.nt+1);
%Ptemp = zeros(length(S0(:)),2*model.nt+2);
Stemp = [];
Ptemp = [];
Utemp = [];
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(:);

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

  %  V = zeros(totlength, model.newton_steps+1);
  V_last = [S(:,t); U(:,t); P(:,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, S(:,t), X, Slength, Ulength+Slength);
  %[(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, tempsols] = gplmm(fsolve_fun, V_last, gplmmopts);

  V_new = real(V_new);
  disp(['t= ', num2str(model.t)]);

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

  %% TODO: This is bad! Isn't it?
  S(:,t+1) = min(V_new(1:Slength),1-eps);
  Stemp = [Stemp, tempsols(1:Slength,:)];
  U(:,t+1) = V_new(Slength+1:Slength+Ulength);
  Utemp = [Utemp, tempsols(Slength+1:Slength+Ulength,:)];
  P(:,t+1) = V_new(Slength+Ulength+1:end);
  Ptemp = [Ptemp, tempsols(Slength+Ulength+1:end,:)];

end

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, Uoff, Poff)

  Stemp = X(1:Uoff);
  %  Stemp(Stemp < 0) = 0+eps;
  %  Stemp(Stemp > 1) = 1-eps;
  arglength = length(X);
  L_sat_arg.S = Stemp;
  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 = (Stemp - Sold)./descr.dt + sat_ret;

  L_vel_arg.S = Stemp;
  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 = X(Uoff+1:Poff);
  L_div_arg.Uoff = Uoff;
  L_div_arg.m    = arglength;
  [div_ret, dummy, div_Jret] = descr.L_divergence.apply(descr, model_data, L_div_arg);

%  prs_ret = X(Poff+1);
  [prs_ret, dummy, prs_Jret] = descr.L_pressure.apply(descr, model_data, L_vel_arg);

  ret = [sat_ret; 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 = size(vel_Jret,1);
    sp_minus_U = sparse(1:n, Uoff+1:Poff, -1, n,arglength);
    J = [ sparse(1:Uoff, 1:Uoff, 1/descr.dt * ones(1,Uoff), Uoff, arglength) ...
          + sat_Jret; ...
          vel_Jret + sp_minus_U; ...
          div_Jret; ...
          prs_Jret ];
%          zeros(1,Poff), 1, zeros(1,Uoff-1) ];

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

function Xproj = box_projection(X, Slength)

  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];