quicwind/mat_mul_v_test.m

   1 function err = mat_mul_v_test
   2 % test Mmul_v_test against mat_gen_wind_flux_div
   3
   4 disp('mat_mul_v_test')
   5
   6 mesh_len=[3,3,3];
   7 h=rand(1,3);
   8 nx = mesh_len(1);
   9 ny = mesh_len(2);
  10 nz = mesh_len(3);
  11
  12 % Test mesh with a hill terrain (unit vectors for now)
  13 X = regular_mesh(mesh_len, h, 1.0);
  14 X = add_terrain_to_mesh(X,'hill','squash',0.4);
  15
  16 % Initialize inputs for testing
  17 testDT_vec = rand(prod(mesh_len),1);
  18 testD_vec = rand((nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1),1);
  19 testA_vec = rand((nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1),1);
  20 testMT_vec = rand((nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1),1);
  21 testM_vec = rand((nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1),1);
  22
  23 % Obtain sparse matrices for testing
  24 D = mat_gen_wind_flux_div(X,'D');
  25 M = mat_gen_wind_flux_div(X,'M');
  26 A = sparse(diag(rand((nx+1)*ny*nz + nx*(ny+1)*nz + nx*ny*(nz+1),1)));
  27
  28 % Generate output vectors
  29 resDT_mul = transpose(D) * testDT_vec;
  30 resDT_loops = Dmul_v(X,'t',testDT_vec);
  31
  32 resD_mul = D * testD_vec;
  33 resD_loops = Dmul_v(X,'n',testD_vec);
  34
  35 resA_mul = A * testA_vec;
  36 resA_loops = Amul_v(X,A,testA_vec);
  37
  38 resM_mul = M * testM_vec;
  39 resM_loops = Mmul_v(X,'n',testM_vec);
  40
  41 resMT_mul = transpose(M) * testMT_vec;
  42 resMT_loops = Mmul_v(X,'t',testMT_vec);
  43
  44 disp('Compare M*v results')
  45 err_Mv = big(resM_mul - resM_loops)
  46
  47 disp('Compare M^T*v results')
  48 err_MTv = big(resMT_mul - resMT_loops)
  49
  50 disp('Compare A*v results')
  51 err_Av = big(resA_mul - resA_loops)
  52
  53 disp('Compare D*v results')
  54 err_Dv = big(resD_mul - resD_loops)
  55
  56 disp('Compare D^T*v results')
  57 err_DTv = big(resDT_mul - resDT_loops)
  58
  59 errs = [err_Mv,err_MTv,err_Av,err_Dv,err_DTv];
  60
  61 err = max(errs);
  62
  63
  64
  65
  66