femwind/fortran/module_multigrid.f90

   1 module module_multigrid
   2    use module_common
   3    use module_utils
   4    ! use module_f_assembly
   5    use module_ndt_assembly
   6    ! use module_w_assembly
   7    use module_ndt_mult
   8    use module_coarsening
   9    use module_boundary_conditions
  10    use module_sweeps
  11
  12 integer::nlevels                                ! number of levels
  13
  14 contains
  15
  16 subroutine multigrid_setup(mg)
  17 implicit none
  18 ! set up the mg_level structure
  19 ! input: mg(1)%X,Y,Z,dx,dy,dz already set
  20 !*** arguments
  21 type(mg_type),intent(inout)::mg(:)  ! multigrid level
  22 !*** local
  23
  24 integer::k,l,m,nzc
  25 real::s
  26 integer::   &
  27     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
  28     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
  29     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  30     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
  31     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
  32     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
  33     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
  34     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte          ! coarse grid tile
  35
  36 !*** executable
  37
  38 ! decide on the grid coarsening and compute scalars and 1D index arrays
  39
  40     mg(1)%nn = mg(1)%nx *  mg(1)%ny * mg(1)%nz
  41     mg(1)%level = 1
  42
  43     print *,'multigrid_setup received'
  44     call print_mg_dims(mg,1)
  45
  46     do l=1,max_levels
  47
  48         ! decide on coarsening
  49
  50         print *,'multigrid level ',l,' grid size ', mg(l)%nx, mg(l)%ny, mg(l)%nz
  51         ! decide about the next level
  52         call coarsening_icl(mg(l)%cr_x,mg(l)%cr_y,mg(l)%icl_z,  &
  53                             mg(l)%dx,mg(l)%dy,mg(l)%dz,         &
  54                             params%A,params%minaspect,params%maxaspect)
  55
  56         ! get horizontal coarsening lists
  57
  58         call coarsening_hzc2icl(mg(l)%icl_x, mg(l)%icl_y, &
  59                                 mg(l)%cr_x,mg(l)%cr_y,&
  60                                 mg(l)%nx, mg(l)%ny)
  61
  62         ! update coarse mesh scalars
  63
  64         mg(l+1)%nx = size(mg(l)%icl_x)
  65         mg(l+1)%ny = size(mg(l)%icl_y)
  66         mg(l+1)%nz = size(mg(l)%icl_z)
  67         mg(l+1)%nn = mg(l+1)%nx *  mg(l+1)%ny * mg(l+1)%nz
  68         mg(l+1)%dx = mg(l)%dx * mg(l)%cr_x
  69         mg(l+1)%dy = mg(l)%dy * mg(l)%cr_y
  70         mg(l+1)%level = l+1
  71
  72         call print_mg_dims(mg,l+1)
  73
  74         if (mg(l+1)%nn >= mg(l)%nn .or. l == max_levels) then
  75             print *,'stopping at ',l,' levels because coarse ndof ',mg(l+1)%nn, &
  76                 '>= fine ', mg(l)%nn,' or l= ',l,' = ',max_levels
  77             nlevels = l
  78             exit ! the loop
  79         endif
  80
  81         ! coarsen dz
  82
  83         allocate(mg(l+1)%dz(mg(l+1)%nz-1))
  84         do k=1,mg(l+1)%nz - 1
  85             ! print *,'computing coarse dz of layer ',k
  86             s = 0.
  87             do m=mg(l)%icl_z(k),mg(l)%icl_z(k+1)-1
  88                 ! print *,'adding fine dz ',m,' equal ',mg(l)%dz(m)
  89                 s = s + mg(l)%dz(m)
  90             enddo
  91             mg(l+1)%dz(k)=s
  92             ! print *,'height dz of coarse layer',k,' is ',s
  93         enddo
  94     enddo
  95     print *,nlevels,' levels'
  96
  97     ! allocate grid variables
  98     do l=1,nlevels
  99         call allocate_mg_level(mg(l))
 100     enddo
 101
 102     ! get coarse grid coordinates
 103     do l=1,nlevels-1
 104         call get_mg_dims(mg(l), &
 105             ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
 106             ifms, ifme, kfms,kfme, jfms, jfme,            &
 107             ifps, ifpe, kfps,kfpe, jfps, jfpe,           & ! fire patch bounds
 108             ifts, ifte, kfts,kfte, jfts,jfte)
 109         call get_mg_dims(mg(l+1), &
 110             ifcds, ifcde, kfcds,kfcde, jfcds, jfcde,            & ! fire grid dimensions
 111             ifcms, ifcme, kfcms,kfcme, jfcms, jfcme,            &
 112             ifcps, ifcpe, kfcps,kfcpe, jfcps, jfcpe,           & ! fire patch bounds
 113             ifcts, ifcte, kfcts,kfcte, jfcts,jfcte)
 114         call coarsening_grid(l, &
 115             ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 116             ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 117             ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 118             ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 119             ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 120             ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 121             ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 122             ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 123             mg(l)%icl_x, mg(l)%icl_y, mg(l)%icl_z,        & ! indices of coarse grid
 124             mg(l)%X, mg(l)%Y, mg(l)%Z,                    & ! fine grid coordinates
 125             mg(l+1)%X, mg(l+1)%Y, mg(l+1)%Z)                ! coarse grid coordinates
 126     enddo
 127
 128     do l=1,nlevels
 129         if(params%debug_level >=l)call write_3d(mg,l,mg(l)%X,'X','v')
 130         if(params%debug_level >=l)call write_3d(mg,l,mg(l)%Y,'Y','v')
 131         if(params%debug_level >=l)call write_3d(mg,l,mg(l)%Z,'Z','v')
 132     enddo
 133
 134
 135     ! assemble the stiffness matrices
 136     do l=1,nlevels
 137         call get_mg_dims(mg(l),                         &
 138             ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 139             ifms, ifme, kfms, kfme, jfms, jfme,         &
 140             ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 141             ifts, ifte, kfts, kfte, jfts, jfte)
 142         allocate(mg(l)%Kglo(ifms:ifme, kfms:kfme, jfms:jfme, msize))
 143         call ndt_assembly(                              &
 144             ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 145             ifms, ifme, kfms, kfme, jfms, jfme,         &
 146             ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 147             ifts, ifte, kfts, kfte, jfts, jfte,         &
 148             params%A, mg(l)%X,mg(l)%Y,mg(l)%Z, 1,              &
 149             mg(l)%Kglo)
 150         call ndt_boundary_conditions(                   &
 151             ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 152             ifms, ifme, kfms, kfme, jfms, jfme,         &
 153             ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 154             ifts, ifte, kfts, kfte, jfts, jfte,         &
 155             mg(l)%Kglo)
 156     enddo
 157
 158 end subroutine multigrid_setup
 159
 160
 161 recursive subroutine multigrid_cycle(mg,l,rate)
 162 implicit none
 163
 164 !*** arguments
 165
 166 type(mg_type),intent(inout)::mg(:)  ! multigrid levels
 167 integer, intent(in)::l           ! level
 168 real, intent(out)::rate
 169
 170 !*** local
 171
 172 integer::   &
 173     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 174     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 175     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 176     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 177     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 178     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 179     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 180     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte          ! coarse grid tile
 181
 182 integer::it, maxit, nit, nsmooth, cycles
 183 logical::coarse,coarsest
 184 real:: restol, norm_f, norm_res, diftol, reldif
 185 character(len=10)::it_kind
 186
 187 !*** executable
 188
 189 call get_mg_dims(mg(l),                         &
 190     ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 191     ifms, ifme, kfms, kfme, jfms, jfme,         &
 192     ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 193     ifts, ifte, kfts, kfte, jfts, jfte)
 194 if(l<nlevels) then
 195     call get_mg_dims(mg(l+1),                         &
 196     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 197     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 198     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 199     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte)         ! coarse grid tile
 200 !*** temporary fix, meshes are in vertices not elements here.
 201 !*** should be fixed using snode in all subroutines called from here
 202 !    ifcte = ifcte + 1
 203 !    jfcte = jfcte + 1
 204 !    kfcte = kfcte + 1
 205 endif
 206
 207 !*** temporary fix, meshes are in vertices not elements here.
 208 !*** should be fixed using snode in all subroutines called from here
 209 !ifte = ifte + 1
 210 !jfte = jfte + 1
 211 !kfte = kfte + 1
 212
 213 ! initialize
 214 cycles = 0
 215
 216 ! decide on number of iterations
 217 coarsest = l.eq.nlevels
 218 diftol = 0. ! default if params%check_reldif is false
 219 if(coarsest)then
 220     ! coarsest level
 221     maxit = params%coarsest_iter
 222     nsmooth = 0
 223     if(params%check_reldif)diftol = params%diftol_coarsest
 224 elseif(l == 1)then
 225     ! top level
 226     maxit = params%maxit
 227     nsmooth = params%nsmooth
 228     if(params%check_reldif)diftol = params%diftol_finest
 229 else
 230     ! some coarse level in between
 231     maxit = params%maxit_coarse
 232     nsmooth = params%nsmooth_coarse
 233     if(params%check_reldif)diftol = params%diftol_coarse
 234 endif
 235
 236 mg(l)%lambda=0.                             ! initial solution 0
 237
 238 if(params%check_relres)then
 239     norm_f = norm2( mg(l)%f(ifts: ifte, kfts: kfte, jfts: jfte))  !
 240     restol = norm_f * params%restol
 241     print *,'starting mg cycle level ',l,' norm_f ',norm_f,' residual tolerance ',restol
 242 endif
 243
 244
 245 do it=1,maxit
 246     coarse = mod(it,nsmooth+1)==0 .and. l < nlevels
 247     if(params%print_level>=l)print *,'level=',l,' it=',it,' coarse=',coarse
 248     if(coarse)then                                      ! coarse correction
 249         if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'cc_lambda_in')
 250         it_kind='coarse correction'
 251         ! compute residual residual = f - Kglo*lambda
 252         call ndt_mult(  &
 253             ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 254             ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 255             ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 256             ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 257             mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res)
 258         if(params%debug_level >=l)call  write_array(mg(l)%res(ifts: ifte, kfts: kfte, jfts:jfte),'cc_res')
 259         ! restriction: f_coarse = R*residual
 260         call restriction(   &
 261             ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 262             ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 263             ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 264             ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile boundss
 265             ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 266             ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 267             ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 268             ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 269             mg(l+1)%f,mg(l)%res,                          &
 270             mg(l)%cr_x,mg(l)%cr_y,mg(l)%icl_z,            &
 271             mg(l)%X,mg(l)%Y,mg(l)%Z)
 272         if(params%debug_level >=l)call  write_array(mg(l+1)%f(ifcts: ifcte, kfcts: kfcte, jfcts:jfcte),'cc_f_c')
 273         if(params%debug_level >=l)call pause
 274         !if (l.eq.2) call crash('stopping here')
 275         ! call self on level l+1
 276         call multigrid_cycle(mg,l+1,rate)
 277
 278         ! prolongation lambda = lambda + P*lambda_coarse
 279
 280         if(params%debug_level >=l)call  write_array(mg(l+1)%lambda(ifcts: ifcte+1, kfcts: kfcte+1, jfcts:jfcte+1),'cc_lambda_c')
 281         call prolongation(   &
 282             ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 283             ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 284             ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 285             ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 286             ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 287             ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 288             ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 289             ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 290             mg(l)%lambda,mg(l+1)%lambda,                  &
 291             mg(l)%cr_x,mg(l)%cr_y,mg(l)%icl_z,            &
 292             mg(l)%X,mg(l)%Y,mg(l)%Z)
 293     else
 294         it_kind = 'smoothing'
 295         if(coarsest)it_kind='coarsest'
 296         if(params%debug_level >=l)call  write_array(mg(l)%f(ifts: ifte, kfts: kfte, jfts:jfte),'sweep_f_in')
 297         if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'sweep_lambda_in')
 298         call sweeps( &
 299             ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 300             ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 301             ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 302             ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 303             mg(l)%Kglo, mg(l)%f, mg(l)%lambda, reldif)
 304         if(reldif < diftol)then
 305             if(params%print_level>=l)print *,'level ',l,' iteration ',it,' rel diff=',reldif,'< diftol=',diftol,' exiting'
 306             exit
 307         endif
 308         if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'sweep_lambda_out')
 309         if((.not.coarsest.or.it==maxit).and.(params%debug_level >=l))call pause
 310     endif
 311
 312     if(params%check_relres)then
 313         call ndt_mult(  &
 314             ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 315             ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 316             ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 317             ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 318             mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res)
 319         if(params%debug_level >=l)call  write_array(mg(l+1)%res(ifcts: ifcte, kfcts: kfcte, jfcts:jfcte),'it_res')
 320         norm_res = norm2( mg(l)%res(ifts: ifte, kfts: kfte, jfts: jfte))   ! residual norm
 321         if (mod(it,params%nsmooth+1)==params%nsmooth)then
 322             cycles=cycles+1;
 323         endif
 324         if(params%print_level>=l)print *,'cycle ',cycles,' level ',l, 'avg rate ',rate
 325         rate = (norm_res/norm_f)**(1d0/cycles);
 326         if(params%print_level>=l)print *,'level ',l,' iteration ',it,' ',it_kind,' residual ',norm_res, &
 327             ' norm_f ',norm_f,' rate ',rate
 328         if(norm_res < restol)exit
 329     endif
 330     if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'lambda_out')
 331
 332
 333 enddo
 334
 335 end subroutine multigrid_cycle
 336
 337 end module  module_multigrid