MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_constrained_transport.t
Go to the documentation of this file.
2  implicit none
3  public
4 
5  !> velocities store for constrained transport
7  double precision, dimension(:^D&,:), allocatable :: vnorm,cbarmin,cbarmax
8  double precision, dimension(:^D&,:,:), allocatable :: vbarc,vbarlc,vbarrc
9  end type ct_velocity
10 
11  type(ct_velocity), save :: vcts
12 
13 contains
14  !> re-calculate the magnetic field from the vector potential in a completely
15  !> divergency free way
16  subroutine recalculateb
19  use mod_physics
21 
22  integer :: igrid,iigrid
23 
24  call init_comm_fix_conserve(1,ndim,^nd)
25 
26  !$OMP PARALLEL DO PRIVATE(igrid)
27  do iigrid=1,igridstail; igrid=igrids(iigrid);
28  ! Make zero the magnetic fluxes
29  ! Fake advance, storing electric fields at edges
30  block=>ps(igrid)
31  call fake_advance(igrid,1,^nd,ps(igrid))
32 
33  end do
34  !$OMP END PARALLEL DO
35 
36  ! Do correction
37  call recvflux(1,ndim)
38  call sendflux(1,ndim)
39  call fix_conserve(ps,1,ndim,1,^nd)
40 
41  call fix_edges(ps,1,^nd)
42 
43  ! Now we fill the centers for the staggered variables
44  !$OMP PARALLEL DO PRIVATE(igrid)
45  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
46  call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
47  ! update cell center magnetic field
48  call phys_face_to_center(ixm^ll,ps(igrid))
49  call phys_to_conserved(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
50  end do
51  !$OMP END PARALLEL DO
52 
53  call getbc(global_time,0.d0,ps,1,nwflux+nwaux)
54 
55  end subroutine recalculateb
56 
57  !> fake advance a step to calculate magnetic field
58  subroutine fake_advance(igrid,idim^LIM,s)
61 
62  integer :: igrid,idim^LIM
63  type(state) :: s
64 
65  double precision :: dx^D
66  double precision :: fC(ixg^t,1:nwflux,1:ndim)
67  double precision :: fE(ixg^t,7-2*ndim:3)
68 
69  dx^d=rnode(rpdx^d_,igrid);
70 
71  call fake_update(ixg^ll,s,fc,fe,dx^d)
72 
73  call store_flux(igrid,fc,idim^lim,^nd)
74  call store_edge(igrid,ixg^ll,fe,idim^lim)
75 
76  end subroutine fake_advance
77 
78  !>fake update magnetic field from vector potential
79  subroutine fake_update(ixI^L,s,fC,fE,dx^D)
81 
82  integer :: ixI^L
83  type(state) :: s
84  double precision :: fC(ixi^s,1:nwflux,1:ndim)
85  double precision :: fE(ixi^s,7-2*ndim:3)
86  double precision :: dx^D
87 
88  integer :: ixIs^L,ixO^L,idir
89  double precision :: A(s%ixgs^s,1:3)
90 
91  associate(ws=>s%ws,x=>s%x)
92 
93  a=zero
94  ws=zero
95 
96  ixis^l=s%ixGs^l;
97  ixo^l=ixi^l^lsubnghostcells;
98 
99  fc=0.d0
100  call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, a)
101 
102  ! This is important only in 3D
103  do idir=7-2*ndim,3
104  fe(ixi^s,idir) =-a(ixi^s,idir)
105  end do
106 
107  end associate
108  end subroutine fake_update
109 
110  !> calculate magnetic field from vector potential A at cell edges
111  subroutine b_from_vector_potentiala(ixIs^L, ixI^L, ixO^L, ws, x, A)
114 
115  integer, intent(in) :: ixIs^L, ixI^L, ixO^L
116  double precision, intent(inout) :: ws(ixis^s,1:nws),A(ixis^s,1:3)
117  double precision, intent(in) :: x(ixi^s,1:ndim)
118 
119  integer :: ixC^L, hxC^L, idim1, idim2, idir
120  double precision :: xC(ixis^s,1:ndim),xCC(ixis^s,1:ndim)
121  double precision :: circ(ixis^s,1:ndim)
122 
123  a=zero
124  ! extend one layer of cell center locations in xCC
125  xcc=0.d0
126  xcc(ixi^s,1:ndim)=x(ixi^s,1:ndim)
127  {
128  xcc(ixismin^d^d%ixI^s,1:ndim)=x(iximin^d^d%ixI^s,1:ndim)
129  xcc(ixismin^d^d%ixIs^s,^d)=x({iximin^dd,},^d)-block%dx({iximin^dd,},^d)
130  \}
131  {^ifthreed
132  xcc(iximin1:iximax1,ixismin2,ixismin3,1)=x(iximin1:iximax1,iximin2,iximin3,1)
133  xcc(ixismin1,iximin2:iximax2,ixismin3,2)=x(iximin1,iximin2:iximax2,iximin3,2)
134  xcc(ixismin1,ixismin2,iximin3:iximax3,3)=x(iximin1,iximin2,iximin3:iximax3,3)
135  }
136 
137  do idir=7-2*ndim,3
138  ixcmax^d=ixomax^d;
139  ixcmin^d=ixomin^d-1+kr(idir,^d);
140  do idim1=1,ndim
141  ! Get edge coordinates
142  if (idim1/=idir) then
143  xc(ixc^s,idim1)=xcc(ixc^s,idim1)+half*block%dx(ixc^s,idim1)
144  else
145  xc(ixc^s,idim1)=xcc(ixc^s,idim1)
146  end if
147  end do
148  ! Initialise vector potential at the edge
149  call usr_init_vector_potential(ixis^l, ixc^l, xc, a(ixis^s,idir), idir)
150  a(ixc^s,idir)=a(ixc^s,idir)*block%dsC(ixc^s,idir)
151  end do
152 
153  ! Take the curl of the vector potential
154  circ=zero
155  ! Calculate circulation on each face
156  do idim1=1,ndim ! Coordinate perpendicular to face
157  ixcmax^d=ixomax^d;
158  ixcmin^d=ixomin^d-kr(idim1,^d);
159  do idim2=1,ndim
160  do idir=7-2*ndim,3 ! Direction of line integral
161  if(lvc(idim1,idim2,idir)==0) cycle
162  ! Assemble indices
163  hxc^l=ixc^l-kr(idim2,^d);
164  ! Add line integrals in direction idir
165  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
166  +lvc(idim1,idim2,idir)* &
167  (a(ixc^s,idir)&
168  -a(hxc^s,idir))
169  end do
170  end do
171  ! Divide by the area of the face to get B
172  where(block%surfaceC(ixc^s,idim1)==0)
173  circ(ixc^s,idim1)=zero
174  elsewhere
175  circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
176  end where
177  ws(ixc^s,idim1) = circ(ixc^s,idim1)
178  end do
179 
180  end subroutine b_from_vector_potentiala
181 
182  !> Reconstruct scalar q within ixO^L to 1/2 dx in direction idir
183  !> Return both left and right reconstructed values
184  subroutine reconstruct(ixI^L,ixC^L,idir,q,qL,qR)
187 
188  integer, intent(in) :: ixI^L, ixC^L, idir
189  double precision, intent(in) :: q(ixi^s)
190  double precision, intent(out) :: qL(ixi^s), qR(ixi^s)
191 
192  double precision :: qC(ixi^s)
193  double precision,dimension(ixI^S) :: dqC,ldq,rdq
194  integer :: ixO^L,jxC^L,gxC^L,hxC^L
195 
196  jxc^l=ixc^l+kr(idir,^d);
197  gxcmin^d=ixcmin^d-kr(idir,^d);gxcmax^d=jxcmax^d;
198  hxc^l=gxc^l+kr(idir,^d);
199 
200  qr(gxc^s) = q(hxc^s)
201  ql(gxc^s) = q(gxc^s)
202 
203  select case (typelimiter)
204 
205  case (limiter_ppm)
206  ! the ordinary grid-index:
207  ixomin^d=ixcmin^d+kr(idir,^d);
208  ixomax^d=ixcmax^d;
209  call ppmlimitervar(ixi^l,ixo^l,idir,q,q,ql,qr)
210 
211  case (limiter_mp5)
212  call mp5limitervar(ixi^l,ixc^l,idir,q,ql,qr)
213 
214  case (limiter_weno5)
215  call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,1)
216 
217  case (limiter_wenozp5)
218  call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,3)
219 
220  case default
221 
222  dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
223  call dwlimiter2(dqc,ixi^l,gxc^l,idir,typelimiter,ldq,rdq)
224  ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
225  qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
226  end select
227 
228  end subroutine reconstruct
229 
230 end module mod_constrained_transport
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
update ghost cells of all blocks including physical boundaries
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
subroutine fake_update(ixIL, s, fC, fE, dxD)
fake update magnetic field from vector potential
integer, parameter limiter_wenozp5
Definition: mod_limiter.t:35
integer, parameter limiter_mp5
Definition: mod_limiter.t:28
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
subroutine recalculateb
re-calculate the magnetic field from the vector potential in a completely divergency free way ...
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:85
Module for flux conservation near refinement boundaries.
velocities store for constrained transport
subroutine, public fix_edges(psuse, idimLIM)
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
integer, parameter limiter_ppm
Definition: mod_limiter.t:27
subroutine, public sendflux(idimLIM)
Module with slope/flux limiters.
Definition: mod_limiter.t:2
Module with all the methods that users can customize in AMRVAC.
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:60
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
double precision global_time
The global simulation time.
subroutine reconstruct(ixIL, ixCL, idir, q, qL, qR)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
procedure(init_vector_potential), pointer usr_init_vector_potential
subroutine, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw, a2max)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
Definition: mod_limiter.t:128
subroutine fake_advance(igrid, idimLIM, s)
fake advance a step to calculate magnetic field
integer, parameter limiter_weno5
Definition: mod_limiter.t:31
double precision, dimension(ndim) dxlevel
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine, public recvflux(idimLIM)
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:61
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
subroutine, public store_edge(igrid, ixIL, fE, idimLIM)