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