1 % demo of a
explicit FV scheme
for advection-diffusion
3 % Simple transport-diffusion problem: square domain with parabolic
4 % velocity profile, flow from left to right (or vice versa, depending on
5 % the sign of the specified velocity)
6 % upper half of the square is neuman-boundary with noflow conditions,
7 % lower half is dirichlet-boundary.
9 % A simulation for a square and a triangular grid is performed. For
10 % the triangular grid a more restrictive time-step is required, so
11 % 4 times as many time-slices are computed.
13 % Bernard Haasdonk 18.4.2006
17 model = lin_evol_model_default;
18 gparams = {
struct(
'gridtype',
'rectgrid'),...
19 struct(
'gridtype',
'triagrid')};
21 %%% 1. common parameters
22 disp(
'Full FV-simulation:');
23 disp(
'-------------------');
24 disp(
'setting parameters');
26 model.xrange = [0,1]; % grid x-coordinate range
27 model.yrange = [0,1]; % grid y-coordinate range
28 model.xnumintervals = 100; % number of grid cells in x-direction
29 model.ynumintervals = 100; % number of grid cells in y-direction
31 % set everything to dirichlet (first column)
32 % then the upper half boundary to neuman (second column)
33 model.bnd_rect_corner1 = [-1, -1; ...
35 model.bnd_rect_corner2 = [ 2 , 2; ...
37 model.bnd_rect_index = [-1, -2]; % first is dirichlet, second neuman
40 model.verbose = 20; % 0 : no output,
41 % >9 : informative output
42 % >19 :
verbose informative output
44 % >99 : desperate debug output, with dbstops
46 model.T = 0.25; % maximum time
48 %model.nt = 250; % number of time steps
53 % alternatively generate
triagrid: smaller dt required
for stability
54 %model.nt = model.nt*4;
56 %
set data functions and simulation information
58 model.neumann_values_ptr = @neumann_values_homogeneous;
62 % parameters inducing
"affine dependency" into the problem
63 model.diffusivity_ptr = @diffusivity_homogeneous;
65 model.k = 0.004; % diffusion coefficient 0 to 0.004 is OK
66 % although CFL condition is not completely
67 % met
for highest value
68 model.laplacian_ptr = @(glob, U, model) U;
69 model.laplacian_derivative_ptr = @(glob, U, model) ones(length(U),1);
74 model.flux_linear = 1;
76 model.c = -2; % maximum velocity 0 to 2 is OK
77 model.flux_quad_degree = 1;
79 model.init_values_ptr = @init_values_blobs;
80 model.beta = 0.25; % weighting of Gauss-blobs 0 to 1 is reasonable
81 model.gamma = 100; % 100-1000 width of Gauss-blob of initial values
82 model.radius = 0.1; % radius of cutoff-function
84 plot_params.show_colorbar = 1;
85 plot_params.no_lines = 1;
86 plot_params.axis_equal = 1;
87 plot_params.plot = @plot_element_data;
91 %%% 2.
"complete" simulation, complexities and plot
93 gridnames = {
'rectangular',
'triangular'};
97 disp([
'starting full simulation for ',gridnames{g},
' grid and ',...
98 ' numerical flux ',func2str(fluxes{g})]);
100 model.num_conv_flux_ptr = fluxes{g};
104 model_data = gen_model_data(model);
106 % use standard unit-square grid ...
108 % ... and update boundary data
109 model_data.grid = set_boundary_types(model_data.grid, model);
111 U = fv_conv_diff(model,model_data);
113 disp([
'===> Computation time of full simulation : ',num2str(t)]);
114 disp([
'===> Memory requirement for full simulation : ', ...
115 num2str(numel(U)*8/(1024*1024)),
' MB']);
116 % U = rand(model.nx*model.ny,model.nt+1);
118 model.title = ['Simulation ',gridnames{g},' grid, ',...
119 func2str(fluxes{g}),' numerical flux.'];