lin_evol_opt_detailed_der_simulation.m

function sim_data = lin_evol_opt_detailed_der_simulation(model,model_data)
%function sim_data = lin_evol_opt_detailed_der_simulation(model,model_data)
%
% Function performing a detailed simulation and calculating the derivative
% solution with respect to all parameters.
% 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 (grid, params)
%                 example: init_values_cog
% operators_algorithm: name of function for computing the
%                 L_E,L_I,b-operators with arguments (grid, params)
%                 example operators_conv_diff
% 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.
% compute_output_functional: flag indicating, whether output
%                 functional is to be computed
%
% compute_derivative_info: flag if derivative information shall be
%                 calculated or not
%
% return fields of sim_data:
%    U : sequence of DOF vectors
%    y : sequence of output functional values (if activated)
% 
% Markus Dihlmann 15.10.10

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

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

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


model.decomp_mode = 1; % == components;
model.t = 0;


if ~isfield(model,'optimization')
    disp('Model has no field optimization.')
    model.optimization =[];
end


if isfield(model.optimization,'params_to_optimize')
    indx_mu_i = find(model.optimization.params_to_optimize); %indexes i of \mu_i for which derivative information is needed
else
    indx_mu_i = 1:length(get_mu(model));
end
nr_indx_mu_i= length(indx_mu_i); %nr of parameters to be optimized
%nr_derivatives = size(model.optimization.params_to_optimize,2); %total nr of parameters in the model



if isfield(model,'save_time_indices')
  save_time_index = zeros(1,model.nt+1);
  save_time_index(model.save_time_indices(:)+1) = 1;
  save_time_indices = find(save_time_index==1)-1;
else
  save_time_index = ones(1,model.nt+1);
  save_time_indices = 1:(model.nt+1);
end;

% initial values by midpoint evaluation
U0_comp = model.init_values_algorithm(model,model_data);

model.decomp_mode=2; %==> coefficients
U0_coeff = model.init_values_algorithm(model, model_data);
if model.compute_derivative_info
    U0_der_coeff = model.rb_derivative_init_values_coefficients(model);
end

%assembling initial values true solution
Ut = lincomb_sequence(U0_comp,U0_coeff);

%assembling initial values derivative solution
if model.compute_derivative_info
Ut_der = zeros(length(Ut),nr_indx_mu_i);
for i=1:nr_indx_mu_i
   Ut_der(:,i) = lincomb_sequence(U0_comp, U0_der_coeff(indx_mu_i(i),:)'); 
end

U_der=cell(nr_indx_mu_i,1);
for i=1:nr_indx_mu_i
    U_der{i} = sparse(length(Ut), length(save_time_indices));
end


end


U = sparse(length(Ut),length(save_time_indices));


time_indices = zeros(1,length(save_time_indices));
time =  zeros(1,length(save_time_indices));
model.dt = model.T/model.nt;
%operators_required = 1;
if ~isfield(model,'data_const_in_time')
  model.data_const_in_time = 0;
end,

% get operator components
if model.affinely_decomposed
  old_decomp_mode = model.decomp_mode;
  model.decomp_mode = 1;
  [L_I_comp, L_E_comp, b_comp] = model.operators_ptr(model, model_data);
  model.decomp_mode = old_decomp_mode;
end;

t_column = 1; % next column index to be filled in output
t_ind = 0; % t_ind between 0 and nt
t = 0; % absolute time between 0 and T

% store init data
if save_time_index(t_ind+1)
  U(:,t_column) = Ut;
  if model.compute_derivative_info
  for i=1:nr_indx_mu_i
      U_der{i}(:,t_column) = Ut_der(:,i);
  end
  end
  time_indices(t_column) = t_ind;
  time(t_column) = t;
  t_column = t_column + 1;  
end;

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

%model.plot(U0,model_data.grid,params);

%allocate memory for derivation matrices:
L_E_derived = cell(nr_indx_mu_i,1);
L_I_derived = cell(nr_indx_mu_i,1);
b_derived = cell(nr_indx_mu_i,1);

% loop over fv-steps
for t_ind = 1:model.nt  
  t = (t_ind-1)*model.dt;
  if model.verbose >= 10
    disp(['entered time-loop step ',num2str(t_ind)]);
  end;
  if model.verbose == 9
    fprintf('.');
  end;
  % get matrices and bias-vector 
  
  if (t_ind==1)||(~model.data_const_in_time)%operators_required
    
    model.t = t;
    
    % new assembly in every iteration
    if ~model.affinely_decomposed
      [L_I, L_E, b] = model.operators_ptr(model, model_data);
    else
      % or simple assembly based on affine parameter decomposition
      % => turns out to be slower for few components
      
      % test affine decomposition
      if model.debug
        disp('test affine decomposition of operators')
        test_affine_decomp(model.operators_ptr,3,1,model,model_data);
        disp('OK')
      end;
      model.decomp_mode = 2;
      [L_I_coeff, L_E_coeff, b_coeff] = model.operators_ptr(model, ...
                                                  model_data);
      model.decomp_mode = old_decomp_mode;
      L_I = lincomb_sequence(L_I_comp, L_I_coeff);
      L_E = lincomb_sequence(L_E_comp, L_E_coeff);    
      b = lincomb_sequence(b_comp, b_coeff);
    end;
    
    %if model.data_const_in_time
    %  operators_required = 0;
    %end;
    
    if model.compute_derivative_info
    % Get Operators for derivative information:
    [sLi,sLe,sb] = model.rb_derivative_operators_coefficients_ptr(model);
    %sLe=-sLe;
    
   if ~isempty(L_I_comp) %only do the linear combination if the coefficients sLi are not empty
        for i=1:nr_indx_mu_i
            L_I_derived{i}=sparse(size(L_I_comp{1},1),size(L_I_comp{1},2));
        end
   end
    
    for i=1:nr_indx_mu_i
        L_E_derived{i} = sparse(size(L_E_comp{1},1),size(L_E_comp{1},2));
        b_derived{i} = sparse(length(b_comp{1}),1);
    end
    
    
    for i=1:nr_indx_mu_i
        if ~isequal(sLi,[]) %only do the linear combination if the coefficients sLi are not empty
            L_I_derived{i} = lincomb_sequence(L_I_comp,sLi(indx_mu_i(i)));
        else
            L_I_derived{i}=0;
        end
        L_E_derived{i} = lincomb_sequence(L_E_comp,sLe(indx_mu_i(i),:));
        b_derived{i} = lincomb_sequence(b_comp,sb(indx_mu_i(i),:));
    end
    end

  end;

  %rhs = L_E * Ut + b;
  rhs = (speye(size(L_E)) + model.dt*L_E) * Ut + model.dt*b;
  Ut_old=Ut;
  
  if isempty(find(L_I, 1))%isequal(L_I, speye(size(L_I)))
    Ut = rhs;
  else
    % solve linear system
    %    disp('check symmetry and choose solver accordingly!');
    %    keyboard;
    % nonsymmetric solvers:
    %  [U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
    %  [U(:,t+1), flag] = cgs(L_I,rhs,[],1000);
    %
    % symmetric solver, non pd:
    % omit symmlq:
    % see bug_symmlq.mat for a very strange bug: cannot solve identity system!
    % reported to matlab central, but no solution up to now.
    % [U(:,t+1), flag] = symmlq(L_I,rhs,[],1000);
    %
    %  [U(:,t+1), flag] = minres(L_I,rhs,[],1000);
    % symmetric solver, pd:
    %[U(:,t+1), flag] = pcg(L_I,rhs,[],1000);
    % bicgstab works also quite well:
    %[U(:,t+1), flag] = bicgstab(L_I,rhs,[],1000);
      Ut = (speye(size(L_I))+model.dt*L_I)\rhs;%Ut + model.dt*(L_I \ rhs);

    %    if flag>0
    %      disp(['error in system solution, solver return flag = ', ...
    %       num2str(flag)]);
    %      keyboard;
    %    end;

  end;
  
  
  if save_time_index(t_ind)
    U(:,t_column) = Ut;
    time_indices(t_column) = t_ind;
    time(t_column) = t;
    t_column = t_column + 1;  
  end;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %calculation derivative solution
  if model.compute_derivative_info
  for i=1:nr_indx_mu_i
    rhs_der = (speye(size(L_E)) + model.dt*L_E) * Ut_der(:,i) +...
        model.dt*L_E_derived{i} * Ut_old + model.dt*L_I_derived{i}*Ut + model.dt*b_derived{i};
    if isequal(L_I,zeros(size(L_I)))
        Ut_der(:,i) = rhs_der;
    else
        Ut_der(:,i) = (speye(size(L_I))+ model.dt*L_I)\rhs_der;
    end
    if save_time_index(t_ind)
        U_der{i}(:,t_column-1) = Ut_der(:,i);
    end;
    
  end
  end
  
end;


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


sim_data.U = U;
if model.compute_derivative_info
sim_data.U_der = U_der;
end

sim_data.time = time;
sim_data.time_indices = time_indices;

if model.compute_output_functional
  % get operators
  model.decomp_mode=1;
  v = model.operators_output(model,model_data);
  sim_data.y = (v{1}') * U;
  sim_data.y_der=[];
  if model.compute_derivative_info
      for i = 1:nr_indx_mu_i
          sim_data.y_der(i,:) = (v{1}')*U_der{i};
      end
  end
end;