quicwind/wind2flux_test.m

   1 function errmax=wind2flux_test
   2 nx=50; ny=30; nz=10;
   3 h=[rand,rand,1];
   4 hh=rand(1,3);
   5 % corner nodes
   6 X = uniform_mesh([nx,ny,nz],hh);
   7 % gradient for sizing
   8 Usize = grad3z(zeros(size(X{1})-1),h,1);
   9
  10 disp('constant wind')
  11 c=[rand,rand,rand];
  12 for i=1:3,
  13     Uconst{i} = 0*Usize{i}+c(i);
  14 end
  15
  16 errmax=0
  17 testing_wind(Uconst)
  18 fprintf('mesh size %g %g %g max error %g\n',nx,ny,nz,errmax)
  19
  20 function testing_wind(U)
  21
  22     fl=wind2flux(U,X);
  23     d=div3(fl);
  24     disp('divergence zero except at the bottom')
  25     err=big(d(:,:,2:end))
  26
  27     maxh = max(X{3}(:))-hh(3);
  28
  29     disp('terrain slope in x direction')
  30     disp('divergence zero except at the bottom')
  31     thx = 0.1*hh(1)*[0:nx]'*ones(1,ny+1);
  32     thx = thx/max(thx(:))*maxh;
  33     test_terrain(thx)
  34
  35     disp('terrain slope in y direction')
  36     disp('divergence zero except at the bottom')
  37     thy = 0.1*hh(2)*ones(nx+1,1)*[0:ny];
  38     thy = thy/max(thy(:))*maxh;
  39     test_terrain(thy)
  40
  41     disp('terrain slope in random constant direction')
  42     thxy = rand*thx + rand*thy;
  43     thxy = thxy/max(thxy(:))*maxh;
  44     test_terrain(thxy)
  45
  46     disp('roof slope in x direction')
  47     half = floor(nx/2);
  48     xs=hh(1)*[0:half,half-1:-1:2*half-nx];
  49     ys=ones(1,ny+1);
  50     th = 0.1*xs'*ys;
  51
  52     test_terrain(th)
  53
  54     disp('roof slope in y direction')
  55     half = floor(ny/2);
  56     xs=ones(1,nx+1);
  57     ys=hh(2)*[0:half,half-1:-1:2*half-ny];
  58     th = 0.1*xs'*ys;
  59     th = th/max(th(:))*maxh;
  60     test_terrain(th)
  61
  62     disp('pyramid roof slope')
  63     half = floor(ny/2);
  64     xs=[0:nx]; ys=[0:ny];
  65     [ii,jj]=ndgrid(xs,ys);
  66     th = nx+ny-abs(ii-nx/2)-abs(jj-ny/2);
  67     th = th/max(th(:))*maxh;
  68     test_terrain(th)
  69
  70     disp('random terrain')
  71     th = min(hh)*0.1*rand(size(X{1}(:,:,1)));
  72     th = th/max(th(:))*maxh;
  73     test_terrain(th)
  74
  75     function test_terrain(t)
  76         XX=add_terrain_to_mesh(X, t, 'shift');
  77         test_divergence
  78         XX=add_terrain_to_mesh(X, t, 'compress');
  79         test_divergence
  80
  81         function test_divergence
  82             fl=wind2flux(U,XX);
  83             d=div3(fl);
  84             disp('divergence zero except at the bottom')
  85             err=big(d(:,:,2:end));
  86             fprintf('err=%g\n',err)
  87             errmax=max(errmax,err)
  88         end
  89     end
  90
  91 end
  92
  93 end
  94
  95
  96
  97
  98