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