detailed_simulation_impes.m

function sim_data = detailed_simulation_impes(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`.
%
% 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_algorithm(model,model_data);

S = zeros(length(S0(:)),model.nt+1);
U = 0.5*ones(model_data.gn_edges, model.nt+1);
P = 0.5*ones(length(S0(:)),model.nt+1);
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);
  UP_last = [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_imp(model, model_data, S(:, t), X, 0, Ulength);

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

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

  U(:,t+1) = UP_new(1:Ulength);
  P(:,t+1) = UP_new(Ulength+1:end);

  S(:,t+1) = S(:,t) - model.dt * sat_fun_es(model, model_data, U(:,t+1), S(:,t));

end

sim_data.S = S;
sim_data.U = U;
sim_data.P = P;
sim_data.exit_flags = exflags;
sim_data.outputs    = outtext;

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

function ret = sat_fun_es(descr, model_data, Uold, Sold)
  L_sat_arg.S = Sold;
  L_sat_arg.U = Uold;
  L_sat_arg.m = length(Sold);
  ret = descr.L_saturation.apply(descr, model_data, L_sat_arg);

function [ret, J] = sat_fun_imp(descr, model_data, Sold, X, Uoff, Poff)

  arglength = length(X);

  L_vel_arg.S = Sold;
  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, dummy, prs_Jret] = descr.L_pressure.apply(descr, model_data, L_vel_arg);

  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 = [ vel_Jret + sp_minus_U;
          div_Jret;
          prs_Jret ];

%    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;
  SX(SX>1) = 1-eps;
  Xproj = [SX; X(Slength+1:end)];