quicwind/mass_cons_flux_test.m

   1 function err = mass_cons_flux_test
   2 disp('mass_cons_flux_test - testing mass consistent flux approximation with terrain')
   3
   4 %% Settings
   5 % options:
   6 %   'constant' : constant wind to x direction,
   7 %   'log' : vertical log profile,
   8 init_wind = 'log';
   9 terrain_following = true; % still under developmentssss
  10 u0 = 1; % initial wind magnitude
  11 h0 = .5; % roughness high: only for log profile cases ('log' or 'terrain')
  12 mesh_len=[20,20,20];
  13 mesh_disp=[20,20,20];
  14 h=ones(1,3);
  15 X = regular_mesh(mesh_len,h,1.2);
  16 X = add_terrain_to_mesh(X,'hill','squash',0.05);
  17 x = X{1};
  18 y = X{2};
  19 z = X{3};
  20 [nx,ny,nz] = size(X{1});
  21
  22 %% Initial wind
  23 % declare initial wind
  24 U0 = grad3z(ones(size(x)-1),mesh_len);
  25 switch init_wind
  26     case 'constant'
  27         U0{1} = u0*ones(size(U0{1}));
  28         U0{2} = 0*ones(size(U0{2}));
  29         U0{3} = 0*ones(size(U0{3}));
  30     case 'log'
  31         hz = z-z(:,:,1);
  32         h_at_u = (hz(:,1:end-1,1:end-1)+hz(:,2:end,1:end-1)+hz(:,1:end-1,2:end)+hz(:,2:end,2:end))/4;
  33         mlogh = log(max(h_at_u(:)));
  34         U0{1} = u0*log(max(h_at_u/h0,1))/mlogh;
  35         U0{2} = 0*ones(size(U0{2}));
  36         U0{3} = 0*ones(size(U0{3}));
  37     otherwise
  38         error('mass_cons_flux_test: unknown init_wind')
  39 end
  40 plot_wind_above_terrain(U0,X,mesh_disp)
  41 if terrain_following
  42     U0 = terrain_correction(U0,X);
  43     plot_wind_above_terrain(U0,X,mesh_disp)
  44 end
  45
  46 %% Direct solve
  47 [U,Lambda_d] = mass_cons_flux(U0,X,'direct','check');
  48 plot_wind_above_terrain(U,X,mesh_disp)
  49
  50 %% PCG solve
  51 [V,Lambda_pcg,err] = mass_cons_flux(U0,X,'pcg','check');
  52 plot_wind_above_terrain(V,X,mesh_disp)
  53
  54 %% Test results
  55 direct_vs_pcg = big(cell2vec(U)-cell2vec(V))
  56 lambda_diff = big(Lambda_d-Lambda_pcg)
  57
  58 end