femwind/fortran/ndt_mult_test.f90

   1 program ndt_mult_test
   2
   3 use module_ndt_mult   ! testing only
   4 use module_utils ! to read and write matrices as text files from matlab
   5
   6 implicit none
   7
   8 real, pointer:: kmat(:,:,:,:), u(:,:,:), y(:,:,:), &  ! fortran is not case sensitive
   9                 kmat_m(:,:,:,:), u_m(:,:,:), y_m(:,:,:),  r(:,:,:)
  10 real, pointer::a(:)
  11 integer :: s(4),n(3)
  12
  13 integer :: msize, &
  14     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  15     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  16     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  17     ifts, ifte, kfts, kfte, jfts,jfte                            ! fire tile bounds
  18 integer :: i,j,k,jx
  19 real:: siz,relres
  20
  21 ! read input arrays in ikj index ordering and tight bounds
  22 call read_array_nd(a,s,'kmat')
  23 allocate(kmat_m(s(1),s(2),s(3),s(4)))
  24 kmat_m = reshape(a,s)
  25 call read_array_nd(a,n,'u')
  26 allocate(u_m(n(1),n(2),n(3)))
  27 u_m = reshape(a,n)
  28
  29 if (s(1).ne.n(1).or.s(2).ne.n(2).or.s(3).ne.n(3))call crash('ndt_mult_test: inconsistent size kmat and u')
  30
  31 ifts = 1
  32 ifte = n(1) - 1
  33 kfts = 1
  34 kfte = n(2) - 1
  35 jfts = 1
  36 jfte = n(3) - 1
  37
  38 ifds = ifts
  39 ifde = ifte
  40 jfds = jfts
  41 jfde = jfte
  42 kfds = kfts
  43 kfde = kfte
  44
  45 msize = s(4)
  46 if(msize.ne.14)call crash('msize must be 14')
  47 ifms = ifts-1
  48 ifme = ifte+2
  49 kfms = kfts-1
  50 kfme = kfte+2
  51 jfms = jfts-1
  52 jfme = jfte+2
  53
  54 ! allocate a little bigger with zeros in extra areas
  55 allocate(Kmat(ifms:ifme,kfms:kfme,jfms:jfme,1:msize))
  56 allocate(   u(ifms:ifme,kfms:kfme,jfms:jfme))
  57 allocate(   y(ifms:ifme,kfms:kfme,jfms:jfme))
  58 allocate(   r(ifms:ifme,kfms:kfme,jfms:jfme))
  59 kmat = 0.
  60 u = 0.
  61 y = 0.
  62
  63 print *,'copying the input data'
  64 Kmat(1:s(1),1:s(2),1:s(3),1:s(4))=Kmat_m(1:s(1),1:s(2),1:s(3),1:s(4))
  65 u(1:s(1),1:s(2),1:s(3))=u_m(1:s(1),1:s(2),1:s(3))
  66
  67 write(*,'(a)')'ntd_mult now computing r = y - Kmat*u, calling with y=0'
  68 call ndt_mult(  &
  69   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  70   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  71   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  72   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
  73   kmat, u, y, r, siz, relres)
  74
  75 write(*,'(a,3i8)')'copying the resulting -y=K*u to array size ',n
  76 allocate(y_m(1:n(1),1:n(2),1:n(3)))
  77 y_m(1:s(1),1:s(2),1:s(3))=-r(1:s(1),1:s(2),1:s(3))
  78
  79 call write_array_nd(reshape(y_m,(/product(n)/)),n,'y')
  80
  81 end program ndt_mult_test