MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_thermal_conduction.t
Go to the documentation of this file.
1 !> Thermal conduction for HD and MHD
2 !> Adaptation of mod_thermal_conduction for the mod_supertimestepping
3 !> In order to use it set use_mhd_tc=1 (for the mhd impl) or 2 (for the hd impl) in mhd_list (for the mhd module both hd and mhd impl can be used)
4 !> or use_new_hd_tc in hd_list parameters to true
5 !> (for the hd module, hd implementation has to be used)
6 !> The TC is set by calling one
7 !> tc_init_hd_for_total_energy and tc_init_mhd_for_total_energy might
8 !> The second argument: ixArray has to be [rho_,e_,mag(1)] for mhd (Be aware that the other components of the mag field are assumed consecutive) and [rho_,e_] for hd
9 !> additionally when internal energy equation is solved, an additional element of this array is eaux_: the index of the internal energy variable.
10 !>
11 !> 10.07.2011 developed by Chun Xia and Rony Keppens
12 !> 01.09.2012 moved to modules folder by Oliver Porth
13 !> 13.10.2013 optimized further by Chun Xia
14 !> 12.03.2014 implemented RKL2 super timestepping scheme to reduce iterations
15 !> and improve stability and accuracy up to second order in time by Chun Xia.
16 !> 23.08.2014 implemented saturation and perpendicular TC by Chun Xia
17 !> 12.01.2017 modulized by Chun Xia
18 !>
19 !> PURPOSE:
20 !> IN MHD ADD THE HEAT CONDUCTION SOURCE TO THE ENERGY EQUATION
21 !> S=DIV(KAPPA_i,j . GRAD_j T)
22 !> where KAPPA_i,j = tc_k_para b_i b_j + tc_k_perp (I - b_i b_j)
23 !> b_i b_j = B_i B_j / B**2, I is the unit matrix, and i, j= 1, 2, 3 for 3D
24 !> IN HD ADD THE HEAT CONDUCTION SOURCE TO THE ENERGY EQUATION
25 !> S=DIV(tc_k_para . GRAD T)
26 !> USAGE:
27 !> 1. in mod_usr.t -> subroutine usr_init(), add
28 !> unit_length=your length unit
29 !> unit_numberdensity=your number density unit
30 !> unit_velocity=your velocity unit
31 !> unit_temperature=your temperature unit
32 !> before call (m)hd_activate()
33 !> 2. to switch on thermal conduction in the (m)hd_list of amrvac.par add:
34 !> (m)hd_thermal_conduction=.true.
35 !> 3. in the tc_list of amrvac.par :
36 !> tc_perpendicular=.true. ! (default .false.) turn on thermal conduction perpendicular to magnetic field
37 !> tc_saturate=.false. ! (default .true. ) turn off thermal conduction saturate effect
38 !> tc_slope_limiter='MC' ! choose limiter for slope-limited anisotropic thermal conduction in MHD
39 
41  use mod_global_parameters, only: std_len
42  use mod_geometry
43  implicit none
44  private
45  !> Coefficient of thermal conductivity (parallel to magnetic field)
46  double precision, public :: tc_k_para
47 
48  !> Coefficient of thermal conductivity perpendicular to magnetic field
49  double precision, public :: tc_k_perp
50 
51  !> The adiabatic index
52  double precision :: tc_gamma_1
53 
54  !> small_e, see if it is worth having this a global var
55  ! or it is fine to recalculate it inside handle_small_e
56  double precision :: small_e
57  !> Name of slope limiter for transverse component of thermal flux
58  character(len=std_len) :: tc_slope_limiter
59 
60 
61  !> Logical switch for test constant conductivity
62  logical :: tc_constant
63  !> Calculate thermal conduction perpendicular to magnetic field (.true.) or not (.false.)
64  logical, public :: tc_perpendicular=.false.
65 
66  !> Consider thermal conduction saturation effect (.true.) or not (.false.)
67  logical :: tc_saturate=.true.
68 
69  !> Indices of the variables
70  integer :: rho_=-1,mag(1:3)=-1,e_=-1,eaux_=-1,tcoff_=-1
71 
73  abstract interface
74  subroutine get_temperature_subr(w,x,ixI^L,ixO^L,res)
76  integer, intent(in) :: ixi^l, ixo^l
77  double precision, intent(in) :: w(ixi^s,nw)
78  double precision, intent(in) :: x(ixi^s,1:ndim)
79  double precision, intent(out):: res(ixi^s)
80  end subroutine get_temperature_subr
81  end interface
82  procedure(get_temperature_subr), pointer :: get_temperature_from_eint => null()
83  procedure(get_temperature_subr), pointer :: get_temperature_from_conserved => null()
84 
85 contains
86 
87  !> Init TC coeffiecients: MHD case
88  subroutine tc_init_mhd_params(phys_gamma, ixArray)
90  double precision, intent(in) :: phys_gamma
91  integer, intent(in) :: ixarray(:)
92  rho_ = ixarray(1)
93  e_ = ixarray(2)
94  mag(1) = ixarray(3)
95  mag(2) = mag(1) + 1
96  mag(3) = mag(2) + 1
97  if(size(ixarray).eq.4) eaux_ = ixarray(4)
98  if(phys_trac) tcoff_=iw_tcoff
99 
100  tc_gamma_1=phys_gamma-1d0
101  small_e = small_pressure/tc_gamma_1
102  tc_slope_limiter='MC'
103 
104  tc_k_para=0.d0
105 
106  tc_k_perp=0.d0
107 
109 
110  if(tc_k_para==0.d0 .and. tc_k_perp==0.d0) then
111  if(si_unit) then
112  ! Spitzer thermal conductivity with SI units
114  ! thermal conductivity perpendicular to magnetic field
116  else
117  ! Spitzer thermal conductivity with cgs units
119  ! thermal conductivity perpendicular to magnetic field
121  end if
122  else
123  tc_constant=.true.
124  end if
125 
126  contains
127 
128  !> Read tc module parameters from par file: MHD case
129  subroutine tc_params_read_mhd(files)
131  character(len=*), intent(in) :: files(:)
132  integer :: n
133 
134  namelist /tc_list/ tc_perpendicular, tc_saturate, tc_slope_limiter, tc_k_para, tc_k_perp
135 
136  do n = 1, size(files)
137  open(unitpar, file=trim(files(n)), status="old")
138  read(unitpar, tc_list, end=111)
139 111 close(unitpar)
140  end do
141 
142  end subroutine tc_params_read_mhd
143 
144  end subroutine tc_init_mhd_params
145 
146  !> Initialize the module
147  !> this adds the term in the STS methods linked lists when the total energy equation is used
148  !> Params:
149  !> gamma
150  !> ixArray : an array with the indices of the variables
151  !> mhd_get_temperature_from_etot, mhd_get_temperature_from_eint subroutines which calculates temperature
152  !> mhd_ei_to_e, mhd_e_to_ei subroutines which calculates e_tot from e_int and e_int from e_tot
153  subroutine tc_init_mhd_for_total_energy(phys_gamma, ixArray, mhd_get_temperature_from_etot, mhd_get_temperature_from_eint)
157  double precision, intent(in) :: phys_gamma
158  integer, intent(in) :: ixArray(:)
159 
160  interface
161  subroutine mhd_get_temperature_from_etot(w,x,ixI^L,ixO^L,res)
163 
164  integer, intent(in) :: ixI^L, ixO^L
165  double precision, intent(in) :: w(ixi^s,nw)
166  double precision, intent(in) :: x(ixi^s,1:ndim)
167  double precision, intent(out):: res(ixi^s)
168  end subroutine mhd_get_temperature_from_etot
169 
170  subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
172  integer, intent(in) :: ixI^L, ixO^L
173  double precision, intent(in) :: w(ixi^s, 1:nw)
174  double precision, intent(in) :: x(ixi^s, 1:ndim)
175  double precision, intent(out):: res(ixi^s)
176  end subroutine mhd_get_temperature_from_eint
177  end interface
178 
179  call tc_init_mhd_params(phys_gamma, ixarray)
180  call sts_init()
181  get_temperature_from_conserved => mhd_get_temperature_from_etot
182  get_temperature_from_eint => mhd_get_temperature_from_eint
183  call add_sts_method(get_tc_dt_mhd,sts_set_source_tc_mhd,e_,1,e_,1,.false.)
185 
187 
188  end subroutine tc_init_mhd_for_total_energy
189 
190  !> Initialize tc module: MHD case
191  !> this adds the term in the STS methods linked lists when the internal energy equation is used
192  !> Params:
193  !> gamma
194  !> ixArray : an array with the indices of the variables
195  !> mhd_get_temperature_from_eint subroutines which calculates temperature
196  subroutine tc_init_mhd_for_internal_energy(phys_gamma, ixArray, mhd_get_temperature_from_eint)
199  double precision, intent(in) :: phys_gamma
200  integer, intent(in) :: ixArray(:)
201  interface
202  subroutine mhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
204  integer, intent(in) :: ixI^L, ixO^L
205  double precision, intent(in) :: w(ixi^s, 1:nw)
206  double precision, intent(in) :: x(ixi^s, 1:ndim)
207  double precision, intent(out):: res(ixi^s)
208  end subroutine mhd_get_temperature_from_eint
209  end interface
210 
211  call tc_init_mhd_params(phys_gamma, ixarray)
212  call sts_init()
213  get_temperature_from_conserved => mhd_get_temperature_from_eint
214  get_temperature_from_eint => mhd_get_temperature_from_eint
215  call add_sts_method(get_tc_dt_mhd,sts_set_source_tc_mhd,e_,1,e_,1,.false.)
216 
218 
219  end subroutine tc_init_mhd_for_internal_energy
220 
221  subroutine tc_init_hd_params(phys_gamma, ixArray)
223  double precision, intent(in) :: phys_gamma
224  integer, intent(in) :: ixArray(:)
225  rho_ = ixarray(1)
226  e_ = ixarray(2)
227  if(size(ixarray).eq.3) eaux_ = ixarray(3)
228  if(phys_trac) tcoff_=iw_tcoff
229 
230  tc_gamma_1=phys_gamma - 1d0
231  small_e = small_pressure/tc_gamma_1
232  tc_k_para=0.d0
234  if(tc_k_para==0.d0 ) then
235  if(si_unit) then
236  ! Spitzer thermal conductivity with SI units
238  else
239  ! Spitzer thermal conductivity with cgs units
241  end if
242  end if
243 
244  contains
245 
246  !> Read tc parameters from par file: HD case
247  subroutine tc_params_read_hd(files)
249  character(len=*), intent(in) :: files(:)
250  integer :: n
251 
252  namelist /tc_list/ tc_saturate, tc_k_para
253 
254  do n = 1, size(files)
255  open(unitpar, file=trim(files(n)), status="old")
256  read(unitpar, tc_list, end=111)
257 111 close(unitpar)
258  end do
259 
260  end subroutine tc_params_read_hd
261  end subroutine tc_init_hd_params
262 
263  !> Initialize the module for the HD
264  !> this adds the term in the STS methods linked lists when the total energy equation is used
265  !> Params:
266  !> gamma
267  !> ixArray : an array with the indices of the variables
268  !> mhd_get_temperature_from_etot, mhd_get_temperature_from_eint subroutines which calculates temperature
269  !> mhd_ei_to_e, mhd_e_to_ei subroutines which calculates e_tot from e_int and e_int from e_tot
270  subroutine tc_init_hd_for_total_energy(phys_gamma,ixArray,hd_get_temperature_from_etot,hd_get_temperature_from_eint)
274 
275  double precision, intent(in) :: phys_gamma
276  integer, intent(in) :: ixArray(:)
277 
278  interface
279  subroutine hd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
281  integer, intent(in) :: ixI^L, ixO^L
282  double precision, intent(in) :: w(ixi^s, 1:nw)
283  double precision, intent(in) :: x(ixi^s, 1:ndim)
284  double precision, intent(out):: res(ixi^s)
285  end subroutine hd_get_temperature_from_etot
286 
287  subroutine hd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
289  integer, intent(in) :: ixI^L, ixO^L
290  double precision, intent(in) :: w(ixi^s, 1:nw)
291  double precision, intent(in) :: x(ixi^s, 1:ndim)
292  double precision, intent(out):: res(ixi^s)
293  end subroutine hd_get_temperature_from_eint
294  end interface
295 
296  call tc_init_hd_params(phys_gamma, ixarray)
297 
298  get_temperature_from_eint => hd_get_temperature_from_eint
299  get_temperature_from_conserved => hd_get_temperature_from_etot
300  call sts_init()
301  call add_sts_method(get_tc_dt_hd,sts_set_source_tc_hd,e_,1,e_,1,.false.)
303 
305 
306  end subroutine tc_init_hd_for_total_energy
307 
308  !> Initialize the module
309  !> this adds the term in the STS methods linked lists when the internal energy equation is used
310  !> Params:
311  !> gamma
312  !> ixArray : an array with the indices of the variables
313  !> mhd_get_temperature_from_eint subroutines which calculates temperature
314  subroutine tc_init_hd_for_internal_energy(phys_gamma,ixArray,hd_get_temperature_from_eint)
318 
319  double precision, intent(in) :: phys_gamma
320  integer, intent(in) :: ixArray(:)
321 
322  interface
323  subroutine hd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
325  integer, intent(in) :: ixI^L, ixO^L
326  double precision, intent(in) :: w(ixi^s, 1:nw)
327  double precision, intent(in) :: x(ixi^s, 1:ndim)
328  double precision, intent(out):: res(ixi^s)
329  end subroutine hd_get_temperature_from_eint
330  end interface
331 
332  call tc_init_hd_params(phys_gamma, ixarray)
333  get_temperature_from_eint => hd_get_temperature_from_eint
334  get_temperature_from_conserved => hd_get_temperature_from_eint
335  call sts_init()
336  call add_sts_method(get_tc_dt_hd,sts_set_source_tc_hd,e_,1,e_,1,.false.)
337 
339 
340  end subroutine tc_init_hd_for_internal_energy
341 
342  !> Get the explicut timestep for the TC (mhd implementation)
343  function get_tc_dt_mhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
344  !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
345  !where tc_k_para_i=tc_k_para*B_i**2/B**2
346  !and T=p/rho
348 
349  integer, intent(in) :: ixI^L, ixO^L
350  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
351  double precision, intent(in) :: w(ixi^s,1:nw)
352  double precision :: dtnew
353 
354  double precision :: dxinv(1:ndim),mf(ixi^s,1:ndir)
355  double precision :: tmp2(ixi^s),tmp(ixi^s),Te(ixi^s),B2(ixi^s)
356  double precision :: dtdiff_tcond,maxtmp2
357  integer :: idim,ix^D
358 
359  ^d&dxinv(^d)=one/dx^d;
360 
361  !temperature
362  call get_temperature_from_conserved(w,x,ixi^l,ixo^l,te)
363 
364  !tc_k_para_i
365  if(tc_constant) then
366  tmp(ixo^s)=tc_k_para
367  else
368  tmp(ixo^s)=tc_k_para*dsqrt(te(ixo^s)**5)/w(ixo^s,rho_)
369  end if
370 
371  ! B
372  if(b0field) then
373  mf(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
374  else
375  mf(ixo^s,:)=w(ixo^s,iw_mag(:))
376  end if
377  ! B^-2
378  b2(ixo^s)=sum(mf(ixo^s,:)**2,dim=ndim+1)
379  ! B_i**2/B**2
380  where(b2(ixo^s)/=0.d0)
381  ^d&mf(ixo^s,^d)=mf(ixo^s,^d)**2/b2(ixo^s);
382  elsewhere
383  ^d&mf(ixo^s,^d)=1.d0;
384  end where
385 
386  if(tc_saturate) b2(ixo^s)=22.d0*dsqrt(te(ixo^s))
387  dtnew=bigdouble
388  do idim=1,ndim
389  tmp2(ixo^s)=tmp(ixo^s)*mf(ixo^s,idim)
390  if(tc_saturate) then
391  where(tmp2(ixo^s)>b2(ixo^s))
392  tmp2(ixo^s)=b2(ixo^s)
393  end where
394  end if
395  maxtmp2=maxval(tmp2(ixo^s))
396  if(maxtmp2==0.d0) maxtmp2=smalldouble
397  ! dt< dx_idim**2/((gamma-1)*tc_k_para_i/rho*B_i**2/B**2)
398  dtdiff_tcond=1.d0/tc_gamma_1/(maxtmp2*dxinv(idim)**2)
399  ! limit the time step
400  dtnew=min(dtnew,dtdiff_tcond)
401  end do
402  dtnew=dtnew/dble(ndim)
403 
404  end function get_tc_dt_mhd
405 
406  !> anisotropic thermal conduction with slope limited symmetric scheme
407  !> Sharma 2007 Journal of Computational Physics 227, 123
408  subroutine sts_set_source_tc_mhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
410  use mod_fix_conserve
411 
412  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
413  double precision, intent(in) :: x(ixi^s,1:ndim)
414  double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
415  double precision, intent(in) :: my_dt
416  logical, intent(in) :: fix_conserve_at_step
417 
418  !! qd store the heat conduction energy changing rate
419  double precision :: qd(ixi^s)
420  double precision :: qvec(ixi^s,1:ndim)
421 
422  double precision, dimension(ixI^S,1:ndir) :: mf,Bc,Bcf
423  double precision, dimension(ixI^S,1:ndim) :: gradT
424  double precision, dimension(ixI^S) :: Te,ka,kaf,ke,kef,qdd,qe,Binv,minq,maxq,Bnorm
425  double precision, allocatable, dimension(:^D&,:,:) :: fluxall
426  double precision :: alpha,dxinv(ndim)
427  integer, dimension(ndim) :: lowindex
428  integer :: idims,idir,ix^D,ix^L,ixC^L,ixA^L,ixB^L
429 
430  ! coefficient of limiting on normal component
431  if(ndim<3) then
432  alpha=0.75d0
433  else
434  alpha=0.85d0
435  end if
436  ix^l=ixo^l^ladd1;
437 
438  dxinv=1.d0/dxlevel
439 
440  call get_temperature_from_eint(w, x, ixi^l, ixi^l, te) !calculate Te in whole domain (+ghosts)
441  ! T gradient at cell faces
442  ! B vector
443  if(b0field) then
444  mf(ixi^s,:)=w(ixi^s,iw_mag(:))+block%B0(ixi^s,:,0)
445  else
446  mf(ixi^s,:)=w(ixi^s,iw_mag(:))
447  end if
448  ! |B|
449  binv(ix^s)=dsqrt(sum(mf(ix^s,:)**2,dim=ndim+1))
450  where(binv(ix^s)/=0.d0)
451  binv(ix^s)=1.d0/binv(ix^s)
452  elsewhere
453  binv(ix^s)=bigdouble
454  end where
455  ! b unit vector: magnetic field direction vector
456  do idims=1,ndim
457  mf(ix^s,idims)=mf(ix^s,idims)*binv(ix^s)
458  end do
459  ! ixC is cell-corner index
460  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
461  ! b unit vector at cell corner
462  bc=0.d0
463  {do ix^db=0,1\}
464  ixamin^d=ixcmin^d+ix^d;
465  ixamax^d=ixcmax^d+ix^d;
466  bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)+mf(ixa^s,1:ndim)
467  {end do\}
468  bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)*0.5d0**ndim
469  ! T gradient at cell faces
470  gradt=0.d0
471  do idims=1,ndim
472  ixbmin^d=ixmin^d;
473  ixbmax^d=ixmax^d-kr(idims,^d);
474  call gradientc(te,ixi^l,ixb^l,idims,minq)
475  gradt(ixb^s,idims)=minq(ixb^s)
476  end do
477  if(tc_constant) then
478  if(tc_perpendicular) then
479  ka(ixc^s)=tc_k_para-tc_k_perp
480  ke(ixc^s)=tc_k_perp
481  else
482  ka(ixc^s)=tc_k_para
483  end if
484  else
485  ! conductivity at cell center
486  if(phys_trac) then
487  minq(ix^s)=te(ix^s)
488  where(minq(ix^s) < w(ix^s,tcoff_))
489  minq(ix^s)=w(ix^s,tcoff_)
490  end where
491  minq(ix^s)=tc_k_para*sqrt(minq(ix^s)**5)
492  else
493  minq(ix^s)=tc_k_para*sqrt(te(ix^s)**5)
494  end if
495  ka=0.d0
496  {do ix^db=0,1\}
497  ixbmin^d=ixcmin^d+ix^d;
498  ixbmax^d=ixcmax^d+ix^d;
499  ka(ixc^s)=ka(ixc^s)+minq(ixb^s)
500  {end do\}
501  ! cell corner conductivity
502  ka(ixc^s)=0.5d0**ndim*ka(ixc^s)
503  ! compensate with perpendicular conductivity
504  if(tc_perpendicular) then
505  minq(ix^s)=tc_k_perp*w(ix^s,rho_)**2*binv(ix^s)**2/dsqrt(te(ix^s))
506  ke=0.d0
507  {do ix^db=0,1\}
508  ixbmin^d=ixcmin^d+ix^d;
509  ixbmax^d=ixcmax^d+ix^d;
510  ke(ixc^s)=ke(ixc^s)+minq(ixb^s)
511  {end do\}
512  ! cell corner conductivity: k_parallel-k_perpendicular
513  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
514  where(ke(ixc^s)<ka(ixc^s))
515  ka(ixc^s)=ka(ixc^s)-ke(ixc^s)
516  elsewhere
517  ka(ixc^s)=0.d0
518  ke(ixc^s)=ka(ixc^s)
519  end where
520  end if
521  end if
522  if(tc_slope_limiter=='no') then
523  ! calculate thermal conduction flux with symmetric scheme
524  do idims=1,ndim
525  !qd corner values
526  qd=0.d0
527  {do ix^db=0,1 \}
528  if({ ix^d==0 .and. ^d==idims | .or.}) then
529  ixbmin^d=ixcmin^d+ix^d;
530  ixbmax^d=ixcmax^d+ix^d;
531  qd(ixc^s)=qd(ixc^s)+gradt(ixb^s,idims)
532  end if
533  {end do\}
534  ! temperature gradient at cell corner
535  qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
536  end do
537  ! b grad T at cell corner
538  qd(ixc^s)=sum(qvec(ixc^s,1:ndim)*bc(ixc^s,1:ndim),dim=ndim+1)
539  do idims=1,ndim
540  ! TC flux at cell corner
541  gradt(ixc^s,idims)=ka(ixc^s)*bc(ixc^s,idims)*qd(ixc^s)
542  if(tc_perpendicular) gradt(ixc^s,idims)=gradt(ixc^s,idims)+ke(ixc^s)*qvec(ixc^s,idims)
543  end do
544  ! TC flux at cell face
545  qvec=0.d0
546  do idims=1,ndim
547  ixb^l=ixo^l-kr(idims,^d);
548  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
549  {do ix^db=0,1 \}
550  if({ ix^d==0 .and. ^d==idims | .or.}) then
551  ixbmin^d=ixamin^d-ix^d;
552  ixbmax^d=ixamax^d-ix^d;
553  qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
554  end if
555  {end do\}
556  qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
557  if(tc_saturate) then
558  ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
559  ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
560  bcf=0.d0
561  {do ix^db=0,1 \}
562  if({ ix^d==0 .and. ^d==idims | .or.}) then
563  ixbmin^d=ixamin^d-ix^d;
564  ixbmax^d=ixamax^d-ix^d;
565  bcf(ixa^s,idims)=bcf(ixa^s,idims)+bc(ixb^s,idims)
566  end if
567  {end do\}
568  ! averaged b at face centers
569  bcf(ixa^s,idims)=bcf(ixa^s,idims)*0.5d0**(ndim-1)
570  ixb^l=ixa^l+kr(idims,^d);
571  qd(ixa^s)=2.75d0*(w(ixa^s,rho_)+w(ixb^s,rho_))*dsqrt(0.5d0*(te(ixa^s)+te(ixb^s)))**3*dabs(bcf(ixa^s,idims))
572  {do ix^db=ixamin^db,ixamax^db\}
573  if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
574  qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
575  end if
576  {end do\}
577  end if
578  end do
579  else
580  ! calculate thermal conduction flux with slope-limited symmetric scheme
581  qvec=0.d0
582  do idims=1,ndim
583  ixb^l=ixo^l-kr(idims,^d);
584  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
585  ! calculate normal of magnetic field
586  ixb^l=ixa^l+kr(idims,^d);
587  bnorm(ixa^s)=0.5d0*(mf(ixa^s,idims)+mf(ixb^s,idims))
588  bcf=0.d0
589  kaf=0.d0
590  kef=0.d0
591  {do ix^db=0,1 \}
592  if({ ix^d==0 .and. ^d==idims | .or.}) then
593  ixbmin^d=ixamin^d-ix^d;
594  ixbmax^d=ixamax^d-ix^d;
595  bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)+bc(ixb^s,1:ndim)
596  kaf(ixa^s)=kaf(ixa^s)+ka(ixb^s)
597  if(tc_perpendicular) kef(ixa^s)=kef(ixa^s)+ke(ixb^s)
598  end if
599  {end do\}
600  ! averaged b at face centers
601  bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)*0.5d0**(ndim-1)
602  ! averaged thermal conductivity at face centers
603  kaf(ixa^s)=kaf(ixa^s)*0.5d0**(ndim-1)
604  if(tc_perpendicular) kef(ixa^s)=kef(ixa^s)*0.5d0**(ndim-1)
605  ! limited normal component
606  minq(ixa^s)=min(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
607  maxq(ixa^s)=max(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
608  ! eq (19)
609  !corner values of gradT
610  qdd=0.d0
611  {do ix^db=0,1 \}
612  if({ ix^d==0 .and. ^d==idims | .or.}) then
613  ixbmin^d=ixcmin^d+ix^d;
614  ixbmax^d=ixcmax^d+ix^d;
615  qdd(ixc^s)=qdd(ixc^s)+gradt(ixb^s,idims)
616  end if
617  {end do\}
618  ! temperature gradient at cell corner
619  qdd(ixc^s)=qdd(ixc^s)*0.5d0**(ndim-1)
620  ! eq (21)
621  qe=0.d0
622  {do ix^db=0,1 \}
623  qd(ixc^s)=qdd(ixc^s)
624  if({ ix^d==0 .and. ^d==idims | .or.}) then
625  ixbmin^d=ixamin^d-ix^d;
626  ixbmax^d=ixamax^d-ix^d;
627  where(qd(ixb^s)<=minq(ixa^s))
628  qd(ixb^s)=minq(ixa^s)
629  elsewhere(qd(ixb^s)>=maxq(ixa^s))
630  qd(ixb^s)=maxq(ixa^s)
631  end where
632  qvec(ixa^s,idims)=qvec(ixa^s,idims)+bc(ixb^s,idims)**2*qd(ixb^s)
633  if(tc_perpendicular) qe(ixa^s)=qe(ixa^s)+qd(ixb^s)
634  end if
635  {end do\}
636  qvec(ixa^s,idims)=kaf(ixa^s)*qvec(ixa^s,idims)*0.5d0**(ndim-1)
637  ! add normal flux from perpendicular conduction
638  if(tc_perpendicular) qvec(ixa^s,idims)=qvec(ixa^s,idims)+kef(ixa^s)*qe(ixa^s)*0.5d0**(ndim-1)
639  ! limited transverse component, eq (17)
640  ixbmin^d=ixamin^d;
641  ixbmax^d=ixamax^d+kr(idims,^d);
642  do idir=1,ndim
643  if(idir==idims) cycle
644  qd(ixi^s)=slope_limiter(gradt(ixi^s,idir),ixi^l,ixb^l,idir,-1)
645  qd(ixi^s)=slope_limiter(qd,ixi^l,ixa^l,idims,1)
646  qvec(ixa^s,idims)=qvec(ixa^s,idims)+kaf(ixa^s)*bnorm(ixa^s)*bcf(ixa^s,idir)*qd(ixa^s)
647  end do
648 
649  ! consider magnetic null point
650  !where(Binv(ixA^S)==0.d0)
651  ! qvec(ixA^S,idims)=tc_k_para*(0.5d0*(Te(ixA^S)+Te(ixB^S)))**2.5d0*gradT(ixA^S,idims)
652  !end where
653 
654  if(tc_saturate) then
655  ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
656  ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
657  ixb^l=ixa^l+kr(idims,^d);
658  qd(ixa^s)=2.75d0*(w(ixa^s,rho_)+w(ixb^s,rho_))*dsqrt(0.5d0*(te(ixa^s)+te(ixb^s)))**3*dabs(bnorm(ixa^s))
659  {do ix^db=ixamin^db,ixamax^db\}
660  if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
661  ! write(*,*) 'it',it,qvec(ix^D,idims),qd(ix^D),' TC saturated at ',&
662  ! x(ix^D,:),' rho',w(ix^D,rho_),' Te',Te(ix^D)
663  qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
664  end if
665  {end do\}
666  end if
667  end do
668  end if
669 
670  qd=0.d0
671  if(slab_uniform) then
672  do idims=1,ndim
673  qvec(ix^s,idims)=dxinv(idims)*qvec(ix^s,idims)
674  ixb^l=ixo^l-kr(idims,^d);
675  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
676  end do
677  else
678  do idims=1,ndim
679  qvec(ix^s,idims)=qvec(ix^s,idims)*block%surfaceC(ix^s,idims)
680  ixb^l=ixo^l-kr(idims,^d);
681  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
682  end do
683  qd(ixo^s)=qd(ixo^s)/block%dvolume(ixo^s)
684  end if
685 
686  if(fix_conserve_at_step) then
687  allocate(fluxall(ixi^s,1,1:ndim))
688  fluxall(ixi^s,1,1:ndim)=my_dt*qvec(ixi^s,1:ndim)
689  call store_flux(igrid,fluxall,1,ndim,nflux)
690  deallocate(fluxall)
691  end if
692 
693  wres(ixo^s,e_)=qd(ixo^s)
694 
695  end subroutine sts_set_source_tc_mhd
696 
697  function slope_limiter(f,ixI^L,ixO^L,idims,pm) result(lf)
699  integer, intent(in) :: ixI^L, ixO^L, idims, pm
700  double precision, intent(in) :: f(ixi^s)
701  double precision :: lf(ixi^s)
702 
703  double precision :: signf(ixi^s)
704  integer :: ixB^L
705 
706  ixb^l=ixo^l+pm*kr(idims,^d);
707  signf(ixo^s)=sign(1.d0,f(ixo^s))
708  select case(tc_slope_limiter)
709  case('minmod')
710  ! minmod limiter
711  lf(ixo^s)=signf(ixo^s)*max(0.d0,min(abs(f(ixo^s)),signf(ixo^s)*f(ixb^s)))
712  case ('MC')
713  ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
714  lf(ixo^s)=two*signf(ixo^s)* &
715  max(zero,min(dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s),&
716  signf(ixo^s)*quarter*(f(ixb^s)+f(ixo^s))))
717  case ('superbee')
718  ! Roes superbee limiter (eq.3.51i)
719  lf(ixo^s)=signf(ixo^s)* &
720  max(zero,min(two*dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s)),&
721  min(dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s)))
722  case ('koren')
723  ! Barry Koren Right variant
724  lf(ixo^s)=signf(ixo^s)* &
725  max(zero,min(two*dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s),&
726  (two*f(ixb^s)*signf(ixo^s)+dabs(f(ixo^s)))*third))
727  case default
728  call mpistop("Unknown slope limiter for thermal conduction")
729  end select
730 
731  end function slope_limiter
732 
733  !> Calculate gradient of a scalar q at cell interfaces in direction idir
734  subroutine gradientc(q,ixI^L,ixO^L,idir,gradq)
736 
737  integer, intent(in) :: ixI^L, ixO^L, idir
738  double precision, intent(in) :: q(ixi^s)
739  double precision, intent(inout) :: gradq(ixi^s)
740  integer :: jxO^L
741 
742  associate(x=>block%x)
743 
744  jxo^l=ixo^l+kr(idir,^d);
745  select case(coordinate)
746  case(cartesian)
747  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/dxlevel(idir)
748  case(cartesian_stretched)
749  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,idir)-x(ixo^s,idir))
750  case(spherical)
751  select case(idir)
752  case(1)
753  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,1)-x(ixo^s,1))
754  {^nooned
755  case(2)
756  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/( (x(jxo^s,2)-x(ixo^s,2))*x(ixo^s,1) )
757  }
758  {^ifthreed
759  case(3)
760  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/( (x(jxo^s,3)-x(ixo^s,3))*x(ixo^s,1)*dsin(x(ixo^s,2)) )
761  }
762  end select
763  case(cylindrical)
764  if(idir==phi_) then
765  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/((x(jxo^s,phi_)-x(ixo^s,phi_))*x(ixo^s,r_))
766  else
767  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,idir)-x(ixo^s,idir))
768  end if
769  case default
770  call mpistop('Unknown geometry')
771  end select
772 
773  end associate
774  end subroutine gradientc
775 
776  function get_tc_dt_hd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
777  ! Check diffusion time limit dt < dx_i**2 / ((gamma-1)*tc_k_para_i/rho)
779 
780  integer, intent(in) :: ixI^L, ixO^L
781  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
782  double precision, intent(in) :: w(ixi^s,1:nw)
783  double precision :: dtnew
784 
785  double precision :: dxinv(1:ndim), tmp(ixi^s), Te(ixi^s)
786  double precision :: dtdiff_tcond,dtdiff_tsat
787  integer :: idim,ix^D
788 
789  ^d&dxinv(^d)=one/dx^d;
790 
791  call get_temperature_from_conserved(w,x,ixi^l,ixo^l,te)
792 
793  tmp(ixo^s)=tc_gamma_1*tc_k_para*dsqrt((te(ixo^s))**5)/w(ixo^s,rho_)
794  dtnew = bigdouble
795 
796  do idim=1,ndim
797  ! dt< dx_idim**2/((gamma-1)*tc_k_para_idim/rho)
798  dtdiff_tcond=1d0/maxval(tmp(ixo^s)*dxinv(idim)**2)
799  if(tc_saturate) then
800  ! dt< dx_idim**2/((gamma-1)*sqrt(Te)*5*phi)
801  dtdiff_tsat=1d0/maxval(tc_gamma_1*dsqrt(te(ixo^s))*&
802  5.d0*dxinv(idim)**2)
803  ! choose the slower flux (bigger time scale) between classic and saturated
804  dtdiff_tcond=max(dtdiff_tcond,dtdiff_tsat)
805  end if
806  ! limit the time step
807  dtnew=min(dtnew,dtdiff_tcond)
808  end do
809  dtnew=dtnew/dble(ndim)
810 
811  end function get_tc_dt_hd
812 
813  subroutine sts_set_source_tc_hd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
815  use mod_fix_conserve
816 
817  integer, intent(in) :: ixI^L, ixO^L, igrid, nflux
818  double precision, intent(in) :: x(ixi^s,1:ndim)
819  double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
820  double precision, intent(in) :: my_dt
821  logical, intent(in) :: fix_conserve_at_step
822 
823  double precision :: gradT(ixi^s,1:ndim),Te(ixi^s),ke(ixi^s)
824  double precision :: qvec(ixi^s,1:ndim),qd(ixi^s)
825  double precision, allocatable, dimension(:^D&,:,:) :: fluxall
826 
827  double precision :: dxinv(ndim)
828  integer, dimension(ndim) :: lowindex
829  integer :: idims,ix^D,ix^L,ixC^L,ixA^L,ixB^L,ixD^L
830 
831  ix^l=ixo^l^ladd1;
832  ! ixC is cell-corner index
833  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
834 
835  dxinv=1.d0/dxlevel
836 
837  !calculate Te in whole domain (+ghosts)
838  call get_temperature_from_eint(w, x, ixi^l, ixi^l, te)
839 
840  ! cell corner temperature in ke
841  ke=0.d0
842  ixamax^d=ixmax^d; ixamin^d=ixmin^d-1;
843  {do ix^db=0,1\}
844  ixbmin^d=ixamin^d+ix^d;
845  ixbmax^d=ixamax^d+ix^d;
846  ke(ixa^s)=ke(ixa^s)+te(ixb^s)
847  {end do\}
848  ke(ixa^s)=0.5d0**ndim*ke(ixa^s)
849  ! T gradient (central difference) at cell corners
850  gradt=0.d0
851  do idims=1,ndim
852  ixbmin^d=ixmin^d;
853  ixbmax^d=ixmax^d-kr(idims,^d);
854  call gradient(ke,ixi^l,ixb^l,idims,qd)
855  gradt(ixb^s,idims)=qd(ixb^s)
856  end do
857  ! transition region adaptive conduction
858  if(phys_trac) then
859  where(ke(ixi^s) < w(ixi^s,tcoff_))
860  ke(ixi^s)=w(ixi^s,tcoff_)
861  end where
862  end if
863  ! cell corner conduction flux
864  do idims=1,ndim
865  gradt(ixc^s,idims)=gradt(ixc^s,idims)*tc_k_para*sqrt(ke(ixc^s)**5)
866  end do
867 
868  if(tc_saturate) then
869  ! consider saturation with unsigned saturated TC flux = 5 phi rho c**3
870  ! saturation flux at cell center
871  qd(ix^s)=5.d0*w(ix^s,rho_)*dsqrt(te(ix^s)**3)
872  !cell corner values of qd in ke
873  ke=0.d0
874  {do ix^db=0,1\}
875  ixbmin^d=ixcmin^d+ix^d;
876  ixbmax^d=ixcmax^d+ix^d;
877  ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
878  {end do\}
879  ! cell corner saturation flux
880  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
881  ! magnitude of cell corner conduction flux
882  qd(ixc^s)=norm2(gradt(ixc^s,:),dim=ndim+1)
883  {do ix^db=ixcmin^db,ixcmax^db\}
884  if(qd(ix^d)>ke(ix^d)) then
885  ke(ix^d)=ke(ix^d)/qd(ix^d)
886  do idims=1,ndim
887  gradt(ix^d,idims)=ke(ix^d)*gradt(ix^d,idims)
888  end do
889  end if
890  {end do\}
891  end if
892 
893  ! conductionflux at cell face
894  !face center values of gradT in qvec
895  qvec=0.d0
896  do idims=1,ndim
897  ixb^l=ixo^l-kr(idims,^d);
898  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
899  {do ix^db=0,1 \}
900  if({ ix^d==0 .and. ^d==idims | .or.}) then
901  ixbmin^d=ixamin^d-ix^d;
902  ixbmax^d=ixamax^d-ix^d;
903  qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
904  end if
905  {end do\}
906  qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
907  end do
908 
909  qd=0.d0
910  if(slab_uniform) then
911  do idims=1,ndim
912  qvec(ix^s,idims)=dxinv(idims)*qvec(ix^s,idims)
913  ixb^l=ixo^l-kr(idims,^d);
914  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
915  end do
916  else
917  do idims=1,ndim
918  qvec(ix^s,idims)=qvec(ix^s,idims)*block%surfaceC(ix^s,idims)
919  ixb^l=ixo^l-kr(idims,^d);
920  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
921  end do
922  qd(ixo^s)=qd(ixo^s)/block%dvolume(ixo^s)
923  end if
924 
925  if(fix_conserve_at_step) then
926  allocate(fluxall(ixi^s,1,1:ndim))
927  fluxall(ixi^s,1,1:ndim)=my_dt*qvec(ixi^s,1:ndim)
928  call store_flux(igrid,fluxall,1,ndim,nflux)
929  deallocate(fluxall)
930  end if
931 
932  wres(ixo^s,e_)=qd(ixo^s)
933 
934  end subroutine sts_set_source_tc_hd
935 
936  subroutine handle_small_e(w, x, ixI^L, ixO^L, step)
938  use mod_small_values
939  integer, intent(in) :: ixI^L,ixO^L
940  double precision, intent(inout) :: w(ixi^s,1:nw)
941  double precision, intent(in) :: x(ixi^s,1:ndim)
942  integer, intent(in) :: step
943 
944  integer :: idir
945  logical :: flag(ixi^s,1:nw)
946  character(len=140) :: error_msg
947 
948  flag=.false.
949  where(w(ixo^s,e_)<small_e) flag(ixo^s,e_)=.true.
950  if(any(flag(ixo^s,e_))) then
951  select case (small_values_method)
952  case ("replace")
953  where(flag(ixo^s,e_)) w(ixo^s,e_)=small_e
954  case ("average")
955  call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
956  case default
957  ! small values error shows primitive variables
958  w(ixo^s,e_)=w(ixo^s,e_)*tc_gamma_1
959  do idir = 1, ndir
960  w(ixo^s, iw_mom(idir)) = w(ixo^s, iw_mom(idir))/w(ixo^s,rho_)
961  end do
962  write(error_msg,"(a,i3)") "Thermal conduction step ", step
963  call small_values_error(w, x, ixi^l, ixo^l, flag, error_msg)
964  end select
965  end if
966  end subroutine handle_small_e
967 
968 end module mod_thermal_conduction
integer, parameter cylindrical
Definition: mod_geometry.t:9
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, parameter cartesian_stretched
Definition: mod_geometry.t:8
subroutine tc_params_read_hd(files)
Read tc parameters from par file: HD case.
subroutine tc_init_hd_params(phys_gamma, ixArray)
double precision function, dimension(ixi^s) slope_limiter(f, ixIL, ixOL, idims, pm)
logical, public tc_perpendicular
Calculate thermal conduction perpendicular to magnetic field (.true.) or not (.false.)
double precision function get_tc_dt_hd(w, ixIL, ixOL, dxD, x)
double precision unit_length
Physical scaling factor for length.
double precision unit_density
Physical scaling factor for density.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
subroutine, public set_error_handling_to_head(sts_error_handling)
Set the hook of error handling in the STS update. This method is called before updating the BC...
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine, public tc_init_hd_for_total_energy(phys_gamma, ixArray, hd_get_temperature_from_etot, hd_get_temperature_from_eint)
Initialize the module for the HD this adds the term in the STS methods linked lists when the total en...
integer coordinate
Definition: mod_geometry.t:6
Module for flux conservation near refinement boundaries.
double precision function get_tc_dt_mhd(w, ixIL, ixOL, dxD, x)
Get the explicut timestep for the TC (mhd implementation)
logical b0field
split magnetic field as background B0 field
subroutine, public set_conversion_methods_to_head(sts_before_first_cycle, sts_after_last_cycle)
Set the hooks called before the first cycle and after the last cycle in the STS update This method sh...
procedure(get_temperature_subr), pointer get_temperature_from_eint
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
subroutine, public sts_init()
Initialize sts module.
subroutine handle_small_e(w, x, ixIL, ixOL, step)
subroutine gradientc(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q at cell interfaces in direction idir.
subroutine sts_set_source_tc_hd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
procedure(sub_convert), pointer phys_e_to_ei
Definition: mod_physics.t:65
double precision small_pressure
integer, parameter cartesian
Definition: mod_geometry.t:7
subroutine, public tc_init_mhd_for_internal_energy(phys_gamma, ixArray, mhd_get_temperature_from_eint)
Initialize tc module: MHD case this adds the term in the STS methods linked lists when the internal e...
Thermal conduction for HD and MHD Adaptation of mod_thermal_conduction for the mod_supertimestepping ...
subroutine tc_params_read_mhd(files)
Read tc module parameters from par file: MHD case.
double precision unit_velocity
Physical scaling factor for velocity.
subroutine sts_set_source_tc_mhd(ixIL, ixOL, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
integer, parameter unitpar
file handle for IO
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
subroutine, public tc_init_mhd_for_total_energy(phys_gamma, ixArray, mhd_get_temperature_from_etot, mhd_get_temperature_from_eint)
Initialize the module this adds the term in the STS methods linked lists when the total energy equati...
double precision, public tc_k_para
Coefficient of thermal conductivity (parallel to magnetic field)
subroutine, public small_values_average(ixIL, ixOL, w, x, w_flag, windex)
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
double precision, public tc_k_perp
Coefficient of thermal conductivity perpendicular to magnetic field.
Module for handling problematic values in simulations, such as negative pressures.
double precision unit_magneticfield
Physical scaling factor for magnetic field.
double precision unit_temperature
Physical scaling factor for temperature.
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:208
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
subroutine, public small_values_error(wprim, x, ixIL, ixOL, w_flag, subname)
double precision, dimension(ndim) dxlevel
type(state), pointer block
Block pointer for using one block and its previous state.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
procedure(sub_convert), pointer phys_ei_to_e
Definition: mod_physics.t:64
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
integer, parameter spherical
Definition: mod_geometry.t:10
subroutine, public add_sts_method(sts_getdt, sts_set_sources, startVar, nflux, startwbc, nwbc, evolve_B)
subroutine which added programatically a term to be calculated using STS Params: sts_getdt function c...
character(len=20), public small_values_method
How to handle small values.
double precision unit_numberdensity
Physical scaling factor for number density.
subroutine, public tc_init_hd_for_internal_energy(phys_gamma, ixArray, hd_get_temperature_from_eint)
Initialize the module this adds the term in the STS methods linked lists when the internal energy equ...