quicwind/masscons.m

   1 function u=masscons(u0,h,w)
   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 % <=>
  11 % min_u max_lambda 1/2 sum_i  <D(u-u0),u-u0> + <lambda,div u>
  12 % <=>     (using <lambda,div u> = - <grad lambda, u>)
  13 %    D(u-u0) - grad lambda = 0
  14 %    div u                 = 0
  15 % <=>
  16 %    - div inv(D) grad lambda = u0
  17 %    u = u0 + inv(D) grad lambda