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