MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_finite_volume.t
Go to the documentation of this file.
1 !> Module with finite volume methods for fluxes
3 #include "amrvac.h"
4  implicit none
5  private
6 
7  public :: finite_volume
8  public :: hancock
9  public :: reconstruct_lr
10 
11 contains
12 
13  !> The non-conservative Hancock predictor for TVDLF
14  !>
15  !> on entry:
16  !> input available on ixI^L=ixG^L asks for output on ixO^L=ixG^L^LSUBnghostcells
17  !> one entry: (predictor): wCT -- w_n wnew -- w_n qdt=dt/2
18  !> on exit : (predictor): wCT -- w_n wnew -- w_n+1/2
19  subroutine hancock(qdt,dtfactor,ixI^L,ixO^L,idims^LIM,qtC,sCT,qt,snew,dxs,x)
20  use mod_physics
22  use mod_source, only: addsource2
23  use mod_comm_lib, only: mpistop
24 
25  integer, intent(in) :: ixi^l, ixo^l, idims^lim
26  double precision, intent(in) :: qdt, dtfactor,qtc, qt, dxs(ndim), x(ixi^s,1:ndim)
27  type(state) :: sct, snew
28 
29  double precision, dimension(ixI^S,1:nw) :: wprim, wlc, wrc
30  ! left and right constructed status in primitive form, needed for better performance
31  double precision, dimension(ixI^S,1:nw) :: wlp, wrp
32  double precision, dimension(ixO^S) :: inv_volume
33  double precision :: flc(ixi^s, nwflux), frc(ixi^s, nwflux)
34  double precision :: dxinv(1:ndim)
35  integer :: idims, iw, ix^l, hxo^l
36  logical :: active
37 
38  associate(wct=>sct%w,wnew=>snew%w)
39  ! Expand limits in each idims direction in which fluxes are added
40  ix^l=ixo^l;
41  do idims= idims^lim
42  ix^l=ix^l^laddkr(idims,^d);
43  end do
44  if (ixi^l^ltix^l|.or.|.or.) &
45  call mpistop("Error in Hancock: Nonconforming input limits")
46 
47  wrp=0.d0
48  wlp=0.d0
49  wprim=wct
50  call phys_to_primitive(ixi^l,ixi^l,wprim,x)
51 
52  dxinv=-qdt/dxs
53  if(.not.slab_uniform) inv_volume = 1.d0/block%dvolume(ixo^s)
54  do idims= idims^lim
55  b0i=idims
56  ! Calculate w_j+g_j/2 and w_j-g_j/2
57  ! First copy all variables, then upwind wLC and wRC.
58  ! wLC is to the left of ixO, wRC is to the right of wCT.
59  hxo^l=ixo^l-kr(idims,^d);
60 
61  wrp(hxo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
62  wlp(ixo^s,1:nwflux)=wprim(ixo^s,1:nwflux)
63 
64  ! apply limited reconstruction for left and right status at cell interfaces
65  call reconstruct_lr(ixi^l,ixo^l,hxo^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
66 
67  ! Calculate the fLC and fRC fluxes
68  call phys_get_flux(wrc,wrp,x,ixi^l,hxo^l,idims,frc)
69  call phys_get_flux(wlc,wlp,x,ixi^l,ixo^l,idims,flc)
70 
71  ! Advect w(iw)
72  if (slab_uniform) then
73  if(local_timestep) then
74  do iw=1,nwflux
75  wnew(ixo^s,iw)=wnew(ixo^s,iw)-block%dt(ixo^s)*dtfactor/dxs(idims)* &
76  (flc(ixo^s, iw)-frc(hxo^s, iw))
77  end do
78  else
79  do iw=1,nwflux
80  wnew(ixo^s,iw)=wnew(ixo^s,iw)+dxinv(idims)* &
81  (flc(ixo^s, iw)-frc(hxo^s, iw))
82  end do
83  endif
84  else
85  if(local_timestep) then
86  do iw=1,nwflux
87  wnew(ixo^s,iw)=wnew(ixo^s,iw) - block%dt(ixo^s)*dtfactor * inv_volume &
88  *(block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
89  -block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
90  end do
91  else
92  do iw=1,nwflux
93  wnew(ixo^s,iw)=wnew(ixo^s,iw) - qdt * inv_volume &
94  *(block%surfaceC(ixo^s,idims)*flc(ixo^s, iw) &
95  -block%surfaceC(hxo^s,idims)*frc(hxo^s, iw))
96  end do
97  end if
98  end if
99  end do ! next idims
100  b0i=0
101 
102  if (.not.slab.and.idimsmin==1) call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wprim,wnew,x)
103 
104  call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim), &
105  dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim),&
106  ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
107 
108  ! check and optionally correct unphysical values
109  if(fix_small_values) then
110  call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'exit hancock finite_volume')
111  endif
112  end associate
113  end subroutine hancock
114 
115  !> finite volume method
116  subroutine finite_volume(method,qdt,dtfactor,ixI^L,ixO^L,idims^LIM, &
117  qtC,sCT,qt,snew,fC,fE,dxs,x)
118  use mod_physics
119  use mod_variables
121  use mod_tvd, only:tvdlimit2
122  use mod_source, only: addsource2
123  use mod_usr_methods
124  use mod_comm_lib, only: mpistop
125 
126  integer, intent(in) :: method
127  double precision, intent(in) :: qdt, dtfactor, qtc, qt, dxs(ndim)
128  integer, intent(in) :: ixi^l, ixo^l, idims^lim
129  double precision, dimension(ixI^S,1:ndim), intent(in) :: x
130  type(state) :: sct, snew
131  double precision, dimension(ixI^S,1:nwflux,1:ndim) :: fc
132  double precision, dimension(ixI^S,sdim:3) :: fe
133 
134  ! primitive w at cell center
135  logical :: active
136  integer :: idims, iw, ix^d, hx^d, ix^l, hxo^l, ixc^l, ixcr^l, kxc^l, kxr^l, ii
137  integer, dimension(ixI^S) :: patchf
138  double precision, dimension(ixI^S,1:nw) :: wprim
139  ! left and right constructed status in conservative form
140  double precision, dimension(ixI^S,1:nw) :: wlc, wrc
141  ! left and right constructed status in primitive form, needed for better performance
142  double precision, dimension(ixI^S,1:nw) :: wlp, wrp
143  double precision, dimension(ixI^S,1:nwflux) :: flc, frc
144  double precision, dimension(ixI^S,1:number_species) :: cmaxc
145  double precision, dimension(ixI^S,1:number_species) :: cminc
146  double precision, dimension(ixI^S) :: hspeed
147  double precision, dimension(ixO^S) :: inv_volume
148  double precision, dimension(1:ndim) :: dxinv
149  type(ct_velocity) :: vcts
150 
151  associate(wct=>sct%w, wnew=>snew%w)
152 
153  fc=0.d0
154  flc=0.d0
155  frc=0.d0
156  wlp=0.d0
157  wrp=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  do idims= idims^lim
172  ! use interface value of w0 at idims
173  b0i=idims
174 
175  hxo^l=ixo^l-kr(idims,^d);
176 
177  kxcmin^d=iximin^d; kxcmax^d=iximax^d-kr(idims,^d);
178  kxr^l=kxc^l+kr(idims,^d);
179 
180  if(stagger_grid) then
181  ! ct needs all transverse cells
182  ixcmax^d=ixomax^d+nghostcells-nghostcells*kr(idims,^d);
183  ixcmin^d=hxomin^d-nghostcells+nghostcells*kr(idims,^d);
184  else
185  ! ixC is centered index in the idims direction from ixOmin-1/2 to ixOmax+1/2
186  ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
187  end if
188 
189  ! wRp and wLp are defined at the same locations, and will correspond to
190  ! the left and right reconstructed values at a cell face. Their indexing
191  ! is similar to cell-centered values, but in direction idims they are
192  ! shifted half a cell towards the 'lower' direction.
193  wrp(kxc^s,1:nw)=wprim(kxr^s,1:nw)
194  wlp(kxc^s,1:nw)=wprim(kxc^s,1:nw)
195 
196  ! Determine stencil size
197  {ixcrmin^d = max(ixcmin^d - phys_wider_stencil,ixglo^d)\}
198  {ixcrmax^d = min(ixcmax^d + phys_wider_stencil,ixghi^d)\}
199  ! apply limited reconstruction for left and right status at cell interfaces
200  call reconstruct_lr(ixi^l,ixcr^l,ixcr^l,idims,wprim,wlc,wrc,wlp,wrp,x,dxs(idims))
201 
202  ! special modification of left and right status before flux evaluation
203  call phys_modify_wlr(ixi^l,ixcr^l,qt,wlc,wrc,wlp,wrp,sct,idims)
204 
205  ! evaluate physical fluxes according to reconstructed status
206  call phys_get_flux(wlc,wlp,x,ixi^l,ixc^l,idims,flc)
207  call phys_get_flux(wrc,wrp,x,ixi^l,ixc^l,idims,frc)
208 
209  if(h_correction) then
210  call phys_get_h_speed(wprim,x,ixi^l,ixo^l,idims,hspeed)
211  end if
212  ! estimating bounds for the minimum and maximum signal velocities
213  if(method==fs_tvdlf.or.method==fs_tvdmu) then
214  call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,hspeed,cmaxc)
215  ! index of var velocity appears in the induction eq.
216  if(stagger_grid) call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc(ixi^s,index_v_mag))
217  else
218  call phys_get_cbounds(wlc,wrc,wlp,wrp,x,ixi^l,ixc^l,idims,hspeed,cmaxc,cminc)
219  if(stagger_grid) call phys_get_ct_velocity(vcts,wlp,wrp,ixi^l,ixc^l,idims,cmaxc(ixi^s,index_v_mag),cminc(ixi^s,index_v_mag))
220  end if
221 
222  ! use approximate Riemann solver to get flux at interfaces
223  select case(method)
224  case(fs_hll)
225  do ii=1,number_species
226  call get_riemann_flux_hll(start_indices(ii),stop_indices(ii))
227  end do
228  case(fs_hllc,fs_hllcd)
229  do ii=1,number_species
230  call get_riemann_flux_hllc(start_indices(ii),stop_indices(ii))
231  end do
232  case(fs_hlld)
233  do ii=1,number_species
234  if(ii==index_v_mag) then
235  call get_riemann_flux_hlld(start_indices(ii),stop_indices(ii))
236  else
237  call get_riemann_flux_hll(start_indices(ii),stop_indices(ii))
238  endif
239  end do
240  case(fs_tvdlf)
241  do ii=1,number_species
242  call get_riemann_flux_tvdlf(start_indices(ii),stop_indices(ii))
243  end do
244  case(fs_tvdmu)
246  case default
247  call mpistop('unkown Riemann flux in finite volume')
248  end select
249 
250  end do ! Next idims
251  b0i=0
252  if(stagger_grid) call phys_update_faces(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,snew,vcts)
253  if(slab_uniform) then
254  if(local_timestep) then
255  dxinv(1:ndim)=-dtfactor/dxs(1:ndim)
256  do idims= idims^lim
257  hxomin^d=ixomin^d-kr(idims,^d)\
258  do iw=iwstart,nwflux
259  {do ix^db=hxomin^db,ixomax^db\}
260  fc(ix^d,iw,idims)=block%dt(ix^d)*dxinv(idims)*fc(ix^d,iw,idims)
261  {end do\}
262  {do ix^db=ixomin^db,ixomax^db\}
263  hx^d=ix^d-kr(idims,^d)\
264  wnew(ix^d,iw)=wnew(ix^d,iw)+fc(ix^d,iw,idims)-fc(hx^d,iw,idims)
265  {end do\}
266  end do
267  ! For the MUSCL scheme apply the characteristic based limiter
268  if(method==fs_tvdmu) &
269  call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
270  end do
271  else
272  dxinv(1:ndim)=-qdt/dxs(1:ndim)
273  do idims= idims^lim
274  hxomin^d=ixomin^d-kr(idims,^d)\
275  do iw=iwstart,nwflux
276  {do ix^db=hxomin^db,ixomax^db\}
277  fc(ix^d,iw,idims)=dxinv(idims)*fc(ix^d,iw,idims)
278  {end do\}
279  {do ix^db=ixomin^db,ixomax^db\}
280  hx^d=ix^d-kr(idims,^d)\
281  wnew(ix^d,iw)=wnew(ix^d,iw)+fc(ix^d,iw,idims)-fc(hx^d,iw,idims)
282  {end do\}
283  end do
284  ! For the MUSCL scheme apply the characteristic based limiter
285  if(method==fs_tvdmu) &
286  call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
287  end do
288  end if
289  else
290  inv_volume(ixo^s) = 1.d0/block%dvolume(ixo^s)
291  if(local_timestep) then
292  do idims= idims^lim
293  hxomin^d=ixomin^d-kr(idims,^d)\
294  do iw=iwstart,nwflux
295  {do ix^db=hxomin^db,ixomax^db\}
296  fc(ix^d,iw,idims)=-block%dt(ix^d)*dtfactor*fc(ix^d,iw,idims)*block%surfaceC(ix^d,idims)
297  {end do\}
298  {do ix^db=ixomin^db,ixomax^db\}
299  hx^d=ix^d-kr(idims,^d)\
300  wnew(ix^d,iw)=wnew(ix^d,iw)+(fc(ix^d,iw,idims)-fc(hx^d,iw,idims))*inv_volume(ix^d)
301  {end do\}
302  end do
303  ! For the MUSCL scheme apply the characteristic based limiter
304  if (method==fs_tvdmu) &
305  call tvdlimit2(method,qdt,ixi^l,ixc^l,ixo^l,idims,wlc,wrc,wnew,x,fc,dxs)
306  end do
307  else
308  do idims= idims^lim
309  hxomin^d=ixomin^d-kr(idims,^d)\
310  do iw=iwstart,nwflux
311  {do ix^db=hxomin^db,ixomax^db\}
312  fc(ix^d,iw,idims)=-qdt*fc(ix^d,iw,idims)*block%surfaceC(ix^d,idims)
313  {end do\}
314  {do ix^db=ixomin^db,ixomax^db\}
315  hx^d=ix^d-kr(idims,^d)\
316  wnew(ix^d,iw)=wnew(ix^d,iw)+(fc(ix^d,iw,idims)-fc(hx^d,iw,idims))*inv_volume(ix^d)
317  {end do\}
318  end do
319  end do
320  end if
321  end if
322 
323  if (.not.slab.and.idimsmin==1) &
324  call phys_add_source_geom(qdt,dtfactor,ixi^l,ixo^l,wct,wprim,wnew,x)
325 
326  if(stagger_grid) call phys_face_to_center(ixo^l,snew)
327 
328  ! check and optionally correct unphysical values
329  if(fix_small_values) then
330  call phys_handle_small_values(.false.,wnew,x,ixi^l,ixo^l,'multi-D finite_volume')
331  if(crash) then
332  ! replace erroneous values with values at previous step
333  wnew=wct
334  if(stagger_grid) snew%ws=sct%ws
335  end if
336  end if
337 
338  call addsource2(qdt*dble(idimsmax-idimsmin+1)/dble(ndim),&
339  dtfactor*dble(idimsmax-idimsmin+1)/dble(ndim),&
340  ixi^l,ixo^l,1,nw,qtc,wct,wprim,qt,wnew,x,.false.,active)
341 
342  end associate
343  contains
344 
346  do iw=iwstart,nwflux
347  fc(ixc^s,iw,idims)=half*(flc(ixc^s,iw)+frc(ixc^s,iw))
348  end do
349  end subroutine get_riemann_flux_tvdmu
350 
351  subroutine get_riemann_flux_tvdlf(iws,iwe)
352  integer, intent(in) :: iws,iwe
353 
354  integer :: ix^D,jx^D
355  double precision :: fac(ixC^S),phi
356 
357  fac(ixc^s) = -0.5d0*tvdlfeps*cmaxc(ixc^s,ii)
358  do iw=iws,iwe
359  fc(ixc^s,iw,idims)=0.5d0*(flc(ixc^s, iw)+frc(ixc^s, iw))
360  ! Add TVDLF dissipation to the flux
361  if(flux_type(idims, iw) /= flux_no_dissipation) then
362  if(flux_adaptive_diffusion) then
363  {do ix^db=ixcmin^db,ixcmax^db\}
364  jx^d=ix^d+kr(idims,^d)\
365  !> adaptive diffusion from Rempel et al. 2009, see also Rempel et al. 2014
366  !> the previous version is adopt
367  if(((wrc(ix^d,iw)-wlc(ix^d,iw))*(sct%w(jx^d,iw)-sct%w(ix^d,iw))) .gt. 1.e-18) then
368  phi=min((wrc(ix^d,iw)-wlc(ix^d,iw))**2/((sct%w(jx^d,iw)-sct%w(ix^d,iw))**2+1.e-18),one)
369  else
370  phi=1.d0
371  end if
372  fc(ix^d,iw,idims)=fc(ix^d,iw,idims)+fac(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw))*phi
373  {end do\}
374  else
375  {do ix^db=ixcmin^db,ixcmax^db\}
376  fc(ix^d,iw,idims)=fc(ix^d,iw,idims)+fac(ix^d)*(wrc(ix^d,iw)-wlc(ix^d,iw))
377  {end do\}
378  end if
379  end if
380  end do
381 
382  end subroutine get_riemann_flux_tvdlf
383 
384  subroutine get_riemann_flux_hll(iws,iwe)
385  integer, intent(in) :: iws,iwe
386  integer :: ix^D
387  double precision :: phi
388 
389  do iw=iws,iwe
390  if(flux_type(idims, iw) == flux_tvdlf) then
391  ! CT MHD does not need normal B flux
392  if(stagger_grid) cycle
393  fc(ixc^s,iw,idims) = -tvdlfeps*half*max(cmaxc(ixc^s,ii),dabs(cminc(ixc^s,ii))) * &
394  (wrc(ixc^s,iw)-wlc(ixc^s,iw))
395  else
396  if(flux_adaptive_diffusion) then
397  {do ix^db=ixcmin^db,ixcmax^db\}
398  if(cminc(ix^d,ii) >= zero) then
399  fc(ix^d,iw,idims)=flc(ix^d,iw)
400  else if(cmaxc(ix^d,ii) <= zero) then
401  fc(ix^d,iw,idims)=frc(ix^d,iw)
402  else
403  !> reduced diffusion is from Wang et al. 2024
404  phi=max(abs(cmaxc(ix^d,ii)),abs(cminc(ix^d,ii)))/(cmaxc(ix^d,ii)-cminc(ix^d,ii))
405  fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
406  +phi*cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
407  /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
408  end if
409  {end do\}
410  else
411  {do ix^db=ixcmin^db,ixcmax^db\}
412  if(cminc(ix^d,ii) >= zero) then
413  fc(ix^d,iw,idims)=flc(ix^d,iw)
414  else if(cmaxc(ix^d,ii) <= zero) then
415  fc(ix^d,iw,idims)=frc(ix^d,iw)
416  else
417  fc(ix^d,iw,idims)=(cmaxc(ix^d,ii)*flc(ix^d, iw)-cminc(ix^d,ii)*frc(ix^d,iw)&
418  +cminc(ix^d,ii)*cmaxc(ix^d,ii)*(wrc(ix^d,iw)-wlc(ix^d,iw)))&
419  /(cmaxc(ix^d,ii)-cminc(ix^d,ii))
420  end if
421  {end do\}
422  end if
423  end if
424  end do
425  end subroutine get_riemann_flux_hll
426 
427  subroutine get_riemann_flux_hllc(iws,iwe)
428  integer, intent(in) :: iws, iwe
429  double precision, dimension(ixI^S,1:nwflux) :: whll, Fhll, fCD
430  double precision, dimension(ixI^S) :: lambdaCD
431 
432  integer :: rho_, p_, e_, mom(1:ndir)
433 
434  rho_ = iw_rho
435  if (allocated(iw_mom)) mom(:) = iw_mom(:)
436  e_ = iw_e
437 
438  if(associated(phys_hllc_init_species)) then
439  call phys_hllc_init_species(ii, rho_, mom(:), e_)
440  endif
441 
442  p_ = e_
443 
444  patchf(ixc^s) = 1
445  where(cminc(ixc^s,1) >= zero)
446  patchf(ixc^s) = -2
447  elsewhere(cmaxc(ixc^s,1) <= zero)
448  patchf(ixc^s) = 2
449  endwhere
450  ! Use more diffusive scheme, is actually TVDLF and selected by patchf=4
451  if(method==fs_hllcd) &
452  call phys_diffuse_hllcd(ixi^l,ixc^l,idims,wlc,wrc,flc,frc,patchf)
453 
454  !---- calculate speed lambda at CD ----!
455  if(any(patchf(ixc^s)==1)) &
456  call phys_get_lcd(wlc,wrc,flc,frc,cminc(ixi^s,ii),cmaxc(ixi^s,ii),idims,ixi^l,ixc^l, &
457  whll,fhll,lambdacd,patchf)
458 
459  ! now patchf may be -1 or 1 due to phys_get_lCD
460  if(any(abs(patchf(ixc^s))== 1))then
461  !======== flux at intermediate state ========!
462  call phys_get_wcd(wlc,wrc,whll,frc,flc,fhll,patchf,lambdacd,&
463  cminc(ixi^s,ii),cmaxc(ixi^s,ii),ixi^l,ixc^l,idims,fcd)
464  endif ! Calculate the CD flux
465 
466  do iw=iws,iwe
467  if (flux_type(idims, iw) == flux_tvdlf) then
468  flc(ixc^s,iw)=-tvdlfeps*half*max(cmaxc(ixc^s,ii),abs(cminc(ixc^s,ii))) * &
469  (wrc(ixc^s,iw) - wlc(ixc^s,iw))
470  else
471  where(patchf(ixc^s)==-2)
472  flc(ixc^s,iw)=flc(ixc^s,iw)
473  elsewhere(abs(patchf(ixc^s))==1)
474  flc(ixc^s,iw)=fcd(ixc^s,iw)
475  elsewhere(patchf(ixc^s)==2)
476  flc(ixc^s,iw)=frc(ixc^s,iw)
477  elsewhere(patchf(ixc^s)==3)
478  ! fallback option, reducing to HLL flux
479  flc(ixc^s,iw)=fhll(ixc^s,iw)
480  elsewhere(patchf(ixc^s)==4)
481  ! fallback option, reducing to TVDLF flux
482  flc(ixc^s,iw) = half*((flc(ixc^s,iw)+frc(ixc^s,iw)) &
483  -tvdlfeps * max(cmaxc(ixc^s,ii), dabs(cminc(ixc^s,ii))) * &
484  (wrc(ixc^s,iw)-wlc(ixc^s,iw)))
485  endwhere
486  end if
487 
488  fc(ixc^s,iw,idims)=flc(ixc^s,iw)
489 
490  end do ! Next iw
491  end subroutine get_riemann_flux_hllc
492 
493  !> HLLD Riemann flux from Miyoshi 2005 JCP, 208, 315 and Guo 2016 JCP, 327, 543
494  subroutine get_riemann_flux_hlld(iws,iwe)
495  integer, intent(in) :: iws, iwe
496  double precision, dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
497  double precision, dimension(ixI^S,1:nwflux) :: w2R,w2L
498  double precision, dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
499  double precision, dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
500  ! velocity from the right and the left reconstruction
501  double precision, dimension(ixI^S,ndir) :: vRC, vLC
502  ! magnetic field from the right and the left reconstruction
503  double precision, dimension(ixI^S,ndir) :: BR, BL
504  integer :: ip1,ip2,ip3,idir,ix^D
505  integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
506 
507  associate(sr=>cmaxc,sl=>cminc)
508 
509  rho_ = iw_rho
510  mom(:) = iw_mom(:)
511  mag(:) = iw_mag(:)
512  e_ = iw_e
513  p_ = e_
514 
515  !f1R=0.d0
516  !f1L=0.d0
517  !f2R=0.d0
518  !f2L=0.d0
519  !w1L=0.d0
520  !w1R=0.d0
521  !w2L=0.d0
522  !w2R=0.d0
523  ip1=idims
524  ip3=3
525  vrc(ixc^s,:)=wrp(ixc^s,mom(:))
526  vlc(ixc^s,:)=wlp(ixc^s,mom(:))
527  if(b0field) then
528  br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
529  bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
530  else
531  br(ixc^s,:)=wrc(ixc^s,mag(:))
532  bl(ixc^s,:)=wlc(ixc^s,mag(:))
533  end if
534  if(stagger_grid) then
535  bx(ixc^s)=block%ws(ixc^s,ip1)
536  else
537  ! HLL estimation of normal magnetic field at cell interfaces
538  ! Li, Shenghai, 2005 JCP, 203, 344, equation (33)
539  bx(ixc^s)=(sr(ixc^s,ii)*br(ixc^s,ip1)-sl(ixc^s,ii)*bl(ixc^s,ip1))/(sr(ixc^s,ii)-sl(ixc^s,ii))
540  end if
541  ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
542  ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
543  if(iw_equi_rho>0) then
544  sur(ixc^s) = wrc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
545  else
546  sur(ixc^s) = wrc(ixc^s,rho_)
547  endif
548  sur(ixc^s)=(sr(ixc^s,ii)-vrc(ixc^s,ip1))*sur(ixc^s)
549  if(iw_equi_rho>0) then
550  sul(ixc^s) = wlc(ixc^s,rho_)+ block%equi_vars(ixc^s,iw_equi_rho,ip1)
551  else
552  sul(ixc^s) = wlc(ixc^s,rho_)
553  endif
554  sul(ixc^s)=(sl(ixc^s,ii)-vlc(ixc^s,ip1))*sul(ixc^s)
555  ! Miyoshi equation (38) and Guo euqation (20)
556  sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
557  ptr(ixc^s)+ptl(ixc^s))/(sur(ixc^s)-sul(ixc^s))
558  ! Miyoshi equation (39) and Guo euqation (28)
559  w1r(ixc^s,mom(ip1))=sm(ixc^s)
560  w1l(ixc^s,mom(ip1))=sm(ixc^s)
561  w2r(ixc^s,mom(ip1))=sm(ixc^s)
562  w2l(ixc^s,mom(ip1))=sm(ixc^s)
563  ! Guo equation (22)
564  w1r(ixc^s,mag(ip1))=bx(ixc^s)
565  w1l(ixc^s,mag(ip1))=bx(ixc^s)
566  if(b0field) then
567  ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
568  ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
569  end if
570 
571  ! Miyoshi equation (43) and Guo equation (27)
572  w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,ii)-sm(ixc^s))
573  w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,ii)-sm(ixc^s))
574 
575  ip2=mod(ip1+1,ndir)
576  if(ip2==0) ip2=ndir
577  r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
578  where(abs(r1r(ixc^s))>smalldouble)
579  r1r(ixc^s)=1.d0/r1r(ixc^s)
580  else where
581  r1r(ixc^s)=0.d0
582  end where
583  r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,ii)-sm(ixc^s))-bx(ixc^s)**2
584  where(abs(r1l(ixc^s))>smalldouble)
585  r1l(ixc^s)=1.d0/r1l(ixc^s)
586  else where
587  r1l(ixc^s)=0.d0
588  end where
589  ! Miyoshi equation (44)
590  w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
591  (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
592  w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
593  (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
594  ! partial solution for later usage
595  w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,ii)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
596  w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,ii)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
597  if(ndir==3) then
598  ip3=mod(ip1+2,ndir)
599  if(ip3==0) ip3=ndir
600  ! Miyoshi equation (46)
601  w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
602  (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
603  w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
604  (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
605  ! Miyoshi equation (47)
606  w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
607  w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
608  end if
609  ! Miyoshi equation (45)
610  w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
611  w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
612  if(b0field) then
613  ! Guo equation (26)
614  w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
615  w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
616  end if
617  ! equation (48)
618  if(phys_energy) then
619  ! Guo equation (25) equivalent to Miyoshi equation (41)
620  w1r(ixc^s,p_)=sur(ixc^s)*(sm(ixc^s)-vrc(ixc^s,ip1))+ptr(ixc^s)
621  !w1L(ixC^S,p_)=suL(ixC^S)*(sm(ixC^S)-vLC(ixC^S,ip1))+ptL(ixC^S)
622  w1l(ixc^s,p_)=w1r(ixc^s,p_)
623  if(b0field) then
624  ! Guo equation (32)
625  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)
626  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)
627  end if
628  ! Miyoshi equation (48) and main part of Guo euqation (31)
629  w1r(ixc^s,e_)=((sr(ixc^s,ii)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
630  w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
631  sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,ii)-sm(ixc^s))
632  w1l(ixc^s,e_)=((sl(ixc^s,ii)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
633  w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
634  sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,ii)-sm(ixc^s))
635  if(b0field) then
636  ! Guo equation (31)
637  w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
638  sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
639  w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
640  sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
641  end if
642  if(iw_equi_p>0) then
643 #if !defined(E_RM_W0) || E_RM_W0 == 1
644  w1r(ixc^s,e_)= w1r(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
645  (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
646  w1l(ixc^s,e_)= w1l(ixc^s,e_) + 1d0/(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
647  (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
648 #else
649  w1r(ixc^s,e_)= w1r(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
650  (sm(ixc^s)-vrc(ixc^s,ip1))/(sr(ixc^s,ii)-sm(ixc^s))
651  w1l(ixc^s,e_)= w1l(ixc^s,e_) + phys_gamma /(phys_gamma - 1) * block%equi_vars(ixc^s,iw_equi_p,ip1) * &
652  (sm(ixc^s)-vlc(ixc^s,ip1))/(sl(ixc^s,ii)-sm(ixc^s))
653 #endif
654  endif
655  end if
656 
657  ! Miyoshi equation (49) and Guo equation (35)
658  w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
659  w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
660  w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
661  w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
662 
663  r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
664  r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
665  tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
666  signbx(ixc^s)=sign(1.d0,bx(ixc^s))
667  ! Miyoshi equation (51) and Guo equation (33)
668  s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
669  s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
670  ! Miyoshi equation (59) and Guo equation (41)
671  w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
672  (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
673  w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
674  ! Miyoshi equation (61) and Guo equation (43)
675  w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
676  r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
677  w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
678  if(ndir==3) then
679  ! Miyoshi equation (60) and Guo equation (42)
680  w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
681  (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
682  w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
683  ! Miyoshi equation (62) and Guo equation (44)
684  w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
685  r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
686  w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
687  end if
688  ! Miyoshi equation (63) and Guo equation (45)
689  if(phys_energy) then
690  w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
691  sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
692  w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
693  sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
694  end if
695 
696  ! convert velocity to momentum
697  do idir=1,ndir
698  w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
699  w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
700  w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
701  w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
702  end do
703  if(iw_equi_rho>0) then
704  w1r(ixc^s,rho_) = w1r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
705  w1l(ixc^s,rho_) = w1l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
706  w2r(ixc^s,rho_) = w2r(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
707  w2l(ixc^s,rho_) = w2l(ixc^s,rho_) - block%equi_vars(ixc^s,iw_equi_rho,ip1)
708  endif
709  ! get fluxes of intermedate states
710  do iw=iws,iwe
711  if(flux_type(idims, iw) == flux_special) then
712  ! known flux (fLC=fRC) for normal B and psi_ in GLM method
713  f1l(ixc^s,iw)=flc(ixc^s,iw)
714  f1r(ixc^s,iw)=f1l(ixc^s,iw)
715  f2l(ixc^s,iw)=f1l(ixc^s,iw)
716  f2r(ixc^s,iw)=f1l(ixc^s,iw)
717  else if(flux_type(idims, iw) == flux_hll) then
718  ! using hll flux for tracers
719  f1l(ixc^s,iw)=(sr(ixc^s,ii)*flc(ixc^s, iw)-sl(ixc^s,ii)*frc(ixc^s, iw) &
720  +sr(ixc^s,ii)*sl(ixc^s,ii)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,ii)-sl(ixc^s,ii))
721  f1r(ixc^s,iw)=f1l(ixc^s,iw)
722  f2l(ixc^s,iw)=f1l(ixc^s,iw)
723  f2r(ixc^s,iw)=f1l(ixc^s,iw)
724  else
725  ! construct hlld flux
726  f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,ii)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
727  f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,ii)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
728  f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
729  f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
730  end if
731  end do
732 
733  ! Miyoshi equation (66) and Guo equation (46)
734  {do ix^db=ixcmin^db,ixcmax^db\}
735  if(sl(ix^d,ii)>0.d0) then
736  fc(ix^d,iws:iwe,ip1)=flc(ix^d,iws:iwe)
737  else if(s1l(ix^d)>=0.d0) then
738  fc(ix^d,iws:iwe,ip1)=f1l(ix^d,iws:iwe)
739  else if(sm(ix^d)>=0.d0) then
740  fc(ix^d,iws:iwe,ip1)=f2l(ix^d,iws:iwe)
741  else if(s1r(ix^d)>=0.d0) then
742  fc(ix^d,iws:iwe,ip1)=f2r(ix^d,iws:iwe)
743  else if(sr(ix^d,ii)>=0.d0) then
744  fc(ix^d,iws:iwe,ip1)=f1r(ix^d,iws:iwe)
745  else if(sr(ix^d,ii)<0.d0) then
746  fc(ix^d,iws:iwe,ip1)=frc(ix^d,iws:iwe)
747  end if
748  {end do\}
749 
750  end associate
751  end subroutine get_riemann_flux_hlld
752 
753  !> HLLD Riemann flux from Miyoshi 2005 JCP, 208, 315 and Guo 2016 JCP, 327, 543
754  !> https://arxiv.org/pdf/2108.04991.pdf
755  subroutine get_riemann_flux_hlld_mag2()
756  !use mod_mhd_phys
757  use mod_variables
758  use mod_physics
759  implicit none
760  double precision, dimension(ixI^S,1:nwflux) :: w1R,w1L,f1R,f1L,f2R,f2L
761  double precision, dimension(ixI^S,1:nwflux) :: w2R,w2L
762  double precision, dimension(ixI^S) :: sm,s1R,s1L,suR,suL,Bx
763  double precision, dimension(ixI^S) :: pts,ptR,ptL,signBx,r1L,r1R,tmp
764  ! velocity from the right and the left reconstruction
765  double precision, dimension(ixI^S,ndir) :: vRC, vLC
766  ! magnetic field from the right and the left reconstruction
767  double precision, dimension(ixI^S,ndir) :: BR, BL
768  integer :: ip1,ip2,ip3,idir,ix^D
769  double precision :: phiPres, thetaSM, du, dv, dw
770  integer :: ixV^L, ixVb^L, ixVc^L, ixVd^L, ixVe^L, ixVf^L
771  integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
772  double precision, parameter :: aParam = 4d0
773 
774  rho_ = iw_rho
775  mom(:) = iw_mom(:)
776  mag(:) = iw_mag(:)
777  p_ = iw_e
778  e_ = iw_e
779 
780  associate(sr=>cmaxc,sl=>cminc)
781 
782  f1r=0.d0
783  f1l=0.d0
784  ip1=idims
785  ip3=3
786 
787  vrc(ixc^s,:)=wrp(ixc^s,mom(:))
788  vlc(ixc^s,:)=wlp(ixc^s,mom(:))
789 
790  ! reuse s1L s1R
791  call get_hlld2_modif_c(wlp,x,ixi^l,ixo^l,s1l)
792  call get_hlld2_modif_c(wrp,x,ixi^l,ixo^l,s1r)
793  !phiPres = min(1, maxval(max(s1L(ixO^S),s1R(ixO^S))/cmaxC(ixO^S,1)))
794  phipres = min(1d0, maxval(max(s1l(ixo^s),s1r(ixo^s)))/maxval(cmaxc(ixo^s,1)))
795  phipres = phipres*(2d0 - phipres)
796 
797  !we use here not reconstructed velocity: wprim?
798  ixv^l=ixo^l;
799  !first dim
800  ixvmin1=ixomin1+1
801  ixvmax1=ixomax1+1
802  du = minval(wprim(ixv^s,mom(1))-wprim(ixo^s,mom(1)))
803  if(du>0d0) du=0d0
804  dv = 0d0
805  dw = 0d0
806 
807  {^nooned
808  !second dim
809  !i,j-1,k
810  ixv^l=ixo^l;
811  ixvmin2=ixomin2-1
812  ixvmax2=ixomax2-1
813 
814  !i,j+1,k
815  ixvb^l=ixo^l;
816  ixvbmin2=ixomin2+1
817  ixvbmax2=ixomax2+1
818 
819  !i+1,j,k
820  ixvc^l=ixo^l;
821  ixvcmin1=ixomin1+1
822  ixvcmax1=ixomax1+1
823 
824  !i+1,j-1,k
825  ixvd^l=ixo^l;
826  ixvdmin1=ixomin1+1
827  ixvdmax1=ixomax1+1
828  ixvdmin2=ixomin2-1
829  ixvdmax2=ixomax2-1
830 
831  !i+1,j+1,k
832  ixve^l=ixo^l;
833  ixvemin1=ixomin1+1
834  ixvemax1=ixomax1+1
835  ixvemin2=ixomin2+1
836  ixvemax2=ixomax2+1
837 
838  dv = minval(min(wprim(ixo^s,mom(2))-wprim(ixv^s,mom(2)),&
839  wprim(ixvb^s,mom(2))-wprim(ixo^s,mom(2)),&
840  wprim(ixvc^s,mom(2))-wprim(ixvd^s,mom(2)),&
841  wprim(ixve^s,mom(2))-wprim(ixvc^s,mom(2))&
842  ))
843  if(dv>0d0) dv=0d0}
844 
845  {^ifthreed
846  !third dim
847  !i,j,k-1
848  ixv^l=ixo^l;
849  ixvmin3=ixomin3-1
850  ixvmax3=ixomax3-1
851 
852  !i,j,k+1
853  ixvb^l=ixo^l;
854  ixvbmin3=ixomin3+1
855  ixvbmax3=ixomax3+1
856 
857  !i+1,j,k
858  ixvc^l=ixo^l;
859  ixvcmin1=ixomin1+1
860  ixvcmax1=ixomax1+1
861 
862  !i+1,j,k-1
863  ixvd^l=ixo^l;
864  ixvdmin1=ixomin1+1
865  ixvdmax1=ixomax1+1
866  ixvdmin3=ixomin3-1
867  ixvdmax3=ixomax3-1
868 
869  !i+1,j,k+1
870  ixve^l=ixo^l;
871  ixvemin1=ixomin1+1
872  ixvemax1=ixomax1+1
873  ixvemin3=ixomin3+1
874  ixvemax3=ixomax3+1
875  dw = minval(min(wprim(ixo^s,mom(3))-wprim(ixv^s,mom(3)),&
876  wprim(ixvb^s,mom(3))-wprim(ixo^s,mom(3)),&
877  wprim(ixvc^s,mom(3))-wprim(ixvd^s,mom(3)),&
878  wprim(ixve^s,mom(3))-wprim(ixvc^s,mom(3))&
879  ))
880  if(dw>0d0) dw=0d0}
881  thetasm = maxval(cmaxc(ixo^s,1))
882 
883  thetasm = (min(1d0, (thetasm-du)/(thetasm-min(dv,dw))))**aparam
884  !print*, "HLLD2 ", du,dv,dw, thetaSM, phiPres
885 
886  if(b0field) then
887  br(ixc^s,:)=wrc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
888  bl(ixc^s,:)=wlc(ixc^s,mag(:))+block%B0(ixc^s,:,ip1)
889  else
890  br(ixc^s,:)=wrc(ixc^s,mag(:))
891  bl(ixc^s,:)=wlc(ixc^s,mag(:))
892  end if
893  ! HLL estimation of normal magnetic field at cell interfaces
894  bx(ixc^s)=(sr(ixc^s,index_v_mag)*br(ixc^s,ip1)-sl(ixc^s,index_v_mag)*bl(ixc^s,ip1)-&
895  flc(ixc^s,mag(ip1))-frc(ixc^s,mag(ip1)))/(sr(ixc^s,index_v_mag)-sl(ixc^s,index_v_mag))
896  ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(br(ixc^s,:)**2,dim=ndim+1)
897  ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(bl(ixc^s,:)**2,dim=ndim+1)
898  sur(ixc^s)=(sr(ixc^s,index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,rho_)
899  sul(ixc^s)=(sl(ixc^s,index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,rho_)
900  ! Miyoshi equation (38) and Guo euqation (20)
901  sm(ixc^s)=(sur(ixc^s)*vrc(ixc^s,ip1)-sul(ixc^s)*vlc(ixc^s,ip1)-&
902  thetasm*(ptr(ixc^s)-ptl(ixc^s)) )/(sur(ixc^s)-sul(ixc^s))
903  ! Miyoshi equation (39) and Guo euqation (28)
904  w1r(ixc^s,mom(ip1))=sm(ixc^s)
905  w1l(ixc^s,mom(ip1))=sm(ixc^s)
906  w2r(ixc^s,mom(ip1))=sm(ixc^s)
907  w2l(ixc^s,mom(ip1))=sm(ixc^s)
908  ! Guo equation (22)
909  w1r(ixc^s,mag(ip1))=bx(ixc^s)
910  w1l(ixc^s,mag(ip1))=bx(ixc^s)
911  if(b0field) then
912  ptr(ixc^s)=wrp(ixc^s,p_)+0.5d0*sum(wrc(ixc^s,mag(:))**2,dim=ndim+1)
913  ptl(ixc^s)=wlp(ixc^s,p_)+0.5d0*sum(wlc(ixc^s,mag(:))**2,dim=ndim+1)
914  end if
915 
916  ! Miyoshi equation (43) and Guo equation (27)
917  w1r(ixc^s,rho_)=sur(ixc^s)/(sr(ixc^s,index_v_mag)-sm(ixc^s))
918  w1l(ixc^s,rho_)=sul(ixc^s)/(sl(ixc^s,index_v_mag)-sm(ixc^s))
919 
920  ip2=mod(ip1+1,ndir)
921  if(ip2==0) ip2=ndir
922  r1r(ixc^s)=sur(ixc^s)*(sr(ixc^s,index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
923  where(r1r(ixc^s)/=0.d0)
924  r1r(ixc^s)=1.d0/r1r(ixc^s)
925  endwhere
926  r1l(ixc^s)=sul(ixc^s)*(sl(ixc^s,index_v_mag)-sm(ixc^s))-bx(ixc^s)**2
927  where(r1l(ixc^s)/=0.d0)
928  r1l(ixc^s)=1.d0/r1l(ixc^s)
929  endwhere
930  ! Miyoshi equation (44)
931  w1r(ixc^s,mom(ip2))=vrc(ixc^s,ip2)-bx(ixc^s)*br(ixc^s,ip2)*&
932  (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
933  w1l(ixc^s,mom(ip2))=vlc(ixc^s,ip2)-bx(ixc^s)*bl(ixc^s,ip2)*&
934  (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
935  ! partial solution for later usage
936  w1r(ixc^s,mag(ip2))=(sur(ixc^s)*(sr(ixc^s,index_v_mag)-vrc(ixc^s,ip1))-bx(ixc^s)**2)*r1r(ixc^s)
937  w1l(ixc^s,mag(ip2))=(sul(ixc^s)*(sl(ixc^s,index_v_mag)-vlc(ixc^s,ip1))-bx(ixc^s)**2)*r1l(ixc^s)
938  if(ndir==3) then
939  ip3=mod(ip1+2,ndir)
940  if(ip3==0) ip3=ndir
941  ! Miyoshi equation (46)
942  w1r(ixc^s,mom(ip3))=vrc(ixc^s,ip3)-bx(ixc^s)*br(ixc^s,ip3)*&
943  (sm(ixc^s)-vrc(ixc^s,ip1))*r1r(ixc^s)
944  w1l(ixc^s,mom(ip3))=vlc(ixc^s,ip3)-bx(ixc^s)*bl(ixc^s,ip3)*&
945  (sm(ixc^s)-vlc(ixc^s,ip1))*r1l(ixc^s)
946  ! Miyoshi equation (47)
947  w1r(ixc^s,mag(ip3))=br(ixc^s,ip3)*w1r(ixc^s,mag(ip2))
948  w1l(ixc^s,mag(ip3))=bl(ixc^s,ip3)*w1l(ixc^s,mag(ip2))
949  end if
950  ! Miyoshi equation (45)
951  w1r(ixc^s,mag(ip2))=br(ixc^s,ip2)*w1r(ixc^s,mag(ip2))
952  w1l(ixc^s,mag(ip2))=bl(ixc^s,ip2)*w1l(ixc^s,mag(ip2))
953  if(b0field) then
954  ! Guo equation (26)
955  w1r(ixc^s,mag(:))=w1r(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
956  w1l(ixc^s,mag(:))=w1l(ixc^s,mag(:))-block%B0(ixc^s,:,ip1)
957  end if
958  ! equation (48)
959  if(phys_energy) then
960  ! Guo equation (25) equivalent to Miyoshi equation (41)
961  w1r(ixc^s,p_)=(sur(ixc^s)*ptl(ixc^s) - sul(ixc^s)*ptr(ixc^s) +&
962  phipres * sur(ixc^s)*sul(ixc^s)*(vrc(ixc^s,ip1)-vlc(ixc^s,ip1)))/&
963  (sur(ixc^s)-sul(ixc^s))
964  w1l(ixc^s,p_)=w1r(ixc^s,p_)
965  if(b0field) then
966  ! Guo equation (32)
967  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)
968  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)
969  end if
970  ! Miyoshi equation (48) and main part of Guo euqation (31)
971  w1r(ixc^s,e_)=((sr(ixc^s,index_v_mag)-vrc(ixc^s,ip1))*wrc(ixc^s,e_)-ptr(ixc^s)*vrc(ixc^s,ip1)+&
972  w1r(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vrc(ixc^s,:)*wrc(ixc^s,mag(:)),dim=ndim+1)-&
973  sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)))/(sr(ixc^s,index_v_mag)-sm(ixc^s))
974  w1l(ixc^s,e_)=((sl(ixc^s,index_v_mag)-vlc(ixc^s,ip1))*wlc(ixc^s,e_)-ptl(ixc^s)*vlc(ixc^s,ip1)+&
975  w1l(ixc^s,p_)*sm(ixc^s)+bx(ixc^s)*(sum(vlc(ixc^s,:)*wlc(ixc^s,mag(:)),dim=ndim+1)-&
976  sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)))/(sl(ixc^s,index_v_mag)-sm(ixc^s))
977  if(b0field) then
978  ! Guo equation (31)
979  w1r(ixc^s,e_)=w1r(ixc^s,e_)+(sum(w1r(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
980  sum(wrc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vrc(ixc^s,ip1))/(sr(ixc^s,index_v_mag)-sm(ixc^s))
981  w1l(ixc^s,e_)=w1l(ixc^s,e_)+(sum(w1l(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*sm(ixc^s)-&
982  sum(wlc(ixc^s,mag(:))*block%B0(ixc^s,:,ip1),dim=ndim+1)*vlc(ixc^s,ip1))/(sl(ixc^s,index_v_mag)-sm(ixc^s))
983  end if
984  end if
985 
986  ! Miyoshi equation (49) and Guo equation (35)
987  w2r(ixc^s,rho_)=w1r(ixc^s,rho_)
988  w2l(ixc^s,rho_)=w1l(ixc^s,rho_)
989  w2r(ixc^s,mag(ip1))=w1r(ixc^s,mag(ip1))
990  w2l(ixc^s,mag(ip1))=w1l(ixc^s,mag(ip1))
991 
992  r1r(ixc^s)=sqrt(w1r(ixc^s,rho_))
993  r1l(ixc^s)=sqrt(w1l(ixc^s,rho_))
994  tmp(ixc^s)=1.d0/(r1r(ixc^s)+r1l(ixc^s))
995  signbx(ixc^s)=sign(1.d0,bx(ixc^s))
996  ! Miyoshi equation (51) and Guo equation (33)
997  s1r(ixc^s)=sm(ixc^s)+abs(bx(ixc^s))/r1r(ixc^s)
998  s1l(ixc^s)=sm(ixc^s)-abs(bx(ixc^s))/r1l(ixc^s)
999  ! Miyoshi equation (59) and Guo equation (41)
1000  w2r(ixc^s,mom(ip2))=(r1l(ixc^s)*w1l(ixc^s,mom(ip2))+r1r(ixc^s)*w1r(ixc^s,mom(ip2))+&
1001  (w1r(ixc^s,mag(ip2))-w1l(ixc^s,mag(ip2)))*signbx(ixc^s))*tmp(ixc^s)
1002  w2l(ixc^s,mom(ip2))=w2r(ixc^s,mom(ip2))
1003  ! Miyoshi equation (61) and Guo equation (43)
1004  w2r(ixc^s,mag(ip2))=(r1l(ixc^s)*w1r(ixc^s,mag(ip2))+r1r(ixc^s)*w1l(ixc^s,mag(ip2))+&
1005  r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip2))-w1l(ixc^s,mom(ip2)))*signbx(ixc^s))*tmp(ixc^s)
1006  w2l(ixc^s,mag(ip2))=w2r(ixc^s,mag(ip2))
1007  if(ndir==3) then
1008  ! Miyoshi equation (60) and Guo equation (42)
1009  w2r(ixc^s,mom(ip3))=(r1l(ixc^s)*w1l(ixc^s,mom(ip3))+r1r(ixc^s)*w1r(ixc^s,mom(ip3))+&
1010  (w1r(ixc^s,mag(ip3))-w1l(ixc^s,mag(ip3)))*signbx(ixc^s))*tmp(ixc^s)
1011  w2l(ixc^s,mom(ip3))=w2r(ixc^s,mom(ip3))
1012  ! Miyoshi equation (62) and Guo equation (44)
1013  w2r(ixc^s,mag(ip3))=(r1l(ixc^s)*w1r(ixc^s,mag(ip3))+r1r(ixc^s)*w1l(ixc^s,mag(ip3))+&
1014  r1l(ixc^s)*r1r(ixc^s)*(w1r(ixc^s,mom(ip3))-w1l(ixc^s,mom(ip3)))*signbx(ixc^s))*tmp(ixc^s)
1015  w2l(ixc^s,mag(ip3))=w2r(ixc^s,mag(ip3))
1016  end if
1017  ! Miyoshi equation (63) and Guo equation (45)
1018  if(phys_energy) then
1019  w2r(ixc^s,e_)=w1r(ixc^s,e_)+r1r(ixc^s)*(sum(w1r(ixc^s,mom(:))*w1r(ixc^s,mag(:)),dim=ndim+1)-&
1020  sum(w2r(ixc^s,mom(:))*w2r(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
1021  w2l(ixc^s,e_)=w1l(ixc^s,e_)-r1l(ixc^s)*(sum(w1l(ixc^s,mom(:))*w1l(ixc^s,mag(:)),dim=ndim+1)-&
1022  sum(w2l(ixc^s,mom(:))*w2l(ixc^s,mag(:)),dim=ndim+1))*signbx(ixc^s)
1023  end if
1024 
1025  ! convert velocity to momentum
1026  do idir=1,ndir
1027  w1r(ixc^s,mom(idir))=w1r(ixc^s,mom(idir))*w1r(ixc^s,rho_)
1028  w1l(ixc^s,mom(idir))=w1l(ixc^s,mom(idir))*w1l(ixc^s,rho_)
1029  w2r(ixc^s,mom(idir))=w2r(ixc^s,mom(idir))*w2r(ixc^s,rho_)
1030  w2l(ixc^s,mom(idir))=w2l(ixc^s,mom(idir))*w2l(ixc^s,rho_)
1031  end do
1032 
1033  ! get fluxes of intermedate states
1034  do iw=1,nwflux
1035  ! CT MHD does not need normal B flux
1036  if(stagger_grid .and. flux_type(idims, iw) == flux_tvdlf) cycle
1037  if(flux_type(idims, iw) == flux_special) then
1038  ! known flux (fLC=fRC) for normal B and psi_ in GLM method
1039  f1l(ixc^s,iw)=flc(ixc^s,iw)
1040  f1r(ixc^s,iw)=f1l(ixc^s,iw)
1041  f2l(ixc^s,iw)=f1l(ixc^s,iw)
1042  f2r(ixc^s,iw)=f1l(ixc^s,iw)
1043  else if(flux_type(idims, iw) == flux_hll) then
1044  ! using hll flux for tracers
1045  f1l(ixc^s,iw)=(sr(ixc^s,index_v_mag)*flc(ixc^s, iw)-sl(ixc^s,index_v_mag)*frc(ixc^s, iw) &
1046  +sr(ixc^s,index_v_mag)*sl(ixc^s,index_v_mag)*(wrc(ixc^s,iw)-wlc(ixc^s,iw)))/(sr(ixc^s,index_v_mag)-sl(ixc^s,index_v_mag))
1047  f1r(ixc^s,iw)=f1l(ixc^s,iw)
1048  f2l(ixc^s,iw)=f1l(ixc^s,iw)
1049  f2r(ixc^s,iw)=f1l(ixc^s,iw)
1050  else
1051  f1l(ixc^s,iw)=flc(ixc^s,iw)+sl(ixc^s,index_v_mag)*(w1l(ixc^s,iw)-wlc(ixc^s,iw))
1052  f1r(ixc^s,iw)=frc(ixc^s,iw)+sr(ixc^s,index_v_mag)*(w1r(ixc^s,iw)-wrc(ixc^s,iw))
1053  f2l(ixc^s,iw)=f1l(ixc^s,iw)+s1l(ixc^s)*(w2l(ixc^s,iw)-w1l(ixc^s,iw))
1054  f2r(ixc^s,iw)=f1r(ixc^s,iw)+s1r(ixc^s)*(w2r(ixc^s,iw)-w1r(ixc^s,iw))
1055  end if
1056  end do
1057 
1058  ! Miyoshi equation (66) and Guo equation (46)
1059  {do ix^db=ixcmin^db,ixcmax^db\}
1060  if(sl(ix^d,index_v_mag)>0.d0) then
1061  fc(ix^d,1:nwflux,ip1)=flc(ix^d,1:nwflux)
1062  else if(s1l(ix^d)>=0.d0) then
1063  fc(ix^d,1:nwflux,ip1)=f1l(ix^d,1:nwflux)
1064  else if(sm(ix^d)>=0.d0) then
1065  fc(ix^d,1:nwflux,ip1)=f2l(ix^d,1:nwflux)
1066  else if(s1r(ix^d)>=0.d0) then
1067  fc(ix^d,1:nwflux,ip1)=f2r(ix^d,1:nwflux)
1068  else if(sr(ix^d,index_v_mag)>=0.d0) then
1069  fc(ix^d,1:nwflux,ip1)=f1r(ix^d,1:nwflux)
1070  else if(sr(ix^d,index_v_mag)<0.d0) then
1071  fc(ix^d,1:nwflux,ip1)=frc(ix^d,1:nwflux)
1072  end if
1073  {end do\}
1074 
1075  end associate
1076  end subroutine get_riemann_flux_hlld_mag2
1077 
1078  !> Calculate fast magnetosonic wave speed
1079  subroutine get_hlld2_modif_c(w,x,ixI^L,ixO^L,csound)
1081 
1082  integer, intent(in) :: ixI^L, ixO^L
1083  double precision, intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
1084  double precision, intent(out):: csound(ixI^S)
1085  double precision :: cfast2(ixI^S), AvMinCs2(ixI^S), b2(ixI^S), kmax
1086  double precision :: inv_rho(ixO^S), gamma_A2(ixO^S)
1087  integer :: rho_, p_, e_, mom(1:ndir), mag(1:ndir)
1088 
1089  rho_ = iw_rho
1090  mom(:) = iw_mom(:)
1091  mag(:) = iw_mag(:)
1092  p_ = iw_e
1093  e_ = iw_e
1094 
1095  inv_rho=1.d0/w(ixo^s,rho_)
1096 
1097  ! store |B|^2 in v
1098 
1099  if (b0field) then
1100  b2(ixo^s) = sum((w(ixo^s, mag(:))+block%B0(ixo^s,:,b0i))**2, dim=ndim+1)
1101  else
1102  b2(ixo^s) = sum(w(ixo^s, mag(:))**2, dim=ndim+1)
1103  end if
1104 
1105 
1106  if (b0field) then
1107  avmincs2= w(ixo^s, mag(idims))+block%B0(ixo^s,idims,b0i)
1108  else
1109  avmincs2= w(ixo^s, mag(idims))
1110  end if
1111 
1112 
1113  csound(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
1114 
1115  cfast2(ixo^s) = b2(ixo^s) * inv_rho+csound(ixo^s)
1116  avmincs2(ixo^s) = cfast2(ixo^s)**2-4.0d0*csound(ixo^s) &
1117  * avmincs2(ixo^s)**2 &
1118  * inv_rho
1119 
1120  where(avmincs2(ixo^s)<zero)
1121  avmincs2(ixo^s)=zero
1122  end where
1123 
1124  avmincs2(ixo^s)=sqrt(avmincs2(ixo^s))
1125 
1126  csound(ixo^s) = sqrt(half*(cfast2(ixo^s)+avmincs2(ixo^s)))
1127 
1128  end subroutine get_hlld2_modif_c
1129 
1130  end subroutine finite_volume
1131 
1132  !> Determine the upwinded wLC(ixL) and wRC(ixR) from w.
1133  !> the wCT is only used when PPM is exploited.
1134  subroutine reconstruct_lr(ixI^L,ixL^L,ixR^L,idims,w,wLC,wRC,wLp,wRp,x,dxdim)
1135  use mod_physics
1137  use mod_limiter
1138  use mod_comm_lib, only: mpistop
1139 
1140  integer, intent(in) :: ixi^l, ixl^l, ixr^l, idims
1141  double precision, intent(in) :: dxdim
1142  ! cell center w in primitive form
1143  double precision, dimension(ixI^S,1:nw) :: w
1144  ! left and right constructed status in conservative form
1145  double precision, dimension(ixI^S,1:nw) :: wlc, wrc
1146  ! left and right constructed status in primitive form
1147  double precision, dimension(ixI^S,1:nw) :: wlp, wrp
1148  double precision, dimension(ixI^S,1:ndim) :: x
1149 
1150  integer :: jxr^l, ixc^l, jxc^l, iw
1151  double precision :: ldw(ixi^s), rdw(ixi^s), dwc(ixi^s)
1152  double precision :: a2max
1153 
1154  select case (type_limiter(block%level))
1155  case (limiter_mp5)
1156  call mp5limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
1157  case (limiter_weno3)
1158  call weno3limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
1159  case (limiter_wenoyc3)
1160  call weno3limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
1161  case (limiter_weno5)
1162  call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
1163  case (limiter_weno5nm)
1164  call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,1)
1165  case (limiter_wenoz5)
1166  call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
1167  case (limiter_wenoz5nm)
1168  call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,2)
1169  case (limiter_wenozp5)
1170  call weno5limiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,3)
1171  case (limiter_wenozp5nm)
1172  call weno5nmlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp,3)
1173  case (limiter_weno5cu6)
1174  call weno5cu6limiter(ixi^l,ixl^l,idims,w,wlp,wrp)
1175  case (limiter_teno5ad)
1176  call teno5adlimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp)
1177  case (limiter_weno7)
1178  call weno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp,1)
1179  case (limiter_mpweno7)
1180  call weno7limiter(ixi^l,ixl^l,idims,w,wlp,wrp,2)
1181  case (limiter_venk)
1182  call venklimiter(ixi^l,ixl^l,idims,dxdim,w,wlp,wrp)
1183  if(fix_small_values) then
1184  call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
1185  call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
1186  end if
1187  case (limiter_ppm)
1188  call ppmlimiter(ixi^l,ixm^ll,idims,w,w,wlp,wrp)
1189  if(fix_small_values) then
1190  call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
1191  call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
1192  end if
1193  case default
1194  jxr^l=ixr^l+kr(idims,^d);
1195  ixcmax^d=jxrmax^d; ixcmin^d=ixlmin^d-kr(idims,^d);
1196  jxc^l=ixc^l+kr(idims,^d);
1197  do iw=1,nwflux
1198  if (loglimit(iw)) then
1199  w(ixcmin^d:jxcmax^d,iw)=dlog10(w(ixcmin^d:jxcmax^d,iw))
1200  wlp(ixl^s,iw)=dlog10(wlp(ixl^s,iw))
1201  wrp(ixr^s,iw)=dlog10(wrp(ixr^s,iw))
1202  end if
1203 
1204  dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1205  if(need_global_a2max) then
1206  a2max=a2max_global(idims)
1207  else
1208  select case(idims)
1209  case(1)
1210  a2max=schmid_rad1
1211  {^iftwod
1212  case(2)
1213  a2max=schmid_rad2}
1214  {^ifthreed
1215  case(2)
1216  a2max=schmid_rad2
1217  case(3)
1218  a2max=schmid_rad3}
1219  case default
1220  call mpistop("idims is wrong in mod_limiter")
1221  end select
1222  end if
1223 
1224  ! limit flux from left and/or right
1225  call dwlimiter2(dwc,ixi^l,ixc^l,idims,type_limiter(block%level),ldw,rdw,a2max=a2max)
1226  wlp(ixl^s,iw)=wlp(ixl^s,iw)+half*ldw(ixl^s)
1227  wrp(ixr^s,iw)=wrp(ixr^s,iw)-half*rdw(jxr^s)
1228 
1229  if (loglimit(iw)) then
1230  w(ixcmin^d:jxcmax^d,iw)=10.0d0**w(ixcmin^d:jxcmax^d,iw)
1231  wlp(ixl^s,iw)=10.0d0**wlp(ixl^s,iw)
1232  wrp(ixr^s,iw)=10.0d0**wrp(ixr^s,iw)
1233  end if
1234  end do
1235  if(fix_small_values) then
1236  call phys_handle_small_values(.true.,wlp,x,ixi^l,ixl^l,'reconstruct left')
1237  call phys_handle_small_values(.true.,wrp,x,ixi^l,ixr^l,'reconstruct right')
1238  end if
1239  end select
1240 
1241  wlc(ixl^s,1:nwflux) = wlp(ixl^s,1:nwflux)
1242  wrc(ixr^s,1:nwflux) = wrp(ixr^s,1:nwflux)
1243  call phys_to_conserved(ixi^l,ixl^l,wlc,x)
1244  call phys_to_conserved(ixi^l,ixr^l,wrc,x)
1245  if(nwaux>0)then
1246  wlp(ixl^s,nwflux+1:nwflux+nwaux) = wlc(ixl^s,nwflux+1:nwflux+nwaux)
1247  wrp(ixr^s,nwflux+1:nwflux+nwaux) = wrc(ixr^s,nwflux+1:nwflux+nwaux)
1248  endif
1249 
1250  end subroutine reconstruct_lr
1251 
1252 end module mod_finite_volume
subroutine get_hlld2_modif_c(w, x, ixIL, ixOL, csound)
Calculate fast magnetosonic wave speed.
subroutine get_riemann_flux_tvdmu()
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module with finite volume methods for fluxes.
subroutine, public finite_volume(method, qdt, dtfactor, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, fC, fE, dxs, x)
finite volume method
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.
subroutine, public hancock(qdt, dtfactor, ixIL, ixOL, idimsLIM, qtC, sCT, qt, snew, dxs, x)
The non-conservative Hancock predictor for TVDLF.
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.
logical h_correction
If true, do H-correction to fix the carbuncle problem at grid-aligned shocks.
integer, parameter fs_tvdlf
integer ixghi
Upper index of grid block arrays.
integer, parameter fs_hlld
integer, dimension(3, 3) kr
Kronecker delta tensor.
logical, dimension(:), allocatable loglimit
integer, parameter fs_hllcd
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
integer, parameter fs_tvdmu
integer b0i
background magnetic field location indicator
integer, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
logical need_global_a2max
global value for schmid scheme
integer ixm
the mesh range of a physical block without ghost cells
logical slab
Cartesian geometry or not.
integer, parameter fs_hll
flux schemes
logical b0field
split magnetic field as background B0 field
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
integer nghostcells
Number of ghost cells surrounding a grid.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, parameter fs_hllc
double precision, dimension(ndim) a2max_global
global largest a2 for schmid scheme
integer, parameter ixglo
Lower index of grid block arrays (always 1)
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_get_ct_velocity), pointer phys_get_ct_velocity
Definition: mod_physics.t:83
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:59
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
Definition: mod_physics.t:48
procedure(sub_small_values), pointer phys_handle_small_values
Definition: mod_physics.t:82
procedure(sub_get_flux), pointer phys_get_flux
Definition: mod_physics.t:68
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition: mod_physics.t:67
integer, parameter flux_hll
Indicates the flux should be treated with hll.
Definition: mod_physics.t:54
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition: mod_physics.t:72
integer, parameter flux_special
Indicates the flux should be specially treated.
Definition: mod_physics.t:52
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
Definition: mod_physics.t:43
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:20
procedure(sub_update_faces), pointer phys_update_faces
Definition: mod_physics.t:84
procedure(sub_convert), pointer phys_to_conserved
Definition: mod_physics.t:58
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition: mod_physics.t:60
procedure(sub_get_h_speed), pointer phys_get_h_speed
Definition: mod_physics.t:66
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:27
Module for handling split source terms (split from the fluxes)
Definition: mod_source.t:2
subroutine, public addsource2(qdt, dtfactor, ixIL, ixOL, iwLIM, qtC, wCT, wCTprim, qt, w, x, qsourcesplit, src_active)
Add source within ixO for iws: w=w+qdt*S[wCT].
Definition: mod_source.t:132
Subroutines for TVD-MUSCL schemes.
Definition: mod_tvd.t:2
subroutine, public tvdlimit2(method, qdt, ixIL, ixICL, ixOL, idims, wL, wR, wnew, x, fC, dxs)
Definition: mod_tvd.t:39
Module with all the methods that users can customize in AMRVAC.
integer nw
Total number of variables.
Definition: mod_variables.t:23
integer iwstart
Definition: mod_variables.t:38
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
Definition: mod_variables.t:75
integer, dimension(:), allocatable iw_mom
Indices of the momentum density.
Definition: mod_variables.t:55
integer, dimension(:), allocatable start_indices
the indices in 1:nwflux array are assumed consecutive for each species this array should be of size n...
Definition: mod_variables.t:84
integer, dimension(:), allocatable stop_indices
the indices in 1:nwflux array are assumed consecutive for each species this array should be of size n...
Definition: mod_variables.t:88
integer, dimension(:), allocatable, protected iw_mag
Indices of the magnetic field components.
Definition: mod_variables.t:67
integer iw_rho
Index of the (gas) density.
Definition: mod_variables.t:52
integer nwflux
Number of flux variables.
Definition: mod_variables.t:8
integer index_v_mag
index of the var whose velocity appears in the induction eq.
Definition: mod_variables.t:79
integer iw_e
Index of the energy density.
Definition: mod_variables.t:58