hmm_micro_local_model.m

function model = hmm_micro_local_model(params)
%function model = hmm_micro_local_model(params)
%
% function generating a RBmatlab model to be used for reduction of
% the micromodel with local operator. Params is currently ignored
%
% prototype is the following  Situation:
%
% Makromodell:  u_t +f(u)_x = 0 auf I=(a,\infty), f(u) = u^3
%
%               u(x,0)  = u_L,U_R (,so dass die Loesung eine
%                                  nichtklassische Welle ist)
%               u(a,t)  = a(t)    (Stoerfunktion, so dass "viele"
%                                  verschiedene Riemannprobleme
%                                  am nichklassischen Schock geloest
%                                  werden muessen)
%               
%
% Mikromodell:  
%       u^\eps_t+f(u^\eps)_x = \eps u^\eps_{xx} + \gamma \eps^2 u^\eps_{xxx}

mm = [];
if isempty(params)
  params = [];
end;

%specify default parameters of Frederikes Model:
mm.Time = 2.24*10^(-5);  % endtime of micromodel
if isfield(params,'Time')
  mm.Time= params.Time;
end;  
mm.ul = 4;            % left state
mm.ur = -2;           % right state
mm.alpha = 4;         % alpha parameter of regularization 
mm.n_x = 300;         % number of space steps
if isfield(params,'n_x')
  mm.n_x = params.n_x;         % number of space steps
end;  
mm.epsilon = 10^(-5); % regularization parameter
mm.dx=2*mm.epsilon/5; % space step width
mm.dt_factor = 0.06;
if isfield(params,'dt_factor')
  mm.dt_factor = params.dt_factor;
end;
%mm.w = max(abs(mm.ul),abs(mm.ur));
w = 4; % maximum to be expected value of u
mm.dt = mm.dt_factor*min([mm.dx/(3*w^2), ...
               mm.dx^2/(2*mm.epsilon), ...
               mm.dx^3/(3*mm.alpha*mm.epsilon^2)]);
% Zeitschrittweite fuer ul<-4 vorne 0.06 durch 0.01 ersetzen!!! 
mm.n_t=round(mm.Time/mm.dt);  % Anzahl Zeitschritte

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% generate RBmatlab model from Frederikes
model = nonlin_evol_model_default;
model.base_model = mm;

% specify parametrization
model.mu_names = {'ul','ur'};
%model.mu_ranges = {[1,4],[-5,-1]};
model.mu_ranges = {[1,4],[-3,-1]};
% explanation Frederike: smaller values of ul: too slow evolution
% other values of ur: Oscillations
% Bernard: changed to {[1,4],[-3,-1]}; in order to prevent
% divergent scheme

% copy some fields from base model to model:
model.epsilon = mm.epsilon;
model.alpha = mm.alpha;
model.nt = mm.n_t;
model.T = mm.Time;
model.dt = mm.dt;
model.nx = mm.n_x;
model.dx = mm.dx;
model.ul = mm.ul;
model.ur = mm.ur;

model.name = ['hmm_micro_local_',num2str(model.nt)];

model.xrange = [-40*model.dx, -40*model.dx+(model.nx-1)* ...
            model.dx];
model.xnumintervals = model.nx-1;

% set further fields of nonlin_evol_model:
model.inner_product_matrix_algorithm = @(model,model_data) ...
     eye(model.nx);
model.verbose = 10;
model.debug = 3;
model.plot = @my_plot;
model.no_lines = 0;
model.axis_equal = 0;
model.axis_tight = 0;
% only used for plotting:
model.range_lim = [-6,6];
model.gridtype = 'onedgrid';

% clip values after local evaluation
model.enable_range_limiting = 0;
model.enable_error_estimator = 0;
model.implicit_nonlinear = 0;
model.newton_solver = 0;

% initial values:
model.gen_model_data = @my_gen_model_data;

model.init_values_algorithm = @my_init_values_algorithm;
%model.init_values_algorithm = @init_values_affine_decomposed;
model.init_values_coefficients_ptr = @my_init_values_coefficients;
model.init_values_components_ptr = @my_init_values_components;

% local time stepping operator:
model.L_E_local_ptr = @my_L_E_local;
model.local_stencil_size = 2;

model.l2_error_sequence_algorithm = @my_l2_error_sequence_algorithm;
model.error_algorithm = @my_linftyl2_error_algorithm;

%@(U1,U2,grid,params) ...
%   sqrt(max(sum((U1-U2).^2),0));
% empirical interpolation data:
%model.ei_numintervals = [4,4]; 
model.ei_numintervals = [1,1]; 
model.Mmax = 300; % should be exact then
%model.Mmax = 160; % take 150 for sim, 10 for error-est
model.Mstrich = 0; % no error estimation

model.ei_detailed_savepath = [model.name,'_detailed'];
model.ei_operator_savepath = [model.name,'ei_operator_'];
model.ei_time_indices = 1:model.nt+1;

model.CRB_generation_mode = 'param-time-space-grid';
model.ei_target_error = 'interpol';

% fields for RB-generation:
model.error_norm = 'l2';
%                   'RB_extension_max_error_snapshot'}; 
model.RB_stop_timeout = 3*60*60; % 2 hours              
model.RB_stop_epsilon = 1e-5; 
model.RB_error_indicator = 'error'; % true error
%model.RB_error_indicator = 'estimator'; % Delta from rb_simulation
model.RB_stop_Nmax = 30;
model.RB_generation_mode = 'greedy_uniform_fixed';
model.RB_numintervals = model.ei_numintervals;
model.RB_detailed_train_savepath = model.ei_detailed_savepath;
model.orthonormalize  = @model_orthonormalize_qr;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function model_data = my_gen_model_data(model)
model_data = nonlin_evol_gen_model_data(model);
%params.X = (-40*model.dx:model.dx:-40*model.dx+(model.nx-1)* ...
%           model.dx);
%model_data.grid = onedgrid(params);
%model_data.grid.nelements = model.nx;
% neighbour indices of 1D grid:
%model_data.grid.A = model.dx * ones(1,model.nx);    

function res = my_init_values_algorithm(model,model_data)
glob = [];
if model.decomp_mode < 2
  glob = model_data.grid.X(:);
end;
res = init_values_affine_decomposed(glob,model);

function ucomp = my_init_values_components(x,model)
% indicator function for left and right state regions
u1 = zeros(length(x),1); 
i = find(x<0);
u1(i) = 1;
u2 = 1-u1;
ucomp = {u1,u2};
%for i=1:length(x)
%    if x(i)<0
%        u(i) = ul;
%    else
%        u(i) = ur;
%    end
%end

function usigma = my_init_values_coefficients(model)
usigma = [model.ul; model.ur];

function [inc,dummy] = my_L_E_local(model,model_data,ulocal_ext,index_local)
% index_local: vector of indeces of ulocal_ext, in which increment
%              values are desired.
if isempty(index_local) % global evaluation
  %inc = Phasen_lokal_space(ulocal,model.epsilon,model.alpha,index,...
  %                      model.dx,model.dt);
  u_new = Phasen_lokal(ulocal_ext,...
                       model.epsilon,model.alpha,...
                       index_local,model.dx,model.dt);
  inc = -(u_new - ulocal_ext)/model.dt;
%  keyboard;
else % local evaluation
  index_ext = model_data.grid_local_ext{1}.global_eind;
  [u_new,u_index] = Phasen_lokal(ulocal_ext,model.epsilon,...
                                 model.alpha,index_ext,model.dx, ...
                                 model.dt);

  [index_local_sorted, sort_index] = sort(index_local);
  % => index_local(sort_index) = index_local_sorted;
  % invert sorting:
%  keyboard;
%  invert_sort_index = ...;
  invert_sort_index = zeros(size(sort_index));
  invert_sort_index(sort_index) = 1:length(sort_index);
  if ~isequal(index_local_sorted(invert_sort_index),index_local)
    error('sorting is buggy!!');
  end;
  index = index_ext(index_local_sorted);
%  index = sort(index);
  if length(u_index)==length(index)
    u_new_M = u_new;
    u_index_M = u_index;
  else
    % u_new and u_index are possibly "too" large, pick out target
    % values: 
    posi = zeros(1,max(u_index));
    posi(u_index) = 1:length(u_index);
    posi2 = posi(index);
    u_new_M = u_new(posi2);
    u_index_M = u_index(posi2);
  end;
  if model.debug
    if ~isequal(u_index_M,index)
      error('u_index_M and index not identical!!!!');
    end;
  end;
  if isfield(model,'enable_range_limiting') && model.enable_range_limiting
    i = find(u_new_M<model.range_lim(1));
    if ~isempty(i)
      u_new_M(i) = model.range_lim(1);
      disp('performed range limiting!');
    end;
    i = find(u_new_M>model.range_lim(2));
    if ~isempty(i)
      u_new_M(i) = model.range_lim(2);
    end;
  end;
  
  inc = -(u_new_M(invert_sort_index) - ulocal_ext(index_local))/model.dt;
%  keyboard;

end;
if ~isempty(find(isnan(inc))) && isempty(find(isnan(ulocal_ext)))  
  warning('nan occuring!');
end;
%[u,index_n]=Phasen_lokal(u_sicher,epsilon,alpha,index,dx,dt);
%X = x(index_n);        
% Kontrolle/ global Loesung zum Zeitpunkt T+dt
%[u_kontrolle]=Phasen_lokal(u_sicher,epsilon,alpha,[],dx,dt);
% Fehler zwischen globaler und lokaler Loesung
%fehler = max(abs(u_kontrolle(index_n)-u));
dummy = [];

function p = my_plot(d,grid,params)
p = plot(grid.X,d,'LineWidth',2);
axis tight
if isfield(params,'range_lim') && isnumeric(params.range_lim)
  set(gca,'Ylim',params.range_lim);
else
  set(gca,'Ylim',[min(d),max(d)]);
end;
xlabel('x');

function res = my_l2_error_sequence_algorithm(U1,U2,grid,params)
res = sqrt(max(sum((U1-U2).^2),0));

%keyboard;

function res = my_linftyl2_error_algorithm(U1,U2,grid,params)
res = max(sqrt(max(sum((U1-U2).^2),0)));