MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_finite_volume.t
Go to the documentation of this file.
1 !> Module with finite volume methods for fluxes
3  implicit none
4  private
5 
6  public :: finite_volume
7  public :: hancock
8  public :: reconstruct_lr
9 
10 contains
11 
12  !> The non-conservative Hancock predictor for TVDLF
13  !>
14  !> on entry:
15  !> input available on ixI^L=ixG^L asks for output on ixO^L=ixG^L^LSUBnghostcells
16  !> one entry: (predictor): wCT -- w_n wnew -- w_n qdt=dt/2
17  !> on exit : (predictor): wCT -- w_n wnew -- w_n+1/2
18  subroutine hancock(qdt,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,snew,dx^D,x)
21  use mod_source, only: addsource2
22 
23  integer, intent(in) :: ixI^L, ixO^L, idims^LIM
24  double precision, intent(in) :: qdt, qtC, qt, dx^D, x(ixi^s,1:ndim)
25  type(state) :: sCT, snew
26 
27  double precision, dimension(ixI^S,1:nw) :: wprim, wLC, wRC
28  ! left and right constructed status in primitive form, needed for better performance
29  double precision, dimension(ixI^S,1:nw) :: wLp, wRp
30  double precision, dimension(ixO^S) :: inv_volume
31  double precision :: fLC(ixi^s, nwflux), fRC(ixi^s, nwflux)
32  double precision :: dxinv(1:ndim),dxdim(1:ndim)
33  integer :: idims, iw, ix^L, hxO^L
34 
35  associate(wct=>sct%w,wnew=>snew%w)
36  ! Expand limits in each idims direction in which fluxes are added
37  ix^l=ixo^l;
38  do idims= idims^lim
39  ix^l=ix^l^laddkr(idims,^d);
40  end do
41  if (ixi^l^ltix^l|.or.|.or.) &
42  call mpistop("Error in Hancock: Nonconforming input limits")
43 
44  wprim=wct
45  call phys_to_primitive(ixi^l,ixi^l,wprim,x)
46 
47  ^d&dxinv(^d)=-qdt/dx^d;
48  ^d&dxdim(^d)=dx^d;
49  do idims= idims^lim
50  b0i=idims
51  ! Calculate w_j+g_j/2 and w_j-g_j/2
52  ! First copy all variables, then upwind wLC and wRC.
53  ! wLC is to the left of ixO, wRC is to the right of wCT.
54  hxo^l=ixo^l-kr(idims,^d);
55 
56  wrp(hxo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
57  wlp(ixo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
58 
59  ! apply limited reconstruction for left and right status at cell interfaces
60  call reconstruct_lr(ixi^l,ixo^l,hxo^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxdim(idims))
61 
62  ! Calculate the fLC and fRC fluxes
63  call phys_get_flux(wrc,wrp,x,ixi^l,hxo^l,idims,frc)
64  call phys_get_flux(wlc,wlp,x,ixi^l,ixo^l,idims,flc)
65 
66  ! Advect w(iw)
67  do iw=1,nwflux
68  if (slab_uniform) then
69  if (associated(phys_iw_methods(iw)%inv_capacity)) then
70  call phys_iw_methods(iw)%inv_capacity(ixi^l, ixo^l, wnew, inv_volume)
71  wnew(ixo^s,iw)=wnew(ixo^s,iw)+dxinv(idims) * inv_volume * &
72  (flc(ixo^s, iw)-frc(hxo^s, iw))
73  else
74  wnew(ixo^s,iw)=wnew(ixo^s,iw)+dxinv(idims)* &
75  (flc(ixo^s, iw)-frc(hxo^s, iw))
76  end if
77  else
78  if (associated(phys_iw_methods(iw)%inv_capacity)) then
79  call phys_iw_methods(iw)%inv_capacity(ixi^l, ixo^l, wnew, inv_volume)
80  else
81  inv_volume = 1.0d0
82  end if
83  inv_volume = inv_volume/block%dvolume(ixo^s)
84 
85  wnew(ixo^s,iw)=wnew(ixo^s,iw) - qdt * inv_volume &
86  *(block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
87  -block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
88  end if
89  end do
90  end do ! next idims
91  b0i=0
92 
93  do iw = 1, nwflux
94  if (associated(phys_iw_methods(iw)%inv_capacity)) then
95  ! Copy state before adding source terms
96  wprim(ixo^s, iw) = wnew(ixo^s, iw)
97  end if
98  end do
99 
100  if (.not.slab.and.idimsmin==1) call phys_add_source_geom(qdt,ixi^l,ixo^l,wct,wnew,x)
101 
102  call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim), &
103  ixi^l,ixo^l,1,nw,qtc,wct,qt,wnew,x,.false.)
104 
105  ! If there are capacity functions, now correct the added source terms
106  do iw = 1, nwflux
107  if (associated(phys_iw_methods(iw)%inv_capacity)) then
108  call phys_iw_methods(iw)%inv_capacity(ixi^l, ixo^l, wnew, inv_volume)
109  wnew(ixo^s, iw) = wprim(ixo^s, iw) + inv_volume * &
110  (wnew(ixo^s, iw) - wprim(ixo^s, iw))
111  end if
112  end do
113 
114  ! check and optionally correct unphysical values
115  if(fix_small_values) then
116  call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'exit hancock finite_volume')
117  endif
118  end associate
119  end subroutine hancock
120 
121  !> finite volume method
122  subroutine finite_volume(method,qdt,ixI^L,ixO^L,idims^LIM, &
123  qtC,sCT,qt,snew,sold,fC,fE,dx^D,x)
126  use mod_tvd, only:tvdlimit2
127  use mod_source, only: addsource2
128  use mod_usr_methods
129 
130  integer, intent(in) :: method
131  double precision, intent(in) :: qdt, qtC, qt, dx^D
132  integer, intent(in) :: ixI^L, ixO^L, idims^LIM
133  double precision, dimension(ixI^S,1:ndim), intent(in) :: x
134  type(state) :: sCT, snew, sold
135  double precision, dimension(ixI^S,1:nwflux,1:ndim) :: fC
136  double precision, dimension(ixI^S,7-2*ndim:3) :: fE
137 
138  ! primitive w at cell center
139  double precision, dimension(ixI^S,1:nw) :: wprim
140  ! left and right constructed status in conservative form
141  double precision, dimension(ixI^S,1:nw) :: wLC, wRC
142  ! left and right constructed status in primitive form, needed for better performance
143  double precision, dimension(ixI^S,1:nw) :: wLp, wRp
144  double precision, dimension(ixI^S,1:nwflux) :: fLC, fRC
145  double precision, dimension(ixI^S) :: cmaxC
146  double precision, dimension(ixI^S) :: cminC
147  double precision, dimension(ixO^S) :: inv_volume
148  double precision, dimension(1:ndim) :: dxinv, dxdim
149  integer, dimension(ixI^S) :: patchf
150  integer :: idims, iw, ix^L, hxO^L, ixC^L, ixCR^L, kxC^L, kxR^L
151  type(ct_velocity) :: vcts
152 
153  associate(wct=>sct%w, wnew=>snew%w, wold=>sold%w)
154 
155  fc=0.d0
156  flc=0.d0
157  frc=0.d0
158 
159  ! The flux calculation contracts by one in the idims direction it is applied.
160  ! The limiter contracts the same directions by one more, so expand ixO by 2.
161  ix^l=ixo^l;
162  do idims= idims^lim
163  ix^l=ix^l^ladd2*kr(idims,^d);
164  end do
165  if (ixi^l^ltix^l|.or.|.or.) &
166  call mpistop("Error in fv : Nonconforming input limits")
167 
168  wprim=wct
169  call phys_to_primitive(ixi^l,ixi^l,wprim,x)
170 
171  ^d&dxinv(^d)=-qdt/dx^d;
172  ^d&dxdim(^d)=dx^d;
173  do idims= idims^lim
174  ! use interface value of w0 at idims
175  b0i=idims
176 
177  hxo^l=ixo^l-kr(idims,^d);
178 
179  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
180  kxr^l=kxc^l+kr(idims,^d);
181 
182  if(stagger_grid) then
183  ! ct needs all transverse cells
184  ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d); ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
185  else
186  ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
187  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
188  end if
189 
190  ! wRp and wLp are defined at the same locations, and will correspond to
191  ! the left and right reconstructed values at a cell face. Their indexing
192  ! is similar to cell-centered values, but in direction idims they are
193  ! shifted half a cell towards the 'lower' direction.
194  wrp(kxc^s,1:nw)=wprim(kxr^s,1:nw)
195  wlp(kxc^s,1:nw)=wprim(kxc^s,1:nw)
196 
197  ! Determine stencil size
198  {ixcrmin^d = max(ixcmin^d - phys_wider_stencil,ixglo^d)\}
199  {ixcrmax^d = min(ixcmax^d + phys_wider_stencil,ixghi^d)\}
200 
201  ! apply limited reconstruction for left and right status at cell interfaces
202  call reconstruct_lr(ixi^l,ixcr^l,ixcr^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxdim(idims))
203 
204  ! special modification of left and right status before flux evaluation
205  call phys_modify_wlr(ixi^l,ixcr^l,qt,wlc,wrc,wlp,wrp,sct,idims)
206 
207  ! evaluate physical fluxes according to reconstructed status
208  call phys_get_flux(wlc,wlp,x,ixi^l,ixc^l,idims,flc)
209  call phys_get_flux(wrc,wrp,x,ixi^l,ixc^l,idims,frc)
210 
211  ! estimating bounds for the minimum and maximum signal velocities
212  if(method==fs_tvdlf.or.method==fs_tvdmu) then
213  call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,cmaxc)
214  if(stagger_grid) call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc)
215  else
216  call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,cmaxc,cminc)
217  if(stagger_grid) call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc,cminc)
218  end if
219 
220  ! use approximate Riemann solver to get flux at interfaces
221  select case(method)
222  case(fs_hll)
223  call get_riemann_flux_hll()
224  case(fs_hllc,fs_hllcd)
225  call get_riemann_flux_hllc()
226  case(fs_hlld)
227  call get_riemann_flux_hlld()
228  case(fs_tvdlf)
229  call get_riemann_flux_tvdlf()
230  case(fs_tvdmu)
232  case default
233  call mpistop('unkown Riemann flux in finite volume')
234  end select
235 
236  end do ! Next idims
237  b0i=0
238 
239  if(stagger_grid) call phys_update_faces(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,snew,vcts)
240 
241  do idims= idims^lim
242  hxo^l=ixo^l-kr(idims,^d);
243 
244  ! Multiply the fluxes by -dt/dx since Flux fixing expects this
245  if (slab_uniform) then
246  fc(ixi^s,1:nwflux,idims)=dxinv(idims)*fc(ixi^s,1:nwflux,idims)
247 
248  do iw = iwstart, nwflux
249  if (associated(phys_iw_methods(iw)%inv_capacity)) then
250  call phys_iw_methods(iw)%inv_capacity(ixi^l, ixo^l, wnew, inv_volume)
251  wnew(ixo^s,iw)=wnew(ixo^s,iw) + inv_volume * &
252  (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))
253  else
254  wnew(ixo^s,iw)=wnew(ixo^s,iw) &
255  + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims))
256  end if
257  end do
258  else
259  if (.not. angmomfix) then ! default case
260  if (associated(phys_iw_methods(iw)%inv_capacity)) then
261  call phys_iw_methods(iw)%inv_capacity(ixi^l, ixo^l, wnew, inv_volume)
262  else
263  inv_volume = 1.0d0
264  end if
265  inv_volume = inv_volume/block%dvolume(ixo^s)
266 
267  do iw=iwstart,nwflux
268  fc(ixi^s,iw,idims)=-qdt*fc(ixi^s,iw,idims)*block%surfaceC(ixi^s,idims)
269  wnew(ixo^s,iw)=wnew(ixo^s,iw) + (fc(ixo^s,iw,idims)-fc(hxo^s,iw,idims)) * &
270  inv_volume
271  enddo
272  else
273  ! If angular momentum conserving way to solve the equations,
274  ! some fluxes additions need to be treated specifically
275  call phys_angmomfix(fc,x,wnew,ixi^l,ixo^l,idims)
276  endif
277  end if
278 
279  ! For the MUSCL scheme apply the characteristic based limiter
280  if (method==fs_tvdmu) &
281  call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dx^d)
282 
283  end do ! Next idims
284 
285  ! check and optionally correct unphysical values
286  if(fix_small_values) then
287  call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'multi-D finite_volume')
288  endif
289 
290  if (.not.slab.and.idimsmin==1) &
291  call phys_add_source_geom(qdt,ixi^l,ixo^l,wct,wnew,x)
292 
293  if(stagger_grid) call phys_face_to_center(ixo^l,snew)
294 
295  call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim), &
296  ixi^l,ixo^l,1,nw,qtc,wct,qt,wnew,x,.false.)
297 
298  if(phys_solve_eaux.and.levmin==levmax) then
299  ! synchronize internal energy for uniform grid
300  call phys_energy_synchro(ixi^l,ixo^l,wnew,x)
301  endif
302 
303  ! check and optionally correct unphysical values
304  if(fix_small_values) then
305  call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'exit finite_volume')
306  endif
307  end associate
308  contains
309 
310  subroutine get_riemann_flux_tvdmu()
312  do iw=iwstart,nwflux
313  ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
314  flc(ixc^s, iw)=half*(flc(ixc^s, iw)+frc(ixc^s, iw))
315  fc(ixc^s,iw,idims)=flc(ixc^s, iw)
316  end do
317  end subroutine get_riemann_flux_tvdmu
318 
319  subroutine get_riemann_flux_tvdlf()
320  double precision :: fac(ixc^s)
321 
322  fac = -0.5d0*tvdlfeps*cmaxc(ixc^s)
323 
324  ! Calculate fLC=f(uL_j+1/2) and fRC=f(uR_j+1/2) for each iw
325  do iw=iwstart,nwflux
326 
327  ! To save memory we use fLC to store (F_L+F_R)/2=half*(fLC+fRC)
328  flc(ixc^s, iw)=0.5d0*(flc(ixc^s, iw)+frc(ixc^s, iw))
329 
330  ! Add TVDLF dissipation to the flux
331  if (flux_type(idims, iw) /= flux_no_dissipation) then
332  flc(ixc^s, iw)=flc(ixc^s, iw) + fac*(wrc(ixc^s,iw)-wlc(ixc^s,iw))
333  end if
334 
335  fc(ixc^s,iw,idims)=flc(ixc^s, iw)
336  end do ! Next iw
337  end subroutine get_riemann_flux_tvdlf
338 
339  subroutine get_riemann_flux_hll()
340  integer :: ix^D
341 
342  do iw=iwstart,nwflux
343  if(flux_type(idims, iw) == flux_tvdlf) then
344  fc(ixc^s,iw,idims) = half*(flc(ixc^s, iw) + frc(ixc^s, iw) &
345  -tvdlfeps*max(cmaxc(ixc^s), dabs(cminc(ixc^s))) * &
346  (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
347  else
348  {do ix^db=ixcmin^db,ixcmax^db\}
349  if(cminc(ix^d) >= zero) then
350  fc(ix^d,iw,idims)=flc(ix^d,iw)
351  else if(cmaxc(ix^d) <= zero) then
352  fc(ix^d,iw,idims)=frc(ix^d,iw)
353  else
354  ! Add hll dissipation to the flux
355  fc(ix^d,iw,idims)=(cmaxc(ix^d)*flc(ix^d, iw)-cminc(ix^d)*frc(ix^d,iw)&
356  +tvdlfeps*cminc(ix^d)*cmaxc(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
357  /(cmaxc(ix^d)-cminc(ix^d))
358  end if
359  {end do\}
360  end if
361  end do
362 
363  end subroutine get_riemann_flux_hll
364 
365  subroutine get_riemann_flux_hllc()
366  use mod_mhd_phys
367  implicit none
368  double precision, dimension(ixI^S,1:nwflux) :: whll, Fhll, fCD
369  double precision, dimension(ixI^S) :: lambdaCD
370 
371  patchf(ixc^s) = 1
372  where(cminc(ixc^s) >= zero)
373  patchf(ixc^s) = -2
374  elsewhere(cmaxc(ixc^s) <= zero)
375  patchf(ixc^s) = 2
376  endwhere
377  ! Use more diffusive scheme, is actually TVDLF and selected by patchf=4
378  if(method==fs_hllcd) &
379  call phys_diffuse_hllcd(ixi^l,ixc^l,idims,wlc,wrc,flc,frc,patchf)
380 
381  !---- calculate speed lambda at CD ----!
382  if(any(patchf(ixc^s)==1)) &
383  call phys_get_lcd(wlc,wrc,flc,frc,cminc,cmaxc,idims,ixi^l,ixc^l, &
384  whll,fhll,lambdacd,patchf)
385 
386  ! now patchf may be -1 or 1 due to phys_get_lCD
387  if(any(abs(patchf(ixc^s))== 1))then
388  !======== flux at intermediate state ========!
389  call phys_get_wcd(wlc,wrc,whll,frc,flc,fhll,patchf,lambdacd,&
390  cminc,cmaxc,ixi^l,ixc^l,idims,fcd)
391  endif ! Calculate the CD flux
392 
393  ! use hll flux for the auxiliary internal e
394  if(mhd_energy.and.mhd_solve_eaux) then
395  iw=eaux_
396  fcd(ixc^s, iw) = (cmaxc(ixc^s)*flc(ixc^s, iw)-cminc(ixc^s) * frc(ixc^s, iw) &
397  +cminc(ixc^s)*cmaxc(ixc^s)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(cmaxc(ixc^s)-cminc(ixc^s))
398  end if
399 
400  do iw=iwstart,nwflux
401  if (flux_type(idims, iw) == flux_tvdlf) then
402  flc(ixc^s,iw) = 0.5d0 * (flc(ixc^s,iw) + frc(ixc^s,iw) - tvdlfeps * &
403  max(cmaxc(ixc^s), abs(cminc(ixc^s))) * &
404  (wrc(ixc^s,iw) - wlc(ixc^s,iw)))
405  else
406  where(patchf(ixc^s)==-2)
407  flc(ixc^s,iw)=flc(ixc^s,iw)
408  elsewhere(abs(patchf(ixc^s))==1)
409  flc(ixc^s,iw)=fcd(ixc^s,iw)
410  elsewhere(patchf(ixc^s)==2)
411  flc(ixc^s,iw)=frc(ixc^s,iw)
412  elsewhere(patchf(ixc^s)==3)
413  ! fallback option, reducing to HLL flux
414  flc(ixc^s,iw)=fhll(ixc^s,iw)
415  elsewhere(patchf(ixc^s)==4)
416  ! fallback option, reducing to TVDLF flux
417  flc(ixc^s,iw) = half*((flc(ixc^s,iw)+frc(ixc^s,iw)) &
418  -tvdlfeps * max(cmaxc(ixc^s), dabs(cminc(ixc^s))) * &
419  (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
420  endwhere
421  end if
422 
423  fc(ixc^s,iw,idims)=flc(ixc^s,iw)
424 
425  end do ! Next iw
426  end subroutine get_riemann_flux_hllc
427 
428  !> HLLD Riemann flux from Miyoshi 2005 JCP, 208, 315 and Guo 2016 JCP, 327, 543
429  subroutine get_riemann_flux_hlld()
431  implicit none
432  double precision, dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
433  double precision, dimension(ixI^S,1:nwflux) :: w2R,w2L
434  double precision, dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
435  double precision, dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
436  ! velocity from the right and the left reconstruction
437  double precision, dimension(ixI^S,ndir) :: vRC, vLC
438  ! magnetic field from the right and the left reconstruction
439  double precision, dimension(ixI^S,ndir) :: BR, BL
440  integer :: ip1,ip2,ip3,idir,ix^D
441 
442  associate(sr=>cmaxc,sl=>cminc)
443 
444  f1r=0.d0
445  f1l=0.d0
446  ip1=idims
447  ip3=3
448  vrc(ixc^s,:)=wrp(ixc^s,mom(:))
449  vlc(ixc^s,:)=wlp(ixc^s,mom(:))
450  if(b0field) then
451  br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
452  bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
453  else
454  br(ixc^s,:)=wrc(ixc^s,mag(:))
455  bl(ixc^s,:)=wlc(ixc^s,mag(:))
456  end if
457  ! HLL estimation of normal magnetic field at cell interfaces
458  bx(ixc^s)=(sr(ixc^s)*br(ixc^s,ip1)-sl(ixc^s)*bl(ixc^s,ip1)-&
459  flc(ixc^s,mag(ip1))-frc(ixc^s,mag(ip1)))/(sr(ixc^s)-sl(ixc^s))
460  ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
461  ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
462  sur(ixc^s)=(sr(ixc^s)-vrc(ixc^s,ip1))*wrc(ixc^s,rho_)
463  sul(ixc^s)=(sl(ixc^s)-vlc(ixc^s,ip1))*wlc(ixc^s,rho_)
464  ! Miyoshi equation (38) and Guo euqation (20)
465  sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
466  ptr(ixc^s)+ptl(ixc^s))/(sur(ixc^s)-sul(ixc^s))
467  ! Miyoshi equation (39) and Guo euqation (28)
468  w1r(ixc^s,mom(ip1))=sm(ixc^s)
469  w1l(ixc^s,mom(ip1))=sm(ixc^s)
470  w2r(ixc^s,mom(ip1))=sm(ixc^s)
471  w2l(ixc^s,mom(ip1))=sm(ixc^s)
472  ! Guo equation (22)
473  w1r(ixc^s,mag(ip1))=bx(ixc^s)
474  w1l(ixc^s,mag(ip1))=bx(ixc^s)
475  if(b0field) then
476  ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
477  ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
478  end if
479 
480  ! Miyoshi equation (43) and Guo equation (27)
481  w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s)-sm(ixc^s))
482  w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s)-sm(ixc^s))
483 
484  ip2=mod(ip1+1,ndir)
485  if(ip2==0) ip2=ndir
486  r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s)-sm(ixc^s))-bx(ixc^s)**2
487  where(r1r(ixc^s)/=0.d0)
488  r1r(ixc^s)=1.d0/r1r(ixc^s)
489  endwhere
490  r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s)-sm(ixc^s))-bx(ixc^s)**2
491  where(r1l(ixc^s)/=0.d0)
492  r1l(ixc^s)=1.d0/r1l(ixc^s)
493  endwhere
494  ! Miyoshi equation (44)
495  w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
496  (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
497  w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
498  (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
499  ! partial solution for later usage
500  w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
501  w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
502  if(ndir==3) then
503  ip3=mod(ip1+2,ndir)
504  if(ip3==0) ip3=ndir
505  ! Miyoshi equation (46)
506  w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
507  (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
508  w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
509  (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
510  ! Miyoshi equation (47)
511  w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
512  w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
513  end if
514  ! Miyoshi equation (45)
515  w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
516  w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
517  if(b0field) then
518  ! Guo equation (26)
519  w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
520  w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
521  end if
522  ! equation (48)
523  if(mhd_energy) then
524  ! Guo equation (25) equivalent to Miyoshi equation (41)
525  w1r(ixc^s,p_)=sur(ixc^s)*(sm(ixc^s)-vrc(ixc^s,ip1))+ptr(ixc^s)
526  w1l(ixc^s,p_)=sul(ixc^s)*(sm(ixc^s)-vlc(ixc^s,ip1))+ptl(ixc^s)
527  !if(mhd_solve_eaux) then
528  ! w1R(ixC^S,eaux_)=(w1R(ixC^S,p_)-half*sum(w1R(ixC^S,mag(:))**2,dim=ndim+1))/(mhd_gamma-one)
529  ! w1L(ixC^S,eaux_)=(w1L(ixC^S,p_)-half*sum(w1L(ixC^S,mag(:))**2,dim=ndim+1))/(mhd_gamma-one)
530  !end if
531  if(b0field) then
532  ! Guo equation (32)
533  w1r(ixc^s,p_)=w1r(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wrc(ixc^s,mag(:))-w1r(ixc^s,mag(:))),dim=ndim+1)
534  w1l(ixc^s,p_)=w1l(ixc^s,p_)+sum(block%B0(ixc^s,:,ip1)*(wlc(ixc^s,mag(:))-w1l(ixc^s,mag(:))),dim=ndim+1)
535  end if
536  ! Miyoshi equation (48) and main part of Guo euqation (31)
537  w1r(ixc^s,e_)=((sr(ixc^s)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
538  w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
539  sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s)-sm(ixc^s))
540  w1l(ixc^s,e_)=((sl(ixc^s)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
541  w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
542  sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s)-sm(ixc^s))
543  if(b0field) then
544  ! Guo equation (31)
545  w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
546  sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s)-sm(ixc^s))
547  w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
548  sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s)-sm(ixc^s))
549  end if
550  end if
551 
552  ! Miyoshi equation (49) and Guo equation (35)
553  w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
554  w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
555  w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
556  w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
557 
558  r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
559  r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
560  tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
561  signbx(ixc^s)=sign(1.d0,bx(ixc^s))
562  ! Miyoshi equation (51) and Guo equation (33)
563  s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
564  s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
565  ! Miyoshi equation (59) and Guo equation (41)
566  w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
567  (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
568  w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
569  ! Miyoshi equation (61) and Guo equation (43)
570  w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
571  r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
572  w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
573  if(ndir==3) then
574  ! Miyoshi equation (60) and Guo equation (42)
575  w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
576  (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
577  w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
578  ! Miyoshi equation (62) and Guo equation (44)
579  w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
580  r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
581  w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
582  end if
583  ! Miyoshi equation (63) and Guo equation (45)
584  if(mhd_energy) then
585  w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
586  sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
587  w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
588  sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
589  !if(mhd_solve_eaux) then
590  ! w2R(ixC^S,eaux_)=w2R(ixC^S,e_)-half*(sum(w2R(ixC^S,mag(:))**2,dim=ndim+1)+&
591  ! sum(w2R(ixC^S,mom(:))**2,dim=ndim+1)*w2R(ixC^S,rho_))
592  ! w2L(ixC^S,eaux_)=w2L(ixC^S,e_)-half*(sum(w2L(ixC^S,mag(:))**2,dim=ndim+1)+&
593  ! sum(w2L(ixC^S,mom(:))**2,dim=ndim+1)*w2L(ixC^S,rho_))
594  !end if
595  end if
596 
597  ! convert velocity to momentum
598  do idir=1,ndir
599  w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
600  w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
601  w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
602  w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
603  end do
604 
605  ! get fluxes of intermedate states
606  do iw=1,nwflux
607  if (flux_type(idims, iw) == flux_tvdlf) then
608  !! hll flux for normal B
609  !f1L(ixC^S,iw)=(sR(ixC^S)*fLC(ixC^S, iw)-sL(ixC^S)*fRC(ixC^S, iw) &
610  ! +sR(ixC^S)*sL(ixC^S)*(wRC(ixC^S,iw)-wLC(ixC^S,iw)))/(sR(ixC^S)-sL(ixC^S))
611  ! tvldf flux for normal B
612  f1l(ixc^s,iw)= half*(flc(ixc^s, iw) + frc(ixc^s, iw) &
613  -tvdlfeps*max(sr(ixc^s), dabs(sl(ixc^s))) * &
614  (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
615  f1r(ixc^s,iw)=f1l(ixc^s,iw)
616  f2l(ixc^s,iw)=f1l(ixc^s,iw)
617  f2r(ixc^s,iw)=f1l(ixc^s,iw)
618  else if(flux_type(idims, iw) == flux_special) then
619  ! known flux (fLC=fRC) for normal B and psi_ in GLM method
620  f1l(ixc^s,iw)=flc(ixc^s,iw)
621  f1r(ixc^s,iw)=f1l(ixc^s,iw)
622  f2l(ixc^s,iw)=f1l(ixc^s,iw)
623  f2r(ixc^s,iw)=f1l(ixc^s,iw)
624  else if(flux_type(idims, iw) == flux_hll) then
625  ! using hll flux for eaux and tracers
626  f1l(ixc^s,iw)=(sr(ixc^s)*flc(ixc^s, iw)-sl(ixc^s)*frc(ixc^s, iw) &
627  +sr(ixc^s)*sl(ixc^s)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s)-sl(ixc^s))
628  f1r(ixc^s,iw)=f1l(ixc^s,iw)
629  f2l(ixc^s,iw)=f1l(ixc^s,iw)
630  f2r(ixc^s,iw)=f1l(ixc^s,iw)
631  else
632  f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
633  f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
634  f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
635  f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
636  end if
637  end do
638 
639  ! Miyoshi equation (66) and Guo equation (46)
640  {do ix^db=ixcmin^db,ixcmax^db\}
641  if(sl(ix^d)>0.d0) then
642  fc(ix^d,1:nwflux,ip1)=flc(ix^d,1:nwflux)
643  else if(s1l(ix^d)>=0.d0) then
644  fc(ix^d,1:nwflux,ip1)=f1l(ix^d,1:nwflux)
645  else if(sm(ix^d)>=0.d0) then
646  fc(ix^d,1:nwflux,ip1)=f2l(ix^d,1:nwflux)
647  else if(s1r(ix^d)>=0.d0) then
648  fc(ix^d,1:nwflux,ip1)=f2r(ix^d,1:nwflux)
649  else if(sr(ix^d)>=0.d0) then
650  fc(ix^d,1:nwflux,ip1)=f1r(ix^d,1:nwflux)
651  else if(sr(ix^d)<0.d0) then
652  fc(ix^d,1:nwflux,ip1)=frc(ix^d,1:nwflux)
653  end if
654  {end do\}
655 
656  end associate
657  end subroutine get_riemann_flux_hlld
658 
659  end subroutine finite_volume
660 
661  !> Determine the upwinded wLC(ixL) and wRC(ixR) from w.
662  !> the wCT is only used when PPM is exploited.
663  subroutine reconstruct_lr(ixI^L,ixL^L,ixR^L,idims,w,wLC,wRC,wLp,wRp,x,dxdim)
666  use mod_limiter
667 
668  integer, intent(in) :: ixI^L, ixL^L, ixR^L, idims
669  double precision, intent(in) :: dxdim
670  double precision, dimension(ixI^S,1:nw) :: w
671  ! left and right constructed status in conservative form
672  double precision, dimension(ixI^S,1:nw) :: wLC, wRC
673  ! left and right constructed status in primitive form
674  double precision, dimension(ixI^S,1:nw) :: wLp, wRp
675  double precision, dimension(ixI^S,1:ndim) :: x
676 
677  integer :: jxR^L, ixC^L, jxC^L, iw
678  double precision :: ldw(ixi^s), rdw(ixi^s), dwC(ixi^s)
679  double precision :: a2max
680 
681  select case (type_limiter(block%level))
682  case (limiter_mp5)
683  call mp5limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
684  case (limiter_weno3)
685  call weno3limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
686  case (limiter_wenoyc3)
687  call weno3limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
688  case (limiter_weno5)
689  call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
690  case (limiter_weno5nm)
691  call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
692  case (limiter_wenoz5)
693  call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
694  case (limiter_wenoz5nm)
695  call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
696  case (limiter_wenozp5)
697  call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,3)
698  case (limiter_wenozp5nm)
699  call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,3)
700  case (limiter_weno5cu6)
701  call weno5cu6limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
702  case (limiter_teno5ad)
703  call teno5adlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp)
704  case (limiter_weno7)
705  call weno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp,1)
706  case (limiter_mpweno7)
707  call weno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp,2)
708  case (limiter_exeno7)
709  call exeno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
710  case (limiter_venk)
711  call venklimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp)
712  if(fix_small_values) then
713  call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
714  call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
715  end if
716  case (limiter_ppm)
717  call ppmlimiter(ixi^l,ixm^ll,idims,w,w,wlp,wrp)
718  if(fix_small_values) then
719  call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
720  call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
721  end if
722  case default
723  jxr^l=ixr^l+kr(idims,^d);
724  ixcmax^d=jxrmax^d; ixcmin^d=ixlmin^d-kr(idims,^d);
725  jxc^l=ixc^l+kr(idims,^d);
726  do iw=1,nwflux
727  if (loglimit(iw)) then
728  w(ixcmin^d:jxcmax^d,iw)=dlog10(w(ixcmin^d:jxcmax^d,iw))
729  wlp(ixl^s,iw)=dlog10(wlp(ixl^s,iw))
730  wrp(ixr^s,iw)=dlog10(wrp(ixr^s,iw))
731  end if
732 
733  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
734  if(need_global_a2max) then
735  a2max=a2max_global(idims)
736  else
737  select case(idims)
738  case(1)
739  a2max=schmid_rad1
740  {^iftwod
741  case(2)
742  a2max=schmid_rad2}
743  {^ifthreed
744  case(2)
745  a2max=schmid_rad2
746  case(3)
747  a2max=schmid_rad3}
748  case default
749  call mpistop("idims is wrong in mod_limiter")
750  end select
751  end if
752 
753  ! limit flux from left and/or right
754  call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),ldw,rdw,a2max=a2max)
755  wlp(ixl^s,iw)=wlp(ixl^s,iw)+half*ldw(ixl^s)
756  wrp(ixr^s,iw)=wrp(ixr^s,iw)-half*rdw(jxr^s)
757 
758  if (loglimit(iw)) then
759  w(ixcmin^d:jxcmax^d,iw)=10.0d0**w(ixcmin^d:jxcmax^d,iw)
760  wlp(ixl^s,iw)=10.0d0**wlp(ixl^s,iw)
761  wrp(ixr^s,iw)=10.0d0**wrp(ixr^s,iw)
762  end if
763  end do
764  if(fix_small_values) then
765  call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
766  call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
767  end if
768  end select
769 
770  wlc(ixl^s,1:nw)=wlp(ixl^s,1:nw)
771  wrc(ixr^s,1:nw)=wrp(ixr^s,1:nw)
772  call phys_to_conserved(ixi^l,ixl^l,wlc,x)
773  call phys_to_conserved(ixi^l,ixr^l,wrc,x)
774 
775  end subroutine reconstruct_lr
776 
777 end module mod_finite_volume
type(iw_methods), dimension(max_nw) phys_iw_methods
Special methods defined per variable.
Definition: mod_physics.t:59
This module contains definitions of global parameters and variables and some generic functions/subrou...
Module for handling split source terms (split from the fluxes)
Definition: mod_source.t:2
logical need_global_a2max
global value for schmid scheme
integer, public, protected e_
Index of the energy density (-1 if not present)
Definition: mod_mhd_phys.t:98
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
Definition: mod_mhd_phys.t:101
procedure(sub_energy_synchro), pointer phys_energy_synchro
Definition: mod_physics.t:72
integer b0i
background magnetic field location indicator
subroutine, public hancock(qdt, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, dxD, x)
The non-conservative Hancock predictor for TVDLF.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, parameter limiter_wenoz5
Definition: mod_limiter.t:33
integer, parameter limiter_wenozp5
Definition: mod_limiter.t:35
integer, parameter limiter_mp5
Definition: mod_limiter.t:28
subroutine get_riemann_flux_tvdmu()
integer, parameter limiter_wenoyc3
Definition: mod_limiter.t:30
integer, parameter limiter_mpweno7
Definition: mod_limiter.t:40
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:88
procedure(sub_small_values), pointer phys_handle_small_values
Definition: mod_physics.t:85
double precision, dimension(ndim) a2max_global
global largest a2 for schmid scheme
integer, parameter limiter_teno5ad
Definition: mod_limiter.t:38
subroutine, public finite_volume(method, qdt, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, sold, fC, fE, dxD, x)
finite volume method
integer, parameter fs_hllc
integer, parameter ixglo
Lower index of grid block arrays (always 1)
integer, parameter limiter_weno7
Definition: mod_limiter.t:39
integer, parameter limiter_ppm
Definition: mod_limiter.t:27
logical, public, protected mhd_energy
Whether an energy equation is used.
Definition: mod_mhd_phys.t:8
Module with finite volume methods for fluxes.
Module with slope/flux limiters.
Definition: mod_limiter.t:2
procedure(sub_update_faces), pointer phys_update_faces
Definition: mod_physics.t:87
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter limiter_wenozp5nm
Definition: mod_limiter.t:36
Module with all the methods that users can customize in AMRVAC.
logical, dimension(:), allocatable loglimit
integer, parameter limiter_weno5nm
Definition: mod_limiter.t:32
integer ixghi
Upper index of grid block arrays.
subroutine, public reconstruct_lr(ixIL, ixLL, ixRL, idims, w, wLC, wRC, wLp, wRp, x, dxdim)
Determine the upwinded wLC(ixL) and wRC(ixR) from w. the wCT is only used when PPM is exploited...
integer, dimension(:), allocatable, public, protected mag
Indices of the magnetic field.
Definition: mod_mhd_phys.t:104
logical angmomfix
Enable to strictly conserve the angular momentum (works both in cylindrical and spherical coordinates...
integer, parameter fs_hlld
logical stagger_grid
True for using stagger grid.
procedure(sub_angmomfix), pointer phys_angmomfix
Definition: mod_physics.t:84
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:62
integer ixm
the mesh range (within a block with ghost cells)
integer, parameter fs_tvdmu
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_mhd_phys.t:95
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
integer, parameter limiter_weno5cu6
Definition: mod_limiter.t:37
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:71
logical, public, protected mhd_solve_eaux
Whether auxiliary internal energy is solved.
Definition: mod_mhd_phys.t:61
subroutine get_riemann_flux_hlld()
HLLD Riemann flux from Miyoshi 2005 JCP, 208, 315 and Guo 2016 JCP, 327, 543.
integer, parameter limiter_venk
Definition: mod_limiter.t:25
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_mhd_phys.t:92
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:75
integer, parameter fs_hll
flux schemes
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:131
integer, parameter limiter_exeno7
Definition: mod_limiter.t:41
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, parameter fs_hllcd
integer, parameter limiter_weno5
Definition: mod_limiter.t:31
procedure(sub_get_ct_velocity), pointer phys_get_ct_velocity
Definition: mod_physics.t:86
type(state), pointer block
Block pointer for using one block and its previous state.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer phys_wider_stencil
To use wider stencils in flux calculations. A value of 1 will extend it by one cell in both direction...
Definition: mod_physics.t:18
logical phys_solve_eaux
Solve internal energy and total energy equations.
Definition: mod_physics.t:34
logical slab
Cartesian geometry or not.
integer, parameter fs_tvdlf
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
subroutine, public tvdlimit2(method, qdt, ixIL, ixICL, ixOL, idims, wL, wR, wnew, x, fC, dxD)
Definition: mod_tvd.t:39
Subroutines for TVD-MUSCL schemes.
Definition: mod_tvd.t:2
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:63
Magneto-hydrodynamics module.
Definition: mod_mhd_phys.t:2
integer, public, protected eaux_
Indices of auxiliary internal energy.
Definition: mod_mhd_phys.t:110
integer, parameter limiter_weno3
Definition: mod_limiter.t:29
subroutine, public addsource2(qdt, ixIL, ixOL, iwLIM, qtC, wCT, qt, w, x, qsourcesplit, src_active)
Add source within ixO for iws: w=w+qdt*S[wCT].
Definition: mod_source.t:121
integer, parameter limiter_wenoz5nm
Definition: mod_limiter.t:34
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:70
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition: mod_physics.t:66