1 function flux = conv_flux_to_be_distributed_in_single_files(glob,U,params)
2 ´%function flux = conv_flux(glob,U,params)
4 % function computing the convective flux of a convection problem.
5 % flux is a 2xnpoints matrix, representing the x/y-coordinates of the
6 % velocity in the edge midpoints.
8 % required fields of params:
9 % t: real time value in case of time-dependent flux
11 % if name_flux == 'burgers_parabola' (reasonable domain [0,1]^2)
12 % then a parabolic velocity field going to the right is used with
13 % maximum velocity params.vmax and angle params.vrot_angle to
14 % the x-axis. The computed flux then is v * u^params.flux_pdeg. The input
15 % quantity U is assumed to be bounded by linfty-norm 1
17 % if name_flux == 'parabola' (reasonable domain [0,1]^2)
18 % c : velocity of parabola-profile in params.yrange, constant in x-direction
20 % if name_flux == 'gdl' (reasonable domain [0,1e-3]x[0,0.25e-3])
21 % v_in_max : max velocity
22 % a velocity profile is used, which is computed from an
23 % elliptic problem: parabolic inflow profile on the left,
24 % parabolic outflow on the right, noflow on bottom and middle
25 % of top. Dirichlet values in pressure on two top domains
27 % if name_flux == 'gdl2' (reasonable domain [0,1e-3]x[0,0.25e-3])
28 % lambda : factor for velocity computation: v = - lambda * grad(p)
29 % a velocity field is used, which is computed from an
30 % elliptic problem for the pressure by compute_pressure_gdl2
31 % (some gdl-elements in a row, Neumann-0 and dirichlet
32 % conditions between 4 and 1 for the pressure, linear
33 % decreasing along the channel)
34 % if name_flux == 'gdl_circ_par_mix' then a mixture of parabolic and
35 % circular hand crafted velocity profile is used.
36 % does not work good, produces singularities in fv simulation.
37 % if name_flux == 'gdl_pgrad_and_parabola' then a mixture of a parabola
38 % and a simulated pressure gradient is used. By this,
39 % symmetric boundary conditions can be used
41 % Function supports affine decomposition, i.e. different operation modes
42 % guided by optional field decomp_mode in params. See also the
43 % contents.txt for general explanation
45 % optional fields of params:
46 % mu_names : names of fields to be regarded as parameters in vector mu
47 % decomp_mode: operation mode of the function
48 % 'none' (default): no parameter dependence or decomposition is
49 % performed. output is as described above.
50 % 'components': For each output argument a cell array of output
51 % arguments is returned representing the q-th component
52 % independent of the parameters given in mu_names
53 % 'coefficients': For each output argument a cell array of output
54 % arguments is returned representing the q-th coefficient
55 % dependent of the parameters given in mu_names
57 % in 'coefficient' mode, the parameters in brackets are empty
59 % Bernard Haasdonk 4.4.2006
63 if ~isempty(glob) && size(glob,1) < size(glob,2)
64 warning('coordinates in variable glob are given row-wise, but expected them to be column-wise');
72 % determine affine_decomposition_mode as integer
73 decomp_mode = params.decomp_mode;
74 % flag indicating whether the computation respected the decomposition
75 respected_decomp_mode = 0;
77 %disp('halt in conv_flux')
80 % keep the following up-to-date with conv_flux_linearization!!!
81 %linear_fluxes = {
'gdl2',
'parabola',
'gdl_circ_par_mix',
'gdl',...
82 %
'gdl_pgrad_and_parabola'};
84 %
if ismember(params.name_flux, linear_fluxes)
85 %
if fluxname is linear fluxes
87 % evaluate linear flux
88 lin_flux = conv_flux_linearization(params,U,X,Y);
89 if decomp_mode == 0 % i.e.
'none'
90 flux.lambda = lin_flux.lambda;
91 if ~isempty(lin_flux.Vx) % otherwise size 0/0 => 0/1 !!!!!
92 flux.Fx = lin_flux.Vx(:).* U(:);
93 flux.Fy = lin_flux.Vy(:).* U(:);
98 elseif decomp_mode == 1 % i.e.
'components'
99 flux = cell(size(lin_flux));
100 Q = length(lin_flux);
102 % f.lambda = lin_flux{q}.lambda;
103 f.Fx = lin_flux{q}.Vx(:).* U(:);
104 f.Fy = lin_flux{q}.Vy(:).* U(:);
107 elseif decomp_mode == 2 % i.e.
'coefficients'
108 flux = lin_flux; % copy of coefficient vector from linear flux
110 respected_decomp_mode = 1;
112 %
if no linear flux, then select appropriate nonlinear flux
113 switch params.name_flux
114 case 'burgers_parabola'
115 % i.e. a parabolic velocity field times u^2, which is rotated by
116 % an arbitrary angle around the midpoint (0.5,0.5)
117 % u is assumed to be bounded by 1, such that lambda can be
118 % specified in advance
120 % NOTE: KEEP THIS DATA IMPLEMENTATION CONSISTENT WITH ITS DERIVATIVE
121 % IN conv_flux_linearization !!!
122 P = [X(:)
'-0.5; Y(:)'-0.5];
123 Rinv = [cos(-params.vrot_angle), - sin(-params.vrot_angle); ...
124 sin(-params.vrot_angle), cos(-params.vrot_angle) ];
125 RP = Rinv*P; % rotated coordinates as columns
126 RP(1,:) = RP(1,:) + 0.5;
127 RP(2,:) = RP(2,:) + 0.5;
129 V(1,:) = (1- RP(2,:)) .* RP(2,:) * params.vmax * 4;
132 flux.Fx = RV(1,:)
' .* U(:).^params.flux_pdeg;
133 flux.Fy = RV(2,:)' .* U(:).^params.flux_pdeg;
136 i = find(abs(U)>1.5);
138 disp(['U not bounded by 1 as assumed by flux, lambda may be ',...
141 flux.lambda = 1/ (params.vmax * 2 * umax^(params.flux_pdeg-1));
142 % NOTE: KEEP THIS DATA IMPLEMENTATION CONSISTENT WITH ITS DERIVATIVE
143 % IN conv_flux_linearization !!!
145 % i.e. constant velocity v = (params.flux_vx, params.flux_vy) times u^pdeg
146 % u is assumed to be bounded by 1, such that lambda can be
147 % specified in advance
148 Upower = real(U(:).^params.flux_pdeg);
149 flux.Fx = params.flux_vx * ones(size(U(:))) .* Upower;
150 flux.Fy = params.flux_vy * ones(size(U(:))) .* Upower;
153 i = find(abs(U)>1.5);
155 disp(['U not bounded by 1 as assumed by flux, lambda may be ',...
158 vmax = sqrt(params.flux_vy^2 + params.flux_vx^2);
159 flux.lambda = 1/ (vmax * 2 * umax^(params.flux_pdeg-1));
161 flux.Fx = zeros(size(U));
162 flux.Fy = zeros(size(U));
164 p_mu = spline_select(params);
165 [ breaks, coeffs, pieces, order ] = unmkpp(p_mu);
166 p_mu_d = mkpp(breaks, coeffs(:,1:order-1) .* repmat([order-1:-1:1],pieces,1));
169 denom = 1 + ppval(p_mu, X);
170 flux.Fx = - 0.001 * ppval(p_mu_d, X) ./ denom .* U';
171 flux.Fy = - 0.001 * ppval(p_mu_d, X).^2 .* Y ./ denom.^2 .* U';
172 flux.lambda = 1/ (max(max([flux.Fx, flux.Fy])));
175 flux.Fx = zeros(size(U));
176 flux.Fy = zeros(size(U));
177 % assumption: vmax = 1, umax = 1;
180 flux.lambda = 1/ (vmax * 2 * umax);
183 error('flux function unknown');
188 if decomp_mode>0 & respected_decomp_mode==0
189 error('function does not support affine decomposition!');