quicwind/mass_cons_flux.m

   1 function [U,varargout]=mass_cons_flux(U0,X,method,check)
   2 % mass consistent flux approximation
   3 % input:
   4 %   U0         cell array, wind vectors u v w on staggered grids
   5 %   X          cell array, nodal coordinatees
   6 %   method     'direct' or 'pcg' method for solving system
   7 %   check      flag for optional timing and error checking
   8 % output:
   9 %   U          cell array, wind vectors u v w
  10 %   Lambda     pressure field (optional)
  11 %
  12 %
  13 % using:
  14 %   D = divergence of flux on cell sides
  15 %   M = transformation of wind to flux through cell sides with
  16 %       zero flux boundary condition at the bottom (the terrain)
  17 %   A = weights of the wind vector components, U'*A*U approx ||U||_L2^2
  18 %
  19 % compute:
  20 %
  21 %    min 1/2 (U - U0)'*A*(U - U0) s.t. D*M*U = 0
  22 %
  23 %    <=>    A*(U - U0) + M'*D'*Lambda = 0
  24 %
  25 %           D*M*U                     = 0
  26 %
  27 %    <=>    U = U0 - inv(A)*M'*D'*Lambda
  28 %
  29 %           D*M*inv(A)*M'*D'*Lambda = D*M*U0
  30
  31 if ~exist('check','var')
  32     check = false;
  33 end
  34
  35 if check
  36     % Timing start
  37     tstart=tic;
  38 end
  39
  40 % Creates the weight vector... diagonal for now?
  41 hvr = 1e4; % ratio of horizontal and vertical penalty: Sherman 1978 page 315
  42 s = cell_sizes(X);
  43 Avec = cell2vec({hvr*s.weight_u,hvr*s.weight_v,s.weight_w});
  44 Ainv = diag(sparse(1./Avec));
  45
  46 %method = 'iterative'
  47 switch method
  48     case 'direct'
  49         DM = mat_gen_wind_flux_div(X); % divergence of u0
  50         U0vec = cell2vec(U0);
  51         rhs = DM * U0vec;
  52         L = DM * Ainv * DM';
  53         Lambda = L \ rhs;
  54         Uvec = U0vec - Ainv * DM' * Lambda;
  55         U = vec2cell(Uvec,U0);
  56     case 'pcg'
  57         % uses pcg and implicit multiplication
  58         DM = mat_gen_wind_flux_div(X); % divergence of u0
  59         U0vec = cell2vec(U0);
  60 %         rhs = Mmul_v(X, 'n', U0vec);
  61 %         rhs = Dmul_v(X, 'n', rhs);
  62         rhs = DM * U0vec;
  63         L = DM * Ainv * DM';
  64 %        PC_L = ichol(PC,struct('michol','on'));
  65 %        lhs_mat = @(x) lhs_apply(X, Ainv, x);
  66 %        Lambda = pcg(lhs_mat, rhs, 1e-10, 100, PC_L, PC_L');
  67 %        Lambda = pcg(lhs_mat, rhs, 1e-8, 400);
  68 %        [Lambda,FLAG,RELRES,ITER,RESVEC] = pcg(lhs_mat, rhs, 1e-6, 400);
  69         [Lambda,FLAG,RELRES,ITER,RESVEC] = pcg(L, rhs, 1e-6, 400);
  70         figure;
  71         semilogy(RESVEC)
  72         title('Residuals using PCG')
  73         xlabel('Iterations'),ylabel('Residuals')
  74         Uvec = U0vec - Ainv * DM' * Lambda;
  75         U = vec2cell(Uvec,U0);
  76     otherwise
  77         error('only direct and iterative methods implemented')
  78 end
  79
  80 if check
  81     % Check that the divergence of the computed U is 0
  82     err_div = big(div3(wind2flux(U,X)))
  83     varargout{2} = err_div;
  84     fprintf('mass_cons_flux %g seconds\n',toc(tstart))
  85 end
  86 if length(varargout) > 0
  87     varargout{1} = Lambda;
  88 end
  89 % check no correction in normal speed at the bottom
  90 % if check, err_corr_w_bottom = big(u{3}(:,:,1)-U0{3}(:,:,1)), end
  91 end