quicwind/mat_wind_flux_div_test.m

   1 function err = mat_wind_flux_div_test
   2 disp('mat_wind_flux_div_test')
   3 % test mesh with a hill terrain
   4 % X = uniform_mesh([2,2,2],[.5,.5,.5]);
   5 X = regular_mesh([3,4,5],[0.1,0.1,0.5],1.2);
   6 X = add_terrain_to_mesh(X,'hill','squash',0.4);
   7 % matrix of wind flux divergence
   8 DM = mat_wind_flux_div(X);
   9 % random test wind
  10 U = grad3z(rand(size(X{1})-1),[1 1 1]);
  11 % its flux divergence on the mesh
  12 L = div3(wind2flux(U,X));
  13 disp('compare with matrix-vector multiplication')
  14 Uvec = cell2vec(U);
  15 Lvec = L(:);
  16 err1 = big(DM * Uvec - Lvec);
  17
  18 wind_template=grad3z(rand(size(X{1})-1),[1 1 1]);  % cell matrix with entries size of u,v,w
  19 n = sum(cellfun(@numel,wind_template));   % size of vector this will act on
  20
  21 disp('compare M^T with matrix-vector multiplication')
  22 M = mat_wind_flux_div(X,'M');
  23 MTu = wind2flux_trans(U,X);
  24 err2 = big(M'*Uvec-cell2vec(MTu));
  25
  26 err = max(err1,err2);
  27
  28