MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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  logical, public, protected :: interp_phy_first_row=.true.
33  logical, public, protected :: bc_phy_first_row=.false.
34  logical, public, protected :: boundary_data_primitive=.false.
35 
36  public :: bc_data_init
37  public :: bc_data_set
38  public :: bc_data_get_2d
39  public :: bc_data_get_3d
40 
41 contains
42 
43  subroutine bc_data_init()
45  use mod_comm_lib, only: mpistop
46 
47  integer :: i, iw, ib, n_files, n_bc
48  double precision :: xmax(3)
49  type(bc_data_t) :: bc
50 
52 
53  allocate(bc_data_ix(nwfluxbc, 2*ndim))
54 
55  bc_data_ix(:, :) = -1
56  n_bc = 0
57 
58  if(any(typeboundary(1:nwfluxbc,1:2*ndim)==bc_data ) .or. any(typeboundary(1:nwfluxbc,1:2*ndim)==bc_icarus)) then
60  xmax = bc%origin + (bc%n_points-1) * bc%dx
61  bc_data_time_varying = (bc%n_points(ndim) > 1)
62 
63  endif
64 
65  do ib = 1, 2 * ndim
66  do iw = 1, nwfluxbc
67  if (typeboundary(iw, ib)==bc_data .or. typeboundary(iw, ib)==bc_icarus) then
68  n_bc = n_bc + 1
69  bc_data_ix(iw, ib) = n_bc
70 
71 
72  {^ifoned
73  call mpistop("bc_data_init: 1D case not supported")
74  }
75  {^iftwod
76  if (bc_data_time_varying) then
77  lt_2d(n_bc) = lt2_create_from_data(bc%origin(1:ndim), &
78  xmax(1:ndim), bc%values(:, :, 1:1, n_bc))
79  else
80  ! Use first point in time
81  lt_1d(n_bc) = lt_create_from_data(bc%origin(1), &
82  xmax(1), bc%values(:, 1, 1:1, n_bc))
83  end if
84  }
85  {^ifthreed
86  if (bc_data_time_varying) then
87  lt_3d(n_bc) = lt3_create_from_data(bc%origin(1:ndim), &
88  xmax(1:ndim), bc%values(:, :, :, n_bc:n_bc))
89  else
90  ! Use first point in time
91  lt_2d(n_bc) = lt2_create_from_data(bc%origin(1:ndim-1), &
92  xmax(1:ndim-1), bc%values(:, :, 1:1, n_bc))
93  end if
94  }
95  end if
96  end do
97  end do
98 
99  end subroutine bc_data_init
100 
101  !> Read this module"s parameters from a file
102  subroutine bc_read_params(files)
104  character(len=*), intent(in) :: files(:)
105  integer :: n
106 
109 
110  do n = 1, size(files)
111  open(unitpar, file=trim(files(n)), status="old")
112  read(unitpar, bd_list, end=111)
113 111 close(unitpar)
114  end do
115 
116  end subroutine bc_read_params
117 
118  elemental function bc_data_get_3d(n_bc, x1, x2, qt) result(val)
119  integer, intent(in) :: n_bc
120  double precision, intent(in) :: x1, x2, qt
121  double precision :: val
122 
123  if (bc_data_time_varying) then
124  val = lt3_get_col(lt_3d(n_bc), 1, x1, x2, qt)
125  else
126  val = lt2_get_col(lt_2d(n_bc), 1, x1, x2)
127  end if
128  end function bc_data_get_3d
129 
130  elemental function bc_data_get_2d(n_bc, x1, qt) result(val)
131  integer, intent(in) :: n_bc
132  double precision, intent(in) :: x1, qt
133  double precision :: val
134 
135  if (bc_data_time_varying) then
136  val = lt2_get_col(lt_2d(n_bc), 1, x1, qt)
137  else
138  val = lt_get_col(lt_1d(n_bc), 1, x1)
139  end if
140  end function bc_data_get_2d
141 
142  !> Set boundary conditions according to user data
143  subroutine bc_data_set(qt,ixI^L,ixO^L,iB,w,x)
146  use mod_comm_lib, only: mpistop
147 
148  integer, intent(in) :: ixi^l, ixo^l, ib
149  double precision, intent(in) :: qt, x(ixi^s,1:ndim)
150  double precision, intent(inout) :: w(ixi^s,1:nw)
151  double precision :: tmp(ixo^s)
152  integer :: i, ix, iw, n_bc
153 
154  integer :: ixoo^l
155 
156  ixoo^l=ixo^l;
157 
158  {^iftwod
159  select case (ib)
160  case (1, 2)
161  if(interp_phy_first_row) then
162  if (ib == 1) then
163  ix = ixomax1+1
164  else
165  ix = ixomin1-1
166  end if
167  if(boundary_data_primitive) then
168  call phys_to_primitive(ixi^l,ix,ixomin2,ix,ixomax2,w,x)
169  endif
170  endif
171  ! the reason for this is that not all bc for all the vars
172  ! might be set with bc_data
173  if(bc_phy_first_row)then
174  if (ib == 1) then
175  ixoomax1=ixoomax1+1
176  else
177  ixoomin1=ixoomin1-1
178  endif
179  if(boundary_data_primitive) then
180  if (ib == 1) then
181  ixoomin1=ixoomax1
182  else
183  ixoomax1=ixoomin1
184  endif
185  call phys_to_primitive(ixi^l,ixoo^l,w,x)
186  if (ib == 1) then
187  ixoomin1=ixomin1 ! to be used later in to_conserved
188  else
189  ixoomax1=ixomax1
190  endif
191  endif
192  endif
193  do iw = 1, nwfluxbc
194  n_bc = bc_data_ix(iw, ib)
195  if (n_bc == -1) cycle
196 
197  tmp(ixomin1, ixomin2:ixomax2) = bc_data_get_2d(n_bc, &
198  x(ixomin1, ixomin2:ixomax2, 2), qt)
199  if(interp_phy_first_row) then
200  ! Approximate boundary value by linear interpolation to first ghost
201  ! cell, rest of ghost cells contains the same value
202  do i = 0, ixoomax1-ixoomin1
203  w(ixoomin1+i, ixomin2:ixomax2, iw) = &
204  2 * tmp(ixomin1, ixomin2:ixomax2) - &
205  w(ix, ixomin2:ixomax2, iw)
206  end do
207  else
208  do i = ixoomin1, ixoomax1
209  w(i, ixomin2:ixomax2, iw) = &
210  tmp(ixomin1, ixomin2:ixomax2)
211  end do
212 
213  endif
214  end do
215  if(interp_phy_first_row) then
216  if(boundary_data_primitive) then
217  call phys_to_conserved(ixi^l,ix,ixomin2,ix,ixomax2,w,x)
218  endif
219  endif
220  case (3, 4)
221  if(interp_phy_first_row) then
222  if (ib == 3) then
223  ix = ixomax2+1
224  else
225  ix = ixomin2-1
226  end if
227  if(boundary_data_primitive) then
228  call phys_to_primitive(ixi^l,ixomin1,ix,ixomax1,ix,w,x)
229  endif
230  endif
231  if(bc_phy_first_row)then
232  if (ib == 3) then
233  ixoomax2=ixoomax2+1
234  else
235  ixoomin2=ixoomin2-1
236  endif
237  if(boundary_data_primitive) then
238  if (ib == 3) then
239  ixoomin2=ixoomax2
240  else
241  ixoomax2=ixoomin2
242  endif
243  call phys_to_primitive(ixi^l,ixoo^l,w,x)
244  if (ib == 3) then
245  ixoomin2=ixomin2 ! to be used later in to_conserved
246  else
247  ixoomax2=ixomax2
248  endif
249  endif
250  endif
251  do iw = 1, nwfluxbc
252  n_bc = bc_data_ix(iw, ib)
253  if (n_bc == -1) cycle
254 
255  tmp(ixomin1:ixomax1, ixomin2) = bc_data_get_2d(n_bc, &
256  x(ixomin1:ixomax1, ixomin2, 1), qt)
257 
258  if(interp_phy_first_row) then
259 
260  ! Approximate boundary value by linear interpolation to first ghost
261  ! cell, rest of ghost cells contains the same value
262  do i = 0, ixoomax2-ixoomin2
263  w(ixomin1:ixomax1, ixoomin2+i, iw) = &
264  2 * tmp(ixomin1:ixomax1, ixomin2) - &
265  w(ixomin1:ixomax1, ix, iw)
266  end do
267  else
268  do i = ixoomin2, ixoomax2
269  w(ixomin1:ixomax1, i, iw) = &
270  tmp(ixomin1:ixomax1, ixomin2)
271  end do
272 
273  endif
274  end do
275  if(interp_phy_first_row) then
276  if(boundary_data_primitive) then
277  call phys_to_conserved(ixi^l,ixomin1,ix,ixomax1,ix,w,x)
278  endif
279  endif
280  case default
281  call mpistop("bc_data_set: unknown iB")
282  end select
283  }
284 
285  {^ifthreed
286  select case (ib)
287  case (1, 2)
288  if(interp_phy_first_row) then
289  if (ib == 1) then
290  ix = ixomax1+1
291  else
292  ix = ixomin1-1
293  end if
294  if(boundary_data_primitive) then
295  call phys_to_primitive(ixi^l,ix,ixomin2,ixomin3,ix,ixomax2,ixomax3,w,x)
296  endif
297  endif
298  if(bc_phy_first_row)then
299  if (ib == 1) then
300  ixoomax1=ixoomax1+1
301  else
302  ixoomin1=ixoomin1-1
303  endif
304  if(boundary_data_primitive) then
305  if (ib == 1) then
306  ixoomin1=ixoomax1
307  else
308  ixoomax1=ixoomin1
309  endif
310  call phys_to_primitive(ixi^l,ixoo^l,w,x)
311  if (ib == 1) then
312  ixoomin1=ixomin1 ! to be used later in to_conserved
313  else
314  ixoomax1=ixomax1
315  endif
316  endif
317  endif
318  do iw = 1, nwfluxbc
319  n_bc = bc_data_ix(iw, ib)
320  if (n_bc == -1) cycle
321 
322  tmp(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3) = bc_data_get_3d(n_bc, &
323  x(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3, 2), &
324  x(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3, 3), qt)
325 
326  if(interp_phy_first_row) then
327 
328  ! Approximate boundary value by linear interpolation to first ghost
329  ! cell, rest of ghost cells contains the same value
330  do i = 0, ixoomax1-ixoomin1
331  w(ixoomin1+i, ixomin2:ixomax2, ixomin3:ixomax3, iw) = &
332  2 * tmp(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3) - &
333  w(ix, ixomin2:ixomax2, ixomin3:ixomax3, iw)
334  end do
335  else
336  do i = ixoomin1,ixoomax1
337  w(i, ixomin2:ixomax2, ixomin3:ixomax3, iw) = &
338  tmp(ixomin1, ixomin2:ixomax2, ixomin3:ixomax3)
339  end do
340  end if
341 
342  end do
343  if(interp_phy_first_row) then
344  if(boundary_data_primitive) then
345  call phys_to_conserved(ixi^l,ix,ixomin2,ixomin3,ix,ixomax2,ixomax3,w,x)
346  endif
347  endif
348  case (3, 4)
349  if(interp_phy_first_row) then
350  if (ib == 3) then
351  ix = ixomax2+1
352  else
353  ix = ixomin2-1
354  end if
355  if(boundary_data_primitive) then
356  call phys_to_primitive(ixi^l,ixomin1,ix,ixomin3,ixomax1,ix,ixomax3,w,x)
357  endif
358  endif
359  if(bc_phy_first_row)then
360  if (ib == 3) then
361  ixoomax2=ixoomax2+1
362  else
363  ixoomin2=ixoomin2-1
364  endif
365  if(boundary_data_primitive) then
366  if (ib == 3) then
367  ixoomin2=ixoomax2
368  else
369  ixoomax2=ixoomin2
370  endif
371  call phys_to_primitive(ixi^l,ixoo^l,w,x)
372  if (ib == 3) then
373  ixoomin2=ixomin2 ! to be used later in to_conserved
374  else
375  ixoomax2=ixomax2
376  endif
377  endif
378  endif
379  do iw = 1, nwfluxbc
380  n_bc = bc_data_ix(iw, ib)
381  if (n_bc == -1) cycle
382 
383  tmp(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3) = bc_data_get_3d(n_bc, &
384  x(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3, 1), &
385  x(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3, 3), qt)
386 
387  if(interp_phy_first_row) then
388 
389  ! Approximate boundary value by linear interpolation to first ghost
390  ! cell, rest of ghost cells contains the same value
391  do i = 0, ixoomax2-ixoomin2
392  w(ixomin1:ixomax1, ixoomin2+i, ixomin3:ixomax3, iw) = &
393  2 * tmp(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3) - &
394  w(ixomin1:ixomax1, ix, ixomin3:ixomax3, iw)
395  end do
396  else
397  do i = ixoomin2,ixoomax2
398  w(ixomin1:ixomax1, i, ixomin3:ixomax3, iw) = &
399  tmp(ixomin1:ixomax1, ixomin2, ixomin3:ixomax3)
400  end do
401  end if
402  end do
403  if(interp_phy_first_row) then
404  if(boundary_data_primitive) then
405  call phys_to_conserved(ixi^l,ixomin1,ix,ixomin3,ixomax1,ix,ixomax3,w,x)
406  endif
407  endif
408  case (5, 6)
409  if(interp_phy_first_row) then
410  if (ib == 5) then
411  ix = ixomax3+1
412  else
413  ix = ixomin3-1
414  end if
415  if(boundary_data_primitive) then
416  call phys_to_primitive(ixi^l,ixomin1,ixomax1,ixomin2,ixomax2,ix,ix,w,x)
417  endif
418  endif
419  if(bc_phy_first_row)then
420  if (ib == 5) then
421  ixoomax3=ixoomax3+1
422  else
423  ixoomin3=ixoomin3-1
424  endif
425  if(boundary_data_primitive) then
426  if (ib == 5) then
427  ixoomin3=ixoomax3
428  else
429  ixoomax3=ixoomin3
430  endif
431  call phys_to_primitive(ixi^l,ixoo^l,w,x)
432  if (ib == 5) then
433  ixoomin3=ixomin3 ! to be used later in to_conserved
434  else
435  ixoomax3=ixomax3
436  endif
437  endif
438  endif
439  do iw = 1, nwfluxbc
440  n_bc = bc_data_ix(iw, ib)
441  if (n_bc == -1) cycle
442 
443  tmp(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3) = bc_data_get_3d(n_bc, &
444  x(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3, 1), &
445  x(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3, 2), qt)
446 
447  if(interp_phy_first_row) then
448 
449  ! Approximate boundary value by linear interpolation to first ghost
450  ! cell, rest of ghost cells contains the same value
451  do i = 0, ixoomax3-ixoomin3
452  w(ixomin1:ixomax1, ixomin2:ixomax2, ixoomin3+i, iw) = &
453  2 * tmp(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3) - &
454  w(ixomin1:ixomax1, ixomin2:ixomax2, ix, iw)
455  end do
456  else
457  do i = ixoomin3,ixoomax3
458  w(ixomin1:ixomax1, ixomin2:ixomax2, i, iw) = &
459  tmp(ixomin1:ixomax1, ixomin2:ixomax2, ixomin3)
460  end do
461  end if
462  end do
463  if(interp_phy_first_row) then
464  if(boundary_data_primitive) then
465  call phys_to_conserved(ixi^l,ixomin1,ixomax1,ixomin2,ixomax2,ix,ix,w,x)
466  endif
467  endif
468  case default
469  call mpistop("bc_data_set: unknown iB")
470  end select
471  }
472 
473  if(boundary_data_primitive) then
474  call phys_to_conserved(ixi^l,ixoo^l,w,x)
475  endif
476 
477  end subroutine bc_data_set
478 
479  subroutine read_vtk_structured_points(fname, bc)
480  character(len=*), intent(in) :: fname
481  type(bc_data_t), intent(out) :: bc
482  double precision, allocatable :: tmp_data(:, :, :, :)
483  integer, parameter :: max_variables = 10
484  character(len=40) :: tmp_names(max_variables)
485  character(len=256) :: line
486  character(len=40) :: word, typename
487  integer, parameter :: my_unit = 123
488  integer :: n, n_points_total
489  integer :: n_components
490 
491  open(my_unit, file=trim(fname), status="old", action="read")
492 
493  ! Header, e.g. # vtk DataFile Version 2.0
494  read(my_unit, "(A)") line
495 
496  ! Dataset name
497  read(my_unit, "(A)") line
498 
499  ! ASCII / BINARY
500  read(my_unit, "(A)") line
501 
502  if (line /= "ASCII") then
503  print *, "line: ", trim(line)
504  error stop "read_vtk: not an ASCII file"
505  end if
506 
507  ! DATASET STRUCTURED_POINTS
508  read(my_unit, "(A)") line
509 
510  if (line /= "DATASET STRUCTURED_POINTS") then
511  print *, "line: ", trim(line)
512  error stop "read_vtk must have: DATASET STRUCTURED_POINTS"
513  end if
514 
515  ! DIMENSIONS NX NY NZ
516  read(my_unit, "(A)") line
517  read(line, *) word, bc%n_points
518 
519  if (word /= "DIMENSIONS") then
520  print *, "line: ", trim(line)
521  error stop "read_vtk expects: DIMENSIONS"
522  end if
523 
524  ! SPACING DX DY DZ
525  read(my_unit, *) word, bc%dx
526 
527  if (word /= "SPACING") then
528  print *, "line: ", trim(line)
529  error stop "read_vtk expects: SPACING"
530  end if
531 
532  ! ORIGIN 0 0 0
533  read(my_unit, *) word, bc%origin
534  if (word /= "ORIGIN") then
535  print *, "line: ", trim(line)
536  error stop "read_vtk expects: ORIGIN"
537  end if
538 
539  ! POINT_DATA NPOINTS
540  read(my_unit, *) word, n_points_total
541 
542  if (word /= "POINT_DATA") then
543  print *, "line: ", trim(line)
544  error stop "read_vtk expects: POINT_DATA n_points"
545  end if
546 
547  if (n_points_total /= product(bc%n_points)) &
548  error stop "read_vtk: n_points not equal to NX*NY*NZ"
549 
550  allocate(tmp_data(bc%n_points(1), bc%n_points(2), bc%n_points(3), &
551  max_variables))
552 
553  ! Read all scalar variables
554  do n = 1, max_variables
555 
556  ! SCALARS name type ncomponents
557  read(my_unit, *, end=900) word, tmp_names(n), typename, n_components
558 
559  if (word /= "SCALARS") then
560  print *, "line: ", trim(line)
561  error stop "read_vtk expects: SCALARS name type ncomponents"
562  end if
563 
564  if (n_components /= 1) error stop "read_vtk: ncomponents should be 1"
565 
566  ! LOOKUP_TABLE default
567  read(my_unit, *) word, typename
568 
569  if (word /= "LOOKUP_TABLE") then
570  print *, "line: ", trim(line)
571  error stop "read_vtk expects: LOOKUP_TABLE name"
572  end if
573 
574  ! Read list of data values
575  read(my_unit, *) tmp_data(:, :, :, n)
576  end do
577 
578  ! Done reading variables
579 900 continue
580 
581  close(my_unit)
582 
583  if (n == max_variables + 1) &
584  error stop "read_vtk: increase max_variables"
585 
586  ! Loop index is one higher than number of variables
587  bc%n_variables = n-1
588  bc%values = tmp_data(:, :, :, 1:n-1)
589  bc%names = tmp_names(1:n-1)
590  end subroutine read_vtk_structured_points
591 
592 end module mod_bc_data
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:480
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:131
logical, public, protected bc_phy_first_row
Definition: mod_bc_data.t:33
subroutine bc_read_params(files)
Read this module"s parameters from a file.
Definition: mod_bc_data.t:103
logical, public, protected boundary_data_primitive
Definition: mod_bc_data.t:34
elemental double precision function, public bc_data_get_3d(n_bc, x1, x2, qt)
Definition: mod_bc_data.t:119
subroutine, public bc_data_set(qt, ixIL, ixOL, iB, w, x)
Set boundary conditions according to user data.
Definition: mod_bc_data.t:144
subroutine, public bc_data_init()
Definition: mod_bc_data.t:44
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
logical, public, protected interp_phy_first_row
Definition: mod_bc_data.t:32
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, 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.
integer, parameter bc_icarus
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.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:56
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:55