1 function flux_lin = conv_flux_linearization(model,U,X,Y)
2 %
function flux_lin = conv_flux_linearization(model,[U],[X],[Y])
4 %
function computing the derivative of a convective flux by forward
5 % difference or exact evaluation. flux_lin consists of the fields
6 % Vx, Vy, representing the x/y-coordinates of the
7 % velocity in the edge midpoints.
8 % lambda : bound such that
9 % lambda * sup_u n_jl * f
'(u) <= 1
10 % e.g. lambda := 1/sup|v(x,y)|
11 % only reasonable in decomp_mode==0, otherwise not returned.
13 % For a linear flux in conv_flux the following is identical
15 % (Fx(U),Fy(U)) = (Vx,Vy) * U
17 % If the flux is linear, exact evaluation is performed, U is not used
18 % If the field use_velocity_matrixfile is set, then a filename in
19 % 'velocity_matrixfile
' is expected, where the point evaluation of
20 % the velocity are assumed to be stored.
21 % The file is expected to be located in [rbmatlabhome,'datafuncs/mat
'].
22 % The file is expected to have
23 % variable-vectors Vx, Vy, X, Y, and a scalar lambda. The name_flux
24 % parameter is in this case ignored (might be used during generation of
25 % the matrixfile by gen_velocity_matrixfile)
26 % As access to this matrix is always depending on
27 % expensive point-correspondence searches, a file-caching of the
28 % velocity extraction can be turned on by setting the field
29 % filecache_velocity_matrixfile_extract = 1
30 % This is useful, if the same point-lists X and Y are requested
31 % multiply, e.g. in timedependent case with constant velocity.
32 % Still this multiple filecaching can be expensive, instead by
33 % filecache_velocity_matrixfile_extract = 2, it is assumed, that
34 % the given file can directly be taken as velocities, i.e. the
35 % points coincide. For this a selection of the filename has to be
36 % done by the calling numerical routines.
38 % required fields of model as desired by conv_flux
40 % Function supports affine decomposition, i.e. different operation modes
41 % guided by optional field decomp_mode in model. See also the
42 % contents.txt for general explanation
44 % optional fields of model:
45 % mu_names : names of fields to be regarded as parameters in vector mu
46 % decomp_mode: operation mode of the function
47 % 'none
' (default): no parameter dependence or decomposition is
48 % performed. output is as described above.
49 % 'components
': For each output argument a cell array of output
50 % arguments is returned representing the q-th component
51 % independent of the parameters given in mu_names
52 % 'coefficients
': For each output argument a cell array of output
53 % arguments is returned representing the q-th coefficient
54 % dependent of the parameters given in mu_names
56 % In 'coefficient
' mode, the parameters in brackets are empty
58 % Bernard Haasdonk 21.4.2006
62 % determine affine_decomposition_mode as integer
63 decomp_mode = model.decomp_mode;
64 % flag indicating whether the computation respected the decomposition
65 respected_decomp_mode = 0;
67 % keep the following up-to-date with conv_flux!!!
68 %linear_fluxes = {'gdl2
','parabola
','gdl_circ_par_mix
','gdl
',...
69 % 'gdl_pgrad_and_parabola
'};
71 %if ismember(model.name_flux, linear_fluxes)
73 % check if flux is to be computed or read from file
74 if model.use_velocity_matrixfile==1
76 % affine decomposition is trivial: no dependency on mu
78 if decomp_mode < 2 % i.e. 'none
' or 'components
'
80 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% loading flux matrix.
81 % load V.X: list of X-Koordinaten
82 % load V.Y: list of Y-Koordinaten
83 % load V.Vy: list of Y-component of velocity
84 % load V.Vx: list of Y-component of velocity
85 % load V.lambda: lambda for CFL computation
87 if ~isfield(model,'velocity_matrixfile
')
88 error('field velocity_matrixfile not set!!
');
91 fullfn = fullfile(rbmatlabhome,'datafunc
','data
',...
92 model.velocity_matrixfile);
94 if ~exist(fullfn,'file
') & ~cache('exist
',fullfn)
95 error(['warning: velocity_matrixfile not existing.
',...
96 'Call gen_velocity_matrixfile with suitable
' ...
100 % the following file extraction can be expensive, so the
101 % extraction can additionally be cached.
102 if ~isfield(model,'filecache_velocity_matrixfile_extract
')
103 model.filecache_velocity_matrixfile_extract = 0;
106 if ~model.filecache_velocity_matrixfile_extract % expensive call
107 [flux_lin.Vx, flux_lin.Vy, lambda] = ...
108 velocity_matrixfile_extract(fullfn,X,Y);
109 elseif model.filecache_velocity_matrixfile_extract==1
110 % cached functioncall
111 [flux_lin.Vx, flux_lin.Vy, lambda] = ...
112 filecache_function('velocity_matrixfile_extract
', ...
114 else % mode 2, i.e assume, that velocityfilename is set correctly to
116 % V = load_cached(fullfn);
117 V = cache('load
',fullfn);
118 % find integer indices such that V.X(j)==X(i) and V.Y(j) = Y(i)
119 % if whole flux-matrix is requested
121 if ~isequal(X,V.X(:)) | ~isequal(Y,V.Y(:))
122 error('matrixfile does not fit to requested points!
');
128 % flux_lin.Vy = NaN*ones(size(X));
129 % flux_lin.Vx(i) = V.Vx(j);
130 % flux_lin.Vy(i) = V.Vy(j);
135 flux_lin.lambda = lambda;
137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of loading flux matrix.
138 if decomp_mode == 1 % put single matrix into cell array
139 flux_lin = {flux_lin};
141 else % decomp_mode = 2;
142 flux_lin = 1; % factor one for single component
144 respected_decomp_mode = 1;
146 else % model.use_velocity_matrixfile == 0
147 %%%%%%%%%%% compute the fluxes instead of loading matrix from file
149 if isequal(model.name_flux,'parabola
')
150 % determine lambda_jl to be globally constant such that
151 % lambda_jl * sup_u n_jl * f'(u) <= 1
153 % e.g. lambda := 1/sup|v(x,y)|
155 flux_lin.lambda = 1/(abs(model.c)/4.0+ 1e-10); % some eps for divby0
158 % flux.precomputable = 1; % velocity field can be computed in advance
159 % flux values : velocity field in the edge midpoints
162 elseif decomp_mode == 1
163 dummy.Vy = zeros(length(X),1);
164 dummy.Vx = Y(:) .* (1-Y(:));
166 elseif decomp_mode == 0 % decomp_mode 0
167 flux_lin.Vy = zeros(length(X),1);
168 flux_lin.Vx = model.c * Y(:) .* (1-Y(:));
170 error('unknown decmposition mode!!
');
172 respected_decomp_mode = 1;
173 elseif isequal(model.name_flux,'gdl2
')
174 % load pressure from file and generate velocity field
175 p = load(fullfile(rbmatlabhome,'datafunc
','data
','gdl_pressure2.mat
'),...
177 [ux, uy] = evaluate_gradient(X, Y,p.p,p.model);
178 flux_lin.Vx = -ux * model.lambda;
179 flux_lin.Vy = -uy * model.lambda;
180 % inverse of maximum absolute velocity value
181 % correct up to a constant factor
182 flux_lin.lambda = max([flux_lin.Vx(:); flux_lin.Vy(:)])^-1;
183 % flux.Fx = ux.* U(:);
184 % flux.Fy = uy.* U(:);
185 elseif isequal(model.name_flux,'gdl_circ_par_mix
')
186 % determine lambda_jl to be globally constant such that
187 % lambda_jl * sup_u n_jl * f'(u) <= 1
189 % e.g. lambda := 1/sup|v(x,y)|
191 flux_lin.lambda = 1/abs(model.v_max); % some eps for preventing divby0
193 % flux values : parabolic and circular profile
194 Vx1 = model.v_max * (100e-6)^(-2) * Y(:) .* (200e-6 -Y(:));
195 Vy1 = zeros(length(X),1);
197 % circular field by squared affine radial distance to M = (500,500)*1e-6
198 e = (X(:)-500e-6).^2 + (2.5*Y(:)-500e-6).^2;
199 % desired velocity profile:
200 % vabs(e) = linearly growing from 0 to vmax for 250<=e<=500
201 % vabs(e) = linearly decreasing to 0 for 500<=e<=750
202 % vabs(e) = zero otherwise
203 vabs = zeros(size(e));
205 i = find( se>= 250e-6 & se <= 500e-6);
206 vabs(i) = model.v_max / 250e-6 * (se(i)-250e-6);
207 i = find( se>= 500e-6 & se <= 750e-6);
208 vabs(i) = -model.v_max / 250e-6 * (se(i)-750e-6);
210 % compute velocity directions by gradient of elliptic field
211 dx_e = 2*(X(:)-500e-6);
212 dy_e = 5*(Y(:)-500e-6);
213 norm_grad_e_inv = (dx_e.^2 + dy_e.^2).^(-1);
214 % rotate by 90ΓΈ and normalization
215 Vx2 = - dy_e .* norm_grad_e_inv;
216 Vy2 = dx_e .*norm_grad_e_inv;
217 % scale to desired absolute velocity
221 flux_lin.Vx = (model.v_weight * Vx1 + (1-model.v_weight) * Vx2);
222 flux_lin.Vy = (model.v_weight * Vy1 + (1-model.v_weight) * Vy2);
223 elseif isequal(model.name_flux,'gdl
')
224 % load pressure from file and generate velocity field
225 p0 = load('gdl_pressure_v0000.mat
','p
','model
');
226 p1000 = load('gdl_pressure_v1000.mat
','p
','model
');
227 factor = 1000; % scaling: maximum inflow now 1000/factor m/s
228 [ux0, uy0] = evaluate_gradient(po0.model,X, Y,p0.p);
231 [ux1000, uy1000] = evaluate_gradient(p1000.model,X, Y,p1000.p);
232 ux1000 = -ux1000/factor;
233 uy1000 = -uy1000/factor;
234 % linear interpolation between velocity profiles yields exact solution
235 % of corresponding elliptic boundary value problem due to
236 % linearity of the problem.
237 ux = model.v_in_max * ux1000 + (1-model.v_in_max) * ux0;
238 uy = model.v_in_max * uy1000 + (1-model.v_in_max) * uy0;
240 % inverse of maximum absolute velocity value
241 flux_lin.lambda = max([ux(:); uy(:)])^-1;
242 % correct up to a constant factor
245 elseif isequal(model.name_flux,'gdl_pgrad_and_parabola
')
246 % load pressure from file and generate velocity field
247 p0 = load('gdl_pressure_dirichlet.mat
','p
','model
');
248 factor = 1000; % scaling: maximum inflow now 1000/factor m/s
249 [ux0, uy0] = evaluate_gradient(p0.model,X, Y,p0.p);
252 uypar = zeros(length(X),1);
253 uxpar = model.v_in_max * (model.yrange(2)-model.yrange(1))^-2 * 4 ...
254 * (Y(:)-model.yrange(1)) .* (model.yrange(2)-Y(:));
256 % linear interpolation between velocity profiles yields exact solution
257 % of corresponding elliptic boundary value problem due to
258 % linearity of the problem.
259 ux = model.v_in_max * uxpar + (1-model.v_in_max) * ux0;
260 uy = model.v_in_max * uypar + (1-model.v_in_max) * uy0;
262 % inverse of maximum absolute velocity value
263 flux_lin.lambda = max([ux(:); uy(:)])^-1; % correct up to a factor
266 elseif isequal(model.name_flux, 'transport
')
268 flux_lin = 1; % single component factor 1
269 elseif decomp_mode == 1 % components:
270 tmp.Vx = model.transport_x * ones(size(U));
271 tmp.Vy = model.transport_y * ones(size(U));
273 else % decomp_mode = 0, complete:
274 flux_lin.Vx = model.transport_x * ones(size(U));
275 flux_lin.Vy = model.transport_y * ones(size(U));
276 flux_lin.lambda = max([flux_lin.Vx(:);
277 flux_lin.Vy(:)])^-1; % correct up to a factor
280 respected_decomp_mode = 1;
281 elseif isequal(model.name_flux, 'richards
')
283 flux_lin.Vx = zeros(size(U));
284 flux_lin.Vy = zeros(size(U));
286 % p_mu = spline_select(model);
287 % [ breaks, coeffs, pieces, order ] = unmkpp(p_mu);
288 % p_mu_d = mkpp(breaks, coeffs(1:order-1) .* [order-1:-1:1]);
290 % denom = 1 + ppval(p_mu, X);
291 [res1, res2] = inv_geo_trans_derivative(model,X,Y,...
292 {[1],[2],[1,1],[1,2]},{[1],[2],[2,1],[2,2]},2);
293 d1 = res1{3} + res2{3};
294 d2 = res1{4} + res2{4};
295 flux_lin.Vx = -model.k .* (res1{1} .* d1 + res1{2} .* d2);
296 flux_lin.Vy = -model.k * (res2{1} .* d1 + res2{2} .* d2);
297 % flux_lin.Vx = model.k * ppval(p_mu_d, X) ./ denom;
298 % flux_lin.Vy = -model.k * ppval(p_mu_d, X).^2 .* Y ./ denom.^2;
299 % inverse of maximum absolute velocity value
300 flux_lin.lambda = max([flux_lin.Vx(:); flux_lin.Vy(:)])^-1;
301 % correct up to a factor
303 error('unknown flux!!
');
304 end; %%% end of flux select
307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308 % Non-linear fluxes %
309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
311 % if not is linear flux: exact computation or forward difference
312 elseif isequal(model.name_flux,'burgers_parabola
')
313 % NOTE: KEEP THIS DATA IMPLEMENTATION CONSISTENT WITH ITS FLUX
315 P = [X(:)'-0.5; Y(:)
'-0.5];
316 Rinv = [cos(-model.vrot_angle), - sin(-model.vrot_angle); ...
317 sin(-model.vrot_angle), cos(-model.vrot_angle) ];
318 RP = Rinv*P; % rotated coordinates as columns
319 RP(1,:) = RP(1,:) + 0.5;
320 RP(2,:) = RP(2,:) + 0.5;
322 V(1,:) = (1- RP(2,:)) .* RP(2,:) * model.vmax * 4;
325 flux_lin.Vx = RV(1,:)' .* model.flux_pdeg .* U(:).^(model.flux_pdeg-1);
326 flux_lin.Vy = RV(2,:)
' .* model.flux_pdeg .* U(:).^(model.flux_pdeg-1);
327 % flux.lambda = 1/model.vmax;
329 else % perform forward difference
330 epsilon = max(abs(U)) * 1e-2;
331 if epsilon < eps; % U completely 0
335 fluxU = conv_flux(model,U,X,Y);
336 fluxU2 = conv_flux(model,U+epsilon,X,Y);
338 flux_lin.Vx = (fluxU2.Fx-fluxU.Fx)/epsilon;
339 flux_lin.Vy = (fluxU2.Fy-fluxU.Fy)/epsilon;
340 flux_lin.lambda = fluxU.lambda;
343 if decomp_mode>0 & respected_decomp_mode==0
344 error('function does not support affine decomposition!
');