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     print *,'allocate grid variables'
  98     do l=1,nlevels
  99         call allocate_mg_level(mg(l))
 100     enddo
 101
 102     print *,'computing  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         print *,'assembling the stiffness matrix on level  ',l
 138         call get_mg_dims(mg(l),                         &
 139             ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 140             ifms, ifme, kfms, kfme, jfms, jfme,         &
 141             ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 142             ifts, ifte, kfts, kfte, jfts, jfte)
 143         allocate(mg(l)%Kglo(ifms:ifme, kfms:kfme, jfms:jfme, msize))
 144         call ndt_assembly(                              &
 145             ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 146             ifms, ifme, kfms, kfme, jfms, jfme,         &
 147             ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 148             ifts, ifte, kfts, kfte, jfts, jfte,         &
 149             params%A, mg(l)%X,mg(l)%Y,mg(l)%Z, 1,              &
 150             mg(l)%Kglo)
 151         print *,'applying the boundary conditions on level ',l
 152         call ndt_boundary_conditions(                   &
 153             ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 154             ifms, ifme, kfms, kfme, jfms, jfme,         &
 155             ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 156             ifts, ifte, kfts, kfte, jfts, jfte,         &
 157             mg(l)%Kglo)
 158     enddo
 159
 160     print *,'end multigrid setup'
 161
 162 end subroutine multigrid_setup
 163
 164
 165 recursive subroutine multigrid_cycle(mg,l,rate)
 166 implicit none
 167
 168 !*** arguments
 169
 170 type(mg_type),intent(inout)::mg(:)  ! multigrid levels
 171 integer, intent(in)::l           ! level
 172 real, intent(out)::rate
 173
 174 !*** local
 175
 176 integer::   &
 177     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 178     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 179     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 180     ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 181     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 182     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 183     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 184     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte          ! coarse grid tile
 185
 186 integer::it, maxit, nit, nsmooth, cycles
 187 logical::coarse,coarsest
 188 real:: restol, norm_f, norm_res, diftol, siz, reldif, relres, rhsiz
 189 character(len=10)::it_kind
 190
 191 !*** executable
 192
 193 call get_mg_dims(mg(l),                         &
 194     ifds, ifde, kfds, kfde, jfds, jfde,         & ! fire grid dimensions
 195     ifms, ifme, kfms, kfme, jfms, jfme,         &
 196     ifps, ifpe, kfps, kfpe, jfps, jfpe,         & ! fire patch bounds
 197     ifts, ifte, kfts, kfte, jfts, jfte)
 198 if(l<nlevels) then
 199     call get_mg_dims(mg(l+1),                         &
 200     ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 201     ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 202     ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 203     ifcts, ifcte, kfcts,kfcte, jfcts,jfcte)         ! coarse grid tile
 204 !*** temporary fix, meshes are in vertices not elements here.
 205 !*** should be fixed using snode in all subroutines called from here
 206 !    ifcte = ifcte + 1
 207 !    jfcte = jfcte + 1
 208 !    kfcte = kfcte + 1
 209 endif
 210
 211 !*** temporary fix, meshes are in vertices not elements here.
 212 !*** should be fixed using snode in all subroutines called from here
 213 !ifte = ifte + 1
 214 !jfte = jfte + 1
 215 !kfte = kfte + 1
 216
 217 ! initialize
 218 cycles = 0
 219
 220 ! decide on number of iterations
 221 coarsest = l.eq.nlevels
 222 diftol = 0. ! default if params%check_reldif is false
 223 if(coarsest)then
 224     ! coarsest level
 225     maxit = params%coarsest_iter
 226     nsmooth = 0
 227     if(params%check_reldif)diftol = params%diftol_coarsest
 228 elseif(l == 1)then
 229     ! top level
 230     maxit = params%maxit
 231     nsmooth = params%nsmooth
 232     if(params%check_reldif)diftol = params%diftol_finest
 233 else
 234     ! some coarse level in between
 235     maxit = params%maxit_coarse
 236     nsmooth = params%nsmooth_coarse
 237     if(params%check_reldif)diftol = params%diftol_coarse
 238 endif
 239
 240 mg(l)%lambda=0.                             ! initial solution 0
 241
 242 if(params%check_relres)then
 243     norm_f = norm2( mg(l)%f(ifts: ifte, kfts: kfte, jfts: jfte))  !
 244     restol = norm_f * params%restol
 245     print *,'starting mg cycle level ',l,' norm_f ',norm_f,' residual tolerance ',restol
 246 endif
 247
 248
 249 do it=1,maxit
 250     coarse = mod(it,nsmooth+1)==0 .and. l < nlevels
 251     if(params%print_level>=l)print *,'level=',l,' it=',it,' coarse=',coarse
 252     if(coarse)then                                      ! coarse correction
 253         if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'cc_lambda_in')
 254         it_kind='coarse correction'
 255         ! compute residual residual = f - Kglo*lambda
 256         call ndt_mult(  &
 257             ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 258             ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 259             ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 260             ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 261             mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res, rhsiz, relres)
 262         if(params%print_level>=l)print *,'level=',l,' it=',it,'mult rhsiz=',siz,' rel res =',reldif
 263         if(params%debug_level >=l)call  write_array(mg(l)%res(ifts: ifte, kfts: kfte, jfts:jfte),'cc_res')
 264         ! restriction: f_coarse = R*residual
 265         call restriction(   &
 266             ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 267             ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 268             ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 269             ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile boundss
 270             ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 271             ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 272             ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 273             ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 274             mg(l+1)%f,mg(l)%res,                          &
 275             mg(l)%cr_x,mg(l)%cr_y,mg(l)%icl_z,            &
 276             mg(l)%X,mg(l)%Y,mg(l)%Z)
 277         if(params%debug_level >=l)call  write_array(mg(l+1)%f(ifcts: ifcte, kfcts: kfcte, jfcts:jfcte),'cc_f_c')
 278         if(params%debug_level >=l)call pause
 279         !if (l.eq.2) call crash('stopping here')
 280         ! call self on level l+1
 281         call multigrid_cycle(mg,l+1,rate)
 282
 283         ! prolongation lambda = lambda + P*lambda_coarse
 284
 285         if(params%debug_level >=l)call  write_array(mg(l+1)%lambda(ifcts: ifcte+1, kfcts: kfcte+1, jfcts:jfcte+1),'cc_lambda_c')
 286         call prolongation(   &
 287             ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 288             ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 289             ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 290             ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 291             ifcds, ifcde, kfcds,kfcde, jfcds,jfcde,       & ! coarse grid domain
 292             ifcms, ifcme, kfcms,kfcme, jfcms,jfcme,       & ! coarse grid dimensions
 293             ifcps, ifcpe, kfcps,kfcpe, jfcps,jfcpe,       & ! coarse grid dimensions
 294             ifcts, ifcte, kfcts,kfcte, jfcts,jfcte,       & ! coarse grid tile
 295             mg(l)%lambda,mg(l+1)%lambda,                  &
 296             mg(l)%cr_x,mg(l)%cr_y,mg(l)%icl_z,            &
 297             mg(l)%X,mg(l)%Y,mg(l)%Z)
 298     else
 299         it_kind = 'smoothing'
 300         if(coarsest)it_kind='coarsest'
 301         if(params%debug_level >=l)call  write_array(mg(l)%f(ifts: ifte, kfts: kfte, jfts:jfte),'sweep_f_in')
 302         if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'sweep_lambda_in')
 303         call sweeps( &
 304             ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 305             ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 306             ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 307             ifts, ifte, kfts, kfte, jfts, jfte,           & ! tile dimensions
 308             mg(l)%Kglo, mg(l)%f, mg(l)%lambda, siz, reldif)
 309         if(params%print_level>=l)print *,'level=',l,' it=',it,'sweeps siz=',siz,' rel diff=',reldif
 310         if(reldif < diftol)then
 311             if(params%print_level>=l)print *,'level ',l,' iteration ',it,' rel diff=',reldif,'< diftol=',diftol,' exiting'
 312             exit
 313         endif
 314         if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'sweep_lambda_out')
 315         if((.not.coarsest.or.it==maxit).and.(params%debug_level >=l))call pause
 316     endif
 317
 318     if(params%check_relres)then
 319         call ndt_mult(  &
 320             ifds, ifde, kfds, kfde, jfds, jfde,                       & ! fire domain bounds
 321             ifms, ifme, kfms, kfme, jfms, jfme,                       & ! fire memory bounds
 322             ifps, ifpe, kfps, kfpe, jfps, jfpe,                       & ! fire patch bounds
 323             ifts, ifte, kfts, kfte, jfts,jfte,                        & ! fire tile bounds
 324             mg(l)%Kglo, mg(l)%lambda, mg(l)%f, mg(l)%res, siz, relres)
 325         if(params%debug_level >=l)call  write_array(mg(l+1)%res(ifcts: ifcte, kfcts: kfcte, jfcts:jfcte),'it_res')
 326         norm_res = norm2( mg(l)%res(ifts: ifte, kfts: kfte, jfts: jfte))   ! residual norm
 327         if (mod(it,params%nsmooth+1)==params%nsmooth)then
 328             cycles=cycles+1;
 329         endif
 330         if(params%print_level>=l)print *,'cycle ',cycles,' level ',l, 'avg rate ',rate
 331         rate = (norm_res/norm_f)**(1d0/cycles);
 332         if(params%print_level>=l)print *,'level ',l,' iteration ',it,' ',it_kind,' residual ',norm_res, &
 333             ' norm_f ',norm_f,' rate ',rate
 334         if(norm_res < restol)exit
 335     endif
 336     if(params%debug_level >=l)call  write_array(mg(l)%lambda(ifts: ifte, kfts: kfte, jfts:jfte),'lambda_out')
 337
 338
 339 enddo
 340
 341 end subroutine multigrid_cycle
 342
 343 end module  module_multigrid