test_linear_convection.m

Go to the documentation of this file.

function OK = test_linear_convection
% function simulating a simple linear convection scheme with source term
% u_t + b*Du + cu = 0,               u(x,0) = u0.
% u0 is 1 inside a box and 0 otherwise. This box gets transported in
% direction b with a slowing term c.
% returns 1, if error to exact solution is "small enough".
% returns 0 otherwise

% Martin Drohmann 26.05.2008

OK = 1;

params = lin_evol_model_default;

% rectgrid or trigrid
params.gridtype = 'rectgrid';
% params.grid_initfile = '2dsquaretriang';

%params.gridtype = 'rectgrid';
params.xnumintervals = 40;
params.ynumintervals = 40;
params.xrange = [0,1];
params.yrange = [0,1];

params.bnd_rect_corner1 = [-1, -1];
params.bnd_rect_corner2 = [2, 2];

%params.bnd_rect_index = (-2);
%params.neumann_values_ptr = @neumann_values_outflow;
params.bnd_rect_index = (-1);
params.dirichlet_values_ptr = @dirichlet_values_homogeneous;
params.c_dir = 0.0;

params.T = 1.0;
params.nt = 400;

params.axis_equal = 1;

params.init_values_ptr = @init_values_gradient_box;
%params.implict_operators_algorithm = 'eye';
params.c_init_up = 1.0;
params.c_init_lo = 1.0;

%params.conv_flux_ptr = @velocity_transport;
params.velocity_ptr = @velocity_constant;
params.transport_x = 0.1;
params.transport_y = 0.1;
params.conv_flux_ptr = @conv_flux_linear;
params.flux_linear = 1;
params.divclean_mode = 0;
params.verbose = 1;
params.debug = false;
params.transport_source = 1;
params.no_lines = true;

params.filecache_velocity_matrixfile_extract = 0;

%%%% for using Lax-Friedrichs Flux:
%params.num_conv_flux_ptr  = @fv_num_conv_flux_lax_friedrichs;
%params.lxf_lambda = [];
%%%params.lxf_lambda = 100;
%%%% for using Engquist-Osher Flux:
params.num_conv_flux_ptr  = @fv_num_conv_flux_engquist_osher;
params.conv_flux_derivative_ptr = @conv_flux_forward_difference;
params.flux_quad_degree = 1;
params.fv_expl_diff_weight = 0;
params.fv_expl_conv_weight = 1;
params.fv_expl_react_weight = 0;

% for a detailed simulation we do not need the affine parameter
% decomposition
params.affine_decomp_mode = 'complete';

grid = construct_grid(params);

model_data.grid = grid;
model_data.W = spdiags(grid.A,0,grid.nelements,grid.nelements);
params.decomp_mode = 0;

% initial values by midpoint evaluation
U0 = fv_init_values(params, model_data);

U = zeros(length(U0(:)),params.nt+1);
U(:,1) = U0(:);
params.dt = params.T/params.nt;
dt = params.dt;

params.data_const_in_time = 1;

% loop over fv-steps
for t = 1:params.nt  
  params.t = t;
  if params.verbose > 19
    disp(['entered time-loop step ',num2str(t), ' for transport equation']);
  elseif params.verbose > 9
    fprintf('.');
  end;

  [inc, fluxes] = fv_explicit_space(params, model_data, U(:,t), []);

  U(:,t+1) = U(:,t) - params.dt * inc;

end

%%%% for plotting:
%plot_params.plot_function = @fv_plot
%plot_data_sequence(params,grid,U,plot_params)

EX = zeros(length(U0(:)), params.nt + 1);
EX(:,1) = U0(:);

for t = 1:params.nt
  EX(:,t+1) = exp(-params.transport_source * t * dt ) * ...
      params.init_values_ptr([grid.CX(:)-params.transport_x * t *dt, ...
      grid.CY(:) - params.transport_y * t * dt], params);
end

err = fv_l2_error(U(:,params.nt+1), EX(:,params.nt+1), model_data.W);
if err > 0.2
    OK = 0;
end