mod_oneblock.t

Go to the documentation of this file.

!=============================================================================
! Module oneblock to hold a datastructure as output by convert_type='oneblock'. 
! Currently only ASCII
! Can be handy to read in initial conditions.
! 2015-02-18 Oliver Porth
!=============================================================================
module mod_oneblock
  use mod_comm_lib, only: mpistop
  implicit none
  save
 
  double precision, dimension({^D&:|,},:), allocatable       :: woneblock
  double precision, dimension({^D&:|,},:), allocatable       :: xoneblock
  integer                                                    :: nc^d
  integer                                                    :: unit=15 !file unit to read on
 
contains
 
  subroutine read_oneblock(filename)
  use mod_global_parameters
 
  double precision                         :: time
  integer                                  :: nctot, ix^D, ixp^D
  integer                                  :: idim
  integer,dimension(^ND)                   :: sendbuff
  character(len=*), intent(in)             :: filename
  character(len=1024)                      :: outfilehead
 
  if (mype == 0) write(*,*) 'mod_oneblock: reading ',nw,' variables on unit:', unit
 
  !----------------------------------------
  ! Root does the reading:
  !----------------------------------------
  if (mype == 0) then 
     open(unit,file=filename,status='unknown')
     ! The header information:
     read(unit,'(A)') outfilehead
     read(unit,*) nctot,nc^d
     read(unit,*) time
     
     ! Allocate and read the grid and variables:
     allocate(xoneblock(nc^d,1:^nd))
     allocate(woneblock(nc^d,1:nw))
     
     {do ix^db=1,nc^db\}  
     
     read(unit,*) xoneblock(ix^d,1:^nd), woneblock(ix^d,1:nw)
     
     {end do\}
     
  ! Close the file
     close(unit)
  end if! mype==0
 
  !----------------------------------------
  ! Broadcast what mype=0 read:
  !----------------------------------------
  if (npe>1) then 
     {sendbuff(^d)=nc^d;}
     call mpi_bcast(sendbuff,^nd,mpi_integer,0,icomm,ierrmpi)
     if (mype .ne. 0) then 
        {nc^d=sendbuff(^d);}
        ! Allocate the grid and variables:
        allocate(xoneblock(nc^d,1:^nd))
        allocate(woneblock(nc^d,1:nw))
     end if
     call mpi_bcast(xoneblock,{nc^d|*}*^nd,mpi_double_precision,0,icomm,ierrmpi)
     call mpi_bcast(woneblock,{nc^d|*}*nw,mpi_double_precision,0,icomm,ierrmpi)
  end if! npe>1
 
  end subroutine read_oneblock
 
  subroutine interpolate_oneblock(x,iw,out)
 
  double precision, dimension(^ND),intent(in)             :: x 
  integer, intent(in)                                     :: iw
  double precision, intent(out)                           :: out
  ! .. local ..
  double precision, dimension(^ND)                        :: xloc
  double precision                                        :: xd^D
  {^iftwod
  double precision                                        :: c00, c10
  }
  {^ifthreed
  double precision                                        :: c0, c1, c00, c10, c01, c11
  }
  integer                                                 :: ipivot^D, idir
  integer                                                 :: ic^D, ic1^D, ic2^D
 
  xloc=x
 
  !--------------------------------------------
  ! Hunt for the index closest to the point
  ! This is a bit slow but allows for stretched grids
  ! (still need to be orthogonal for interpolation though)
  !--------------------------------------------
  ipivot^d=1;
  {^ifoned
  ic1 = minloc(dabs(xloc(1)-xoneblock(:,1)),1,mask=.true.)
  }
  {^iftwod
  do idir = 1, ^nd
     select case (idir)
     case (1)
        ic1 = minloc(dabs(xloc(1)-xoneblock(:,ipivot2,idir)),1,mask=.true.)
     case (2)
        ic2 = minloc(dabs(xloc(2)-xoneblock(ipivot1,:,idir)),1,mask=.true.)
     case default
        call mpistop("error1 in interpolate_oneblock")
     end select
  end do
  }
  {^ifthreed
  do idir = 1, ^nd
     select case (idir)
     case (1)
        ic1 = minloc(dabs(xloc(1)-xoneblock(:,ipivot2,ipivot3,idir)),1, &
             mask=.true.)
     case (2)
        ic2 = minloc(dabs(xloc(2)-xoneblock(ipivot1,:,ipivot3,idir)),1, &
             mask=.true.)
     case (3)
        ic3 = minloc(dabs(xloc(3)-xoneblock(ipivot1,ipivot2,:,idir)),1, &
             mask=.true.)
     case default
        call mpistop("error1 in interpolate_oneblock")
     end select
  end do
  }
 
  ! flat interpolation would simply be:
  !out = woneblock(ic^D,iw)
  !return
 
  !-------------------------------------------
  ! Get the left and right indices
  !-------------------------------------------
  {
  if (xoneblock({ic^dd},^d) .lt. xloc(^d)) then
     ic1^d = ic^d
  else
     ic1^d = ic^d -1
  end if
  ic2^d = ic1^d + 1
  \}
 
  !--------------------------------------------
  ! apply flat interpolation if outside of range, 
  ! change point-location to make this easy!
  !--------------------------------------------
  {
  if (ic1^d .lt. 1) then
     ic1^d = 1
     ic2^d = ic1^d + 1
     xloc(^d) = xoneblock(ic^dd,^d)
  end if
  \}
  {
  if (ic2^d .gt. nc^d) then
     ic2^d = nc^d
     ic1^d = ic2^d - 1
     xloc(^d) = xoneblock(ic^dd,^d)
  end if
  \}
 
 
  !-------------------------------------------
  ! linear, bi- and tri- linear interpolations
  !-------------------------------------------
  {^ifoned
  xd1 = (xloc(1)-xoneblock(ic11,1)) / (xoneblock(ic21,1) - xoneblock(ic11,1))
  out = woneblock(ic11,iw) * (1.0d0 - xd1) + woneblock(ic21,iw) * xd1
  }
  {^iftwod
  xd1 = (xloc(1)-xoneblock(ic11,ic12,1)) / (xoneblock(ic21,ic12,1) - xoneblock(ic11,ic12,1))      
  xd2 = (xloc(2)-xoneblock(ic11,ic12,2)) / (xoneblock(ic11,ic22,2) - xoneblock(ic11,ic12,2))
  c00 = woneblock(ic11,ic12,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,iw) * xd1
  c10 = woneblock(ic11,ic22,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,iw) * xd1
  out = c00 * (1.0d0 - xd2) + c10 * xd2
  }
  {^ifthreed
  xd1 = (xloc(1)-xoneblock(ic11,ic12,ic13,1)) / (xoneblock(ic21,ic12,ic13,1) - xoneblock(ic11,ic12,ic13,1))      
  xd2 = (xloc(2)-xoneblock(ic11,ic12,ic13,2)) / (xoneblock(ic11,ic22,ic13,2) - xoneblock(ic11,ic12,ic13,2))      
  xd3 = (xloc(3)-xoneblock(ic11,ic12,ic13,3)) / (xoneblock(ic11,ic12,ic23,3) - xoneblock(ic11,ic12,ic13,3))    
 
  c00 = woneblock(ic11,ic12,ic13,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,ic13,iw) * xd1
  c10 = woneblock(ic11,ic22,ic13,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,ic13,iw) * xd1
  c01 = woneblock(ic11,ic12,ic23,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,ic23,iw) * xd1
  c11 = woneblock(ic11,ic22,ic23,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,ic23,iw) * xd1
 
  c0  = c00 * (1.0d0 - xd2) + c10 * xd2
  c1  = c01 * (1.0d0 - xd2) + c11 * xd2
 
  out = c0 * (1.0d0 - xd3) + c1 * xd3
  }
 
  end subroutine interpolate_oneblock
 
end module mod_oneblock