femwind/fortran/module_common.f90

   1 module module_common
   2
   3 use module_utils
   4
   5 integer, parameter::msize=14
   6 integer, parameter::max_levels=5
   7
   8 ! method parameters
   9 type params_type
  10     real:: minaspect=1./3.,maxaspect=3.
  11     real:: A(3,3)=reshape((/1., 0., 0.,  0., 1., 0.,  0., 0., 1./),(/3, 3/))
  12     integer:: coarsest_iter=100 ! to solve the coarse problem
  13     integer:: maxit=50     ! total iterations
  14     integer:: nsmooth=3    ! smoothing iterations before correction
  15     integer:: nsmooth_coarse=2 ! smoothing iterations on coarse levels
  16     integer:: maxit_coarse=8 ! on levels>1: 2 smoothing, coarse, 2 smoothing, coarse, 2 smoothing
  17     integer:: debug_level=-1 ! multigrid level to debug up to
  18     integer:: msglevel=1 ! type of messages to print
  19     integer:: print_level=2 ! multigrid level to print messages up to
  20     !logical:: check_relres=.true.  ! check residual after every base iteration if to continue
  21     logical:: check_relres=.false.
  22     real:: restol=1e-5 ! residual tolerance to stop iterating
  23     ! real:: restol=1e-6 ! residual tolerance to stop iterating
  24     !logical:: check_reldif=.false.  ! check difference after every base iteration if to continue
  25     logical:: check_reldif=.true.  ! check difference after every base iteration if to continue
  26     ! real:: diftol_finest=1e-6 ! keep iterating while iterates change by this
  27     real:: diftol_finest=1e-5 ! keep iterating while iterates change by this
  28     real:: diftol_coarse=1e-1 ! keep iterating while iterates change by this on corse levels
  29     real:: diftol_coarsest=1e-2 ! keep iterating while iterates change by this on corsest level
  30 end type
  31
  32 type(params_type)::params
  33
  34 ! multigrid structure
  35
  36 type mg_type
  37     real, dimension(3,3):: A                        ! penalty weight matrix
  38     real, pointer, dimension(:,:,:):: X, Y, Z       ! grid vertices
  39     real, pointer, dimension(:,:,:)::        &      !
  40          f, lambda, res                             ! multigrid variables
  41
  42     real, pointer, dimension(:,:,:,:):: Kglo        ! global stiffness matrix
  43
  44     real:: dx,dy                                    ! horizontal spacing, scalar
  45     real, pointer, dimension(:)::dz                 ! vertical spacing of the layers
  46     integer::nx,ny,nz,nn                            ! mesh size in terms of vertices
  47
  48     integer::level
  49
  50     integer:: cr_x, cr_y                            ! coarsening factors
  51     integer, pointer, dimension(:):: icl_x, icl_y, icl_z  ! coarsening indices
  52
  53 end type
  54
  55 type(mg_type):: mg(max_levels+1)  ! the main multigrid structure
  56
  57 contains
  58
  59 subroutine read_params
  60     real:: minaspect=1./3.,maxaspect=3.
  61     real:: A(3,3)=reshape((/1., 0., 0.,  0., 1., 0.,  0., 0., 1./),(/3, 3/))
  62     integer:: coarsest_iter=100 ! to solve the coarse problem
  63     integer:: maxit=50     ! total iterations
  64     integer:: nsmooth=3    ! smoothing iterations before correction
  65     integer:: nsmooth_coarse=2 ! smoothing iterations on coarse levels
  66     integer:: maxit_coarse=8 ! on levels>1: 2 smoothing, coarse, 2 smoothing, coarse, 2 smoothing
  67     integer:: debug_level=-1 ! multigrid level to debug up to
  68     integer:: msglevel=1 ! type of messages to print
  69     integer:: print_level=2 ! multigrid level to print messages up to
  70     !logical:: check_relres=.true.  ! check residual after every base iteration if to continue
  71     logical:: check_relres=.false.
  72     real:: restol=1e-6 ! residual tolerance to stop iterating
  73     !logical:: check_reldif=.false.  ! check difference after every base iteration if to continue
  74     logical:: check_reldif=.true.  ! check difference after every base iteration if to continue
  75     real:: diftol_finest=1e-6 ! keep iterating while iterates change by this
  76     real:: diftol_coarse=1e-1 ! keep iterating while iterates change by this on corse levels
  77     real:: diftol_coarsest=1e-2 ! keep iterating while iterates change by this on corsest level
  78
  79     namelist/params_nml/ &
  80       minaspect, &
  81       A, &
  82       coarsest_iter, &
  83       maxit, &
  84       nsmooth, &
  85       nsmooth_coarse, &
  86       maxit_coarse, &
  87       debug_level, &
  88       msglevel, &
  89       print_level, &
  90       check_relres, &
  91       restol, &
  92       check_reldif, &
  93       diftol_finest, &
  94       diftol_coarse, &
  95       diftol_coarsest
  96
  97     write(*,'(a)')'Default parameters:'
  98     write(*,nml=params_nml)
  99     open(8,file='params.nml',status='old',err=9)
 100     read(8,nml=params_nml)
 101     close(8)
 102      params%minaspect=minaspect
 103      params%A=A
 104      params%coarsest_iter=coarsest_iter
 105      params%maxit=maxit
 106      params%nsmooth=nsmooth
 107      params%nsmooth_coarse=nsmooth_coarse
 108      params%maxit_coarse=maxit_coarse
 109      params%debug_level=debug_level
 110      params%msglevel=msglevel
 111      params%print_level=print_level
 112      params%check_relres=check_relres
 113      params%restol=restol
 114      params%check_reldif=check_reldif
 115      params%diftol_finest=diftol_finest
 116      params%diftol_coarse=diftol_coarse
 117      params%diftol_coarsest=diftol_coarsest
 118     write(*,'(a)')'Using modified parameters:'
 119     write(*,nml=params_nml)
 120     return
 121 9   write(*,'(a)')'Cannot open file params.nml, using defaults'
 122     end subroutine read_params
 123
 124
 125 subroutine allocate_mg_level(mg)    !allocate one multigrid level
 126 type(mg_type),intent(inout)::mg     !multigrid structure
 127 !*** local
 128 integer:: &
 129     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 130     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 131     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 132     ifts, ifte, kfts, kfte, jfts, jfte              ! tile dimensions
 133
 134 !*** executable
 135 call get_mg_dims(mg, &
 136     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
 137     ifms, ifme, kfms,kfme, jfms, jfme,            &
 138     ifps, ifpe, kfps,kfpe, jfps, jfpe,            & ! fire patch bounds
 139     ifts, ifte, kfts,kfte, jfts,jfte)
 140
 141 if (.not.associated(mg%X))allocate(mg%X(ifms: ifme, kfms: kfme, jfms: jfme))
 142 if (.not.associated(mg%Y))allocate(mg%Y(ifms: ifme, kfms: kfme, jfms: jfme))
 143 if (.not.associated(mg%Z))allocate(mg%Z(ifms: ifme, kfms: kfme, jfms: jfme))
 144 allocate(mg%F(ifms: ifme, kfms: kfme, jfms: jfme))
 145 allocate(mg%lambda(ifms: ifme, kfms: kfme, jfms: jfme))
 146 allocate(mg%res(ifms: ifme, kfms: kfme, jfms: jfme))
 147
 148 end subroutine allocate_mg_level
 149
 150 subroutine get_mg_dims(mg_level, &
 151     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 152     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 153     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 154     ifts, ifte, kfts, kfte, jfts, jfte)              ! tile dimensions
 155 implicit none
 156 type(mg_type),intent(in)::mg_level
 157 integer, intent(out)::  &
 158     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 159     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 160     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 161     ifts, ifte, kfts, kfte, jfts, jfte              ! tile dimensions
 162 character(len=256)::msg
 163
 164     ifds = 1
 165     ifde = mg_level%nx-1  ! dimension in elements
 166     jfds = 1
 167     jfde = mg_level%ny-1
 168     kfds = 1
 169     kfde = mg_level%nz-1
 170
 171     ifts = ifds
 172     ifte = ifde
 173     jfts = jfds
 174     jfte = jfde
 175     kfts = kfds
 176     kfte = kfde
 177
 178     ifps = ifds
 179     ifpe = ifde
 180     jfps = jfds
 181     jfpe = jfde
 182     kfps = kfds
 183     kfpe = kfde
 184
 185     ifms = ifps - 1
 186     ifme = ifpe + 2  ! need +1 for vertex grid, and +1 for zero strip
 187     jfms = jfps - 1
 188     jfme = jfpe + 2
 189     kfms = kfps - 1
 190     kfme = kfpe + 2
 191 write(msg,*)"get_mg_dims: domain",ifds, ifde, kfds, kfde, jfds, jfde
 192 call message(msg)
 193 write(msg,*)"get_mg_dims: memory",ifms, ifme, kfms, kfme, jfms, jfme
 194 call message(msg)
 195 write(msg,*)"get_mg_dims: patch ",ifps, ifde, kfps, kfpe, jfps, jfpe
 196 call message(msg)
 197 write(msg,*)"get_mg_dims: tile  ",ifts, ifte, kfts, kfte, jfts, jfte
 198 call message(msg)
 199
 200 end subroutine get_mg_dims
 201
 202 subroutine write_3d(mg,l,var,name,v)
 203 ! write array in text file for debugging
 204 implicit none
 205
 206 !*** arguments
 207 type(mg_type),intent(in)::mg(:)    !multigrid structure
 208 integer,intent(in)::l              !level
 209 real, intent(in)::var(:,:,:)       !the array
 210 character(len=*),intent(in)::name
 211 character(len=1),intent(in)::v     !'v' or 'e'
 212
 213 !*** local
 214 integer:: &
 215     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 216     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 217     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 218     ifts, ifte, kfts, kfte, jfts, jfte              ! tile dimensions
 219 integer::ie,ke,je
 220 character(len=2)::lc
 221
 222 !*** executable
 223
 224 call get_mg_dims(mg(l), &
 225     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
 226     ifms, ifme, kfms,kfme, jfms, jfme,            &
 227     ifps, ifpe, kfps,kfpe, jfps, jfpe,           & ! fire patch bounds
 228     ifts, ifte, kfts,kfte, jfts,jfte)
 229
 230 write(lc,'(i2.2)')l
 231
 232 if(v.eq.'v')then
 233     ie = snode(ifte,ifde,+1)  ! vertex based
 234     je = snode(jfte,jfde,+1)
 235     ke = snode(kfte,kfde,+1)
 236 elseif(v.eq.'e')then
 237     ie=ifte
 238     je=jfte
 239     ke=kfte
 240 else
 241     call crash('write_3d: bad v')
 242 endif
 243
 244 call write_tight(var) ! write dereferenced
 245
 246 contains
 247
 248 subroutine write_tight(varr)
 249     real, intent(in)::varr(ifms:ifme,kfms:kfme,jfms:jfme)
 250     ! for element/cell/midpoint based arrays, like fmw winds
 251     ! open(iu,file='varr.dat',form='unformatted',status='unknown')
 252     ! write(iu)varr(ifts:ifte,kfts:kfte,jfts:jfte)
 253     ! close(iu)
 254     call write_array(varr(ifts:ie,kfts:ke,jfts:je),name//lc)
 255 end subroutine write_tight
 256
 257 end subroutine write_3d
 258
 259 subroutine print_mg_dims(mg,l)
 260 type(mg_type),intent(in)::mg(:)  ! multigrid level
 261     print *,'level ',mg(l)%level
 262     if(associated(mg(l)%X))print *,'X memory shape ',shape(mg(l)%X)
 263     if(associated(mg(l)%Y))print *,'Y memory shape ',shape(mg(l)%Y)
 264     if(associated(mg(l)%Z))print *,'Z memory shape ',shape(mg(l)%Z)
 265     print *,'grid size nx=',mg(l)%nx,' ny=',mg(l)%ny,' nz=',mg(l)%nz, &
 266        ' total ndof=',mg(l)%nn
 267     print *,'dx=',mg(l)%dx,' dy=',mg(l)%dy
 268     if(associated(mg(l)%dz))then
 269         print *,'dz shape',shape(mg(l)%dz)
 270         if(product(shape(mg(l)%dz))>0) then
 271             print *,'size ',size(mg(l)%dz)
 272             print *,'dz=',mg(l)%dz
 273         endif
 274     endif
 275 end subroutine print_mg_dims
 276
 277 end module module_common