demos_parabolic.m

function res = demos_parabolic(step)
%
% step 1: simple simulation with implicit Euler and plot of results
% step 2: instability of expl. Euler
% step 3: smoothing effect of impl. Euler
% step 4: error convergence table for theta-scheme

% 6.1.2015 B. Haasdonk

if nargin < 1
  step = 1;
end;

switch step
 case 1 % example of stable implicit Euler
  disp('Example of Implicit Euler for Heat-Equation')
  % initialize model:
  params = [];
  params.theta = 1.0;  params.K = 100;
  %  params.theta = 0;  params.K = 1300; params.T = 1; % stable!
  %  params.theta = 0;  params.K = 1200; params.T = 1; % instable!
  %  params.theta = 0.5;  params.K = 100;  
  
  model = parabolic_model(params);
%  disp('temporary setting rhs to 0, comment this line out for real model')
%  model.source_function = @(X,model) zeros(size(X,1),1);

% generate grid and fem matrices:
  model_data = gen_model_data(model);
  % perform full simulation 
  sim_data = detailed_simulation(model, model_data); 
  % plot results and provide slider
  plot_params = [];
  plot_params.height_factor = 0.5;
  plot_sim_data(model,model_data,sim_data,plot_params);

 case 2 % demo of instability with explicit Euler

  disp('Example of Instability by Explicit Euler, CFL violation')
  % initialize model:
  params = [];
  params.theta = 0.0;  params.K = 100;
  %  params.theta = 0;  params.K = 1300; params.T = 1; % stable!
  %  params.theta = 0;  params.K = 1200; params.T = 1; % instable!
  %  params.theta = 0.5;  params.K = 100;  
  
  model = parabolic_model(params);
%  disp('temporary setting rhs to 0, comment this line out for real model')
%  model.source_function = @(X,model) zeros(size(X,1),1);

% generate grid and fem matrices:
  model_data = gen_model_data(model);
  % perform full simulation 
  sim_data = detailed_simulation(model, model_data); 
  % plot results and provide slider
  plot_params = [];
  plot_params.clim_tight = 1;
  plot_params.height_factor = 0.5;
  plot_sim_data(model,model_data,sim_data,plot_params);
  
 case 3 % demo of smoothing property of noisy initial data

  disp('Example of smoothing property of heat equation')
  % initialize model:
  params = [];
  params.theta = 1.0;  params.K = 500;
  params.xnumintervals = 40;
  params.ynumintervals = 40;
  %  params.theta = 0;  params.K = 1300; params.T = 1; % stable!
  %  params.theta = 0;  params.K = 1200; params.T = 1; % instable!
  %  params.theta = 0.5;  params.K = 100;  
  
  model = parabolic_model(params);
  model.source_function = @(X,model) zeros(size(X,1),1);
  model.init_value_function = @(X,model) sin(2*pi*X(:,1)) .* sin(2*pi*X(:,2)) ...
      + 0.2 * sin(20*pi*X(:,1)) .* sin(20*pi*X(:,2)) ...;
  
% generate grid and fem matrices:
  model_data = gen_model_data(model);
  % perform full simulation 
  sim_data = detailed_simulation(model, model_data); 
  % plot results and provide slider
  plot_params = [];
  plot_params.clim_tight = 1;
  plot_params.height_factor = 0.5;
  plot_sim_data(model,model_data,sim_data,plot_params);
  
 case 4 % error convergence tables 
  disp('Error convergence table for theta-scheme')

  params = [];
  thetas = [1.0,0.5,0.0,0.0];
  Ks = [10,10,10,640]; % coarse number of timesteps
  
  i_s = [0,1,2,3,4];  % refinement steps

  for mi = 1:length(thetas);
    theta = thetas(mi);
    K = Ks(mi);
    disp('===================================================================')
    disp(['convergence study for theta = ',num2str(theta),', K0 = ',num2str(K)]);
    disp('-----------------------------------------------------------------')  
    disp('    delta_t     | delta_x = delta_y |  ||u(t^K)-u_h^K||_L2 ')
    disp('-----------------------------------------------------------------')  
    for ni = 1:length(i_s);    
      params.theta = theta;  
      params.K = Ks(mi)*2^i_s(ni);
      params.xnumintervals = 4*2^i_s(ni);
      params.ynumintervals = 4*2^i_s(ni);
      %  params.theta = 0;  params.K = 1300; params.T = 1; % stable!
      %  params.theta = 0;  params.K = 1200; params.T = 1; % instable!
      %  params.theta = 0.5;  params.K = 100;  
      model = parabolic_model(params);
      % generate grid and fem matrices:
      model_data = gen_model_data(model);
      % perform full simulation 
      sim_data = detailed_simulation(model, model_data); 
      % plot results and provide slider
      %plot_sim_data(model,model_data,sim_data,[]);
      err = l2_error_at_end_time(model,model_data,sim_data);
      disp(['   ',num2str(model.delta_t,'%10.5f'),'      |      ',...
            num2str(1/params.xnumintervals,'%10.5f'),'      |      ', ...
            num2str(err,'%10.5f')]); 
    end;
  end;  
  
 otherwise
  error('step number unknown');
end;
    
function model = parabolic_model(params);
%default settings, can be overridden by params
model.source_function = @(X,model) sin(pi*X(:,1)).*sin(2*pi*X(:,2)).* ...
    (-5*pi*sin(5*pi*model.t) + 5 * pi*pi *cos(5*pi*model.t)); 
model.init_value_function = @(X,model) sin(pi*X(:,1)) .* sin(2*pi*X(:,2));
model.exact_solution = @(X,model) sin(pi*X(:,1)).*sin(2*pi*X(:,2)).* cos(5*pi*model.t);

%model.source_function = @(X,model) 100*sin(pi*X(:,1)).*sin(2*pi*X(:,1)).*sin(5*pi*model.t); 
%model.init_value_function = @(X,model) sin(2.5*pi*X(:,1)).*sin(2.5*pi*X(:,1)) ...
%    .* (X(:,1)<0.6) .*(X(:,1)>0.4) ...    
%    .* (X(:,2)<0.6) .*(X(:,2)>0.4);
model.xnumintervals = 10;
model.ynumintervals = 10;
model.xrange = [0,1];
model.yrange = [0,1];
model.bnd_rect_corner1 = [-1,-1];
model.bnd_rect_corner2 = [2,2];
model.bnd_rect_index = [-1];
model.t = 0;
model.T = 1;
model.K = 100; % number of timesteps
model.gen_model_data = @my_gen_model_data;
model.detailed_simulation = @my_detailed_simulation;
model.plot_sim_data = @my_plot_sim_data;

%copy fields from params
fs = fields(params);
for i=1:length(fs)
  model = setfield(model,fs{i},getfield(params,fs{i}));
end;
% compute dependent parameters
model.delta_t = model.T/model.K;

% prepare temporary rbmatlab information for matrices/grid generation;
tmp_model = lin_stat_model_default;
tmp_model.has_diffusivity = 1;
tmp_model.has_source = 1;
tmp_model.has_dirichlet_values = 1;
tmp_model.pdeg = 1; % linear finite elements
tmp_model.dimrange = 1; % scalar finite elements;
tmp_model.qdeg = 2; % quadrature order
%copy fields from model
fs = {'xnumintervals','ynumintervals','xrange','yrange'};
for i=1:length(fs)
  tmp_model = setfield(tmp_model,fs{i},getfield(model,fs{i}));
end;
tmp_model.gridtype = 'triagrid';
tmp_model.c_dir = 0;
tmp_model.diffusivity_tensor = @(glob,params) ...
    [ones(1,size(glob,1));...
     zeros(1,size(glob,1));...
     zeros(1,size(glob,1));...
     ones(1,size(glob,1))]';
tmp_model.source = @(glob,model) ...
    model.source_function(glob,model);
tmp_model.dirichlet_values = @(glob,model) ...
    zeros(1,size(glob,1));

tmp_model = elliptic_discrete_model(tmp_model);

model.rbmatlab_model = tmp_model;

function model_data = my_gen_model_data(model);
% function generating the grid and the FEM-matrices
model_data = [];

% construct grid and finite element information
tmp_model_data = gen_model_data(model.rbmatlab_model);
tmp_model = model.rbmatlab_model;
% extract information required for our problem
[A_diff , A_adv, A_reac, A_dirichlet, A_neumann, A_robin] = ...
    fem_matrix_parts_assembly(tmp_model,tmp_model_data.df_info);
model_data.A = spalloc(tmp_model_data.df_info.ndofs,tmp_model_data.df_info.ndofs,10); 
model_data.A = A_diff + A_adv + A_reac + A_neumann + A_robin + A_dirichlet;
model_data.M = tmp_model_data.df_info.l2_inner_product_matrix;
dirichlet_gids = tmp_model_data.df_info.dirichlet_gids;

model_data.M(dirichlet_gids,:)= 0;
model_data.M(:,dirichlet_gids)= 0;
model_data.M(dirichlet_gids,dirichlet_gids)= eye(length(dirichlet_gids));
model_data.grid = tmp_model_data.grid;
model_data.ndofs = tmp_model_data.df_info.ndofs;
model_data.df_info = tmp_model_data.df_info;

function sim_data = my_detailed_simulation(model,model_data);
% function realizing theta-scheme. Result is stored matrix sim_data.U
sim_data = [];
sim_data.U = NaN * ones(model_data.ndofs,model.K+1);
delta_t = model.delta_t;
theta = model.theta;
M = model_data.M;
A = model_data.A;
% initial data projection as first solution column:
init_value_vector = project_function(model,model_data, ...
                                   model.init_value_function);
uk = M\init_value_vector;
sim_data.U(:,1) = uk;
t = 0;
model.t = t;
model.rbmatlab_model.t = t;
f_function = @(X,tmp_model) model.source_function(X,tmp_model);
fk = project_function(model,model_data, f_function);
for k = 0:model.K-1;
  t = t+ delta_t;
  model.t = t;
  model.rbmatlab_model.t = t;
  fkplus1 = project_function(model,model_data,f_function);
  
  % solve theta scheme LGS system  
  S = M + theta* delta_t * A;
  rhs = (M - (1-theta)*delta_t * A) * uk ...
        + delta_t * theta * fkplus1 + delta_t * (1-theta)*fk;
  ukplus1 = S\rhs;
  
  % store result and prepare next iteration
  sim_data.U(:,k+2) = ukplus1;
  fk = fkplus1;
  uk = ukplus1;
  if isnan(norm(uk))
    return;
  end;
end;

function p = my_plot_sim_data(model,model_data,sim_data,plot_params);
plot_params.axis_equal = 1;
plot_params.axis_tight = 1;
df = femdiscfunc([],model_data.df_info);
plot_params.plot = @(grid,data, plot_params) my_plot_slice(grid,data,plot_params,df);
p = plot_sequence(sim_data.U, model_data.grid, plot_params);
%p = plot(model_data.grid,plot_params);

function p = my_plot_slice(grid,data,plot_params,df)
df.dofs = data;
plot(df,plot_params);
if isfield(plot_params,'height_factor')  
  view(-53,22)
%  V =  [0.5736   -0.8192    0.0000    0.1228
%       0.5481    0.3838    1.4735   -1.2027
%       0.6087    0.4263   -1.3267   12.3230
%       0         0         0    1.0000]
%  view(V);
  set(gca,'Zlim',[-0.5,0.5]);
xlabel('x_1');
ylabel('x_2');
end;

function f_vector = project_function(model,model_data,f_function);
% function providing the vector of l2 inner products of basis
% functions and the given function f_function
tmp_model = model.rbmatlab_model;
tmp_model.source = @(varargin) ...
      discrete_volume_values(f_function,varargin{:});
[r_source, r_dirichlet, r_neumann, r_robin] = ...
    fem_rhs_parts_assembly(tmp_model,model_data.df_info);
f_vector = r_source + r_neumann + r_robin + r_dirichlet;
f_vector(model_data.df_info.dirichlet_gids) = 0;

function err = l2_error_at_end_time(model,model_data,sim_data)
model.t = model.T;
model.rbmatlab_model.t = model.T;
exact_solution_projection = project_function(model,model_data, ...
                                   model.exact_solution);
exact_solution_vector = model_data.M\exact_solution_projection;
errsqr = (exact_solution_vector- sim_data.U(:,end))' * model_data.M * ...
      (exact_solution_vector - sim_data.U(:,end)); 
err = sqrt(errsqr);