quicwind/mass_cons_int.m

   1 function [u,varargout]=mass_cons_int(u0,h,w,check)
   2 % mass consistent approximation
   3 % given
   4 %   u0  3 wind vectors on staggered grids
   5 %   h   step, vector length 3
   6 %   w   weights, vector length 3
   7
   8 % solve
   9 % min_u 1/2 <D(u-u0),u-u0> s.t. div u = 0, D=diag[w(1)I,w(2)I,w(3)I]
  10 % <=> min_u max_lambda 1/2 sum_i  <D(u-u0),u-u0> + <lambda,div u>
  11 % (using <lambda,div u> = - <grad lambda, u>)
  12 % <=> D(u-u0) - grad lambda = 0, div u = 0
  13 % <=> inv(D) grad lambda = u - u0, div u = 0
  14 % <=> div inv(D) grad lambda = div u - div u0
  15 % <=>
  16 % ______________________________________________________
  17 %    - div inv(D) grad lambda = div u0,
  18 %     u = u0 + inv(D) grad lambda
  19 % with boundary condidions lambda=0 on the sides and top
  20 % and dlambda/dz = 0 on the bottom
  21 %-------------------------------------------------------
  22
  23 if ~exist('check','var')
  24     check = false;
  25 end
  26 tstart=tic;
  27 % divergence of u0
  28 f = div3(u0,h);
  29 % reflect the right hand side abound the bottom in 3rd coordinate
  30 n = size(f);
  31 fprintf('mass_cons_int mesh size %i %i %i\n',n)
  32 r = zeros(n(1),n(2),2*n(3));
  33 r(:,:,1:n(3))=f(:,:,n(3):-1:1);
  34 r(:,:,n(3)+1:2*n(3))=f;
  35 % solve with zero boundary conditions
  36 lambda=poisson_fft3z(r,h,1./w);
  37 % check symmetry about the bottom
  38 if check, err_sym = big(lambda(:,:,n(3):-1:1)-lambda(:,:,n(3)+1:2*n(3))), end
  39 % extract the upper half of the solution
  40 lambda = lambda(:,:,n(3)+1:2*n(3));
  41 % grad with reflection about the bottom
  42 g = grad3z(lambda,h,true);
  43 % update u
  44 for i=1:3
  45     u{i} = u0{i} + (1/w(i))*g{i};
  46 end
  47 if check,
  48     % check divergence
  49     err_div = big(div3(u,h))
  50     % check no correction in vertical speed at the bottom
  51     err_corr_w_bottom = big(u{3}(:,:,1)-u0{3}(:,:,1))
  52     varargout={max(err_div,err_corr_w_bottom)};
  53 end
  54
  55 fprintf('mass_cons_int %g seconds\n',toc(tstart))
  56 end