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 jfts = 1
  34 jfte = n(2) - 1
  35 kfts = 1
  36 kfte = 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+1
  49 jfms = jfts-1
  50 jfme = jfte+1
  51 kfms = kfts-1
  52 kfme = kfte+1
  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 ! copy the input data
  64 do j=jfts,jfte
  65   do k=kfts,kfte
  66     do i=ifts,ifte
  67       do jx = 1,msize
  68         kmat(i,k,j,jx) = kmat_m(i,j,k,jx)
  69       enddo
  70       u(i,k,j)=u_m(i,j,k)
  71     enddo
  72   enddo
  73 enddo
  74
  75 write(*,'(a)')'ntd_mult now computing r = y - Kmat*u, calling with y=0'
  76 call ndt_mult(  &
  77   ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
  78   ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
  79   ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
  80   ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
  81   kmat, u, y, r, siz, relres)
  82
  83 write(*,'(a,3i8)')'copying -y to array size ',n
  84 allocate(y_m(1:n(1),1:n(2),1:n(3)))
  85 do j=jfts,jfte
  86   do k=kfts,kfte
  87     do i=ifts,ifte
  88       y_m(i,j,k)=-r(i,k,j)
  89     enddo
  90   enddo
  91 enddo
  92
  93 call write_array_nd(reshape(y_m,(/product(n)/)),n,'y')
  94
  95 end program ndt_mult_test