newton_steps.m

Go to the documentation of this file.

function outfile = newton_steps(step)
% small script demonstrating a buckley leverett problem
% discretization and RB model reduction
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Select here, what is to be performed with the richards data %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%step = 1 % single detailed simulation with given data and plot. Run
          % this with varying parameters mu until sure that scheme
          % is stable. Modify dt or the data-functions accordingly,
          % until a nice parameter-domain with uniformly stable
          % detailed scheme is obtained.
%step = 2 % generate colateral reduced basis of L_E operator
%step = 3 % compare single detailed simulation with and without
          % interpolated space operator
%step = 4 % generate dummy reduced basis from single trajectory and check, if
          % ei_interpolation with projection on this space maintains
          % result. A simple reduced simulation can also be
          % performed. All results should be visually identical
%step = 5 % generate reduced basis
%step = 6 % time measurement of reduced simulation and
          % use reduced basis in rb_demo_gui

warning('off', 'RBMatlab:rb_simulation:runtime');

imsavepath = fullfile(rbmatlabresult, 'images');

params.model_type    = 'buckley_leverett';
%params.model_type    = 'exponential_diffusion';
%params.model_type     = 'eoc_nonlinear';
%params.model_size    = 'big';
%params.model_size    = 'medium';
% params.model_size    = 'small';
%  params.model_type    = 'nonlinear';
params.separate_CRBs = false;

% test parameters
params.RB_detailed_test_savepath = 'test_data_100';
params.RB_test_size              = 100;

params.verbose = 2;
params.debug   = 0;
extrapar.Mstrich = 20;
extrapar.N       = 50;
extrapar.M       = 150;

% email to which status messages about processing mu vectors is sent
% params.email   = 'mdrohmann@uni-muenster.de';

%% parameters for visualization
plot_params.show_colorbar = 1;
plot_params.colorbar_mode = 'EastOutside';
plot_params.no_lines      = true;
plot_params.plot_type     = 'patch';
% plot_params.clim          = [0, 0.5];
% plot_params.clim_tight    = true;
plot_params.plot          = @fv_plot;
plot_params.bind_to_model = true;
plot_params.yscale_uicontrols = 0.6;

% output-filenames in rbmatlabresult
CRBfname      = ['newton_', params.model_type, '_CRB_interpol.mat'];
if exist('model', 'var') && isfield(model, 'RB_error_indicator') ...
    && ~strcmp(model.RB_error_indicator, 'error')
  infix = model.RB_error_indicator;
else
  infix = '';
end
detailedfname = ['newton_', params.model_type, infix, '_detailed_interpol.mat'];
params.step7_outputfile = fullfile(rbmatlabresult, [params.model_type, infix, '_step7_output']);
params.step8_outputfile = fullfile(rbmatlabresult, [params.model_type, infix, '_step8_output']);
params.step0file       = [params.model_type, infix, 'step0'];
params.step1file       = [params.model_type, infix, 'step1'];

outfile = '';

switch step
 case 0 % initialize model data
   basis_gen      = newton_oo_model(params);

   model_data     = gen_model_data(basis_gen);
   save(fullfile(rbmatlabtemp, params.step0file), 'model','model_data')
   outfile = fullfile(rbmatlabtemp, params.step0file);
 case 1 % single detailed simulation and plot
   step1_detailed_simulation;
   outfile = fullfile(rbmatlabtemp, params.step1file);
 case 1.5
   load(fullfile(rbmatlabtemp, params.step0file));
   tmppar.model_type         = 'eoc';

   model = newton_model(tmppar);
   model.debug = 1;
   model.newton_solver = 0;
   model.T = 1;
   model.dt = model.T/model.nt;
   model.verbose = 21;
   err = zeros(model.nt+1,7);
   eoc = zeros(4, 1);
   for i = 1:7
     model_data   = gen_model_data(model);

     sim_data     = detailed_simulation(model, model_data);

     ffe = @exact_function_heat_equation;

     exact_data.U = zeros(model_data.grid.nelements, model.nt+1);

     for tt = 1:model.nt+1
       model.t = (tt-1) * model.dt;
       exact_data.U(:,tt) = ffe([model_data.grid.CX(:), model_data.grid.CY(:)], model);
     end
     model.t = 0;

     err(:,i)       = fv_l2_error(sim_data.U, exact_data.U, model_data.W);
     if(i>1)
       eoc(i-1)   = log(max(err(:,i-1))/max(err(:,i)))/log(2);
       disp(['Error: ', num2str(err(end,i)), 'EOC: ', num2str(eoc(i-1))]);
       pause;
     end

     model.xnumintervals = round(model.xnumintervals.*2);
     model.ynumintervals = round(model.ynumintervals.*2);
   end
   outfile = fullfile(rbmatlabtemp, [params.model_type, '_step1_5']);
   save(outfile, 'err', 'eoc');
 case 1.75
   load(fullfile(rbmatlabtemp, params.step0file));
   tmppar.model_type         = 'eoc_nonlinear';
   tmppar.local_mass_matrix  = @fv_local_mass_matrix_rect;
   tmppar.element_quadrature = @cubequadrature;
   tmppar.evaluate_basis     = @fv_evaluate_basis;
   tmppar.pdeg = 0;

   model = newton_model(tmppar);
   %model.debug = 1;
   % model.newton_regularisation = 0;
   %model.verbose         = 21;
   model.dt              = model.T/model.nt;
   model.ynumintervals = 1;
   err = zeros(model.nt+1,5);
   eoc = zeros(4, 1);
   for i = 1:6
     disp(['Computing EOC step no. ', num2str(i)]);
     model_data   = gen_model_data(model);

     sim_data     = detailed_simulation(model, model_data);

     ffe = @exact_function_plaplace;

     width = model_data.grid.X(2) - model_data.grid.X(1);

     Fe = @(einds, loc, grid, params) ffe([grid.X(einds) + loc(:,1)*width,...
                                           grid.Y(einds) + loc(:,2)*width], ...
                                           params);
     tmppar.diff_k0 = model.diff_k0;
     exact_data.U = zeros(model_data.grid.nelements, model.nt+1);
     for tt = 1:model.nt+1
       model.t = (tt-1) * model.dt;
       exact_data.U(:,tt) = ffe([model_data.grid.CX(:), ...
                                 model_data.grid.CY(:)], model);
%       exact_data.U(:,tt) = l2project(Fe, 2, model_data.grid, tmppar);
     end
     model.t = 0;

     err(:,i) = fv_l2_error(sim_data.U, exact_data.U, model_data.W);
     if(i>1)
       eoc(i-1)   = log(max(err(:,i-1))/max(err(:,i)))/log(2);
       disp(['Error: ', num2str(err(end,i)), ' EOC: ', num2str(eoc(i-1))]);
     end

     model.xnumintervals = model.xnumintervals.*2;
     model.nt = model.nt;
     model.dt = model.T / model.nt;
     %model.ynumintervals = model.ynumintervals.*2;
   end

   outfile = fullfile(rbmatlabtemp, [params.model_type, '_step1_75']);
   save(outfile, 'err', 'eoc');
 case 2 % empirical interpolation of space operator
   load(fullfile(rbmatlabresult, params.step0file));
   step2_empirical_interpolation;
   outfile = fullfile(rbmatlabresult, CRBfname);
 case 3 % compare detailed simulation with and without interpolation
   model.detailed_ei_simulation = @nonlin_evol_detailed_simulation;
   step3_detailed_ei_simulation;
 case 4 % construct dummy reduced basis by single trajectory and simulate
   model.detailed_ei_simulation = @nonlin_evol_detailed_simulation;
   detailed_data = gen_detailed_data(basis_gen, model_data);
   save(fullfile(rbmatlabresult, detailedfname), 'detailed_data', 'basis_gen');
%   step4_dummy_reduced_basis;
 case 5 % reduced basis
%   tmp=load(fullfile(rbmatlabresult,CRBfname));
%   detailed_data = tmp.detailed_data;
   step5_rb_generation;
   outfile = fullfile(rbmatlabresult, detailedfname);
 case 6
%   load(fullfile(rbmatlabresult,detailedfname));
   step6_demo_rb_gui;
 case 7 % training-error landscape & time comparisons
  tmp=load(fullfile(rbmatlabresult,detailedfname));
  detailed_data = tmp.detailed_data;


  model.Mstrich   = 0;
  % maximum number of reduced basis functions
  model.Nmax      = size(detailed_data.RB,2);
  % maximum number of collateral reduced basis functions
  model.Mmax      = size(detailed_data.BM{1},2);
  % number of reduced basis sizes to be tested.
  Nsize           = 7;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 10;
  % sample of values for the coupling variable 'c' for which errors are
  % computed.
  csample     = 0.1:0.1:1;
  % axis description for plots with coupling variable 'c'
  cdescr      = '(N,M) = c * (Nmax, Mmax)';
  % size of mu vector test set.
  mu_set_size = 20;

  error_plots = { 'train_set', 'train_set_coupled', ...
                  'error', 'error_coupled' };

  model.M_by_N_ratio = model.Mmax/model.Nmax;

  step7_error_landscape;
  outfile = params.step7_outputfile;

case 7.5
  load(params.step7_outputfile);
  for i = 1:length(output)
    plot_error_landscape(output{i});
  end
case 8
  disp('warning: takes a few hours!');
  tmp = load(fullfile(rbmatlabresult,detailedfname));
  detailed_data = tmp.detailed_data;

  % maximum number of reduced basis vectors
  model.Nmax = model.RB_stop_Nmax;
  % maximum number of collateral reduced basis vectors used for
  model.Mmax = size( detailed_data.BM{1}, 2 ) - extrapar.Mstrich;

  model.M_by_N_ratio = model.Mmax/model.Nmax;

  % number of reduced basis sizes to be tested.
  Nsize           = 7;
  % number of collateral reduced basis sizes to be tested.
  Msize           = 12;
  % Mstrich values to be tested
  Mstrich_samples = 1:extrapar.Mstrich;
  % sample of values for the coupling variable 'c' for which
  % estimators are computed.
  csample         = 0.1:0.1:1;
  % size of mu vector test set
  mu_set_size     = 20;
  % axis description for plots with coupling variable 'c'
  cdescr        = '(N,M) = c * (Nmax, Mmax)';

  step8_estimator_landscape;
  outfile = params.step7_outputfile;
case 8.5
  load(params.step8_outputfile);
  for i = 1:length(output)
    output{i}.bound = 1500;
    output{i}.stab_limit = 10000;
    plot_error_landscape(output{i});
  end
case 9
  tmp = load(fullfile(rbmatlabresult,detailedfname));
  detailed_data = tmp.detailed_data;
  model_data = detailed_data;
  num_samples = 20;
  testdir = [params.model_type, '_', num2str(num_samples), '_test'];
  dirnames = testdir;
  model.Nmax = size(detailed_data.RB,2);

  % generate test data if not existent
  rand('state',123456);
  M = rand_uniform(num_samples,model.mu_ranges);
  save_detailed_simulations(model,model_data,M,testdir)

  % load demo_nonlin_evol_detailed_data_LE_on_RB;
  reduced_data = gen_reduced_data(model,detailed_data);

  settings = load(fullfile(rbmatlabtemp,...
                  testdir,'settings.mat'));

  Msamples = 12;
  Nsamples = 7;
  % storage for maximum l2-errors (i.e. linfty([0,T],l2)-norm)
  errs   = zeros(Msamples,Nsamples);
  model.RB_error_indicator='ei_estimator_test';

  Mstrich_set = floor(1:9.9:100);%[1,2,3,4,5,10,15,20];

  deltas  = cell(size(M,2),length(Mstrich_set));
  lambda  = zeros(size(M,2),length(Mstrich_set));
  lambdas = cell(size(M,2),length(Mstrich_set));
  errors  = cell(size(M,2),length(Mstrich_set));
  parfor mu_ind = 1:size(M,2);
    disp(['processing parameter vector ',num2str(mu_ind),...
          '/',num2str(size(M,2))]);
    deltas_tmp  = cell(1, length(Mstrich_set));
    errors_tmp  = cell(1, length(Mstrich_set));
    lambdas_tmp = cell(1, length(Mstrich_set));
    lambda_tmp  = zeros(1, length(Mstrich_set));
    tsettings = settings;
    textrapar = extrapar;
    tmodel    = model;
    for msi = 1:length(Mstrich_set)

      mtemp = tmodel.set_mu(tmodel,tsettings.M(:,mu_ind));
      sim_data = load_detailed_simulation(mu_ind,tsettings.savepath,tmodel);
      U_H = sim_data.U;

      mtemp.M       = textrapar.M;
      mtemp.Mstrich = Mstrich_set(msi);
      mtemp.N       = textrapar.N;
      disp(['with CRB size ',num2str(mtemp.M),...
            ' Mstrich'' = ',num2str(mtemp.Mstrich),...
            ' RB size ', num2str(mtemp.N)]);
      mtemp.C_I = 1 - mtemp.diff_m;
      simulation_data = rb_simulation(mtemp,reduced_data);
      disp('.');

      deltas_tmp{msi} = simulation_data.Delta;
      simulation_data = rb_reconstruction(mtemp,detailed_data,simulation_data);
      errors_tmp{msi} =  tmodel.error_algorithm(simulation_data.U, ...
                                                  sim_data.U, ...
                                                  detailed_data.W, mtemp);
      lambdas_tmp{msi} = deltas_tmp{msi} ./ errors_tmp{msi};
      lambda_tmp(msi) = max(deltas_tmp{msi}) / max(errors_tmp{msi});
      disp(['true error: ', num2str(max(errors_tmp{msi})), ...
            ' estimator: ', num2str(max(deltas_tmp{msi}))]);
      disp(['effectivity: ', num2str(lambda_tmp(msi))]);

%      plot_sim_data(model,detailed_data,simulation_data,plot_params);

    end
    deltas(mu_ind,:)  = deltas_tmp;
    errors(mu_ind,:)  = errors_tmp;
    lambdas(mu_ind,:) = lambdas_tmp;
    lambda(mu_ind,:)  = lambda_tmp;
%    plot(Mstrich_set,deltas(mu_ind,:));
%    keyboard;
  end
  save(fullfile(rbmatlabresult, 'step9'), 'deltas','errors','lambdas',...
                                          'lambda','M','Mstrich_set');

  outfile = fullfile(rbmatlabresult, 'step9');
case 9.8
  tmp = load(fullfile(rbmatlabresult, detailedfname));
  detailed_data = tmp.detailed_data;

  %% initial data plot (as patch)
  mu_set = {[1,0.00,0.2]};
  plot_params.plot_type = 'patch';
  plot_params.line_width = 1e-10;
  plot_params.no_lines = true;
  timeinstants = 0;
  step10_plot_trajectories;

  %% some sample plots (as contour plot)
  mu_set = {[1,0.01,0.2], [2,0.01,0.2], [4,0.01,0.2], [1,0.01,0], [2,0.01,0],[4,0.02,0]};
%  mu_set = {[1,0.01,0], [1,0.01,0.2]};
%  mu_set = {[0,0.01], [0.1,0.01],[0.2,1.3]};
%  clim   = [0.1, 0.7];
  plot_params.plot_type = 'contour';
  plot_params.contour_lines = 5;
  plot_params.line_width = 2;
  timeinstants = [0:1/4:1]*model.T;
  boxed_snapshots = true;
  %step10_plot_trajectories;

  %% interpolation dofs
  plot_params.plot_type = 'patch';
  plot_params.colormap  = 1/length(detailed_data.TM{1})*repmat((length(detailed_data.TM{1}):-1:1)',1,3);
  plot_params.show_colorbar = 0;
  u = zeros(detailed_data.grid.nelements,1);
  u(detailed_data.TM{1}) = 1:length(detailed_data.TM{1});
  u(1)=50;
  u(end)=50;
  plot_params.plot(detailed_data.grid, u, plot_params);
  title(['Interpolation DOFS for implicit and explicit operator']);
  set(gcf, 'Color', 'none');
  set(gcf, 'InvertHardCopy', 'off');
  showaxes('hide');
  tikzparams.filepath      = fullfile(imsavepath, params.model_type);
  export_fig(fullfile(tikzparams.filepath, 'ei_dofs'),'-pdf','-a0',gcf);
  system(['convert -density 1200x1200 -trim ', ...
      fullfile(tikzparams.filepath, 'ei_dofs.pdf '), ...
      fullfile(tikzparams.filepath, 'ei_dofs.png')]);
  system(['rm ', fullfile(tikzparams.filepath, 'ei_dofs.pdf')]);
  close(gcf);

  %% error convergence of crb generation
  plot(detailed_data.ei_info{1}.max_err_sequence);
  set(gca,'Yscale','log');
  title('EI-interpolation error decrease');
  matlab2tikz(fullfile(tikzparams.filepath, 'ei_error_decrease.tikz'));

  %% error convergence of POD greedy
  plot(detailed_data.RB_info.max_err_sequence);
  set(gca,'Yscale','log');
  title('POD greedy error decrease');
  matlab2tikz(fullfile(tikzparams.filepath, 'pod_greedy_error_decrease.tikz'));
  close(gcf);

  %% ei snapshots
  Msize  = size(detailed_data.QM{1},2);
%  colorm = colormap(hot(Msize));
%  colorm = colorm(Msize:-1:1,:);
  emus = detailed_data.ei_info{1}.extension_mus([1,3],1:Msize);
  midmu = detailed_data.ei_info{1}.extension_mus(2,1:Msize);
  coord1 = emus(:,midmu==0.01);
  coord2 = emus(:,midmu==0.005);
  [coord1u,coord1I,coordJ] = unique(sortrows(coord1'),'rows');
  count1 = [coord1I(1); coord1I(2:end) - coord1I(1:end-1)]*10;
  [coord2u,coord2I,coordJ] = unique(sortrows(coord2'),'rows');
  count2 = [coord2I(1); coord2I(2:end) - coord2I(1:end-1)]*10;

  scatter(coord1u(:,1),coord1u(:,2), count1, 'black', 'filled');
  hold on;
  scatter(coord2u(:,1),coord2u(:,2), count2, [0.7 0.7 0.7], 'd','filled');
  hold off;

  set(gcf, 'Color', 'none');
  set(gcf, 'InvertHardCopy', 'off');

  showaxes('hide');

  im = export_fig(fullfile(tikzparams.filepath, 'ei_mu_snapshots'), '-png','-a3',gcf);
  close(gcf);

  %% RB snapshots
  Nsize  = size(detailed_data.RB,2);
%  colorm = colormap(hot(Msize));
%  colorm = colorm(Msize:-1:1,:);
  emus = detailed_data.RB_info.mu_sequence([1,3],1:Nsize);
  midmu = detailed_data.RB_info.mu_sequence(2,1:Nsize);
  coord1 = emus(:,midmu==0.01);
  coord2 = emus(:,midmu==0.005);
  [coord1u,coord1I,coordJ] = unique(sortrows(coord1'),'rows');
  count1 = [coord1I(1); coord1I(2:end) - coord1I(1:end-1)]*100;
  [coord2u,coord2I,coordJ] = unique(sortrows(coord2'),'rows');
  count2 = [coord2I(1); coord2I(2:end) - coord2I(1:end-1)]*100;

  scatter(coord1u(:,1),coord1u(:,2), count1, 'black', 'filled');
  hold on;
  scatter(coord2u(:,1),coord2u(:,2), count2, [0.7 0.7 0.7], 'd','filled');
  hold off;

  set(gcf, 'Color', 'none');
  set(gcf, 'InvertHardCopy', 'off');

  showaxes('hide');

  im = export_fig(fullfile(tikzparams.filepath, 'rb_mu_snapshots'), '-png','-a3',gcf);
  close(gcf);

  % error convergence of random samples


%  scatter3(coord(1,:),coord(2,:),coord(3,:),20,colorm,'filled');
%  xlabel(['\mu_1 = ',model.mu_names{1}]);
%  ylabel(['\mu_3 = ',model.mu_names{3}]);
  %ylabel(['\mu_2 = c_{init}']);%,model.mu_names{2}]);
%  zlabel('time index k');
%  title(['snapshots selected for collateral basis for operator no.'])
%  set(gca,'Box','on');
%  save(fullfile(tikzparams.filepath, 'ei_mu_snapshots.dat'), 'coord', '-ascii', '-double', '-tabs');
%  matlab2tikz(fullfile(tikzparams.filepath, 'ei_mu_snapshots.tikz'));
%  close(gcf);


case 10
  %tmp = load(fullfile(rbmatlabresult,'strange.mat'));
  %Cmap = tmp.Cmap;
  tmp                  = load(fullfile(rbmatlabresult,detailedfname));
  detailed_data        = tmp.detailed_data;
  model                = tmp.model;
  plot_params.no_lines = 1;

  plot_params.no_axes                = 1;
  plot_params.transparent_background = 1;
  plot_params.show_colorbar          = 0;
  %plot_params.shrink_factor         = 1.13;
  %plot_params.colormap               = Cmap;
  plot_params.clim                   = [0.1,0.6];

  reduced_data     = gen_reduced_data(model, detailed_data);
  model.N = size(detailed_data.RB,2);
  model.M = size(detailed_data.QM,2);
  reduced_data = extract_reduced_data_subset(model, reduced_data);

  cparams.range = cell(length(model.mu_ranges)+1,1);
  cparams.range{1}  = [1, model.nt];
  cparams.range(2:end) = model.mu_ranges(:);
  cparams.numintervals = 2*ones(length(model.mu_ranges)+1,1);
  cgrid = cubegrid(cparams);

  parammat = cgrid.vertex;
  parammat(:,1) = round(parammat(:,1));

  for i=1:length(parammat)
    tstep = parammat(i,1);
    newmodel = model.set_mu(model, parammat(i,2:end));
    rb_sim_data = rb_simulation(newmodel, reduced_data);
    rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);

    % transformed version
    plot_element_data(model_data.grid, rb_sim_data.U(:,tstep), plot_params);

    mustring = strrep(strrep(strrep(strrep(mat2str(parammat(i,2:end)),' ', '_'),'[',''),']',''),'.','p');
    fp = fullfile(imsavepath, params.model_type);
    if ~exist(fp,'dir')
      mkdir(fp);
    end
    fp = fullfile(fp, ['t_', num2str(tstep)]);
    if ~exist(fp,'dir')
      mkdir(fp);
    end
    fn = ['sample_mu_', mustring];
    disp(['Processing ', fn, ' ...']);

    tikzparams.filename      = fn;
    tikzparams.filepath      = fp;
    tikzparams.save_colorbar = 1;

    plot_as_tikzfile(model, tikzparams);

    close(gcf);
  end

case 11
  load('richards_affine_detailed_interpol');

  rdata = gen_reduced_data(model, detailed_data);

  model.N = size(detailed_data,2);
  model.M = 10;
  model.Mstrich = 1;
  rdata2 = extract_reduced_data_subset(model, rdata);
  model.Mstrich = 10;
  rdata1 = extract_reduced_data_subset(model, rdata);

  TM1 = rdata1.TM_local{1}(1:10);
  TM2 = rdata2.TM_local{1}(1:10);

  gl1 = rdata1.grid_local_ext{1};
  gl2 = rdata2.grid_local_ext{1};

  if(any(gl1.CX(TM1) ~= gl2.CX(TM2)))
    disp('CX differs');
  end
  if(any(gl1.CY(TM1) ~= gl2.CY(TM2)))
    disp('CY differs');
  end
  if(any(any(gl1.X(gl1.VI(TM1,:)) ~= gl2.X(gl2.VI(TM2,:)))))
    disp('X differs');
  end
  if(any(any(gl1.Y(gl1.VI(TM1,:)) ~= gl2.Y(gl2.VI(TM2,:)))))
    disp('Y differs');
  end
  if(any(any(gl1.SX(TM1,:) ~= gl2.SX(TM2,:))))
    disp('SX differs');
  end
  if(any(any(gl1.SY(TM1,:) ~= gl2.SY(TM2,:))))
    disp('SY differs');
  end
    if(any(gl1.EL(TM1,:) ~= gl2.EL(TM2,:)))
      disp('EL differs');
    end
    if(any(gl1.DC(TM1,:) ~= gl2.DC(TM2,:)))
      disp('DC differs');
    end
    if(any(gl1.NX(TM1,:) ~= gl2.NX(TM2,:)))
      disp('NX differs');
    end
    if(any(gl1.NY(TM1,:) ~= gl2.NY(TM2,:)))
      disp('NY differs');
    end
    if(any(gl1.ECX(TM1,:) ~= gl2.ECX(TM2,:)))
      disp('ECX differs');
    end
    if(any(gl1.ECY(TM1,:) ~= gl2.ECY(TM2,:)))
      disp('ECY differs');
    end
  %% neighbours
  NBI1 = gl1.NBI(TM1,:);
  NBI2 = gl2.NBI(TM2,:);
  if(any(gl1.CX(NBI1) ~= gl2.CX(NBI2)))
    disp('CX differs');
  end
  if(any(gl1.CY(NBI1) ~= gl2.CY(NBI2)))
    disp('CY differs');
  end
  if(any(any(gl1.X(gl1.VI(NBI1,:)) ~= gl2.X(gl2.VI(NBI2,:)))))
    disp('X differs');
  end
  if(any(any(gl1.Y(gl1.VI(NBI1,:)) ~= gl2.Y(gl2.VI(NBI2,:)))))
    disp('Y differs');
  end
  if(any(any(gl1.SX(NBI1,:) ~= gl2.SX(NBI2,:))))
    disp('SX differs');
  end
  if(any(any(gl1.SY(NBI1,:) ~= gl2.SY(NBI2,:))))
    disp('SY differs');
  end
  for i = 1:length(TM1)
  end

otherwise
  error('step-number is unknown!');
end

end

%| \docupdate