femwind/fortran/module_coarsening.f90

   1 module module_coarsening
   2
   3 use module_utils
   4
   5 contains
   6
   7 subroutine prolongation(   &
   8     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
   9     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
  10     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  11     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
  12     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
  13     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
  14     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
  15     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
  16     u,uc,cr_x,cr_y,icl_z,X,Y,Z)
  17
  18 ! Multiply by the prolongation matrix: u = u + P*uc
  19 ! In:
  20 !   uc      coarse grid vector
  21 !   cr_x, cr_y  coarsening factor in horizontal directions x and y
  22 !   icl_z    1D array, indices of coarse grid in the z directions
  23 !   X,Y,Z   grid coordinates
  24 ! Out:
  25 !   u      fine grid vector
  26
  27 implicit none
  28 !*** arguments
  29
  30 integer, intent(in)::                             &
  31     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
  32     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
  33     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  34     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
  35     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
  36     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
  37     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
  38     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte          ! coarse grid tile
  39
  40 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: X,Y,Z !spatial grid
  41 real, intent(in), dimension(ifcms:ifcme,kfcms:kfcme,jfcms:jfcme):: uc ! coarse vector
  42
  43 integer, intent(in):: cr_x, cr_y, &       ! coarsening factors in the horizonal directions
  44     icl_z(kfcts:kfcte)                      ! indices of coarse grid in the vertical direction
  45
  46 real, intent(inout), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: u ! fine grid interpolant
  47
  48 !*** local
  49 integer :: i,j,k,ic,jc,kc,ifc,jfc,kfc,ifs,ife,jfs,jfe,kfs,kfe
  50 real:: qi,qj,qk
  51
  52 !*** executable
  53
  54 ! zero the output to ready for contributions
  55 ! u(ifts:ifte,kfts:kfte,jfts:jfte) = 0.
  56
  57 if ((icl_z(kfcts) .ne. kfts) .or. (icl_z(kfcte) .ne. kfte)) then
  58     call crash('vertical corsening must include all domain')
  59 endif
  60
  61 do kc=kfcts,kfcte           ! loop over coarse layers
  62     kfc=icl_z(kc);          ! the fine grid number of the coarse layer
  63     if (kfc>kfts) then
  64         kfs=icl_z(kc-1)+1   ! from above previous coarse
  65     else
  66         kfs=kfc            ! itself, there is no previous fine layer
  67     endif
  68     if (kfc<kfde) then
  69         kfe=icl_z(kc+1)-1   ! up to under next layer
  70     else
  71         kfe=kfc            ! itself, there is no next layer
  72     endif
  73     !print *,'vertical coarse layer ',kc,' at ',kfc,' contributes to layers ',kfs,':',kfe
  74     do k=kfs,kfe
  75         ! really needs to be made fine point oriented and draw from coarse points
  76         do jc=jfcts,jfcte
  77             jfc=cr_y*(jc-jfcds)+jfds ! fine grid index of the coarse point
  78             jfs=max(jfc-cr_y+1,jfds) ! start support
  79             jfe=min(jfc+cr_y-1,jfde) ! end support
  80             if (jfc > jfde) then  ! after end of domain not divisible by coarse ratio
  81                jfc = jfde
  82             elseif (jfc > jfde - cr_y .and. jfc < jfde)then ! just before end of domain not divisible by coarse ratio
  83                jfe = jfde -1
  84             endif
  85             !print *,'coarse y ',jc,' at j ',jfc,' contributes to ',jfs,':',jfe
  86                 do ic=ifcts,ifcte
  87                 ifc=cr_y*(ic-ifcds)+ifds ! fine grid index of the coarse point
  88                 ifs=max(ifc-cr_y+1,ifds) ! start support
  89                 ife=min(ifc+cr_y-1,ifde) ! end support
  90                 if (ifc > ifde) then  ! after end of domain not divisible by coarse ratio
  91                    ifc = ifde
  92                 elseif (ifc > ifde - cr_x .and. ifc < ifde)then ! just before end of domain not divisible by coarse ratio
  93                    ife = ifde -1
  94                 endif
  95                 !print *,'coarse x ',ic,' at i ',ifc,' contributes to ',ifs,':',ife
  96                 do j=jfs,jfe
  97                     do i=ifs,ife
  98                         if (i>ifc) then
  99                             qi=(X(i,k,j)-X(ife+1,k,j))/(X(ifc,k,j)-X(ife+1,k,j))
 100                         elseif (i<ifc) then
 101                             qi=(X(i,k,j)-X(ifs-1,k,j))/(X(ifc,k,j)-X(ifs-1,k,j))
 102                         else
 103                             qi=1.
 104                         endif
 105                         if (j>jfc) then
 106                             qj=(Y(i,k,j)-Y(i,k,jfe+1))/(Y(i,k,jfc)-Y(i,k,jfe+1))
 107                         elseif (j<jfc) then
 108                             qj=(Y(i,k,j)-Y(i,k,jfs-1))/(Y(i,k,jfc)-Y(i,k,jfs-1))
 109                         else
 110                             qj=1.
 111                         endif
 112                         if (k>kfc) then
 113                             qk=(Z(i,k,j)-Z(i,kfe+1,j))/(Z(i,kfc,j)-Z(i,kfe+1,j))
 114                         elseif (k<kfc) then
 115                             qk=(Z(i,k,j)-Z(i,kfs-1,j))/(Z(i,kfc,j)-Z(i,kfs-1,j))
 116                         else
 117                             qk=1.
 118                         endif
 119                         u(i,k,j) = u(i,k,j) + qi*qk*qj*uc(ic,kc,jc);
 120                     enddo
 121                 enddo
 122             enddo
 123         enddo
 124     enddo
 125 enddo
 126
 127 end subroutine prolongation
 128
 129 subroutine restriction(   &
 130     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 131     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 132     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 133     ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 134     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 135     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 136     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 137     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 138     uc,u,cr_x,cr_y,icl_z,X,Y,Z)
 139
 140 ! Multiply by the prolongation matrix transpose
 141 ! In:
 142 !   u      fine grid vector
 143 !   cr_x, cr_y  coarsening factor in horizontal directions x and y
 144 !   icl_z    1D array, indices of coarse grid in the z directions
 145 !   X,Y,Z   grid coordinates
 146 ! Out:
 147 !   uc      coarse grid vector
 148
 149 implicit none
 150 !*** arguments
 151
 152 integer, intent(in)::                             &
 153     ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 154     ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 155     ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 156     ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile boundss                ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 157     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 158     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 159     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 160     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte          ! coarse grid tile
 161
 162 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: X,Y,Z !spatial grid
 163 real, intent(out), dimension(ifcms:ifcme,kfcms:kfcme,jfcms:jfcme):: uc ! coarse vector
 164
 165 integer, intent(in):: cr_x, cr_y, &       ! coarsening factors in the horizonal directions
 166     icl_z(kfcts:kfcte)                      ! indices of coarse grid in the vertical direction
 167
 168 real, intent(in), dimension(ifms:ifme,kfms:kfme,jfms:jfme):: u ! fine grid interpolant
 169
 170 !*** local
 171 integer :: i,j,k,ic,jc,kc,ifc,jfc,kfc,ifs,ife,jfs,jfe,kfs,kfe
 172 real:: qi,qj,qk
 173
 174 !*** executable
 175
 176 ! zero the output to ready for contributions
 177 uc(ifcts:ifcte,kfcts:kfcte,jfcts:jfcte) = 0.
 178
 179 if ((icl_z(kfcts) .ne. kfts) .or. (icl_z(kfcte) .ne. kfte)) then
 180     call crash('vertical corsening must include all domain')
 181 endif
 182
 183 do kc=kfcts,kfcte           ! loop over coarse layers
 184     kfc=icl_z(kc);          ! the fine grid number of the coarse layer
 185     if (kfc>kfts) then
 186         kfs=icl_z(kc-1)+1   ! from above previous coarse
 187     else
 188         kfs=kfc            ! itself, there is no previous fine layer
 189     endif
 190     if (kfc<kfde) then
 191         kfe=icl_z(kc+1)-1   ! up to under next layer
 192     else
 193         kfe=kfc            ! itself, there is no next layer
 194     endif
 195     !print *,'vertical coarse layer ',kc,' at ',kfc,' contributes to layers ',kfs,':',kfe
 196     do k=kfs,kfe
 197         ! really needs to be made fine point oriented and draw from coarse points
 198         do jc=jfcts,jfcte
 199             jfc=cr_y*(jc-jfcds)+jfds ! fine grid index of the coarse point
 200             jfs=max(jfc-cr_y+1,jfds) ! start support
 201             jfe=min(jfc+cr_y-1,jfde) ! end support
 202             if (jfc > jfde) then  ! after end of domain not divisible by coarse ratio
 203                jfc = jfde
 204             elseif (jfc > jfde - cr_y .and. jfc < jfde)then ! just before end of domain not divisible by coarse ratio
 205                jfe = jfde -1
 206             endif
 207             !print *,'coarse y ',jc,' at j ',jfc,' contributes to ',jfs,':',jfe
 208                 do ic=ifcts,ifcte
 209                 ifc=cr_y*(ic-ifcds)+ifds ! fine grid index of the coarse point
 210                 ifs=max(ifc-cr_y+1,ifds) ! start support
 211                 ife=min(ifc+cr_y-1,ifde) ! end support
 212                 if (ifc > ifde) then  ! after end of domain not divisible by coarse ratio
 213                    ifc = ifde
 214                 elseif (ifc > ifde - cr_x .and. ifc < ifde)then ! just before end of domain not divisible by coarse ratio
 215                    ife = ifde -1
 216                 endif
 217                 !print *,'coarse x ',ic,' at i ',ifc,' contributes to ',ifs,':',ife
 218                 do j=jfs,jfe
 219                     do i=ifs,ife
 220                         if (i>ifc) then
 221                             qi=(X(i,k,j)-X(ife+1,k,j))/(X(ifc,k,j)-X(ife+1,k,j))
 222                         elseif (i<ifc) then
 223                             qi=(X(i,k,j)-X(ifs-1,k,j))/(X(ifc,k,j)-X(ifs-1,k,j))
 224                         else
 225                             qi=1.
 226                         endif
 227                         if (j>jfc) then
 228                             qj=(Y(i,k,j)-Y(i,k,jfe+1))/(Y(i,k,jfc)-Y(i,k,jfe+1))
 229                         elseif (j<jfc) then
 230                             qj=(Y(i,k,j)-Y(i,k,jfs-1))/(Y(i,k,jfc)-Y(i,k,jfs-1))
 231                         else
 232                             qj=1.
 233                         endif
 234                         if (k>kfc) then
 235                             qk=(Z(i,k,j)-Z(i,kfe+1,j))/(Z(i,kfc,j)-Z(i,kfe+1,j))
 236                         elseif (k<kfc) then
 237                             qk=(Z(i,k,j)-Z(i,kfs-1,j))/(Z(i,kfc,j)-Z(i,kfs-1,j))
 238                         else
 239                             qk=1.
 240                         endif
 241                         uc(ic,kc,jc) = uc(ic,kc,jc) + u(i,k,j)*qi*qk*qj
 242                     enddo
 243                 enddo
 244             enddo
 245         enddo
 246     enddo
 247 enddo
 248
 249 end subroutine restriction
 250
 251 subroutine coarsening_icl(cr_x,cr_y,icl_z,dx,dy,dz,A,minaspect,maxaspect)
 252 ! decide on coarsening
 253 ! in:
 254 !   dx,dy       mesh spacings, scalar
 255 !   dz          verticl element size, vector
 256 !   A           matrix size (3,3), only diagonal used
 257 ! out:
 258 !   cr_x,cr_y   horizontal coarsening factors in directions x and y
 259 !   icl_z       coarse indices in direction z, allocated here
 260 implicit none
 261
 262 !*** arguments
 263 real,intent(in):: dx,dy,dz(:),A(:,:)
 264 real,intent(in):: minaspect,maxaspect
 265 integer,intent(out):: cr_x,cr_y
 266 integer,pointer,intent(out):: icl_z(:)
 267
 268 !*** local
 269 real:: dxy,crit,hzcavg,arat
 270 integer, allocatable:: icl3(:)
 271 integer:: nc3,newlcl,lcl,i,n3
 272
 273
 274 !*** executable
 275     dxy=min(dx,dy)  ! horizontal step
 276     n3 = size(dz)+1 !
 277     print *,'coarsening_icl in: dx=',dx,' dy=',dy,' n3=',n3,' dz=',dz
 278     print *,'coarsening_icl in: A=',A,' minaspect=',minaspect,' maxaspect=',maxaspect
 279     arat = A(3,3)/min(A(1,1),A(2,2))  ! scaled vertical penalty
 280     ! decide on horizontal coarsening factors
 281     crit=(dz(1)/dxy)/arat
 282
 283     if (crit > minaspect) then
 284         cr_x = 2
 285         cr_y = 2
 286     else
 287         cr_x = 1
 288         cr_y = 1
 289     endif
 290     hzcavg=sqrt(real(cr_x*cr_y));
 291     print *,'horizontal coarsening factors ',cr_x,cr_y, ' because weighted height is ', crit
 292     allocate(icl3(n3+1))
 293     icl3=0
 294     lcl=1 ! last coarse level
 295     icl3(1)=lcl
 296     nc3=0
 297     do i=1,n3
 298         newlcl=lcl+1  ! next coarse level by 1
 299         if (lcl+2 <= n3) then
 300             crit = ((dz(lcl)+dz(lcl+1))/(dxy*hzcavg))/arat
 301             if (crit < maxaspect ) then
 302                 newlcl=lcl+2 ! next coarse level by 2
 303             endif
 304         endif
 305         lcl = newlcl;
 306         if (lcl <= n3) then
 307             icl3(i+1)=lcl
 308         else ! at the top already
 309             nc3=i
 310             allocate(icl_z(i))
 311             icl_z = icl3(1:i)
 312             exit
 313         endif
 314     enddo
 315     deallocate(icl3)
 316     if (nc3==0) then
 317         call crash('coarsening_icl: number of coarse layers is 0')
 318     endif
 319     print *,'vertical coarse layers ',icl_z
 320
 321 end subroutine coarsening_icl
 322
 323 subroutine coarsening_hzc2icl(icl_x, icl_y, cr_x, cr_y, n_x, n_y)
 324 ! translate horizontal coarsening factors to index vectors
 325 !
 326 ! In:
 327 !   cr_x, cr_y coarsening factors in x and y
 328 !   n_x, n_y   fine mesh size in x and y
 329 ! Out:
 330 !   icl_x, icl_y    fine indices of coarse grid lines
 331
 332 implicit none
 333
 334 !*** arguments
 335 integer, intent(in)::cr_x, cr_y, n_x, n_y
 336 integer, intent(out), pointer, dimension(:) :: icl_x, icl_y
 337
 338 !*** local
 339 integer :: nc_x, nc_y
 340 integer::i
 341
 342 nc_x = min(n_x,n_x/cr_x + 1)
 343 nc_y = min(n_y,n_y/cr_y + 1)
 344
 345 allocate(icl_x(nc_x))
 346 allocate(icl_y(nc_y))
 347
 348 do i=1,nc_x-1
 349     icl_x(i)=1+cr_x*(i-1)  ! every second node
 350 enddo
 351 icl_x(nc_x)=n_x ! last coarse is last fine
 352
 353 do i=1,nc_y-1
 354     icl_y(i)=1+cr_y*(i-1)  ! every second node
 355 enddo
 356 icl_y(nc_y)=n_y ! last coarse is last fine
 357
 358 end subroutine coarsening_hzc2icl
 359
 360 subroutine coarsening_grid(l, &
 361     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 362     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 363     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 364     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 365     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 366     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 367     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 368     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 369     icl_x, icl_y, icl_z,                          &
 370     X, Y, Z, X_coarse, Y_coarse, Z_coarse)
 371 implicit none
 372 ! compute coarse grid coordinates arrays on one tile
 373
 374 !*** arguments
 375 integer,intent(in)::l                               ! fine level, coarse is l+1
 376 integer, intent(in):: &
 377     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 378     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 379     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 380     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 381     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 382     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 383     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 384     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte          ! coarse grid tile
 385
 386 integer, intent(in):: icl_x(ifcts:),icl_y(jfcts:),icl_z(kfcts:)  ! leaving upper bound to be passed on
 387 real, dimension(ifms:ifme,kfms:kfme,jfms:jfme), intent(in):: X,Y,Z !spatial grid
 388 real, intent(out), dimension(ifcms:ifcme,kfcms:kfcme,jfcms:jfcme):: X_coarse, Y_coarse, Z_coarse ! coarse vector
 389
 390 !*** local
 391 integer::ic,jc,kc,ie,je,ke
 392 character(len=2)::lc
 393
 394 !*** executable
 395 ie = snode(ifcte,ifcde,+1)
 396 je = snode(jfcte,jfcde,+1)
 397 ke = snode(kfcte,kfcde,+1)
 398
 399 do jc=ifcts,je
 400     do kc=kfcts,ke
 401         do ic=ifcts,ie
 402             X_coarse(ic,kc,jc) = X(icl_x(ic),icl_z(kc),icl_y(jc))
 403             Y_coarse(ic,kc,jc) = Y(icl_x(ic),icl_z(kc),icl_y(jc))
 404             Z_coarse(ic,kc,jc) = Z(icl_x(ic),icl_z(kc),icl_y(jc))
 405         enddo
 406     enddo
 407 enddo
 408
 409 end subroutine coarsening_grid
 410
 411 end module module_coarsening
 412