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