lin_evol_detailed_simulation_primal_dual.m

Go to the documentation of this file.

function sim_data = lin_evol_detailed_simulation_primal_dual(model, model_data)
%sim_data = lin_evol_detailed_simulation_primal_dual(model, model_data)
%
% copy of lin_evol_detailed_simulation but extended to a primal/dual
% setting.
%
% function which performs a detailed simulation for the given mu, which has
% to be set in model via model.set_mu. Can either perform a pimal or a dual
% detailed simulation controlled via the field model.want_dual (= 1 or 0).
%
% see lin_ds_detailed_simulation for a description on the time indices
%
% 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
%
% return fields of sim_data:
%    U : sequence of DOF vectors
%    y : sequence of output functional values (if activated)
%
% optional fields of model:
%   starting_time_step: starting time step for simulation (for example in
%                   t-partition). Default value is 0.
%   stopping_tim_step: stopping time step for simulation
%
% special Note: in the european_option_pricing model a functional named
% 'theta' includes a partial timederivative and therefore produces a
% slightly different dual problem. Therefore the 'theta'-case is sometimes
% treated differently. Don't care about this case if you don't use the
% european_option_pricing model.
%
% Bernard Haasdonk 27.8.2009
% 
% primal/dual formulation by Dominik Garmater 20.07 2012



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

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

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

model.decomp_mode = 0; % == complete;

model.dt = model.T/model.nt;

%Fixing values for t-partition
if ~isfield(model,'starting_time_step')
    model.t = 0;
    t_ind_start =0;
    t_ind_stop = model.nt;
else
    t_ind_start = model.starting_time_step;
    t_ind_stop = model.stopping_time_step;
    model.t = (t_ind_start)*model.dt;
end

% see lin_ds_detailed_simulation for a description on the time indices
if isfield(model,'save_time_indices')
  valid_save_time_indices1 = model.save_time_indices>=t_ind_start;
  valid_save_time_indices2 = model.save_time_indices<=t_ind_stop;
  valid_save_time_indices = valid_save_time_indices1.*valid_save_time_indices2;
  valid_save_time_indices = find(valid_save_time_indices);
  save_time_index = zeros(1,model.nt+1);
  save_time_index(model.save_time_indices(valid_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

time_indices = zeros(1,length(save_time_indices));
time =  zeros(1,length(save_time_indices));
t_column = 1; % next column index to be filled in output    
t_ind = t_ind_start; % t_ind between 0 and nt
if model.want_dual % dual case
    t = model.T-model.t; % backwards time between T and 0
else % primal case
    t = model.t; % absolute time between 0 and T
end

model.nt = t_ind_stop-t_ind_start+1;

% initial values and preallocation
if model.want_dual
    % in the dual case the riesz-representant of the functional is part of
    % the initail values (this will be modified in the later first timestep)
    [Ut, ~] = model.operators_output(model,model_data);
    U = zeros(length(Ut),length(save_time_indices));
else % primal case
    Ut = model.init_values_algorithm(model,model_data);
    U = zeros(length(Ut),length(save_time_indices));
end

operators_required = 1;

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

% store init data
if save_time_index(t_ind+1)
  time_indices(t_column) = t_ind;
  time(t_column) = t;      
  % in the dual case it is Ut = -L_I^(-1)*v_l.
  % this is performed later, when L_I is available
  if model.want_dual == 0 %primal
    U(:,t_column) = Ut;
    t_column = t_column + 1;
  end
end

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


% loop over fv-steps
for t_ind = (t_ind_start):(t_ind_stop-1)
  
  % set the time - relevant if your operators are dependant on time.
  if model.want_dual == 0
      model.t = (t_ind)*model.dt; % forward time in primal case
  else
      model.t = model.T - (t_ind)*model.dt; % the backwards time in dual case
  end
  
  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 operators_required
    
    model.t = t;
    
    % new assembly in every iteration
    if ~model.affinely_decomposed
      % model.decomp_mode is 0 here
      if model.want_dual %dual
          [~, ~, ~, L_I, L_E] = model.operators_ptr(model, model_data);
      else %primal
          [L_I, L_E, b] = model.operators_ptr(model, model_data);
      end
      
    else
      % or simple assembly based on affine parameter decomposition 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;
      if model.want_dual %dual
          [~, ~, ~, L_I_coeff, L_E_coeff] = ...
              model.operators_ptr(model, model_data);
          L_I = lincomb_sequence(L_I_comp, L_I_coeff);
          L_E = lincomb_sequence(L_E_comp, L_E_coeff);
      else %primal
          [L_I_coeff, L_E_coeff, b_coeff] = ...
              model.operators_ptr(model, model_data);
          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
      model.decomp_mode = old_decomp_mode;
    end

    if model.data_const_in_time
      operators_required = 0;
    end

  end % if operators required
  
% in the dual case it is Ut = -L_I^(-1)*v_l in the first timestep
if t_ind == t_ind_start && model.want_dual
  if strcmp(model.name_output_functional, 'theta') == 1
    v_l = Ut; % save v_l the riesz-representant in the theta case because it is required in the next timestep
  end
  Ut = - L_I \ Ut;
  U(:,t_column) = Ut;
  t_column = t_column + 1;
end
  
%computing the right-hand-side
if model.want_dual % dual
  if strcmp(model.name_output_functional, 'theta') == 1 && t_ind == t_ind_start
    rhs = L_E * Ut + v_l;
  else
    rhs = L_E * Ut;
  end
else % primal
  rhs = L_E * Ut + b;
end
%computing the solution
if isequal(L_I, speye(size(L_I))) || (model.want_dual && t_ind == t_ind_stop-1) %should fix the last timestep in the dual case
    Ut = rhs;
else
    Ut = L_I \ rhs;
end
  
if save_time_index(t_ind+2)
  U(:,t_column) = Ut;
  time_indices(t_column) = t_ind;
  time(t_column) = t;
  t_column = t_column + 1; 
end
  
end %end of for-loop

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

% output
sim_data.U = U;
sim_data.time = time;
sim_data.time_indices = time_indices;

% only in the primal case an output does make sense
if isfield(model,'name_output_functional') && model.want_dual == 0
  
% compute the riesz-representant of the functional
  [v, ~] = model.operators_output(model, model_data);
  
% the theta case is treated seperatly. See eop_theta_functional.
  if strcmp(model.name_output_functional, 'theta') == 1
    y = eop_theta_functional(sim_data.U, v);
    sim_data.y = y(:);
  else
    sim_data.y = (v(:)') * U;
  end
end