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 :: xC(ixi^s,1:ndim), A(ixi^s,1:3)
90  double precision :: circ(ixi^s,1:ndim), dxidir
91 
92  associate(ws=>s%ws,x=>s%x)
93 
94  a(:^d&,:)=zero
95  ws(:^d&,:)=zero
96 
97  ixis^l=s%ixGs^l;
98  ixo^l=ixi^l^lsubnghostcells;
99 
100  fc=0.d0
101  call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, a)
102 
103  ! This is important only in 3D
104  do idir=7-2*ndim,3
105  fe(ixi^s,idir) =-a(ixi^s,idir)
106  end do
107 
108  end associate
109  end subroutine fake_update
110 
111  !> calculate magnetic field from vector potential A at cell edges
112  subroutine b_from_vector_potentiala(ixIs^L, ixI^L, ixO^L, ws, x, A)
115 
116  integer, intent(in) :: ixIs^L, ixI^L, ixO^L
117  double precision, intent(inout) :: ws(ixis^s,1:nws),A(ixi^s,1:3)
118  double precision, intent(in) :: x(ixi^s,1:ndim)
119 
120  integer :: ixC^L, hxC^L, ixCp^L, ixCm^L, hxO^L, idim, idim1, idim2, idir
121  double precision :: xC(ixi^s,1:ndim)
122  double precision :: circ(ixi^s,1:ndim)
123 
124  a=zero
125  ws=zero
126 
127  {ixcmax^d=ixomax^d;}
128  {ixcmin^d=ixomin^d-1;} ! Extend range by one
129  do idir=7-2*ndim,3
130  do idim=1,ndim
131  ! Get edge coordinates
132  if (idim/=idir) then
133  xc(ixc^s,idim)=x(ixc^s,idim)+half*block%dx(ixc^s,idim)
134  else
135  xc(ixc^s,idim)=x(ixc^s,idim)
136  end if
137  end do
138  ! Initialise vector potential at the edge
139  call usr_init_vector_potential(ixi^l, ixc^l, xc, a(ixi^s,idir), idir)
140  end do
141 
142  ! Set NaN to zero (can happen e.g. on axis):
143  where(a(ixi^s,:)/=a(ixi^s,:))
144  a(ixi^s,:)=zero
145  end where
146 
147  ! sub integrals A ds
148  a(ixc^s,1:ndir)=a(ixc^s,1:ndir)*block%dsC(ixc^s,1:ndir)
149 
150  ! Take the curl of the vector potential
151  circ(:^d&,:) = zero
152 
153  ! Calculate circulation on each face
154  do idim1=1,ndim ! Coordinate perpendicular to face
155  ixcmax^d=ixomax^d;
156  ixcmin^d=ixomin^d-kr(idim1,^d);
157  do idim2=1,ndim
158  do idir=7-2*ndim,3 ! Direction of line integral
159  if(lvc(idim1,idim2,idir)==0) cycle
160  ! Assemble indices
161  hxc^l=ixc^l-kr(idim2,^d);
162  ! Add line integrals in direction idir
163  circ(ixc^s,idim1)=circ(ixc^s,idim1)&
164  +lvc(idim1,idim2,idir)* &
165  (a(ixc^s,idir)&
166  -a(hxc^s,idir))
167  end do
168  end do
169  ! Divide by the area of the face to get B
170  where(block%surfaceC(ixc^s,idim1)==0)
171  circ(ixc^s,idim1)=zero
172  elsewhere
173  circ(ixc^s,idim1)=circ(ixc^s,idim1)/block%surfaceC(ixc^s,idim1)
174  end where
175  ws(ixc^s,idim1) = circ(ixc^s,idim1)
176  end do
177 
178  end subroutine b_from_vector_potentiala
179 
180  !> Reconstruct scalar q within ixO^L to 1/2 dx in direction idir
181  !> Return both left and right reconstructed values
182  subroutine reconstruct(ixI^L,ixC^L,idir,q,qL,qR)
185 
186  integer, intent(in) :: ixI^L, ixC^L, idir
187  double precision, intent(in) :: q(ixi^s)
188  double precision, intent(out) :: qL(ixi^s), qR(ixi^s)
189 
190  double precision :: qC(ixi^s)
191  double precision,dimension(ixI^S) :: dqC,ldq,rdq
192  integer :: ixO^L,jxC^L,gxC^L,hxC^L
193 
194  jxc^l=ixc^l+kr(idir,^d);
195  gxcmin^d=ixcmin^d-kr(idir,^d);gxcmax^d=jxcmax^d;
196  hxc^l=gxc^l+kr(idir,^d);
197 
198  qr(gxc^s) = q(hxc^s)
199  ql(gxc^s) = q(gxc^s)
200 
201  select case (typelimiter)
202 
203  case (limiter_ppm)
204  ! the ordinary grid-index:
205  ixomin^d=ixcmin^d+kr(idir,^d);
206  ixomax^d=ixcmax^d;
207  call ppmlimitervar(ixi^l,ixo^l,idir,q,q,ql,qr)
208 
209  case (limiter_mp5)
210  call mp5limitervar(ixi^l,ixc^l,idir,q,ql,qr)
211 
212  case (limiter_weno5)
213  call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,1)
214 
215  case (limiter_wenozp5)
216  call weno5limiter(ixi^l,ixc^l,idir,dxlevel(idir),q,ql,qr,3)
217 
218  case default
219 
220  dqc(gxc^s)= qr(gxc^s)-ql(gxc^s)
221  call dwlimiter2(dqc,ixi^l,gxc^l,idir,typelimiter,ldq,rdq)
222  ql(ixc^s) = ql(ixc^s) + half*ldq(ixc^s)
223  qr(ixc^s) = qr(ixc^s) - half*rdq(jxc^s)
224  end select
225 
226  end subroutine reconstruct
227 
228 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:33
integer, parameter limiter_mp5
Definition: mod_limiter.t:26
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 ...
integer ndir
Number of spatial dimensions (components) for vector variables.
procedure(sub_face_to_center), pointer phys_face_to_center
Definition: mod_physics.t:65
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:25
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:46
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 fake_advance(igrid, idimLIM, s)
fake advance a step to calculate magnetic field
integer, parameter limiter_weno5
Definition: mod_limiter.t:29
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:47
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
Definition: mod_limiter.t:119
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
subroutine, public store_edge(igrid, ixIL, fE, idimLIM)