femwind/fortran/module_sweeps.f90

   1 module module_sweeps
   2 use module_utils
   3
   4 contains
   5
   6 subroutine sweeps(ifds, ifde, kfds,kfde, jfds, jfde,               &     ! fire grid dimensions
   7                   ifms, ifme, kfms,kfme, jfms, jfme,               &
   8                   ifps, ifpe, kfps, kfpe, jfps, jfpe,              &     ! fire patch bounds
   9                   ifts, ifte, kfts, kfte, jfts,jfte,               &
  10                   Kmat, F, x, siz, reldif)                             !input and output matrices/vectors
  11
  12 implicit none
  13
  14 !*** arguments
  15
  16 integer, intent(in)::                                              &
  17     ifds, ifde, kfds, kfde, jfds, jfde,                            &     ! fire domain bounds
  18     ifms, ifme, kfms, kfme, jfms, jfme,                            &     ! fire memory bounds
  19     ifps, ifpe, kfps, kfpe, jfps, jfpe,                            &     ! fire patch bounds
  20     ifts, ifte, kfts, kfte, jfts,jfte                                    ! fire tile bounds
  21
  22
  23 integer, parameter:: msize = 14
  24 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme,msize):: Kmat  ! global stiffness matrix
  25 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: F           ! global load vector
  26 real,intent(out), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: x           ! output vector
  27 real, intent(out):: siz, reldif    ! size of input, relative difference betweet input and output, in max norm
  28
  29 !*** local
  30
  31 integer:: i, j, k, r1, r2, ie,je, ke
  32 real:: t, dif, old, new
  33
  34 !*** executable
  35
  36 !call write_array(F(ifts: ifte, kfts: kfte, jfts:jfte),'F_sweeps_in')
  37 !call write_array(x(ifts: ifte, kfts: kfte, jfts:jfte),'x_sweeps_in')
  38
  39 dif = 0.
  40 siz = 0.
  41 ie = snode(ifte, ifde, +1)
  42 je = snode(jfte, jfde, +1)
  43 ke = snode(kfte, kfde, +1)
  44
  45 do r1 = 1,2
  46     do r2 = 1,2
  47         do j =r2,je,2
  48             do k = 1,ke
  49                do i = r1,ie,2
  50                   old = x(i,k,j)
  51                   new =     (1/Kmat(i  ,k  ,j  , 1))*               &
  52                             ( F(i,k,j) -                            &
  53                             (                                       &
  54                              Kmat(i-1,k-1,j-1,14)*x(i-1,k-1,j-1) +  &
  55                              Kmat(i  ,k-1,j-1,13)*x(i  ,k-1,j-1) +  &
  56                              Kmat(i+1,k-1,j-1,12)*x(i+1,k-1,j-1) +  &
  57                              Kmat(i-1,k-1,j  ,11)*x(i-1,k-1,j  ) +  &
  58                              Kmat(i  ,k-1,j  ,10)*x(i  ,k-1,j  ) +  &
  59                              Kmat(i+1,k-1,j  , 9)*x(i+1,k-1,j  ) +  &
  60                              Kmat(i-1,k-1,j+1, 8)*x(i-1,k-1,j+1) +  &
  61                              Kmat(i  ,k-1,j+1, 7)*x(i  ,k-1,j+1) +  &
  62                              Kmat(i+1,k-1,j+1, 6)*x(i+1,k-1,j+1) +  &
  63                              Kmat(i-1,k  ,j-1, 5)*x(i-1,k  ,j-1) +  &
  64                              Kmat(i  ,k  ,j-1, 4)*x(i  ,k  ,j-1) +  &
  65                              Kmat(i+1,k  ,j-1, 3)*x(i+1,k  ,j-1) +  &
  66                              Kmat(i-1,k  ,j  , 2)*x(i-1,k  ,j  ) +  &
  67                              Kmat(i  ,k  ,j  , 2)*x(i+1,k  ,j  ) +  &
  68                              Kmat(i  ,k  ,j  , 3)*x(i-1,k  ,j+1) +  &
  69                              Kmat(i  ,k  ,j  , 4)*x(i  ,k  ,j+1) +  &
  70                              Kmat(i  ,k  ,j  , 5)*x(i+1,k  ,j+1) +  &
  71                              Kmat(i  ,k  ,j  , 6)*x(i-1,k+1,j-1) +  &
  72                              Kmat(i  ,k  ,j  , 7)*x(i  ,k+1,j-1) +  &
  73                              Kmat(i  ,k  ,j  , 8)*x(i+1,k+1,j-1) +  &
  74                              Kmat(i  ,k  ,j  , 9)*x(i-1,k+1,j  ) +  &
  75                              Kmat(i  ,k  ,j  ,10)*x(i  ,k+1,j  ) +  &
  76                              Kmat(i  ,k  ,j  ,11)*x(i+1,k+1,j  ) +  &
  77                              Kmat(i  ,k  ,j  ,12)*x(i-1,k+1,j+1) +  &
  78                              Kmat(i  ,k  ,j  ,13)*x(i  ,k+1,j+1) +  &
  79                              Kmat(i  ,k  ,j  ,14)*x(i+1,k+1,j+1)    &
  80                             )                                       &
  81                             )
  82                   x(i,k,j) = new
  83                   ! accumulate squared differeces and size
  84                   t = new - old
  85                   dif = max(dif,abs(t))
  86                   siz = max(siz,abs(old),abs(new))
  87               end do
  88            end do
  89         end do
  90     end do
  91 end do
  92
  93 reldif = dif/max(tiny(siz),siz)
  94
  95 !call write_array(x(ifts: ifte, kfts: kfte, jfts:jfte),'x_sweeps_out')
  96
  97 end subroutine sweeps
  98
  99 end module module_sweeps