lin_evol_opt_detailed_der_simulation_old.m

function sim_data = lin_evol_opt_detailed_der_simulation_old(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

disp('OLD VERSION of derivative detailed simulation !!!!!!!!!!!!!!!!!!!!!!!')

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:size(U0_der_coeff,1)
   Ut_der(:,i) = lincomb_sequence(U0_comp, U0_der_coeff(i,:)'); %%% FRAGE: PASST DAS MIT Ut_der? Hat 3 dim statt 2?!
end
%%%%%%%%%%%%%%%%%%%%%% HIER SPARSE MATRITZEN VERWENDEN!!!???
% original:
U_der = zeros(length(Ut), length(save_time_indices),nr_indx_mu_i);

%new:
%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: HIER SPARSE MATRITZEN VERWENDEN!!!

end

% preallocate U
%%%%% DEBUG D, 06.12.10
%original:
U = zeros(length(Ut),length(save_time_indices));
%new:
%U = sparse(length(Ut),length(save_time_indices));
%%%%% END DEBUG D, 06.12.10

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
      %%%%%%% DEBUG E, 06.12.10
      %original:
      U_der(:,t_column,i)=Ut_der(:,i);
      %new:
      %U_der{i}(:,t_column) = Ut_der(:,i); %%% FRAGE: PASST DAS MIT Ut_der? Hat 3 dim statt 2?!      
      %%%%%%% END DEBUG E, 06.12.10
  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);
    %debug 28.10.10
   if ~isempty(L_I_comp) %only do the linear combination if the coefficients sLi are not empty
        %%%%%%%%%%%%%%%%%%%%%%% HIER AUCH NOCH ÄNDERN FALLS ES
        %%%%%%%%%%%%%%%%%%%%%%% LÄUFT!!!!!!!!!! 07.12.10
        %old:
        L_I_derived=zeros(size(L_I_comp{1},1),size(L_I_comp{1},2),nr_indx_mu_i);
        %new:
        %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: HIER AUCH NOCH ÄNDERN FALLS ES
        %%%%%%%%%%%%%%%%%%%%%%% LÄUFT!!!!!!!!!! 07.12.10
   end
    
    % DEBUG 06.12.10
    %%% original:
     L_E_derived=zeros(size(L_E_comp{1},1),size(L_E_comp{1},2),nr_indx_mu_i);
     b_derived=zeros(length(b_comp{1}),nr_indx_mu_i);
    %%% end original
    %%%%%%new try:
    %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
    
    %%%%%
    %End DEBUG 06.12.10
    
    for i=1:nr_indx_mu_i
        if ~isequal(sLi,[]) %only do the linear combination if the coefficients sLi are not empty
            %%%%%%%%%%% DEBUG G, 07.12.10
            %old:
            L_I_derived(:,:,i) = lincomb_sequence(L_I_comp,sLi(indx_mu_i(i),:));
            %new:
            %L_I_derived{i} = lincomb_sequence(L_I_comp,sLi(indx_mu_i(i)));
            %%%%%%%%%%% End: DEBUG G, 07.12.10
        else
            %%%%%%% DEBUG H, 07.12.10
            %old:
            L_I_derived(1,1,i); % ?????????????? ÄNDERN IN L_I_derived{1}=0 ??
            %new:
            %L_I_derived{i}=0;
            %%%%%%% END DEBUG H, 07.12.10
        end
%%%%%% DEBUG A, 06.12.10
%%%%%% original:        
         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),:));
%%%%%% new
%        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 DEBUG A, 06.12.10
    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
    %%%% DEBUG C, 06.12.10
    %%%% original:
     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
    %%%% new:
%    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
    %%%% end DEBUG C, 06.12.10
    
    
    %%%%% DEBUG F, 06.12.10
    %original:
    if save_time_index(t_ind)
        U_der(:,t_column-1,i) = Ut_der(:,i);
    end
    %new:
    %if save_time_index(t_ind)
    %    U_der{i}(:,t_column-1) = Ut_der(:,i);
    %end;
    
    %%%%% END DEBUG F, 06.12.10
    
  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
          %%%%% DEBUG G, 06.12.10
          % original:
          sim_data.y_der(:,i) = -(v{1}')*U_der(:,:,i);
          %new:
          %sim_data.y_der(:,i) = (-1).*(v{1}')*U_der{i};
          %%%%% END DEBUG G, 06.12.10
      end
  end
end;

%| \docupdate