1 function flux_lin = velocity_to_be_distributed_in_single_files(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 affine_decomp_mode=='none
', 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 affine_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 % affine_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 = get_affine_decomp_mode(model);
64 decomp_mode = model.decomp_mode;
65 % flag indicating whether the computation respected the decomposition
66 respected_decomp_mode = 0;
68 % keep the following up-to-date with conv_flux!!!
69 %linear_fluxes = {'gdl2
','parabola
','gdl_circ_par_mix
','gdl
',...
70 % 'gdl_pgrad_and_parabola
'};
72 %if ismember(model.name_flux, linear_fluxes)
74 % check if flux is to be computed or read from file
75 if isfield(model,'use_velocity_matrixfile
') & ...
76 (model.use_velocity_matrixfile==1)
78 % affine decomposition is trivial: no dependency on mu
80 if decomp_mode < 2 % i.e. 'none
' or 'components
'
82 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% loading flux matrix.
83 % load V.X: list of X-Koordinaten
84 % load V.Y: list of Y-Koordinaten
85 % load V.Vy: list of Y-component of velocity
86 % load V.Vx: list of Y-component of velocity
87 % load V.lambda: lambda for CFL computation
89 if ~isfield(model,'velocity_matrixfile
')
90 error('field velocity_matrixfile not set!!
');
93 fullfn = fullfile(rbmatlabhome,'datafunc
','mat
',...
94 model.velocity_matrixfile);
96 if ~exist(fullfn,'file
') & ~cache('exist
',fullfn)
97 error(['warning: velocity_matrixfile not existing.
',...
98 'Call gen_velocity_matrixfile with suitable
' ...
102 % the following file extraction can be expensive, so the
103 % extraction can additionally be cached.
104 if ~isfield(model,'filecache_velocity_matrixfile_extract
')
105 model.filecache_velocity_matrixfile_extract = 0;
108 if ~model.filecache_velocity_matrixfile_extract % expensive call
109 [flux_lin.Vx, flux_lin.Vy, lambda] = ...
110 velocity_matrixfile_extract(fullfn,X,Y);
111 elseif model.filecache_velocity_matrixfile_extract==1
112 % cached functioncall
113 [flux_lin.Vx, flux_lin.Vy, lambda] = ...
114 filecache_function('velocity_matrixfile_extract
', ...
116 else % mode 2, i.e assume, that velocityfilename is set correctly to
118 % V = load_cached(fullfn);
119 V = cache('load
',fullfn);
120 % find integer indices such that V.X(j)==X(i) and V.Y(j) = Y(i)
121 % if whole flux-matrix is requested
123 if ~isequal(X,V.X(:)) | ~isequal(Y,V.Y(:))
124 error('matrixfile does not fit to requested points!
');
130 % flux_lin.Vy = NaN*ones(size(X));
131 % flux_lin.Vx(i) = V.Vx(j);
132 % flux_lin.Vy(i) = V.Vy(j);
137 flux_lin.lambda = lambda;
139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of loading flux matrix.
140 if decomp_mode == 1 % put single matrix into cell array
141 flux_lin = {flux_lin};
143 else % decomp_mode = 2;
144 flux_lin = 1; % factor one for single component
146 respected_decomp_mode = 1;
148 else %%%%%%%%%%% compute the fluxes instead of loading matrix from file
150 if isequal(model.name_flux,'parabola
')
151 % determine lambda_jl to be globally constant such that
152 % lambda_jl * sup_u n_jl * f'(u) <= 1
154 % e.g. lambda := 1/sup|v(x,y)|
156 flux_lin.lambda = 1/(abs(model.c)/4.0+ 1e-10); % some eps for divby0
159 % flux.precomputable = 1; % velocity field can be computed in advance
160 % flux values : velocity field in the edge midpoints
163 elseif decomp_mode == 1
164 dummy.Vy = zeros(length(X),1);
165 dummy.Vx = Y(:) .* (1-Y(:));
167 elseif decomp_mode == 0 % decomp_mode 0
168 flux_lin.Vy = zeros(length(X),1);
169 flux_lin.Vx = model.c * Y(:) .* (1-Y(:));
171 error('unknown decmposition mode!!
');
173 respected_decomp_mode = 1;
174 elseif isequal(model.name_flux,'gdl2
')
175 % load pressure from file and generate velocity field
176 p = load(fullfile(rbmatlabhome,'datafunc
','mat
','gdl_pressure2.mat
'),...
178 [ux, uy] = evaluate_gradient(X, Y,p.p,p.model);
179 flux_lin.Vx = -ux * model.lambda;
180 flux_lin.Vy = -uy * model.lambda;
181 % inverse of maximum absolute velocity value
182 % correct up to a constant factor
183 flux_lin.lambda = max([flux_lin.Vx(:); flux_lin.Vy(:)])^-1;
184 % flux.Fx = ux.* U(:);
185 % flux.Fy = uy.* U(:);
186 elseif isequal(model.name_flux,'gdl_circ_par_mix
')
187 % determine lambda_jl to be globally constant such that
188 % lambda_jl * sup_u n_jl * f'(u) <= 1
190 % e.g. lambda := 1/sup|v(x,y)|
192 flux_lin.lambda = 1/abs(model.v_max); % some eps for preventing divby0
194 % flux values : parabolic and circular profile
195 Vx1 = model.v_max * (100e-6)^(-2) * Y(:) .* (200e-6 -Y(:));
196 Vy1 = zeros(length(X),1);
198 % circular field by squared affine radial distance to M = (500,500)*1e-6
199 e = (X(:)-500e-6).^2 + (2.5*Y(:)-500e-6).^2;
200 % desired velocity profile:
201 % vabs(e) = linearly growing from 0 to vmax for 250<=e<=500
202 % vabs(e) = linearly decreasing to 0 for 500<=e<=750
203 % vabs(e) = zero otherwise
204 vabs = zeros(size(e));
206 i = find( se>= 250e-6 & se <= 500e-6);
207 vabs(i) = model.v_max / 250e-6 * (se(i)-250e-6);
208 i = find( se>= 500e-6 & se <= 750e-6);
209 vabs(i) = -model.v_max / 250e-6 * (se(i)-750e-6);
211 % compute velocity directions by gradient of elliptic field
212 dx_e = 2*(X(:)-500e-6);
213 dy_e = 5*(Y(:)-500e-6);
214 norm_grad_e_inv = (dx_e.^2 + dy_e.^2).^(-1);
215 % rotate by 90ΓΈ and normalization
216 Vx2 = - dy_e .* norm_grad_e_inv;
217 Vy2 = dx_e .*norm_grad_e_inv;
218 % scale to desired absolute velocity
222 flux_lin.Vx = (model.v_weight * Vx1 + (1-model.v_weight) * Vx2);
223 flux_lin.Vy = (model.v_weight * Vy1 + (1-model.v_weight) * Vy2);
224 elseif isequal(model.name_flux,'gdl
')
225 % load pressure from file and generate velocity field
226 p0 = load('gdl_pressure_v0000.mat
','p
','model
');
227 p1000 = load('gdl_pressure_v1000.mat
','p
','model
');
228 factor = 1000; % scaling: maximum inflow now 1000/factor m/s
229 [ux0, uy0] = evaluate_gradient(po0.model,X, Y,p0.p);
232 [ux1000, uy1000] = evaluate_gradient(p1000.model,X, Y,p1000.p);
233 ux1000 = -ux1000/factor;
234 uy1000 = -uy1000/factor;
235 % linear interpolation between velocity profiles yields exact solution
236 % of corresponding elliptic boundary value problem due to
237 % linearity of the problem.
238 ux = model.v_in_max * ux1000 + (1-model.v_in_max) * ux0;
239 uy = model.v_in_max * uy1000 + (1-model.v_in_max) * uy0;
241 % inverse of maximum absolute velocity value
242 flux_lin.lambda = max([ux(:); uy(:)])^-1;
243 % correct up to a constant factor
246 elseif isequal(model.name_flux,'gdl_pgrad_and_parabola
')
247 % load pressure from file and generate velocity field
248 p0 = load('gdl_pressure_dirichlet.mat
','p
','model
');
249 factor = 1000; % scaling: maximum inflow now 1000/factor m/s
250 [ux0, uy0] = evaluate_gradient(p0.model,X, Y,p0.p);
253 uypar = zeros(length(X),1);
254 uxpar = model.v_in_max * (model.yrange(2)-model.yrange(1))^-2 * 4 ...
255 * (Y(:)-model.yrange(1)) .* (model.yrange(2)-Y(:));
257 % linear interpolation between velocity profiles yields exact solution
258 % of corresponding elliptic boundary value problem due to
259 % linearity of the problem.
260 ux = model.v_in_max * uxpar + (1-model.v_in_max) * ux0;
261 uy = model.v_in_max * uypar + (1-model.v_in_max) * uy0;
263 % inverse of maximum absolute velocity value
264 flux_lin.lambda = max([ux(:); uy(:)])^-1; % correct up to a factor
267 elseif isequal(model.name_flux, 'transport
')
269 flux_lin = 1; % single component factor 1
270 elseif decomp_mode == 1 % components:
271 tmp.Vx = model.transport_x * ones(size(U));
272 tmp.Vy = model.transport_y * ones(size(U));
274 else % decomp_mode = 0, complete:
275 flux_lin.Vx = model.transport_x * ones(size(U));
276 flux_lin.Vy = model.transport_y * ones(size(U));
277 flux_lin.lambda = max([flux_lin.Vx(:);
278 flux_lin.Vy(:)])^-1; % correct up to a factor
281 respected_decomp_mode = 1;
282 elseif isequal(model.name_flux, 'richards
')
284 flux_lin.Vx = zeros(size(U));
285 flux_lin.Vy = zeros(size(U));
287 % p_mu = spline_select(model);
288 % [ breaks, coeffs, pieces, order ] = unmkpp(p_mu);
289 % p_mu_d = mkpp(breaks, coeffs(1:order-1) .* [order-1:-1:1]);
291 % denom = 1 + ppval(p_mu, X);
292 [res1, res2] = inv_geo_trans_derivative(model,X,Y,{[1],[2],[1,1],[1,2]},{[1],[2],[2,1],[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; % correct up to a factor
302 error('unknown flux!!
');
303 end; %%% end of flux select
308 % if not is linear flux: exact computation or forward difference
309 elseif isequal(model.name_flux,'burgers_parabola
')
310 % NOTE: KEEP THIS DATA IMPLEMENTATION CONSISTENT WITH ITS FLUX
312 P = [X(:)'-0.5; Y(:)
'-0.5];
313 Rinv = [cos(-model.vrot_angle), - sin(-model.vrot_angle); ...
314 sin(-model.vrot_angle), cos(-model.vrot_angle) ];
315 RP = Rinv*P; % rotated coordinates as columns
316 RP(1,:) = RP(1,:) + 0.5;
317 RP(2,:) = RP(2,:) + 0.5;
319 V(1,:) = (1- RP(2,:)) .* RP(2,:) * model.vmax * 4;
322 flux_lin.Vx = RV(1,:)' .* model.flux_pdeg .* U(:).^(model.flux_pdeg-1);
323 flux_lin.Vy = RV(2,:)
' .* model.flux_pdeg .* U(:).^(model.flux_pdeg-1);
324 % flux.lambda = 1/model.vmax;
326 else % perform forward difference
327 epsilon = max(abs(U)) * 1e-2;
328 if epsilon < eps; % U completely 0
332 fluxU = conv_flux(model,U,X,Y);
333 fluxU2 = conv_flux(model,U+epsilon,X,Y);
335 flux_lin.Vx = (fluxU2.Fx-fluxU.Fx)/epsilon;
336 flux_lin.Vy = (fluxU2.Fy-fluxU.Fy)/epsilon;
337 flux_lin.lambda = fluxU.lambda;
340 if decomp_mode>0 & respected_decomp_mode==0
341 error('function does not support affine decomposition!
');