test_local_operator_evaluation.m

Go to the documentation of this file.

% small script testing, whether detailed and local operator 
% evaluation are the same for detailed_nonlin_evol_simulation

% Bernard Haasdonk 21.5.2007

load demo_nonlin_evol_params;
% only perform one step
params.T = params.T/params.nt;
params.nt = 1;
params.c_init = 1; % hereby nontrivial structure is visible
params.k = 0; % damit nur explizit

params.filecache_velocity_matrixfile_extract = 1;       
%filecache_clear;

grid = construct_grid(params);

% construct a subgrid with some cells

eindmask = zeros(grid.nelements,1);
% select randomly 10 points
rand('state',1234);
eindmask(ceil(rand(1,10)*grid.nelements)) = 1; 
for i = 1:5;
  % grow region
  eind2 = find(eindmask);
  NBI = grid.NBI(eind2,:);
  inner_edge = find(NBI>0);
  eindmask(NBI(inner_edge))= 1;
end;

% generate subgrid
eind = find(eindmask);
subgrid = gridpart(grid,eind);

% and determine inner elements
inner_els_mask = ones(subgrid.nelements,1);
[i,j] = find(subgrid.NBI==-10); % i.e.cut edge
inner_els_mask(i) = 0;
eind_local = find(inner_els_mask);

%subplot(2,1,1), plot(grid);
%subplot(2,1,2), plot(subgrid);

% generate initial data
params.affine_decomp_mode = 'complete';
U0 = feval(params.init_values_algorithm,grid,params);
U0sub = U0(eind);
params.dt = params.T/params.nt;

%detailed explicit step
inc = feval(params.L_E_local_name, ... 
            U0, [], ...
            grid, params);
U = U0 - params.dt * inc;

% local explicit step on subgrid
incsub = feval(params.L_E_local_name, ... 
                U0sub, eind_local, ...
                subgrid, params);
Usub_local = U0sub(eind_local) - params.dt * incsub;
Usub = zeros(subgrid.nelements,1);
Usub(eind_local) = Usub_local;

subplot(3,1,1), plot_element_data(grid,U,params);
title('complete solution after 1 step')
subplot(3,1,2), plot_element_data(subgrid,Usub,params);
title('local solution after 1 step')
set(gca,'Xlim',params.xrange,'Ylim',params.yrange);

err = zeros(subgrid.nelements,1);
err(eind_local) = Usub(eind_local) - U(eind(eind_local));
subplot(3,1,3), plot_element_data(subgrid,err,params);
set(gca,'Xlim',params.xrange,'Ylim',params.yrange);
maxerr = abs(max(err));
title(['local error after 1 step: max(err) =',num2str(maxerr)]);


% das folgende geht natuerlich nicht, weil 'cut'-fluxes nicht
% spezifiziert sind
%Usub = detailed_simulation(subgrid,params);

disp('detailed_simulation und local simulation give');
if maxerr ~= 0 
  disp('different results! Test failed');
else
  disp('identical result! Test OK');
end;



% TO BE ADJUSTED TO NEW SYNTAX
%| \docupdate