Merge branch 'fixf'
[wrf-fire-matlab.git] / femwind / ndt_test.m
blob9263cce5939ebb6c124fb360d881ff630a1963e5
1 format compact
2 % test for ndt_assembly and ndt_mult
3 nel = [5,4,3]
4 %nel=[1 1 1]
5 n=nel+1
6 h = [1,1,1]
7 expand=1.3 
8 A = diag([1 1 1])
9 lambda=[]
10 params=[]
11 u0=[];
13 X = regular_mesh(nel,h,expand);
14 % X = add_terrain_to_mesh(X, 'hill', 'shift', 0.1)
15 X = add_terrain_to_mesh(X, 'hill', 'squash', 0.1)  % more thorough testing
17 plot_mesh_3d(X)
19 K=nd_assembly(A,X,lambda,params);
20 [n1,n2,n3,m1,m2,m3]=size(K);
21 K27 = reshape(K,n1,n2,n3,m1*m2*m3);  % make to n x 27 nd format
23 disp('converting to reduced storage format with 14 numbers per row')
24 K14 = ndt_convert(K27,14); 
25 disp('multiplying st 14 by all ones')
26 x1=ones(n1,n2,n3);
27 y1=ndt_mult(K14,x1,3);
28 K14_err_zero=big(y1)  % should be zero
29 if abs(K14_err_zero)>1e-10, error('should be zero'),end
31 disp('testing multiplication by ones should give zero')
32 x1=ones(n1,n2,n3);
33 y27=ndt_mult(K27,x1);
34 K27_err_zero=big(y27)  % should be zero
37 disp('compare matrix vector multiply nd and ndt 27 numbers per row')
38 xr=rand(n1,n2,n3);
39 y27=ndt_mult(K27,xr);
40 yn=nd_mult(K,xr);
41 Kn_27_err=big(yn-y27)
43 y14=ndt_mult(K14,xr);
44 K14_err_rand=big(y14-y27)  % should be zero
45 if abs(K14_err_rand)>1e-10, error('should be zero'),end
47 disp('convert to sparse compare matrix-vector multiply')
48 Ks = ndt_convert(K,'sparse');  
49 ys=Ks*xr(:);
50 K_Ks_rand_err=big(ys(:)-y27(:))
51 if abs(K_Ks_rand_err)>1e-10, error('should be zero'),end
53 if exist('fortran/ndt_boundary_conditions_test.exe')
54     disp('testing if ndt_boundary_conditions gives same result in fortran')
55     ndt_boundary_conditions_fortran(K14);
56 else
57     warning('fortran/ndt_boundary_conditions_test.exe not available')
58 end
60 disp('testing multiply vec by st 14 vs st 27')
62 % test same results for ndt_mult from matlab and fortran
63 if exist('fortran/ndt_mult_test.exe')
64     disp('testing if ndt_mult same result in fortran')
65     err=ndt_mult_fortran(K14,xr);
66     if abs(err)<1e-6
67         fprintf('error %g OK\n',err)
68     else
69         error(sprintf('error %g too large',err))
70     end
71 else
72     warning('fortran/ndt_mult_test.exe not available')
73 end
75 disp('sparse assembly vs ndt_assembly test')
76 K_sparse=sparse_assembly(A,X,u0,lambda,params);
77 err_mat_sparse=big(Ks-K_sparse);
78 if abs(err_mat_sparse)>1e-10, error('should be zero'),end
80 K27a = ndt_assembly(A,X,u0,lambda,params,27);
81 K27a_err = big(K27 - K27a)
82 if abs(K27a_err)>1e-10, error('should be zero'),end
85 K14a = ndt_assembly(A,X,u0,lambda,params,14);
86 disp('multiplying ndt_assembly st 14 by all ones')
87 x=ones(n1,n2,n3);
88 y1=ndt_mult(K14,x1);
89 K14a_err_zero=big(y1)  % should be zero
90 if K14a_err_zero > 1e-6
91     error('residual should be zero')
92 end
94 disp('comparing st 14 with converged st 27')
95 K14a_err = big(K14 - K14a)
96 if abs(K14a_err)>1e-10, error('should be zero'),end
98 disp('ndt_assembly_fortran')
99 err = ndt_assembly_fortran(A,X,u0,lambda,params,14)
100 if err > 1e-6
101     error('ndt_assembly_fortran error too large')