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
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
 
character(len=*), intent(in)             :: filename
! .. local ..
character(len=1024)                      :: outfilehead
integer                                  :: nctot, ix^D, ixp^D
double precision                         :: time
integer                                  :: idim
integer,dimension(^ND)                   :: sendbuff
!-----------------------------------------------------------------------------
 
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
integer                                                 :: ic^D, ic1^D, ic2^D
double precision                                        :: xd^D
{^iftwod
double precision                                        :: c00, c10
}
{^ifthreed
double precision                                        :: c0, c1, c00, c10, c01, c11
}
integer                                                 :: ipivot^D, idir
!-----------------------------------------------------------------------------
 
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
!=============================================================================