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