MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
boundary_conditions.t
Go to the documentation of this file.
1 !> fill ghost cells at a physical boundary
2 subroutine bc_phys(iside,idims,time,qdt,s,ixG^L,ixB^L)
4  use mod_bc_data, only: bc_data_set
6  use mod_physics
7  use mod_mhd_phys, only: get_divb
8 
9  integer, intent(in) :: iside, idims, ixG^L,ixB^L
10  double precision, intent(in) :: time,qdt
11  type(state), intent(inout) :: s
12  double precision :: wtmp(ixg^s,1:nwflux)
13 
14  integer :: idir, is
15  integer :: ixOs^L,hxO^L,jxO^L
16  double precision :: Q(ixg^s),Qp(ixg^s)
17  integer :: iw, iB, ix^D, ixO^L, ixM^L, nghostcellsi,iib^D
18  logical :: isphysbound
19 
20  associate(x=>s%x,w=>s%w,ws=>s%ws)
21  select case (idims)
22  {case (^d)
23  if (iside==2) then
24  ! maximal boundary
25  ib=ismax^d
26  ixomin^dd=ixbmax^d+1-nghostcells^d%ixOmin^dd=ixbmin^dd;
27  ixomax^dd=ixbmax^dd;
28  ! cont/symm/asymm types
29  do iw=1,nwflux+nwaux
30  select case (typeboundary(iw,ib))
31  case ("symm")
32  w(ixo^s,iw) = w(ixomin^d-1:ixomin^d-nghostcells:-1^d%ixO^s,iw)
33  case ("asymm")
34  w(ixo^s,iw) =-w(ixomin^d-1:ixomin^d-nghostcells:-1^d%ixO^s,iw)
35  case ("cont")
36  do ix^d=ixomin^d,ixomax^d
37  w(ix^d^d%ixO^s,iw) = w(ixomin^d-1^d%ixO^s,iw)
38  end do
39  case("noinflow")
40  if (iw==1+^d)then
41  do ix^d=ixomin^d,ixomax^d
42  w(ix^d^d%ixO^s,iw) = max(w(ixomin^d-1^d%ixO^s,iw),zero)
43  end do
44  else
45  do ix^d=ixomin^d,ixomax^d
46  w(ix^d^d%ixO^s,iw) = w(ixomin^d-1^d%ixO^s,iw)
47  end do
48  end if
49  case ("special", "bc_data")
50  ! skip it here, do AFTER all normal type boundaries are set
51  case ("character")
52  ! skip it here, do AFTER all normal type boundaries are set
53  case ("aperiodic")
54  !this just multiplies the variables with (-), they have been set from neighbors just like periodic.
55  w(ixo^s,iw) = - w(ixo^s,iw)
56  case ("periodic")
57  ! call mpistop("periodic bc info should come from neighbors")
58  case default
59  write (unitterm,*) "Undefined boundarytype ",typeboundary(iw,ib), &
60  "for variable iw=",iw," and side iB=",ib
61  end select
62  end do
63  if(stagger_grid) then
64  do idir=1,nws
65  ! At this stage, extrapolation is applied only to the tangential components
66  if(idir==^d) cycle
67  ixosmax^dd=ixomax^dd;
68  ixosmin^dd=ixomin^dd-kr(^dd,idir);
69  select case(typeboundary(iw_mag(idir),ib))
70  case ("symm")
71  ws(ixos^s,idir) = ws(ixosmin^d-1:ixosmin^d-nghostcells:-1^d%ixOs^s,idir)
72  case ("asymm")
73  ws(ixos^s,idir) =-ws(ixosmin^d-1:ixosmin^d-nghostcells:-1^d%ixOs^s,idir)
74  case ("cont")
75  do ix^d=ixosmin^d,ixosmax^d
76  ws(ix^d^d%ixOs^s,idir) = ws(ixosmin^d-1^d%ixOs^s,idir)
77  end do
78  case ("periodic")
79  case ("special")
80  ! skip it here, do AFTER all normal type boundaries are set
81  case default
82  write (unitterm,*) "Undefined boundarytype ",typeboundary(iw,ib), &
83  "for variable iw=",iw," and side iB=",ib
84  end select
85  end do
86  ! Now that the tangential components are set,
87  ! fill the normal components using a prescription for the divergence.
88  ! This prescription is given by the typeB for the normal component.
89  do idir=1,nws
90  ! Consider only normal direction
91  if (idir/=^d) cycle
92  ixos^l=ixo^l;
93  hxo^l=ixo^l-nghostcells*kr(^dd,^d);
94  ! Calculate divergence and partial divergence
95  call get_divb(s%w,ixg^l,hxo^l,q)
96  select case(typeboundary(iw_mag(idir),ib))
97  case("symm")
98  ws(ixos^s,idir)=zero
99  do ix^d=0,nghostcells-1
100  call get_divb(s%w,ixg^l,ixo^l,qp)
101  ws(ixosmin^d+ix^d^d%ixOs^s,idir)=&
102  (q(hxomax^d-ix^d^d%hxO^s)*s%dvolume(hxomax^d-ix^d^d%hxO^s)&
103  -qp(ixomin^d+ix^d^d%ixO^s)*s%dvolume(ixomin^d+ix^d^d%ixO^s))&
104  /s%surfaceC(ixosmin^d+ix^d^d%ixOs^s,^d)
105  end do
106  case("asymm")
107  ws(ixos^s,idir)=zero
108  do ix^d=0,nghostcells-1
109  call get_divb(s%w,ixg^l,ixo^l,qp)
110  ws(ixosmin^d+ix^d^d%ixOs^s,idir)=&
111  (-q(hxomax^d-ix^d^d%hxO^s)*s%dvolume(hxomax^d-ix^d^d%hxO^s)&
112  -qp(ixomin^d+ix^d^d%ixO^s)*s%dvolume(ixomin^d+ix^d^d%ixO^s))&
113  /s%surfaceC(ixosmin^d+ix^d^d%ixOs^s,^d)
114  end do
115  case("cont")
116  ws(ixos^s,idir)=zero
117  do ix^d=0,nghostcells-1
118  call get_divb(s%w,ixg^l,ixo^l,qp)
119  ws(ixosmin^d+ix^d^d%ixOs^s,idir)=&
120  (q(hxomax^d^d%hxO^s)*s%dvolume(hxomax^d^d%hxO^s)&
121  -qp(ixomin^d+ix^d^d%ixO^s)*s%dvolume(ixomin^d+ix^d^d%ixO^s))&
122  /s%surfaceC(ixosmin^d+ix^d^d%ixOs^s,^d)
123  end do
124  case("periodic")
125  end select
126  end do
127  ! Fill cell averages
128  call phys_face_to_center(ixo^l,s)
129  end if
130  else
131  ! minimal boundary
132  ib=ismin^d
133  ixomin^dd=ixbmin^dd;
134  ixomax^dd=ixbmin^d-1+nghostcells^d%ixOmax^dd=ixbmax^dd;
135  ! cont/symm/asymm types
136  do iw=1,nwflux+nwaux
137  select case (typeboundary(iw,ib))
138  case ("symm")
139  w(ixo^s,iw) = w(ixomax^d+nghostcells:ixomax^d+1:-1^d%ixO^s,iw)
140  case ("asymm")
141  w(ixo^s,iw) =-w(ixomax^d+nghostcells:ixomax^d+1:-1^d%ixO^s,iw)
142  case ("cont")
143  do ix^d=ixomin^d,ixomax^d
144  w(ix^d^d%ixO^s,iw) = w(ixomax^d+1^d%ixO^s,iw)
145  end do
146  case("noinflow")
147  if (iw==1+^d)then
148  do ix^d=ixomin^d,ixomax^d
149  w(ix^d^d%ixO^s,iw) = min(w(ixomax^d+1^d%ixO^s,iw),zero)
150  end do
151  else
152  do ix^d=ixomin^d,ixomax^d
153  w(ix^d^d%ixO^s,iw) = w(ixomax^d+1^d%ixO^s,iw)
154  end do
155  end if
156  case ("special", "bc_data")
157  ! skip it here, do AFTER all normal type boundaries are set
158  case ("character")
159  ! skip it here, do AFTER all normal type boundaries are set
160  case ("aperiodic")
161  !this just multiplies the variables with (-), they have been set from neighbors just like periodic.
162  w(ixo^s,iw) = - w(ixo^s,iw)
163  case ("periodic")
164  ! call mpistop("periodic bc info should come from neighbors")
165  case default
166  write (unitterm,*) "Undefined boundarytype ",typeboundary(iw,ib), &
167  "for variable iw=",iw," and side iB=",ib
168  end select
169  end do
170  if(stagger_grid) then
171  do idir=1,nws
172  ! At this stage, extrapolation is applied only to the tangential components
173  if(idir==^d) cycle
174  ixosmax^dd=ixomax^dd;
175  ixosmin^dd=ixomin^dd-kr(^dd,idir);
176  select case(typeboundary(iw_mag(idir),ib))
177  case ("symm")
178  ws(ixos^s,idir) = ws(ixosmax^d+nghostcells:ixosmax^d+1:-1^d%ixOs^s,idir)
179  case ("asymm")
180  ws(ixos^s,idir) =-ws(ixosmax^d+nghostcells:ixosmax^d+1:-1^d%ixOs^s,idir)
181  case ("cont")
182  do ix^d=ixosmin^d,ixosmax^d
183  ws(ix^d^d%ixOs^s,idir) = ws(ixosmax^d+1^d%ixOs^s,idir)
184  end do
185  case ("periodic")
186  case ("special")
187  ! skip it here, do AFTER all normal type boundaries are set
188  case default
189  write (unitterm,*) "Undefined boundarytype ",typeboundary(iw,ib), &
190  "for variable iw=",iw," and side iB=",ib
191  end select
192  end do
193  ! Now that the tangential components are set,
194  ! fill the normal components using a prescription for the divergence.
195  ! This prescription is given by the typeB for the normal component.
196  do idir=1,nws
197  ! Consider only normal direction
198  if (idir/=^d) cycle
199  ixos^l=ixo^l-kr(^dd,^d);
200  jxo^l=ixo^l+nghostcells*kr(^dd,^d);
201  ! Calculate divergence and partial divergence
202  call get_divb(s%w,ixg^l,jxo^l,q)
203  select case(typeboundary(iw_mag(idir),ib))
204  case("symm")
205  ws(ixos^s,idir)=zero
206  do ix^d=0,nghostcells-1
207  call get_divb(s%w,ixg^l,ixo^l,qp)
208  ws(ixosmax^d-ix^d^d%ixOs^s,idir)=&
209  -(q(jxomin^d+ix^d^d%jxO^s)*s%dvolume(jxomin^d+ix^d^d%jxO^s)&
210  -qp(ixomax^d-ix^d^d%ixO^s)*s%dvolume(ixomax^d-ix^d^d%ixO^s))&
211  /s%surfaceC(ixosmax^d-ix^d^d%ixOs^s,^d)
212  end do
213  case("asymm")
214  ws(ixos^s,idir)=zero
215  do ix^d=0,nghostcells-1
216  call get_divb(s%w,ixg^l,ixo^l,qp)
217  ws(ixosmax^d-ix^d^d%ixOs^s,idir)=&
218  -(-q(jxomin^d+ix^d^d%jxO^s)*s%dvolume(jxomin^d+ix^d^d%jxO^s)&
219  -qp(ixomax^d-ix^d^d%ixO^s)*s%dvolume(ixomax^d-ix^d^d%ixO^s))&
220  /s%surfaceC(ixosmax^d-ix^d^d%ixOs^s,^d)
221  end do
222  case("cont")
223  ws(ixos^s,idir)=zero
224  do ix^d=0,nghostcells-1
225  call get_divb(s%w,ixg^l,ixo^l,qp)
226  ws(ixosmax^d-ix^d^d%ixOs^s,idir)=&
227  -(q(jxomin^d^d%jxO^s)*s%dvolume(jxomin^d^d%jxO^s)&
228  -qp(ixomax^d-ix^d^d%ixO^s)*s%dvolume(ixomax^d-ix^d^d%ixO^s))&
229  /s%surfaceC(ixosmax^d-ix^d^d%ixOs^s,^d)
230  end do
231  case("periodic")
232  end select
233  end do
234  ! Fill cell averages
235  call phys_face_to_center(ixo^l,s)
236  end if
237  end if \}
238  end select
239 
240  ! do user defined special boundary conditions
241  if (any(typeboundary(1:nwflux+nwaux,ib)=="special")) then
242  if (.not. associated(usr_special_bc)) &
243  call mpistop("usr_special_bc not defined")
244  call usr_special_bc(time,ixg^l,ixo^l,ib,w,x)
245  end if
246 
247  ! fill boundary conditions from external data vtk files
248  if (any(typeboundary(1:nwflux+nwaux,ib)=="bc_data")) then
249  call bc_data_set(time,ixg^l,ixo^l,ib,w,x)
250  end if
251 
252  {#IFDEF EVOLVINGBOUNDARY
253  if (any(typeboundary(1:nwflux,ib)=="character")) then
254  ixm^l=ixm^ll;
255  if(ixgmax1==ixghi1) then
256  nghostcellsi=nghostcells
257  else
258  nghostcellsi=ceiling(nghostcells*0.5d0)
259  end if
260  select case (idims)
261  {case (^d)
262  if (iside==2) then
263  ! maximal boundary
264  ixomin^dd=ixgmax^d+1-nghostcellsi^d%ixOmin^dd=ixbmin^dd;
265  ixomax^dd=ixbmax^dd;
266  if(all(w(ixo^s,1:nwflux)==0.d0)) then
267  do ix^d=ixomin^d,ixomax^d
268  w(ix^d^d%ixO^s,1:nwflux) = w(ixomin^d-1^d%ixO^s,1:nwflux)
269  end do
270  end if
271  if(qdt>0.d0.and.ixgmax^d==ixghi^d) then
272  ixomin^dd=ixomin^d^d%ixOmin^dd=ixmmin^dd;
273  ixomax^dd=ixomax^d^d%ixOmax^dd=ixmmax^dd;
274  wtmp(ixg^s,1:nw)=pso(saveigrid)%w(ixg^s,1:nw)
275  call characteristic_project(idims,iside,ixg^l,ixo^l,wtmp,x,dxlevel,qdt)
276  w(ixo^s,1:nwflux)=wtmp(ixo^s,1:nwflux)
277  end if
278  else
279  ! minimal boundary
280  ixomin^dd=ixbmin^dd;
281  ixomax^dd=ixgmin^d-1+nghostcellsi^d%ixOmax^dd=ixbmax^dd;
282  if(all(w(ixo^s,1:nwflux)==0.d0)) then
283  do ix^d=ixomin^d,ixomax^d
284  w(ix^d^d%ixO^s,1:nwflux) = w(ixomax^d+1^d%ixO^s,1:nwflux)
285  end do
286  end if
287  if(qdt>0.d0.and.ixgmax^d==ixghi^d) then
288  ixomin^dd=ixomin^d^d%ixOmin^dd=ixmmin^dd;
289  ixomax^dd=ixomax^d^d%ixOmax^dd=ixmmax^dd;
290  wtmp(ixg^s,1:nw)=pso(saveigrid)%w(ixg^s,1:nw)
291  call characteristic_project(idims,iside,ixg^l,ixo^l,wtmp,x,dxlevel,qdt)
292  w(ixo^s,1:nwflux)=wtmp(ixo^s,1:nwflux)
293  end if
294  end if \}
295  end select
296  if(ixgmax1==ixghi1) then
297  call identifyphysbound(saveigrid,isphysbound,iib^d)
298  if(iib1==-1.and.iib2==-1) then
299  do ix2=nghostcells,1,-1
300  do ix1=nghostcells,1,-1
301  w(ix^d,1:nwflux)=(w(ix1+1,ix2+1,1:nwflux)+w(ix1+1,ix2,1:nwflux)+w(ix1,ix2+1,1:nwflux))/3.d0
302  end do
303  end do
304  end if
305  if(iib1== 1.and.iib2==-1) then
306  do ix2=nghostcells,1,-1
307  do ix1=ixmmax1+1,ixgmax1
308  w(ix^d,1:nwflux)=(w(ix1-1,ix2+1,1:nwflux)+w(ix1-1,ix2,1:nwflux)+w(ix1,ix2+1,1:nwflux))/3.d0
309  end do
310  end do
311  end if
312  if(iib1==-1.and.iib2== 1) then
313  do ix2=ixmmax2+1,ixgmax2
314  do ix1=nghostcells,1,-1
315  w(ix^d,1:nwflux)=(w(ix1+1,ix2-1,1:nwflux)+w(ix1+1,ix2,1:nwflux)+w(ix1,ix2-1,1:nwflux))/3.d0
316  end do
317  end do
318  end if
319  if(iib1== 1.and.iib2== 1) then
320  do ix2=ixmmax2+1,ixgmax2
321  do ix1=ixmmax1+1,ixgmax1
322  w(ix^d,1:nwflux)=(w(ix1-1,ix2-1,1:nwflux)+w(ix1-1,ix2,1:nwflux)+w(ix1,ix2-1,1:nwflux))/3.d0
323  end do
324  end do
325  end if
326  end if
327  end if
328  }
329  !end do
330 
331  end associate
332 end subroutine bc_phys
333 
334 !> fill inner boundary values
335 subroutine getintbc(time,ixG^L)
338 
339  double precision, intent(in) :: time
340  integer, intent(in) :: ixG^L
341 
342  integer :: iigrid, igrid, ixO^L,level
343 
344  ixo^l=ixg^l^lsubnghostcells;
345 
346  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
347  !do iigrid=1,igridstail; igrid=igrids(iigrid);
348  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
349  block=>ps(igrid)
352  level=node(plevel_,igrid)
353  saveigrid=igrid
354 
355  if (associated(usr_internal_bc)) then
356  call usr_internal_bc(level,time,ixg^l,ixo^l,ps(igrid)%w,ps(igrid)%x)
357  end if
358  end do
359 
360 end subroutine getintbc
This module contains definitions of global parameters and variables and some generic functions/subrou...
procedure(internal_bc), pointer usr_internal_bc
integer, dimension(:), allocatable type_gradient_limiter
Type of slope limiter used for computing gradients or divergences, when typegrad or typediv are set t...
integer, parameter plevel_
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:65
subroutine, public bc_data_set(qt, ixIL, ixOL, iB, w, x)
Set boundary conditions according to user data.
Definition: mod_bc_data.t:125
Module to set boundary conditions from user data.
Definition: mod_bc_data.t:2
integer nghostcells
Number of ghost cells surrounding a grid.
Module with all the methods that users can customize in AMRVAC.
integer ixghi
Upper index of grid block arrays.
logical stagger_grid
True for using stagger grid.
character(len=std_len), dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
integer, parameter unitterm
Unit for standard output.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
subroutine bc_phys(iside, idims, time, qdt, s, ixGL, ixBL)
fill ghost cells at a physical boundary
subroutine characteristic_project(idims, iside, ixIL, ixOL, w, x, dxndim, qdt)
Definition: eigensystem.t:182
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
double precision, dimension(ndim) dxlevel
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
integer, dimension(:,:), allocatable node
Magneto-hydrodynamics module.
Definition: mod_mhd_phys.t:2
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
subroutine getintbc(time, ixGL)
fill inner boundary values
procedure(special_bc), pointer usr_special_bc