MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_oneblock.t
Go to the documentation of this file.
1 !=============================================================================
2 ! Module oneblock to hold a datastructure as output by convert_type='oneblock'.
3 ! Currently only ASCII
4 ! Can be handy to read in initial conditions.
5 ! 2015-02-18 Oliver Porth
6 !=============================================================================
8  use mod_comm_lib, only: mpistop
9  implicit none
10  save
11 
12  double precision, dimension({^D&:|,},:), allocatable :: woneblock
13  double precision, dimension({^D&:|,},:), allocatable :: xoneblock
14  integer :: nc^d
15  integer :: unit=15 !file unit to read on
16 
17 contains
18 
19  subroutine read_oneblock(filename)
21 
22  character(len=*), intent(in) :: filename
23  character(len=1024) :: outfilehead
24  integer :: nctot, ix^D, ixp^D
25  double precision :: time
26  integer :: idim
27  integer,dimension(^ND) :: sendbuff
28 
29  if (mype == 0) write(*,*) 'mod_oneblock: reading ',nw,' variables on unit:', unit
30 
31  !----------------------------------------
32  ! Root does the reading:
33  !----------------------------------------
34  if (mype == 0) then
35  open(unit,file=filename,status='unknown')
36  ! The header information:
37  read(unit,'(A)') outfilehead
38  read(unit,*) nctot,nc^d
39  read(unit,*) time
40 
41  ! Allocate and read the grid and variables:
42  allocate(xoneblock(nc^d,1:^nd))
43  allocate(woneblock(nc^d,1:nw))
44 
45  {do ix^db=1,nc^db\}
46 
47  read(unit,*) xoneblock(ix^d,1:^nd), woneblock(ix^d,1:nw)
48 
49  {end do\}
50 
51  ! Close the file
52  close(unit)
53  end if! mype==0
54 
55  !----------------------------------------
56  ! Broadcast what mype=0 read:
57  !----------------------------------------
58  if (npe>1) then
59  {sendbuff(^d)=nc^d;}
60  call mpi_bcast(sendbuff,^nd,mpi_integer,0,icomm,ierrmpi)
61  if (mype .ne. 0) then
62  {nc^d=sendbuff(^d);}
63  ! Allocate the grid and variables:
64  allocate(xoneblock(nc^d,1:^nd))
65  allocate(woneblock(nc^d,1:nw))
66  end if
67  call mpi_bcast(xoneblock,{nc^d|*}*^nd,mpi_double_precision,0,icomm,ierrmpi)
68  call mpi_bcast(woneblock,{nc^d|*}*nw,mpi_double_precision,0,icomm,ierrmpi)
69  end if! npe>1
70 
71  end subroutine read_oneblock
72 
73  subroutine interpolate_oneblock(x,iw,out)
74 
75  double precision, dimension(^ND),intent(in) :: x
76  integer, intent(in) :: iw
77  double precision, intent(out) :: out
78  ! .. local ..
79  double precision, dimension(^ND) :: xloc
80  integer :: ic^D, ic1^D, ic2^D
81  double precision :: xd^D
82  {^iftwod
83  double precision :: c00, c10
84  }
85  {^ifthreed
86  double precision :: c0, c1, c00, c10, c01, c11
87  }
88  integer :: ipivot^D, idir
89 
90  xloc=x
91 
92  !--------------------------------------------
93  ! Hunt for the index closest to the point
94  ! This is a bit slow but allows for stretched grids
95  ! (still need to be orthogonal for interpolation though)
96  !--------------------------------------------
97  ipivot^d=1;
98  {^ifoned
99  ic1 = minloc(dabs(xloc(1)-xoneblock(:,1)),1,mask=.true.)
100  }
101  {^iftwod
102  do idir = 1, ^nd
103  select case (idir)
104  case (1)
105  ic1 = minloc(dabs(xloc(1)-xoneblock(:,ipivot2,idir)),1,mask=.true.)
106  case (2)
107  ic2 = minloc(dabs(xloc(2)-xoneblock(ipivot1,:,idir)),1,mask=.true.)
108  case default
109  call mpistop("error1 in interpolate_oneblock")
110  end select
111  end do
112  }
113  {^ifthreed
114  do idir = 1, ^nd
115  select case (idir)
116  case (1)
117  ic1 = minloc(dabs(xloc(1)-xoneblock(:,ipivot2,ipivot3,idir)),1, &
118  mask=.true.)
119  case (2)
120  ic2 = minloc(dabs(xloc(2)-xoneblock(ipivot1,:,ipivot3,idir)),1, &
121  mask=.true.)
122  case (3)
123  ic3 = minloc(dabs(xloc(3)-xoneblock(ipivot1,ipivot2,:,idir)),1, &
124  mask=.true.)
125  case default
126  call mpistop("error1 in interpolate_oneblock")
127  end select
128  end do
129  }
130 
131  ! flat interpolation would simply be:
132  !out = woneblock(ic^D,iw)
133  !return
134 
135  !-------------------------------------------
136  ! Get the left and right indices
137  !-------------------------------------------
138  {
139  if (xoneblock({ic^dd},^d) .lt. xloc(^d)) then
140  ic1^d = ic^d
141  else
142  ic1^d = ic^d -1
143  end if
144  ic2^d = ic1^d + 1
145  \}
146 
147  !--------------------------------------------
148  ! apply flat interpolation if outside of range,
149  ! change point-location to make this easy!
150  !--------------------------------------------
151  {
152  if (ic1^d .lt. 1) then
153  ic1^d = 1
154  ic2^d = ic1^d + 1
155  xloc(^d) = xoneblock(ic^dd,^d)
156  end if
157  \}
158  {
159  if (ic2^d .gt. nc^d) then
160  ic2^d = nc^d
161  ic1^d = ic2^d - 1
162  xloc(^d) = xoneblock(ic^dd,^d)
163  end if
164  \}
165 
166 
167  !-------------------------------------------
168  ! linear, bi- and tri- linear interpolations
169  !-------------------------------------------
170  {^ifoned
171  xd1 = (xloc(1)-xoneblock(ic11,1)) / (xoneblock(ic21,1) - xoneblock(ic11,1))
172  out = woneblock(ic11,iw) * (1.0d0 - xd1) + woneblock(ic21,iw) * xd1
173  }
174  {^iftwod
175  xd1 = (xloc(1)-xoneblock(ic11,ic12,1)) / (xoneblock(ic21,ic12,1) - xoneblock(ic11,ic12,1))
176  xd2 = (xloc(2)-xoneblock(ic11,ic12,2)) / (xoneblock(ic11,ic22,2) - xoneblock(ic11,ic12,2))
177  c00 = woneblock(ic11,ic12,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,iw) * xd1
178  c10 = woneblock(ic11,ic22,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,iw) * xd1
179  out = c00 * (1.0d0 - xd2) + c10 * xd2
180  }
181  {^ifthreed
182  xd1 = (xloc(1)-xoneblock(ic11,ic12,ic13,1)) / (xoneblock(ic21,ic12,ic13,1) - xoneblock(ic11,ic12,ic13,1))
183  xd2 = (xloc(2)-xoneblock(ic11,ic12,ic13,2)) / (xoneblock(ic11,ic22,ic13,2) - xoneblock(ic11,ic12,ic13,2))
184  xd3 = (xloc(3)-xoneblock(ic11,ic12,ic13,3)) / (xoneblock(ic11,ic12,ic23,3) - xoneblock(ic11,ic12,ic13,3))
185 
186  c00 = woneblock(ic11,ic12,ic13,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,ic13,iw) * xd1
187  c10 = woneblock(ic11,ic22,ic13,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,ic13,iw) * xd1
188  c01 = woneblock(ic11,ic12,ic23,iw) * (1.0d0 - xd1) + woneblock(ic21,ic12,ic23,iw) * xd1
189  c11 = woneblock(ic11,ic22,ic23,iw) * (1.0d0 - xd1) + woneblock(ic21,ic22,ic23,iw) * xd1
190 
191  c0 = c00 * (1.0d0 - xd2) + c10 * xd2
192  c1 = c01 * (1.0d0 - xd2) + c11 * xd2
193 
194  out = c0 * (1.0d0 - xd3) + c1 * xd3
195  }
196 
197  end subroutine interpolate_oneblock
198 
199 end module mod_oneblock
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer mype
The rank of the current MPI task.
integer unit
Definition: mod_oneblock.t:15
integer nc
Definition: mod_oneblock.t:14
subroutine read_oneblock(filename)
Definition: mod_oneblock.t:20
double precision, dimension({^d &:|,},:), allocatable xoneblock
Definition: mod_oneblock.t:13
subroutine interpolate_oneblock(x, iw, out)
Definition: mod_oneblock.t:74
double precision, dimension({^d &:|,},:), allocatable woneblock
Definition: mod_oneblock.t:12