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, allocatable :: woneblock({^d&:|,},:)
  double precision, allocatable :: xoneblock({^d&:|,},:)
 
  integer, private :: nc^d
  integer, private :: unit=15 !file unit to read on
 
contains
 
  subroutine read_oneblock(filename)
    use mod_global_parameters
 
    character(len=*), intent(in) :: filename
 
    double precision     :: time
    integer              :: ix^D, ixp^D, idim
    integer              :: sendbuff(^ND)
    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,*) 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 /= 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, intent(in)  :: x(^ND)
    integer,          intent(in)  :: iw
    double precision, intent(out) :: out
    ! .. local ..
    double precision :: xloc(^ND)
    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("error 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("error 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) < 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 < 1) then
      ic1^d = 1
      ic2^d = ic1^d + 1
      xloc(^d) = xoneblock(ic^dd,^d)
    end if
    \}
    {
    if (ic2^d > 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