rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
conv_flux_linearization.m
1 function flux_lin = conv_flux_linearization(model,U,X,Y)
2 %function flux_lin = conv_flux_linearization(model,[U],[X],[Y])
3 %
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.
12 %
13 % For a linear flux in conv_flux the following is identical
14 %
15 % (Fx(U),Fy(U)) = (Vx,Vy) * U
16 %
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.
37 %
38 % required fields of model as desired by conv_flux
39 %
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
43 %
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
55 %
56 % In 'coefficient' mode, the parameters in brackets are empty
57 
58 % Bernard Haasdonk 21.4.2006
59 
60 flux_lin = [];
61 
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;
66 
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'};
70 
71 %if ismember(model.name_flux, linear_fluxes)
72 if model.flux_linear
73  % check if flux is to be computed or read from file
74  if model.use_velocity_matrixfile==1
75 
76  % affine decomposition is trivial: no dependency on mu
77 
78  if decomp_mode < 2 % i.e. 'none' or 'components'
79 
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
86 
87  if ~isfield(model,'velocity_matrixfile')
88  error('field velocity_matrixfile not set!!');
89  end;
90 
91  fullfn = fullfile(rbmatlabhome,'datafunc','data',...
92  model.velocity_matrixfile);
93 
94  if ~exist(fullfn,'file') & ~cache('exist',fullfn)
95  error(['warning: velocity_matrixfile not existing. ',...
96  'Call gen_velocity_matrixfile with suitable' ...
97  ' divclean_mode.']);
98  end;
99 
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;
104  end;
105 
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', ...
113  fullfn,X,Y);
114  else % mode 2, i.e assume, that velocityfilename is set correctly to
115  % requested points
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
120  X = X(:); Y = Y(:);
121  if ~isequal(X,V.X(:)) | ~isequal(Y,V.Y(:))
122  error('matrixfile does not fit to requested points!');
123  end;
124  i = 1:length(X);
125  flux_lin.Vx = V.Vx;
126  flux_lin.Vy = V.Vy;
127  %NaN*ones(size(X));
128  % flux_lin.Vy = NaN*ones(size(X));
129  % flux_lin.Vx(i) = V.Vx(j);
130  % flux_lin.Vy(i) = V.Vy(j);
131  lambda = V.lambda;
132  end;
133 
134  if decomp_mode == 0
135  flux_lin.lambda = lambda;
136  end;
137  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of loading flux matrix.
138  if decomp_mode == 1 % put single matrix into cell array
139  flux_lin = {flux_lin};
140  end;
141  else % decomp_mode = 2;
142  flux_lin = 1; % factor one for single component
143  end;
144  respected_decomp_mode = 1;
145 
146  else % model.use_velocity_matrixfile == 0
147  %%%%%%%%%%% compute the fluxes instead of loading matrix from file
148 
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
152  % f'(u) = v(x,y)
153  % e.g. lambda := 1/sup|v(x,y)|
154  if decomp_mode == 0
155  flux_lin.lambda = 1/(abs(model.c)/4.0+ 1e-10); % some eps for divby0
156  end;
157 
158  % flux.precomputable = 1; % velocity field can be computed in advance
159  % flux values : velocity field in the edge midpoints
160  if decomp_mode==2
161  flux_lin = model.c;
162  elseif decomp_mode == 1
163  dummy.Vy = zeros(length(X),1);
164  dummy.Vx = Y(:) .* (1-Y(:));
165  flux_lin = {dummy};
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(:));
169  else
170  error('unknown decmposition mode!!');
171  end;
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'),...
176  'p','model');
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
188  % f'(u) = v(x,y)
189  % e.g. lambda := 1/sup|v(x,y)|
190 
191  flux_lin.lambda = 1/abs(model.v_max); % some eps for preventing divby0
192 
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);
196 
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));
204  se = sqrt(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);
209 
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
218  Vx2 = Vx2 .* vabs;
219  Vy2 = Vy2 .* vabs;
220 
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);
229  ux0 = -ux0/factor;
230  uy0 = -uy0/factor;
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;
239 
240  % inverse of maximum absolute velocity value
241  flux_lin.lambda = max([ux(:); uy(:)])^-1;
242  % correct up to a constant factor
243  flux_lin.Vx = ux;
244  flux_lin.Vy = uy;
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);
250  ux0 = -ux0/factor;
251  uy0 = -uy0/factor;
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(:));
255 
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;
261 
262  % inverse of maximum absolute velocity value
263  flux_lin.lambda = max([ux(:); uy(:)])^-1; % correct up to a factor
264  flux_lin.Vx = ux;
265  flux_lin.Vy = uy;
266  elseif isequal(model.name_flux, 'transport')
267  if decomp_mode == 2
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));
272  flux_lin = {tmp};
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
278 
279  end;
280  respected_decomp_mode = 1;
281  elseif isequal(model.name_flux, 'richards')
282 
283  flux_lin.Vx = zeros(size(U));
284  flux_lin.Vy = zeros(size(U));
285 
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]);
289 
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
302  else
303  error('unknown flux!!');
304  end; %%% end of flux select
305  end;
306 
307  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308  % Non-linear fluxes %
309  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
310  %
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
314  % IN conv_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;
321  V = zeros(size(RP));
322  V(1,:) = (1- RP(2,:)) .* RP(2,:) * model.vmax * 4;
323  V(2,:) = 0;
324  RV = Rinv^(-1) * V;
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;
328 
329 else % perform forward difference
330  epsilon = max(abs(U)) * 1e-2;
331  if epsilon < eps; % U completely 0
332  epsilon = 1e-2;
333  end;
334 
335  fluxU = conv_flux(model,U,X,Y);
336  fluxU2 = conv_flux(model,U+epsilon,X,Y);
337 
338  flux_lin.Vx = (fluxU2.Fx-fluxU.Fx)/epsilon;
339  flux_lin.Vy = (fluxU2.Fy-fluxU.Fy)/epsilon;
340  flux_lin.lambda = fluxU.lambda;
341 end;
342 
343 if decomp_mode>0 & respected_decomp_mode==0
344  error('function does not support affine decomposition!');
345 end;
346 
347 
348 %| \docupdate