MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_bc_data.t
Go to the documentation of this file.
1 !> Module to set boundary conditions from user data
2 module mod_bc_data
4 
5  implicit none
6  private
7 
8  integer, parameter :: max_boundary_conds = 10
9 
10  type bc_data_t
11  integer :: n_variables
12  double precision :: origin(3)
13  double precision :: dx(3)
14  integer :: n_points(3)
15  character(len=40), allocatable :: names(:)
16  double precision, allocatable :: values(:, :, :, :)
17  end type bc_data_t
18 
19  type(lt_t) :: lt_1d(max_boundary_conds)
20  type(lt2_t) :: lt_2d(max_boundary_conds)
21  type(lt3_t) :: lt_3d(max_boundary_conds)
22 
23  !> Whether boundary condition data is time varying
24  logical, public, protected :: bc_data_time_varying = .false.
25 
26  !> Integer array for indexing lookup tables per variable per direction
27  integer, public, protected, allocatable :: bc_data_ix(:, :)
28 
29  public :: bc_data_init
30  public :: bc_data_set
31  public :: bc_data_get_2d
32  public :: bc_data_get_3d
33 
34 contains
35 
36  subroutine bc_data_init()
38 
39  integer :: i, iw, ib, n_files, n_bc
40  character(len=std_len) :: bc_name, fname
41  double precision :: xmax(3)
42  type(bc_data_t) :: bc
43 
44  allocate(bc_data_ix(nwfluxbc, 2*ndim))
45 
46  bc_data_ix(:, :) = -1
47  n_bc = 0
48 
49  do ib = 1, 2 * ndim
50  do iw = 1, nwfluxbc
51  bc_name = typeboundary(iw, ib)
52  if (bc_name(1:4) == "vtk:") then
53 
54  n_bc = n_bc + 1
55  fname = bc_name(5:)
56  bc_data_ix(iw, ib) = n_bc
57 
58  ! Other routines don't have to parse the full name
59  typeboundary(iw, ib) = "bc_data"
60 
61  call read_vtk_structured_points(trim(fname), bc)
62  xmax = bc%origin + (bc%n_points-1) * bc%dx
63 
64  if (n_bc == 1) then
65  bc_data_time_varying = (bc%n_points(ndim) > 1)
66  else if (bc_data_time_varying .neqv. (bc%n_points(ndim) > 1)) then
67  call mpistop("bc_data_init: only some files are time varying")
68  end if
69 
70  {^ifoned
71  call mpistop("bc_data_init: 1D case not supported")
72  }
73  {^iftwod
74  if (bc_data_time_varying) then
75  lt_2d(n_bc) = lt2_create_from_data(bc%origin(1:ndim), &
76  xmax(1:ndim), bc%values(:, :, 1:1, 1))
77  else
78  ! Use first point in time
79  lt_1d(n_bc) = lt_create_from_data(bc%origin(1), &
80  xmax(1), bc%values(:, 1, 1:1, 1))
81  end if
82  }
83  {^ifthreed
84  if (bc_data_time_varying) then
85  lt_3d(n_bc) = lt3_create_from_data(bc%origin(1:ndim), &
86  xmax(1:ndim), bc%values(:, :, :, 1:1))
87  else
88  ! Use first point in time
89  lt_2d(n_bc) = lt2_create_from_data(bc%origin(1:ndim-1), &
90  xmax(1:ndim-1), bc%values(:, :, 1:1, 1))
91  end if
92  }
93  end if
94  end do
95  end do
96 
97  end subroutine bc_data_init
98 
99  elemental function bc_data_get_3d(n_bc, x1, x2, qt) result(val)
100  integer, intent(in) :: n_bc
101  double precision, intent(in) :: x1, x2, qt
102  double precision :: val
103 
104  if (bc_data_time_varying) then
105  val = lt3_get_col(lt_3d(n_bc), 1, x1, x2, qt)
106  else
107  val = lt2_get_col(lt_2d(n_bc), 1, x1, x2)
108  end if
109  end function bc_data_get_3d
110 
111  elemental function bc_data_get_2d(n_bc, x1, qt) result(val)
112  integer, intent(in) :: n_bc
113  double precision, intent(in) :: x1, qt
114  double precision :: val
115 
116  if (bc_data_time_varying) then
117  val = lt2_get_col(lt_2d(n_bc), 1, x1, qt)
118  else
119  val = lt_get_col(lt_1d(n_bc), 1, x1)
120  end if
121  end function bc_data_get_2d
122 
123  !> Set boundary conditions according to user data
124  subroutine bc_data_set(qt,ixI^L,ixO^L,iB,w,x)
126  integer, intent(in) :: ixI^L, ixO^L, iB
127  double precision, intent(in) :: qt, x(ixi^s,1:ndim)
128  double precision, intent(inout) :: w(ixi^s,1:nw)
129  double precision :: tmp(ixo^s)
130  integer :: i, ix, iw, n_bc
131 
132  {^iftwod
133  select case (ib)
134  case (1, 2)
135  do iw = 1, nwfluxbc
136  n_bc = bc_data_ix(iw, ib)
137  if (n_bc == -1) cycle
138 
139  tmp(ixomin1, ixomin2:ixomax2) = bc_data_get_2d(n_bc, &
140  x(ixomin1, ixomin2:ixomax2, 2), qt)
141 
142  if (ib == 1) then
143  ix = ixomax1+1
144  else
145  ix = ixomin1-1
146  end if
147 
148  ! Approximate boundary value by linear interpolation to first ghost
149  ! cell, rest of ghost cells contains the same value
150  do i = 0, ixomax1-ixomin1
151  w(ixomin1+i, ixomin2:ixomax2, iw) = &
152  2 * tmp(ixomin1, ixomin2:ixomax2) - &
153  w(ix, ixomin2:ixomax2, iw)
154  end do
155  end do
156  case (3, 4)
157  do iw = 1, nwfluxbc
158  n_bc = bc_data_ix(iw, ib)
159  if (n_bc == -1) cycle
160 
161  tmp(ixomin1:ixomax1, ixomin2) = bc_data_get_2d(n_bc, &
162  x(ixomin1:ixomax1, ixomin2, 1), qt)
163 
164  if (ib == 3) then
165  ix = ixomax2+1
166  else
167  ix = ixomin2-1
168  end if
169 
170  ! Approximate boundary value by linear interpolation to first ghost
171  ! cell, rest of ghost cells contains the same value
172  do i = 0, ixomax2-ixomin2
173  w(ixomin1:ixomax1, ixomin2+i, iw) = &
174  2 * tmp(ixomin1:ixomax1, ixomin2) - &
175  w(ixomin1:ixomax1, ix, iw)
176  end do
177  end do
178  case default
179  call mpistop("bc_data_set: unknown iB")
180  end select
181  }
182 
183  {^ifthreed
184  select case (ib)
185  case (1, 2)
186  do iw = 1, nwfluxbc
187  n_bc = bc_data_ix(iw, ib)
188  if (n_bc == -1) cycle
189 
190  tmp(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3) = bc_data_get_3d(n_bc, &
191  x(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3, 2), &
192  x(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3, 3), qt)
193 
194  if (ib == 1) then
195  ix = ixomax1+1
196  else
197  ix = ixomin1-1
198  end if
199 
200  ! Approximate boundary value by linear interpolation to first ghost
201  ! cell, rest of ghost cells contains the same value
202  do i = 0, ixomax1-ixomin1
203  w(ixomin1+i, ixomin2:ixomax2, ixomin3:ixomax3, iw) = &
204  2 * tmp(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3) - &
205  w(ix, ixomin2:ixomax2, ixomin3:ixomax3, iw)
206  end do
207  end do
208  case (3, 4)
209  do iw = 1, nwfluxbc
210  n_bc = bc_data_ix(iw, ib)
211  if (n_bc == -1) cycle
212 
213  tmp(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3) = bc_data_get_3d(n_bc, &
214  x(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3, 1), &
215  x(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3, 3), qt)
216 
217  if (ib == 3) then
218  ix = ixomax2+1
219  else
220  ix = ixomin2-1
221  end if
222 
223  ! Approximate boundary value by linear interpolation to first ghost
224  ! cell, rest of ghost cells contains the same value
225  do i = 0, ixomax2-ixomin2
226  w(ixomin1:ixomax1, ixomin2+i, ixomin3:ixomax3, iw) = &
227  2 * tmp(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3) - &
228  w(ixomin1:ixomax1, ix, ixomin3:ixomax3, iw)
229  end do
230  end do
231  case (5, 6)
232  do iw = 1, nwfluxbc
233  n_bc = bc_data_ix(iw, ib)
234  if (n_bc == -1) cycle
235 
236  tmp(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3) = bc_data_get_3d(n_bc, &
237  x(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3, 1), &
238  x(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3, 2), qt)
239 
240  if (ib == 5) then
241  ix = ixomax3+1
242  else
243  ix = ixomin3-1
244  end if
245 
246  ! Approximate boundary value by linear interpolation to first ghost
247  ! cell, rest of ghost cells contains the same value
248  do i = 0, ixomax3-ixomin3
249  w(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3+i, iw) = &
250  2 * tmp(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3) - &
251  w(ixomin1:ixomax1, ixomin2:ixomax2, ix, iw)
252  end do
253  end do
254  case default
255  call mpistop("bc_data_set: unknown iB")
256  end select
257  }
258 
259  end subroutine bc_data_set
260 
261  subroutine read_vtk_structured_points(fname, bc)
262  character(len=*), intent(in) :: fname
263  type(bc_data_t), intent(out) :: bc
264  double precision, allocatable :: tmp_data(:, :, :, :)
265  integer, parameter :: max_variables = 10
266  character(len=40) :: tmp_names(max_variables)
267  character(len=256) :: line
268  character(len=40) :: word, typename
269  integer, parameter :: my_unit = 123
270  integer :: n, n_points_total
271  integer :: n_components
272 
273  open(my_unit, file=trim(fname), status="old", action="read")
274 
275  ! Header, e.g. # vtk DataFile Version 2.0
276  read(my_unit, "(A)") line
277 
278  ! Dataset name
279  read(my_unit, "(A)") line
280 
281  ! ASCII / BINARY
282  read(my_unit, "(A)") line
283 
284  if (line /= "ASCII") then
285  print *, "line: ", trim(line)
286  error stop "read_vtk: not an ASCII file"
287  end if
288 
289  ! DATASET STRUCTURED_POINTS
290  read(my_unit, "(A)") line
291 
292  if (line /= "DATASET STRUCTURED_POINTS") then
293  print *, "line: ", trim(line)
294  error stop "read_vtk must have: DATASET STRUCTURED_POINTS"
295  end if
296 
297  ! DIMENSIONS NX NY NZ
298  read(my_unit, "(A)") line
299  read(line, *) word, bc%n_points
300 
301  if (word /= "DIMENSIONS") then
302  print *, "line: ", trim(line)
303  error stop "read_vtk expects: DIMENSIONS"
304  end if
305 
306  ! SPACING DX DY DZ
307  read(my_unit, *) word, bc%dx
308 
309  if (word /= "SPACING") then
310  print *, "line: ", trim(line)
311  error stop "read_vtk expects: SPACING"
312  end if
313 
314  ! ORIGIN 0 0 0
315  read(my_unit, *) word, bc%origin
316  if (word /= "ORIGIN") then
317  print *, "line: ", trim(line)
318  error stop "read_vtk expects: ORIGIN"
319  end if
320 
321  ! POINT_DATA NPOINTS
322  read(my_unit, *) word, n_points_total
323 
324  if (word /= "POINT_DATA") then
325  print *, "line: ", trim(line)
326  error stop "read_vtk expects: POINT_DATA n_points"
327  end if
328 
329  if (n_points_total /= product(bc%n_points)) &
330  error stop "read_vtk: n_points not equal to NX*NY*NZ"
331 
332  allocate(tmp_data(bc%n_points(1), bc%n_points(2), bc%n_points(3), &
333  max_variables))
334 
335  ! Read all scalar variables
336  do n = 1, max_variables
337 
338  ! SCALARS name type ncomponents
339  read(my_unit, *, end=900) word, tmp_names(n), typename, n_components
340 
341  if (word /= "SCALARS") then
342  print *, "line: ", trim(line)
343  error stop "read_vtk expects: SCALARS name type ncomponents"
344  end if
345 
346  if (n_components /= 1) error stop "read_vtk: ncomponents should be 1"
347 
348  ! LOOKUP_TABLE default
349  read(my_unit, *) word, typename
350 
351  if (word /= "LOOKUP_TABLE") then
352  print *, "line: ", trim(line)
353  error stop "read_vtk expects: LOOKUP_TABLE name"
354  end if
355 
356  ! Read list of data values
357  read(my_unit, *) tmp_data(:, :, :, n)
358  end do
359 
360  ! Done reading variables
361 900 continue
362 
363  close(my_unit)
364 
365  if (n == max_variables + 1) &
366  error stop "read_vtk: increase max_variables"
367 
368  ! Loop index is one higher than number of variables
369  bc%n_variables = n-1
370  bc%values = tmp_data(:, :, :, 1:n-1)
371  bc%names = tmp_names(1:n-1)
372  end subroutine read_vtk_structured_points
373 
374 end module mod_bc_data
This module contains definitions of global parameters and variables and some generic functions/subrou...
elemental real(dp) function, public lt_get_col(my_lt, col_ix, x)
Get the value of a single column at x.
elemental double precision function, public bc_data_get_3d(n_bc, x1, x2, qt)
Definition: mod_bc_data.t:100
type(lt3_t) function, public lt3_create_from_data(x_min, x_max, spaced_data)
This function returns a new lookup table from regularly spaced data.
subroutine read_vtk_structured_points(fname, bc)
Definition: mod_bc_data.t:262
subroutine, public bc_data_set(qt, ixIL, ixOL, iB, w, x)
Set boundary conditions according to user data.
Definition: mod_bc_data.t:125
A Fortran 90 module for creating 1D and 2D lookup tables. These tables can be used to efficiently int...
Module to set boundary conditions from user data.
Definition: mod_bc_data.t:2
type(lt_t) function, public lt_create_from_data(x_min, x_max, spaced_data)
This function returns a new lookup table from regularly spaced data.
subroutine, public bc_data_init()
Definition: mod_bc_data.t:37
character(len=std_len), dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
elemental real(dp) function, public lt2_get_col(my_lt, col_ix, x1, x2)
Get the value of a single column at x.
logical, public, protected bc_data_time_varying
Whether boundary condition data is time varying.
Definition: mod_bc_data.t:24
elemental real(dp) function, public lt3_get_col(my_lt, col_ix, x1, x2, x3)
Get the value of a single column at x.
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, parameter ndim
Number of spatial dimensions for grid variables.
elemental double precision function, public bc_data_get_2d(n_bc, x1, qt)
Definition: mod_bc_data.t:112
integer, dimension(:, :), allocatable, public, protected bc_data_ix
Integer array for indexing lookup tables per variable per direction.
Definition: mod_bc_data.t:27
type(lt2_t) function, public lt2_create_from_data(x_min, x_max, spaced_data)
This function returns a new lookup table from regularly spaced data.