rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
velocity_to_be_distributed_in_single_files.m
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])
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 affine_decomp_mode=='none', 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 affine_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 % 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
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 = 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;
67 
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'};
71 
72 %if ismember(model.name_flux, linear_fluxes)
73 if model.flux_linear
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)
77 
78  % affine decomposition is trivial: no dependency on mu
79 
80  if decomp_mode < 2 % i.e. 'none' or 'components'
81 
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
88 
89  if ~isfield(model,'velocity_matrixfile')
90  error('field velocity_matrixfile not set!!');
91  end;
92 
93  fullfn = fullfile(rbmatlabhome,'datafunc','mat',...
94  model.velocity_matrixfile);
95 
96  if ~exist(fullfn,'file') & ~cache('exist',fullfn)
97  error(['warning: velocity_matrixfile not existing. ',...
98  'Call gen_velocity_matrixfile with suitable' ...
99  ' divclean_mode.']);
100  end;
101 
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;
106  end;
107 
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', ...
115  fullfn,X,Y);
116  else % mode 2, i.e assume, that velocityfilename is set correctly to
117  % requested points
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
122  X = X(:); Y = Y(:);
123  if ~isequal(X,V.X(:)) | ~isequal(Y,V.Y(:))
124  error('matrixfile does not fit to requested points!');
125  end;
126  i = 1:length(X);
127  flux_lin.Vx = V.Vx;
128  flux_lin.Vy = V.Vy;
129  %NaN*ones(size(X));
130 % flux_lin.Vy = NaN*ones(size(X));
131 % flux_lin.Vx(i) = V.Vx(j);
132 % flux_lin.Vy(i) = V.Vy(j);
133  lambda = V.lambda;
134  end;
135 
136  if decomp_mode == 0
137  flux_lin.lambda = lambda;
138  end;
139  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of loading flux matrix.
140  if decomp_mode == 1 % put single matrix into cell array
141  flux_lin = {flux_lin};
142  end;
143  else % decomp_mode = 2;
144  flux_lin = 1; % factor one for single component
145  end;
146  respected_decomp_mode = 1;
147 
148  else %%%%%%%%%%% compute the fluxes instead of loading matrix from file
149 
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
153  % f'(u) = v(x,y)
154  % e.g. lambda := 1/sup|v(x,y)|
155  if decomp_mode == 0
156  flux_lin.lambda = 1/(abs(model.c)/4.0+ 1e-10); % some eps for divby0
157  end;
158 
159  % flux.precomputable = 1; % velocity field can be computed in advance
160  % flux values : velocity field in the edge midpoints
161  if decomp_mode==2
162  flux_lin = model.c;
163  elseif decomp_mode == 1
164  dummy.Vy = zeros(length(X),1);
165  dummy.Vx = Y(:) .* (1-Y(:));
166  flux_lin = {dummy};
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(:));
170  else
171  error('unknown decmposition mode!!');
172  end;
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'),...
177  'p','model');
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
189  % f'(u) = v(x,y)
190  % e.g. lambda := 1/sup|v(x,y)|
191 
192  flux_lin.lambda = 1/abs(model.v_max); % some eps for preventing divby0
193 
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);
197 
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));
205  se = sqrt(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);
210 
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
219  Vx2 = Vx2 .* vabs;
220  Vy2 = Vy2 .* vabs;
221 
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);
230  ux0 = -ux0/factor;
231  uy0 = -uy0/factor;
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;
240 
241  % inverse of maximum absolute velocity value
242  flux_lin.lambda = max([ux(:); uy(:)])^-1;
243  % correct up to a constant factor
244  flux_lin.Vx = ux;
245  flux_lin.Vy = uy;
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);
251  ux0 = -ux0/factor;
252  uy0 = -uy0/factor;
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(:));
256 
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;
262 
263  % inverse of maximum absolute velocity value
264  flux_lin.lambda = max([ux(:); uy(:)])^-1; % correct up to a factor
265  flux_lin.Vx = ux;
266  flux_lin.Vy = uy;
267  elseif isequal(model.name_flux, 'transport')
268  if decomp_mode == 2
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));
273  flux_lin = {tmp};
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
279 
280  end;
281  respected_decomp_mode = 1;
282  elseif isequal(model.name_flux, 'richards')
283 
284  flux_lin.Vx = zeros(size(U));
285  flux_lin.Vy = zeros(size(U));
286 
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]);
290 
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
301  else
302  error('unknown flux!!');
303  end; %%% end of flux select
304  end;
305 
306 
307 
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
311  % IN conv_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;
318  V = zeros(size(RP));
319  V(1,:) = (1- RP(2,:)) .* RP(2,:) * model.vmax * 4;
320  V(2,:) = 0;
321  RV = Rinv^(-1) * V;
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;
325 
326 else % perform forward difference
327  epsilon = max(abs(U)) * 1e-2;
328  if epsilon < eps; % U completely 0
329  epsilon = 1e-2;
330  end;
331 
332  fluxU = conv_flux(model,U,X,Y);
333  fluxU2 = conv_flux(model,U+epsilon,X,Y);
334 
335  flux_lin.Vx = (fluxU2.Fx-fluxU.Fx)/epsilon;
336  flux_lin.Vy = (fluxU2.Fy-fluxU.Fy)/epsilon;
337  flux_lin.lambda = fluxU.lambda;
338 end;
339 
340 if decomp_mode>0 & respected_decomp_mode==0
341  error('function does not support affine decomposition!');
342 end;
343 
344 
345 %| \docupdate