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-6 ! residual tolerance to stop iterating
  23     !logical:: check_reldif=.false.  ! check difference after every base iteration if to continue
  24     logical:: check_reldif=.true.  ! check difference after every base iteration if to continue
  25     real:: diftol_finest=1e-6 ! keep iterating while iterates change by this
  26     real:: diftol_coarse=1e-1 ! keep iterating while iterates change by this on corse levels
  27     real:: diftol_coarsest=1e-2 ! keep iterating while iterates change by this on corsest level
  28 end type
  29
  30 type(params_type)::params
  31
  32 ! multigrid structure
  33
  34 type mg_type
  35     real, dimension(3,3):: A                        ! penalty weight matrix
  36     real, pointer, dimension(:,:,:):: X, Y, Z       ! grid vertices
  37     real, pointer, dimension(:,:,:)::        &      !
  38          f, lambda, res                             ! multigrid variables
  39
  40     real, pointer, dimension(:,:,:,:):: Kglo        ! global stiffness matrix
  41
  42     real:: dx,dy                                    ! horizontal spacing, scalar
  43     real, pointer, dimension(:)::dz                 ! vertical spacing of the layers
  44     integer::nx,ny,nz,nn                            ! mesh size in terms of vertices
  45
  46     integer::level
  47
  48     integer:: cr_x, cr_y                            ! coarsening factors
  49     integer, pointer, dimension(:):: icl_x, icl_y, icl_z  ! coarsening indices
  50
  51 end type
  52
  53 type(mg_type):: mg(max_levels+1)  ! the main multigrid structure
  54
  55 contains
  56
  57 subroutine allocate_mg_level(mg)    !allocate one multigrid level
  58 type(mg_type),intent(inout)::mg     !multigrid structure
  59 !*** local
  60 integer:: &
  61     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
  62     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
  63     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  64     ifts, ifte, kfts, kfte, jfts, jfte              ! tile dimensions
  65
  66 !*** executable
  67 call get_mg_dims(mg, &
  68     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
  69     ifms, ifme, kfms,kfme, jfms, jfme,            &
  70     ifps, ifpe, kfps,kfpe, jfps, jfpe,            & ! fire patch bounds
  71     ifts, ifte, kfts,kfte, jfts,jfte)
  72
  73 if (.not.associated(mg%X))allocate(mg%X(ifms: ifme, kfms: kfme, jfms: jfme))
  74 if (.not.associated(mg%Y))allocate(mg%Y(ifms: ifme, kfms: kfme, jfms: jfme))
  75 if (.not.associated(mg%Z))allocate(mg%Z(ifms: ifme, kfms: kfme, jfms: jfme))
  76 allocate(mg%F(ifms: ifme, kfms: kfme, jfms: jfme))
  77 allocate(mg%lambda(ifms: ifme, kfms: kfme, jfms: jfme))
  78 allocate(mg%res(ifms: ifme, kfms: kfme, jfms: jfme))
  79
  80 end subroutine allocate_mg_level
  81
  82 subroutine get_mg_dims(mg_level, &
  83     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
  84     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
  85     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  86     ifts, ifte, kfts, kfte, jfts, jfte)              ! tile dimensions
  87 implicit none
  88 type(mg_type),intent(in)::mg_level
  89 integer, intent(out)::  &
  90     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
  91     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
  92     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
  93     ifts, ifte, kfts, kfte, jfts, jfte              ! tile dimensions
  94 character(len=256)::msg
  95
  96     ifds = 1
  97     ifde = mg_level%nx-1  ! dimension in elements
  98     jfds = 1
  99     jfde = mg_level%ny-1
 100     kfds = 1
 101     kfde = mg_level%nz-1
 102
 103     ifts = ifds
 104     ifte = ifde
 105     jfts = jfds
 106     jfte = jfde
 107     kfts = kfds
 108     kfte = kfde
 109
 110     ifps = ifds
 111     ifpe = ifde
 112     jfps = jfds
 113     jfpe = jfde
 114     kfps = kfds
 115     kfpe = kfde
 116
 117     ifms = ifps - 1
 118     ifme = ifpe + 2  ! need +1 for vertex grid, and +1 for zero strip
 119     jfms = jfps - 1
 120     jfme = jfpe + 2
 121     kfms = kfps - 1
 122     kfme = kfpe + 2
 123 write(msg,*)"get_mg_dims: domain",ifds, ifde, kfds, kfde, jfds, jfde
 124 call message(msg)
 125 write(msg,*)"get_mg_dims: memory",ifms, ifme, kfms, kfme, jfms, jfme
 126 call message(msg)
 127 write(msg,*)"get_mg_dims: patch ",ifps, ifde, kfps, kfpe, jfps, jfpe
 128 call message(msg)
 129 write(msg,*)"get_mg_dims: tile  ",ifts, ifte, kfts, kfte, jfts, jfte
 130 call message(msg)
 131
 132 end subroutine get_mg_dims
 133
 134 subroutine write_3d(mg,l,var,name,v)
 135 ! write array in text file for debugging
 136 implicit none
 137
 138 !*** arguments
 139 type(mg_type),intent(in)::mg(:)    !multigrid structure
 140 integer,intent(in)::l              !level
 141 real, intent(in)::var(:,:,:)       !the array
 142 character(len=*),intent(in)::name
 143 character(len=1),intent(in)::v     !'v' or 'e'
 144
 145 !*** local
 146 integer:: &
 147     ifds, ifde, kfds, kfde, jfds, jfde,           & ! fire grid dimensions
 148     ifms, ifme, kfms, kfme, jfms, jfme,           & ! memory dimensions
 149     ifps, ifpe, kfps, kfpe, jfps, jfpe,           & ! fire patch bounds
 150     ifts, ifte, kfts, kfte, jfts, jfte              ! tile dimensions
 151 integer::ie,ke,je
 152 character(len=2)::lc
 153
 154 !*** executable
 155
 156 call get_mg_dims(mg(l), &
 157     ifds, ifde, kfds,kfde, jfds, jfde,            & ! fire grid dimensions
 158     ifms, ifme, kfms,kfme, jfms, jfme,            &
 159     ifps, ifpe, kfps,kfpe, jfps, jfpe,           & ! fire patch bounds
 160     ifts, ifte, kfts,kfte, jfts,jfte)
 161
 162 write(lc,'(i2.2)')l
 163
 164 if(v.eq.'v')then
 165     ie = snode(ifte,ifde,+1)  ! vertex based
 166     je = snode(jfte,jfde,+1)
 167     ke = snode(kfte,kfde,+1)
 168 elseif(v.eq.'e')then
 169     ie=ifte
 170     je=jfte
 171     ke=kfte
 172 else
 173     call crash('write_3d: bad v')
 174 endif
 175
 176 call write_tight(var) ! write dereferenced
 177
 178 contains
 179
 180 subroutine write_tight(varr)
 181     real, intent(in)::varr(ifms:ifme,kfms:kfme,jfms:jfme)
 182     ! for element/cell/midpoint based arrays, like fmw winds
 183     ! open(iu,file='varr.dat',form='unformatted',status='unknown')
 184     ! write(iu)varr(ifts:ifte,kfts:kfte,jfts:jfte)
 185     ! close(iu)
 186     call write_array(varr(ifts:ie,kfts:ke,jfts:je),name//lc)
 187 end subroutine write_tight
 188
 189 end subroutine write_3d
 190
 191 subroutine print_mg_dims(mg,l)
 192 type(mg_type),intent(in)::mg(:)  ! multigrid level
 193     print *,'level ',mg(l)%level
 194     if(associated(mg(l)%X))print *,'X memory shape ',shape(mg(l)%X)
 195     if(associated(mg(l)%Y))print *,'Y memory shape ',shape(mg(l)%Y)
 196     if(associated(mg(l)%Z))print *,'Z memory shape ',shape(mg(l)%Z)
 197     print *,'grid size nx=',mg(l)%nx,' ny=',mg(l)%ny,' nz=',mg(l)%nz, &
 198        ' total ndof=',mg(l)%nn
 199     print *,'dx=',mg(l)%dx,' dy=',mg(l)%dy
 200     if(associated(mg(l)%dz))then
 201         print *,'dz shape',shape(mg(l)%dz)
 202         if(product(shape(mg(l)%dz))>0) then
 203             print *,'size ',size(mg(l)%dz)
 204             print *,'dz=',mg(l)%dz
 205         endif
 206     endif
 207 end subroutine print_mg_dims
 208
 209 end module module_common