MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_boundary_conditions.t
Go to the documentation of this file.
2 
3  implicit none
4  private
5 
6  public :: bc_phys
7  public :: getintbc
8 
9 contains
10 
11  !> fill ghost cells at a physical boundary
12  subroutine bc_phys(iside,idims,time,qdt,s,ixG^L,ixB^L)
14  use mod_bc_data, only: bc_data_set
16  use mod_physics
18 
19  integer, intent(in) :: iside, idims, ixg^l,ixb^l
20  double precision, intent(in) :: time,qdt
21  type(state), intent(inout) :: s
22  double precision :: wtmp(ixg^s,1:nwflux)
23 
24  double precision :: q(ixg^s),qp(ixg^s)
25  integer :: idir, is
26  integer :: ixos^l,hxo^l,jxo^l
27  integer :: iw, ib, ix^d, ixo^l, ixm^l, nghostcellsi,iib^d
28  logical :: isphysbound
29 
30  associate(x=>s%x,w=>s%w,ws=>s%ws)
31  select case (idims)
32  {case (^d)
33  if (iside==2) then
34  ! maximal boundary
35  ib=2*^d
36  ixomin^dd=ixbmax^d+1-nghostcells^d%ixOmin^dd=ixbmin^dd;
37  ixomax^dd=ixbmax^dd;
38  ! cont/symm/asymm types
39  do iw=1,nwflux+nwaux
40  select case (typeboundary(iw,ib))
41  case (bc_special)
42  ! skip it here, do AFTER all normal type boundaries are set
43  case (bc_cont)
44  do ix^d=ixomin^d,ixomax^d
45  w(ix^d^d%ixO^s,iw) = w(ixomin^d-1^d%ixO^s,iw)
46  end do
47  case (bc_symm)
48  w(ixo^s,iw) = w(ixomin^d-1:ixomin^d-nghostcells:-1^d%ixO^s,iw)
49  case (bc_asymm)
50  w(ixo^s,iw) =-w(ixomin^d-1:ixomin^d-nghostcells:-1^d%ixO^s,iw)
51  case (bc_periodic)
52  ! skip it here, periodic bc info should come from neighbors
53  case(bc_noinflow)
54  if (iw==1+^d)then
55  do ix^d=ixomin^d,ixomax^d
56  w(ix^d^d%ixO^s,iw) = max(w(ixomin^d-1^d%ixO^s,iw),zero)
57  end do
58  else
59  do ix^d=ixomin^d,ixomax^d
60  w(ix^d^d%ixO^s,iw) = w(ixomin^d-1^d%ixO^s,iw)
61  end do
62  end if
63  case (bc_aperiodic)
64  !this just multiplies the variables with (-), they have been set from neighbors just like periodic.
65  w(ixo^s,iw) = - w(ixo^s,iw)
66  case (bc_data)
67  ! skip it here, do AFTER all normal type boundaries are set
68  case (bc_icarus)
69  ! skip it here, do AFTER all normal type boundaries are set
70  case (bc_character)
71  ! skip it here, do AFTER all normal type boundaries are set
72  case default
73  write (unitterm,*) "Undefined boundarytype found in bc_phys", &
74  "for variable iw=",iw," and side iB=",ib
75  end select
76  end do
77  if(stagger_grid) then
78  do idir=1,nws
79  ! At this stage, extrapolation is applied only to the tangential components
80  if(idir==^d) cycle
81  ixosmax^dd=ixomax^dd;
82  ixosmin^dd=ixomin^dd-kr(^dd,idir);
83  select case(typeboundary(iw_mag(idir),ib))
84  case (bc_special)
85  ! skip it here, do AFTER all normal type boundaries are set
86  case (bc_cont)
87  do ix^d=ixosmin^d,ixosmax^d
88  ws(ix^d^d%ixOs^s,idir) = ws(ixosmin^d-1^d%ixOs^s,idir)
89  end do
90  case (bc_symm)
91  ws(ixos^s,idir) = ws(ixosmin^d-1:ixosmin^d-nghostcells:-1^d%ixOs^s,idir)
92  case (bc_asymm)
93  ws(ixos^s,idir) =-ws(ixosmin^d-1:ixosmin^d-nghostcells:-1^d%ixOs^s,idir)
94  case (bc_periodic)
95  case default
96  write (unitterm,*) "Undefined boundarytype found in bc_phys", &
97  "for variable iw=",iw," and side iB=",ib
98  end select
99  end do
100  ! Now that the tangential components are set,
101  ! fill the normal components using a prescription for the divergence.
102  ! This prescription is given by the typeB for the normal component.
103  do idir=1,nws
104  ! Consider only normal direction
105  if (idir/=^d) cycle
106  ixos^l=ixo^l;
107  select case(typeboundary(iw_mag(idir),ib))
108  case(bc_cont)
109  hxo^l=ixo^l-nghostcells*kr(^dd,^d);
110  ! Calculate divergence and partial divergence
111  call get_divb(s%w,ixg^l,hxo^l,q)
112  ws(ixos^s,idir)=zero
113  do ix^d=0,nghostcells-1
114  call get_divb(s%w,ixg^l,ixo^l,qp)
115  ws(ixosmin^d+ix^d^d%ixOs^s,idir)=&
116  (q(hxomax^d^d%hxO^s)*s%dvolume(hxomax^d^d%hxO^s)&
117  -qp(ixomin^d+ix^d^d%ixO^s)*s%dvolume(ixomin^d+ix^d^d%ixO^s))&
118  /s%surfaceC(ixosmin^d+ix^d^d%ixOs^s,^d)
119  end do
120  ! Fill cell averages
121  call phys_face_to_center(ixo^l,s)
122  case(bc_symm)
123  ! (a)symmetric normal B ensures symmetric divB
124  ws(ixos^s,idir)= ws(ixosmin^d-2:ixosmin^d-nghostcells-1:-1^d%ixOs^s,idir)
125  ! Fill cell averages
126  call phys_face_to_center(ixo^l,s)
127  case(bc_asymm)
128  ! (a)symmetric normal B ensures symmetric divB
129  ws(ixos^s,idir)=-ws(ixosmin^d-2:ixosmin^d-nghostcells-1:-1^d%ixOs^s,idir)
130  ! Fill cell averages
131  call phys_face_to_center(ixo^l,s)
132  case(bc_periodic)
133  ! Fill cell averages
134  call phys_face_to_center(ixo^l,s)
135  end select
136  end do
137  end if
138  else
139  ! minimal boundary
140  ib=2*^d-1
141  ixomin^dd=ixbmin^dd;
142  ixomax^dd=ixbmin^d-1+nghostcells^d%ixOmax^dd=ixbmax^dd;
143  ! cont/symm/asymm types
144  do iw=1,nwflux+nwaux
145  select case (typeboundary(iw,ib))
146  case (bc_special)
147  ! skip it here, do AFTER all normal type boundaries are set
148  case (bc_cont)
149  do ix^d=ixomin^d,ixomax^d
150  w(ix^d^d%ixO^s,iw) = w(ixomax^d+1^d%ixO^s,iw)
151  end do
152  case (bc_symm)
153  w(ixo^s,iw) = w(ixomax^d+nghostcells:ixomax^d+1:-1^d%ixO^s,iw)
154  case (bc_asymm)
155  w(ixo^s,iw) =-w(ixomax^d+nghostcells:ixomax^d+1:-1^d%ixO^s,iw)
156  case (bc_periodic)
157  ! skip it here, periodic bc info should come from neighbors
158  case(bc_noinflow)
159  if (iw==1+^d)then
160  do ix^d=ixomin^d,ixomax^d
161  w(ix^d^d%ixO^s,iw) = min(w(ixomax^d+1^d%ixO^s,iw),zero)
162  end do
163  else
164  do ix^d=ixomin^d,ixomax^d
165  w(ix^d^d%ixO^s,iw) = w(ixomax^d+1^d%ixO^s,iw)
166  end do
167  end if
168  case (bc_aperiodic)
169  !this just multiplies the variables with (-), they have been set from neighbors just like periodic.
170  w(ixo^s,iw) = - w(ixo^s,iw)
171  case (bc_data)
172  ! skip it here, do AFTER all normal type boundaries are set
173  case (bc_icarus)
174  ! skip it here, do AFTER all normal type boundaries are set
175  case (bc_character)
176  ! skip it here, do AFTER all normal type boundaries are set
177  case default
178  write (unitterm,*) "Undefined boundarytype found in bc_phys", &
179  "for variable iw=",iw," and side iB=",ib
180  end select
181  end do
182  if(stagger_grid) then
183  do idir=1,nws
184  ! At this stage, extrapolation is applied only to the tangential components
185  if(idir==^d) cycle
186  ixosmax^dd=ixomax^dd;
187  ixosmin^dd=ixomin^dd-kr(^dd,idir);
188  select case(typeboundary(iw_mag(idir),ib))
189  case (bc_special)
190  ! skip it here, periodic bc info should come from neighbors
191  case (bc_cont)
192  do ix^d=ixosmin^d,ixosmax^d
193  ws(ix^d^d%ixOs^s,idir) = ws(ixosmax^d+1^d%ixOs^s,idir)
194  end do
195  case (bc_symm)
196  ws(ixos^s,idir) = ws(ixosmax^d+nghostcells:ixosmax^d+1:-1^d%ixOs^s,idir)
197  case (bc_asymm)
198  ws(ixos^s,idir) =-ws(ixosmax^d+nghostcells:ixosmax^d+1:-1^d%ixOs^s,idir)
199  case (bc_periodic)
200  case default
201  write (unitterm,*) "Undefined boundarytype in bc_phys", &
202  "for variable iw=",iw," and side iB=",ib
203  end select
204  end do
205  ! Now that the tangential components are set,
206  ! fill the normal components using a prescription for the divergence.
207  ! This prescription is given by the typeB for the normal component.
208  do idir=1,nws
209  ! Consider only normal direction
210  if (idir/=^d) cycle
211  ixos^l=ixo^l-kr(^dd,^d);
212  select case(typeboundary(iw_mag(idir),ib))
213  case(bc_cont)
214  jxo^l=ixo^l+nghostcells*kr(^dd,^d);
215  ! Calculate divergence and partial divergence
216  call get_divb(s%w,ixg^l,jxo^l,q)
217  ws(ixos^s,idir)=zero
218  do ix^d=0,nghostcells-1
219  call get_divb(s%w,ixg^l,ixo^l,qp)
220  ws(ixosmax^d-ix^d^d%ixOs^s,idir)=&
221  -(q(jxomin^d^d%jxO^s)*s%dvolume(jxomin^d^d%jxO^s)&
222  -qp(ixomax^d-ix^d^d%ixO^s)*s%dvolume(ixomax^d-ix^d^d%ixO^s))&
223  /s%surfaceC(ixosmax^d-ix^d^d%ixOs^s,^d)
224  end do
225  ! Fill cell averages
226  call phys_face_to_center(ixo^l,s)
227  case(bc_symm)
228  ! (a)symmetric normal B ensures symmetric divB
229  ws(ixos^s,idir)= ws(ixosmax^d+nghostcells+1:ixosmax^d+2:-1^d%ixOs^s,idir)
230  ! Fill cell averages
231  call phys_face_to_center(ixo^l,s)
232  case(bc_asymm)
233  ! (a)symmetric normal B ensures symmetric divB
234  ws(ixos^s,idir)=-ws(ixosmax^d+nghostcells+1:ixosmax^d+2:-1^d%ixOs^s,idir)
235  ! Fill cell averages
236  call phys_face_to_center(ixo^l,s)
237  case(bc_periodic)
238  ! Fill cell averages
239  call phys_face_to_center(ixo^l,s)
240  end select
241  end do
242  end if
243  end if \}
244  end select
245 
246  ! do user defined special boundary conditions
247  if (any(typeboundary(1:nwflux+nwaux,ib)==bc_special)) then
248  if (.not. associated(usr_special_bc)) &
249  call mpistop("usr_special_bc not defined")
250  call usr_special_bc(time,ixg^l,ixo^l,ib,w,x)
251  end if
252 
253  ! fill boundary conditions from external data vtk files
254  if (any(typeboundary(1:nwflux+nwaux,ib)==bc_data)) then
255  call bc_data_set(time,ixg^l,ixo^l,ib,w,x)
256  end if
257 
258  ! fill boundary conditions from external data vtk files and do user defined special boundary conditions
259  if (any(typeboundary(1:nwflux+nwaux,ib)==bc_icarus)) then
260  call bc_data_set(time,ixg^l,ixo^l,ib,w,x)
261  if (.not. associated(usr_special_bc)) &
262  call mpistop("usr_special_bc not defined")
263  call usr_special_bc(time,ixg^l,ixo^l,ib,w,x)
264  end if
265 
266  end associate
267  end subroutine bc_phys
268 
269  !> fill inner boundary values
270  subroutine getintbc(time,ixG^L)
273 
274  double precision, intent(in) :: time
275  integer, intent(in) :: ixg^l
276 
277  integer :: iigrid, igrid, ixo^l
278 
279  ixo^l=ixg^l^lsubnghostcells;
280 
281  !$OMP PARALLEL DO SCHEDULE(dynamic) PRIVATE(igrid)
282  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
283  block=>ps(igrid)
284  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
285 
286  if (associated(usr_internal_bc)) then
287  call usr_internal_bc(node(plevel_,igrid),time,ixg^l,ixo^l,ps(igrid)%w,ps(igrid)%x)
288  end if
289  end do
290  !$OMP END PARALLEL DO
291 
292  end subroutine getintbc
293 
294 end module mod_boundary_conditions
Module to set boundary conditions from user data.
Definition: mod_bc_data.t:2
subroutine, public bc_data_set(qt, ixIL, ixOL, iB, w, x)
Set boundary conditions according to user data.
Definition: mod_bc_data.t:146
subroutine, public getintbc(time, ixGL)
fill inner boundary values
subroutine, public bc_phys(iside, idims, time, qdt, s, ixGL, ixBL)
fill ghost cells at a physical boundary
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer, parameter bc_noinflow
integer, parameter bc_asymm
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter bc_character
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
logical stagger_grid
True for using stagger grid.
integer, parameter bc_data
integer, parameter plevel_
double precision, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
integer, parameter unitterm
Unit for standard output.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer, parameter bc_cont
integer, parameter bc_icarus
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
integer, parameter bc_aperiodic
integer, dimension(:,:), allocatable node
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:85
Module with all the methods that users can customize in AMRVAC.
procedure(special_bc), pointer usr_special_bc
procedure(internal_bc), pointer usr_internal_bc