MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_hd_hllc.t
Go to the documentation of this file.
1 !> Hydrodynamics HLLC module
2 module mod_hd_hllc
3  use mod_hd_phys
4 
5  implicit none
6  private
7 
8  public :: hd_hllc_init
9 
10 contains
11 
12  subroutine hd_hllc_init()
14 
18  end subroutine hd_hllc_init
19 
20  subroutine hd_diffuse_hllcd(ixI^L,ixO^L,idim,wLC,wRC,fLC,fRC,patchf)
21 
22  ! when method is hllcd or hllcd1 then:
23 
24  ! this subroutine is to enforce regions where we AVOID HLLC
25  ! and use TVDLF instead: this is achieved by setting patchf to 4 in
26  ! certain regions. An additional input parameter is nxdiffusehllc
27  ! which sets the size of the fallback region.
28 
30 
31  integer, intent(in) :: ixI^L,ixO^L,idim
32  double precision, dimension(ixI^S,1:nw), intent(in) :: wRC,wLC
33  double precision, dimension(ixI^S,1:nwflux),intent(in) :: fLC, fRC
34  integer , dimension(ixI^S), intent(inout) :: patchf
35 
36  integer :: ixOO^D,TxOO^L
37 
38  ! In a user-controlled region around any point with flux sign change between
39  ! left and right, ensure fallback to TVDLF
40  {do ixoo^d= ixo^lim^d\}
41  {
42  txoomin^d= max(ixoo^d - nxdiffusehllc*kr(idim,^d), ixomin^d);
43  txoomax^d= min(ixoo^d + nxdiffusehllc*kr(idim,^d), ixomax^d);
44  \}
45  if(abs(patchf(ixoo^d)) == 1 .or. abs(patchf(ixoo^d)) == 4)Then
46  if(any(frc(ixoo^d,1:nwflux)*flc(ixoo^d,1:nwflux)<-smalldouble))Then
47  where(abs(patchf(txoo^s))==1)
48  patchf(txoo^s) = 4
49  endwhere
50  endif
51  endif
52  {enddo^d&\}
53 
54  end subroutine hd_diffuse_hllcd
55 
56  subroutine hd_get_lcd(wLC,wRC,fLC,fRC,cmin,cmax,idim,ixI^L,ixO^L, &
57  whll,Fhll,lambdaCD,patchf)
58  ! Calculate lambda at CD and set the patchf to know the orientation
59  ! of the riemann fan and decide on the flux choice
60  ! We also compute here the HLL flux and w value, for fallback strategy
61 
62  ! In this nul version, we simply compute nothing and ensure TVDLF fallback
64 
65  integer, intent(in) :: ixI^L,ixO^L,idim
66  double precision, dimension(ixI^S,1:nw), intent(in) :: wLC,wRC
67  double precision, dimension(ixI^S,1:nwflux), intent(in) :: fLC,fRC
68  double precision, dimension(ixI^S), intent(in) :: cmax,cmin
69  integer , dimension(ixI^S), intent(inout) :: patchf
70  double precision, dimension(ixI^S,1:nwflux), intent(out) :: Fhll,whll
71  double precision, dimension(ixI^S), intent(out) :: lambdaCD
72 
73  logical , dimension(ixI^S) :: Cond_patchf
74  double precision :: Epsilon
75  integer :: iw,ix^D
76 
77  !--------------------------------------------
78  ! on entry, patch is preset to contain values from -2,1,2,4
79  ! -2: take left flux, no computation here
80  ! +2: take right flux, no computation here
81  ! +4: take TVDLF flux, no computation here
82  ! 1: compute the characteristic speed for the CD
83 
84  cond_patchf(ixo^s)=(abs(patchf(ixo^s))==1)
85  lambdacd=0.d0
86 
87  do iw=1,nwflux
88  where(cond_patchf(ixo^s))
89  !============= compute HLL flux ==============!
90  fhll(ixo^s,iw)= (cmax(ixo^s)*flc(ixo^s,iw)-cmin(ixo^s)*frc(ixo^s,iw) &
91  + cmin(ixo^s)*cmax(ixo^s)*(wrc(ixo^s,iw)-wlc(ixo^s,iw)))&
92  /(cmax(ixo^s)-cmin(ixo^s))
93  !======== compute intermediate HLL state =======!
94  whll(ixo^s,iw) = (cmax(ixo^s)*wrc(ixo^s,iw)-cmin(ixo^s)*wlc(ixo^s,iw)&
95  +flc(ixo^s,iw)-frc(ixo^s,iw))/(cmax(ixo^s)-cmin(ixo^s))
96  end where
97  end do
98 
99  ! deduce the characteristic speed at the CD
100  where(cond_patchf(ixo^s))
101  lambdacd(ixo^s)=whll(ixo^s,mom(idim))/whll(ixo^s,rho_)
102  end where
103 
104  {do ix^db=ixomin^db,ixomax^db\}
105  if(cond_patchf(ix^d)) then
106  ! double check whether obtained speed is in between min and max speeds given
107  ! and identify in which part of the Riemann fan the time-axis is
108  if(cmin(ix^d) < zero .and. lambdacd(ix^d)>zero&
109  .and.lambdacd(ix^d)<cmax(ix^d)) then
110  patchf(ix^d) = -1
111  else if(cmax(ix^d) > zero .and. lambdacd(ix^d) < zero&
112  .and.lambdacd(ix^d)>cmin(ix^d)) then
113  patchf(ix^d) = 1
114  else if(lambdacd(ix^d) >= cmax(ix^d) .or. &
115  lambdacd(ix^d) <= cmin(ix^d)) then
116  lambdacd(ix^d) = zero
117  ! we will fall back to HLL flux case in this degeneracy
118  patchf(ix^d) = 3
119  end if
120  end if
121  {end do\}
122 
123  where(patchf(ixo^s)== 3)
124  cond_patchf(ixo^s)=.false.
125  end where
126 
127  ! handle the specific case where the time axis is exactly on the CD
128  if(any(cond_patchf(ixo^s).and.lambdacd(ixo^s)==zero))then
129  ! determine which sector (forward or backward) of the Riemann fan is smallest
130  ! and select left or right flux accordingly
131  {do ix^db=ixomin^db,ixomax^db\}
132  if(lambdacd(ix^d)==zero.and.cond_patchf(ix^d)) then
133  if(-cmin(ix^d)>=cmax(ix^d)) then
134  patchf(ix^d) = 1
135  else
136  patchf(ix^d) = -1
137  end if
138  end if
139  {end do\}
140  end if
141 
142  ! eigenvalue lambda for contact is near zero: decrease noise by this trick
143  if(flathllc)then
144  epsilon=1.0d-6
145  where(cond_patchf(ixo^s).and. &
146  dabs(lambdacd(ixo^s))/max(cmax(ixo^s),epsilon)< epsilon .and. &
147  dabs(lambdacd(ixo^s))/max(dabs(cmin(ixo^s)),epsilon)< epsilon)
148  lambdacd(ixo^s) = zero
149  end where
150  end if
151 
152  end subroutine hd_get_lcd
153 
154  subroutine hd_get_wcd(wLC,wRC,whll,fRC,fLC,Fhll,patchf,lambdaCD,cmin,cmax,&
155  ixI^L,ixO^L,idim,f)
156  ! compute the intermediate state U*
157  ! only needed where patchf=-1/1
158 
159  ! reference Li S., JCP, 203, 2005, 344-357
160  ! reference T. Miyoski, Kusano JCP, 2008, 2005.
161 
163  use mod_dust
164 
165  integer, intent(in) :: ixI^L,ixO^L,idim
166  double precision, dimension(ixI^S,1:nw), intent(in) :: wRC,wLC
167  double precision, dimension(ixI^S,1:nwflux), intent(in):: whll, Fhll
168  double precision, dimension(ixI^S), intent(in) :: lambdaCD
169  double precision, dimension(ixI^S), intent(in) :: cmax,cmin
170  double precision, dimension(ixI^S,1:nwflux), intent(in):: fRC,fLC
171  double precision, dimension(ixI^S,1:nwflux),intent(out):: f
172  double precision, dimension(ixI^S,1:nw) :: wCD,wSub
173  double precision, dimension(ixI^S,1:nwflux) :: fSub
174  double precision, dimension(ixI^S) :: vSub,cspeed,pCD
175  integer , dimension(ixI^S), intent(in) :: patchf
176 
177  double precision :: csmls
178  integer :: n, iw, ix^D
179 
180  !-------------- auxiliary Speed and array-------------!
181  {do ix^db=ixomin^db,ixomax^db\}
182  if(patchf(ix^d)==1) then
183  cspeed(ix^d)=cmax(ix^d)
184  vsub(ix^d)=wrc(ix^d,mom(idim))/wrc(ix^d,rho_)
185  wsub(ix^d,:)=wrc(ix^d,:)
186  fsub(ix^d,:)=frc(ix^d,:)
187  else if(patchf(ix^d)==-1) then
188  cspeed(ix^d)=cmin(ix^d)
189  vsub(ix^d)=wlc(ix^d,mom(idim))/wlc(ix^d,rho_)
190  wsub(ix^d,:)=wlc(ix^d,:)
191  fsub(ix^d,:)=flc(ix^d,:)
192  end if
193  {end do\}
194 
195  {do ix^db=ixomin^db,ixomax^db\}
196  if(abs(patchf(ix^d))==1) then
197  csmls=one/(cspeed(ix^d)-lambdacd(ix^d))
198  wcd(ix^d,rho_) = wsub(ix^d,rho_)&
199  *(cspeed(ix^d)-vsub(ix^d))*csmls
200  do n=1,hd_n_tracer
201  iw = tracer(n)
202  wcd(ix^d,iw) = wsub(ix^d,iw)*(cspeed(ix^d)-vsub(ix^d))*csmls
203  end do
204 
205  !------- Momentum ------!
206  do iw=1, ndir
207  if(iw /= idim)then
208  ! eq. 21 22
209  wcd(ix^d,mom(iw))=(cspeed(ix^d)*wsub(ix^d,mom(iw))-fsub(ix^d,mom(iw)))*csmls
210  else
211  ! eq. 20
212  wcd(ix^d,mom(iw)) = wcd(ix^d,rho_) * lambdacd(ix^d)
213  endif
214  enddo
215  if(hd_energy) then
216  pcd(ix^d) = wsub(ix^d,rho_)*(cspeed(ix^d)-vsub(ix^d))&
217  *(lambdacd(ix^d)-vsub(ix^d))&
218  +fsub(ix^d,mom(idim))-wsub(ix^d,mom(idim))*vsub(ix^d)
219  ! Eq 31
220  wcd(ix^d,e_) = (cspeed(ix^d)*wsub(ix^d,e_) &
221  -fsub(ix^d,e_)+lambdacd(ix^d)*pcd(ix^d))*csmls
222  end if
223  !if(hd_dust) then
224  ! do n=1,dust_n_species
225  ! wCD(ix^D,dust_rho(n)) = wSub(ix^D,dust_rho(n))*(cspeed(ix^D)-vSub(ix^D))*csmls
226  ! do iw=1,ndir
227  ! if(iw /= idim)then
228  ! ! eq. 21 22
229  ! wCD(ix^D,dust_mom(iw,n))=wSub(ix^D,dust_mom(iw,n))*(cspeed(ix^D)-vSub(ix^D))*csmls
230  ! else
231  ! ! eq. 20
232  ! wCD(ix^D,dust_mom(iw,n)) = wCD(ix^D,dust_rho(n)) * lambdaCD(ix^D)
233  ! endif
234  ! end do
235  ! end do
236  !end if
237 
238  if(hd_dust) then
239  do iw=1,nwflux
240  if(iw>=dust_rho(1).and.iw<=dust_mom(ndir,dust_n_species)) then
241  ! use HLL flux for dust
242  f(ix^d,iw)=fhll(ix^d,iw)
243  else
244  ! f_i=fsub+lambda (wCD-wSub)
245  f(ix^d,iw)=fsub(ix^d,iw)+cspeed(ix^d)*(wcd(ix^d,iw)-wsub(ix^d,iw))
246  end if
247  end do
248  else
249  do iw=1,nwflux
250  ! f_i=fsub+lambda (wCD-wSub)
251  f(ix^d,iw)=fsub(ix^d,iw)+cspeed(ix^d)*(wcd(ix^d,iw)-wsub(ix^d,iw))
252  end do
253  end if
254  end if
255  {end do\}
256 
257  end subroutine hd_get_wcd
258 
259 end module mod_hd_hllc
Module for including dust species, which interact with the gas through a drag force.
Definition: mod_dust.t:3
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(3, 3) kr
Kronecker delta tensor.
Hydrodynamics HLLC module.
Definition: mod_hd_hllc.t:2
subroutine, public hd_hllc_init()
Definition: mod_hd_hllc.t:13
subroutine hd_get_wcd(wLC, wRC, whll, fRC, fLC, Fhll, patchf, lambdaCD, cmin, cmax, ixIL, ixOL, idim, f)
Definition: mod_hd_hllc.t:156
subroutine hd_get_lcd(wLC, wRC, fLC, fRC, cmin, cmax, idim, ixIL, ixOL, whll, Fhll, lambdaCD, patchf)
Definition: mod_hd_hllc.t:58
subroutine hd_diffuse_hllcd(ixIL, ixOL, idim, wLC, wRC, fLC, fRC, patchf)
Definition: mod_hd_hllc.t:21
Hydrodynamics physics module.
Definition: mod_hd_phys.t:2
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition: mod_hd_phys.t:51
integer, public, protected rho_
Index of the density (in the w array)
Definition: mod_hd_phys.t:48
procedure(sub_get_wcd), pointer phys_get_wcd
procedure(sub_get_lcd), pointer phys_get_lcd
procedure(sub_diffuse_hllcd), pointer phys_diffuse_hllcd