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 !>
3 !> 10.07.2011 developed by Chun Xia and Rony Keppens
4 !> 01.09.2012 moved to modules folder by Oliver Porth
5 !> 13.10.2013 optimized further by Chun Xia
6 !> 12.03.2014 implemented RKL2 super timestepping scheme to reduce iterations
7 !> and improve stability and accuracy up to second order in time by Chun Xia.
8 !> 23.08.2014 implemented saturation and perpendicular TC by Chun Xia
9 !> 12.01.2017 modulized by Chun Xia
10 !>
11 !> PURPOSE:
12 !> IN MHD ADD THE HEAT CONDUCTION SOURCE TO THE ENERGY EQUATION
13 !> S=DIV(KAPPA_i,j . GRAD_j T)
14 !> where KAPPA_i,j = tc_k_para b_i b_j + tc_k_perp (I - b_i b_j)
15 !> b_i b_j = B_i B_j / B**2, I is the unit matrix, and i, j= 1, 2, 3 for 3D
16 !> IN HD ADD THE HEAT CONDUCTION SOURCE TO THE ENERGY EQUATION
17 !> S=DIV(tc_k_para . GRAD T)
18 !> USAGE:
19 !> 1. in mod_usr.t -> subroutine usr_init(), add
20 !> unit_length=<your length unit>
21 !> unit_numberdensity=<your number density unit>
22 !> unit_velocity=<your velocity unit>
23 !> unit_temperature=<your temperature unit>
24 !> before call (m)hd_activate()
25 !> 2. to switch on thermal conduction in the (m)hd_list of amrvac.par add:
26 !> (m)hd_thermal_conduction=.true.
27 !> 3. in the tc_list of amrvac.par :
28 !> tc_perpendicular=.true. ! (default .false.) turn on thermal conduction perpendicular to magnetic field
29 !> tc_saturate=.false. ! (default .true. ) turn off thermal conduction saturate effect
30 !> tc_dtpar=0.9/0.45/0.3 ! stable time step coefficient for 1D/2D/3D, decrease it for more stable run
31 !> tc_slope_limiter='MC' ! choose limiter for slope-limited anisotropic thermal conduction in MHD
32 
34  use mod_global_parameters, only: std_len
35  implicit none
36  !> Coefficient of thermal conductivity (parallel to magnetic field)
37  double precision, public :: tc_k_para
38 
39  !> Coefficient of thermal conductivity perpendicular to magnetic field
40  double precision, public :: tc_k_perp
41 
42  !> Time step of thermal conduction
43  double precision :: dt_tc
44 
45  !> Number of sub-steps of supertime stepping
46  integer, public :: s
47 
48  !> Index of the density (in the w array)
49  integer, private :: rho_
50 
51  !> Index of the energy density (-1 if not present)
52  integer, private, protected :: e_
53 
54  !> The adiabatic index
55  double precision, private :: tc_gamma
56 
57  !> The small_est allowed energy
58  double precision, private :: small_e
59 
60  !> Time step coefficient
61  double precision, private :: tc_dtpar=0.9d0
62 
63  !> Maximal substeps of TC within one fluid time step to limit fluid time step
64  integer :: tc_ncycles=1000
65 
66  !> Calculate thermal conduction perpendicular to magnetic field (.true.) or not (.false.)
67  logical, private :: tc_perpendicular=.false.
68 
69  !> Consider thermal conduction saturation effect (.true.) or not (.false.)
70  logical, private :: tc_saturate=.true.
71 
72  !> Logical switch for prepare mpi datatype only once
73  logical, private :: first=.true.
74 
75  !> Logical switch for test constant conductivity
76  logical, private :: tc_constant=.false.
77 
78  !> Whether to conserve fluxes at the current partial step
79  logical :: fix_conserve_at_step = .true.
80 
81  !> Name of slope limiter for transverse component of thermal flux
82  character(len=std_len), private :: tc_slope_limiter
83 
84  procedure(thermal_conduction), pointer :: phys_thermal_conduction => null()
85  procedure(get_heatconduct), pointer :: phys_get_heatconduct => null()
86  procedure(getdt_heatconduct), pointer :: phys_getdt_heatconduct => null()
87 
88  abstract interface
89  subroutine thermal_conduction
90  ! Meyer 2012 MNRAS 422,2102
93  end subroutine thermal_conduction
94 
95  subroutine get_heatconduct(tmp,tmp1,tmp2,ixI^L,ixO^L,w,x,qvec)
97 
98  integer, intent(in) :: ixI^L, ixO^L
99  double precision, intent(in) :: x(ixi^s,1:ndim), w(ixi^s,1:nw)
100  !! tmp store the heat conduction energy changing rate
101  double precision, intent(out) :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
102  double precision :: qvec(ixi^s,1:ndim)
103  end subroutine get_heatconduct
104 
105  subroutine getdt_heatconduct(w,ixG^L,ix^L,dtnew,dx^D,x)
107 
108  integer, intent(in) :: ixG^L, ix^L
109  double precision, intent(in) :: dx^D, x(ixg^s,1:ndim)
110  ! note that depending on small_values_method=='error' etc, w values may change
111  ! through call to getpthermal
112  double precision, intent(inout) :: w(ixg^s,1:nw), dtnew
113  end subroutine getdt_heatconduct
114  end interface
115 
116 contains
117  !> Read this module"s parameters from a file
118  subroutine tc_params_read(files)
120  character(len=*), intent(in) :: files(:)
121  integer :: n
122 
123  namelist /tc_list/ tc_perpendicular, tc_saturate, tc_dtpar, tc_slope_limiter, tc_k_para, tc_k_perp, tc_ncycles
124 
125  do n = 1, size(files)
126  open(unitpar, file=trim(files(n)), status="old")
127  read(unitpar, tc_list, end=111)
128 111 close(unitpar)
129  end do
130 
131  end subroutine tc_params_read
132 
133  !> Initialize the module
134  subroutine thermal_conduction_init(phys_gamma)
136  use mod_physics
137 
138  double precision, intent(in) :: phys_gamma
139 
140  tc_gamma=phys_gamma
141 
142  tc_dtpar=tc_dtpar/dble(ndim)
143 
144  tc_slope_limiter='MC'
145 
146  tc_k_para=0.d0
147 
148  tc_k_perp=0.d0
149 
151 
153  if(physics_type=='hd') then
156  else if(physics_type=='mhd') then
159  end if
160 
161  rho_ = iw_rho
162  e_ = iw_e
163 
164  small_e = small_pressure/(tc_gamma - 1.0d0)
165 
166  if(tc_k_para==0.d0 .and. tc_k_perp==0.d0) then
167  if(si_unit) then
168  ! Spitzer thermal conductivity with SI units
170  ! thermal conductivity perpendicular to magnetic field
172  else
173  ! Spitzer thermal conductivity with cgs units
175  ! thermal conductivity perpendicular to magnetic field
177  end if
178  else
179  tc_constant=.true.
180  end if
181 
182  end subroutine thermal_conduction_init
183 
184  subroutine do_thermal_conduction
185  ! Meyer 2012 MNRAS 422,2102
188  use mod_fix_conserve
189 
190  double precision :: omega1,cmu,cmut,cnu,cnut
191  double precision, allocatable :: bj(:)
192  integer:: iigrid, igrid,j
193  logical :: evenstep, stagger_flag
194 
195  ! not do fix conserve and getbc for staggered values if stagger is used
196  stagger_flag=stagger_grid
197  stagger_grid=.false.
198 
199  ! point bc mpi datatype to partial type for thermalconduction
206  ! create bc mpi datatype for ghostcells update
207  if(first) then
208  call create_bc_mpi_datatype(e_,1)
209  first=.false.
210  end if
211 
212  call init_comm_fix_conserve(1,ndim,1)
214 
215  do iigrid=1,igridstail; igrid=igrids(iigrid);
216  if(.not. allocated(ps2(igrid)%w)) allocate(ps2(igrid)%w(ixg^t,1:nw))
217  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
218  ps1(igrid)%w=ps(igrid)%w
219  ps2(igrid)%w=ps(igrid)%w
220  ps3(igrid)%w=ps(igrid)%w
221  end do
222 
223  allocate(bj(0:s))
224  bj(0)=1.d0/3.d0
225  bj(1)=bj(0)
226  if(s>1) then
227  omega1=4.d0/dble(s**2+s-2)
228  cmut=omega1/3.d0
229  else
230  omega1=0.d0
231  cmut=1.d0
232  endif
233 
234  !$OMP PARALLEL DO PRIVATE(igrid)
235  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
236  block=>ps(igrid)
239  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
240  call evolve_step1(igrid,cmut,dt_tc,ixg^ll,ixm^ll,ps1(igrid)%w,ps(igrid)%w,&
241  ps(igrid)%x,ps3(igrid)%w)
242  end do
243  !$OMP END PARALLEL DO
244  ! fix conservation of AMR grid by replacing flux from finer neighbors
245  if (fix_conserve_at_step) then
246  call recvflux(1,ndim)
247  call sendflux(1,ndim)
248  call fix_conserve(ps1,1,ndim,e_,1)
249  end if
250  bcphys=.false.
251  call getbc(global_time,0.d0,ps1,e_,1)
252  if(s==1) then
253  do iigrid=1,igridstail; igrid=igrids(iigrid);
254  ps(igrid)%w(ixg^t,e_)=ps1(igrid)%w(ixg^t,e_)
255  end do
256  ! point bc mpi data type back to full type for (M)HD
263  bcphys=.true.
264  stagger_grid=stagger_flag
265  deallocate(bj)
266  return
267  endif
268  evenstep=.true.
269  do j=2,s
270  bj(j)=dble(j**2+j-2)/dble(2*j*(j+1))
271  cmu=dble(2*j-1)/dble(j)*bj(j)/bj(j-1)
272  cmut=omega1*cmu
273  cnu=dble(1-j)/dble(j)*bj(j)/bj(j-2)
274  cnut=(bj(j-1)-1.d0)*cmut
275  if(evenstep) then
276  !$OMP PARALLEL DO PRIVATE(igrid)
277  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
278  block=>ps(igrid)
281  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
282  call evolve_stepj(igrid,cmu,cmut,cnu,cnut,dt_tc,ixg^ll,ixm^ll,ps1(igrid)%w,&
283  ps2(igrid)%w,ps(igrid)%w,ps(igrid)%x,ps3(igrid)%w)
284  end do
285  !$OMP END PARALLEL DO
286  ! fix conservation of AMR grid by replacing flux from finer neighbors
287  if (fix_conserve_at_step) then
288  call recvflux(1,ndim)
289  call sendflux(1,ndim)
290  call fix_conserve(ps2,1,ndim,e_,1)
291  end if
292  call getbc(global_time,0.d0,ps2,e_,1)
293  evenstep=.false.
294  else
295  !$OMP PARALLEL DO PRIVATE(igrid)
296  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
297  block=>ps(igrid)
300  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
301  call evolve_stepj(igrid,cmu,cmut,cnu,cnut,dt_tc,ixg^ll,ixm^ll,ps2(igrid)%w,&
302  ps1(igrid)%w,ps(igrid)%w,ps(igrid)%x,ps3(igrid)%w)
303  end do
304  !$OMP END PARALLEL DO
305  ! fix conservation of AMR grid by replacing flux from finer neighbors
306  if (fix_conserve_at_step) then
307  call recvflux(1,ndim)
308  call sendflux(1,ndim)
309  call fix_conserve(ps1,1,ndim,e_,1)
310  end if
311  call getbc(global_time,0.d0,ps1,e_,1)
312  evenstep=.true.
313  end if
314  end do
315  if(evenstep) then
316  do iigrid=1,igridstail; igrid=igrids(iigrid);
317  ps(igrid)%w(ixg^t,e_)=ps1(igrid)%w(ixg^t,e_)
318  end do
319  else
320  do iigrid=1,igridstail; igrid=igrids(iigrid);
321  ps(igrid)%w(ixg^t,e_)=ps2(igrid)%w(ixg^t,e_)
322  end do
323  end if
324  deallocate(bj)
325  ! point bc mpi data type back to full type for (M)HD
332  bcphys=.true.
333 
334  ! restore stagger_grid value
335  stagger_grid=stagger_flag
336 
337  end subroutine do_thermal_conduction
338 
339  subroutine evolve_stepj(igrid,qcmu,qcmut,qcnu,qcnut,qdt,ixI^L,ixO^L,w1,w2,w,x,wold)
341  use mod_fix_conserve
342 
343  integer, intent(in) :: igrid,ixI^L,ixO^L
344  double precision, intent(in) :: qcmu,qcmut,qcnu,qcnut,qdt
345  double precision, intent(in) :: w1(ixi^s,1:nw),w(ixi^s,1:nw),wold(ixi^s,1:nw)
346  double precision, intent(in) :: x(ixi^s,1:ndim)
347  double precision, intent(inout) :: w2(ixi^s,1:nw)
348 
349  double precision :: fC(ixi^s,1:1,1:ndim)
350  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
351 
352  call phys_get_heatconduct(tmp,tmp1,tmp2,ixi^l,ixo^l,w1,x,fc)
353 
354  w2(ixo^s,e_)=qcmu*w1(ixo^s,e_)+qcnu*w2(ixo^s,e_)+(1.d0-qcmu-qcnu)*w(ixo^s,e_)&
355  +qcmut*qdt*tmp(ixo^s)+qcnut*wold(ixo^s,e_)
356 
357  if (fix_conserve_at_step) then
358  fc=qcmut*qdt*fc
359  call store_flux(igrid,fc,1,ndim,1)
360  end if
361 
362  end subroutine evolve_stepj
363 
364  subroutine evolve_step1(igrid,qcmut,qdt,ixI^L,ixO^L,w1,w,x,wold)
366  use mod_fix_conserve
368 
369  integer, intent(in) :: igrid,ixI^L,ixO^L
370  double precision, intent(in) :: qcmut, qdt, w(ixi^s,1:nw), x(ixi^s,1:ndim)
371  double precision, intent(out) ::w1(ixi^s,1:nw),wold(ixi^s,1:nw)
372 
373  double precision :: fC(ixi^s,1:1,1:ndim)
374  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
375  integer :: lowindex(ndim), ix^D
376 
377  call phys_get_heatconduct(tmp,tmp1,tmp2,ixi^l,ixo^l,w,x,fc)
378 
379  wold(ixo^s,e_)=qdt*tmp(ixo^s)
380  ! update internal energy
381  tmp1(ixo^s) = tmp1(ixo^s) + qcmut*wold(ixo^s,e_)
382 
383  ! ensure you never trigger negative pressure
384  ! hence code up energy change with respect to kinetic and magnetic
385  ! part(nonthermal)
386  if(small_values_method=='error') then
387  if(any(tmp1(ixo^s)<small_e).and. .not.crash) then
388  lowindex=minloc(tmp1(ixo^s))
389  ^d&lowindex(^d)=lowindex(^d)+ixomin^d-1;
390  write(*,*)'too small internal energy = ',minval(tmp1(ixo^s)),'at x=',&
391  x(^d&lowindex(^d),1:ndim),lowindex,' with limit=',small_e,&
392  ' on time=',global_time,' step=',it, 'where w(1:nwflux)=',w(^d&lowindex(^d),1:nwflux)
393  crash=.true.
394  else
395  if(solve_internal_e) then
396  w1(ixo^s,e_) = tmp1(ixo^s)
397  else
398  w1(ixo^s,e_) = tmp2(ixo^s)+tmp1(ixo^s)
399  end if
400  end if
401  else
402  where(tmp1(ixo^s)<small_e)
403  tmp1(ixo^s)=small_e
404  endwhere
405  if(solve_internal_e) then
406  w1(ixo^s,e_) = tmp1(ixo^s)
407  else
408  w1(ixo^s,e_) = tmp2(ixo^s)+tmp1(ixo^s)
409  end if
410  end if
411 
412  if (fix_conserve_at_step) then
413  fc=qcmut*qdt*fc
414  call store_flux(igrid,fc,1,ndim,1)
415  end if
416 
417  end subroutine evolve_step1
418 
419  !> anisotropic thermal conduction with slope limited symmetric scheme
420  !> Sharma 2007 Journal of Computational Physics 227, 123
421  subroutine mhd_get_heatconduct(qd,tmp1,tmp2,ixI^L,ixO^L,w,x,qvec)
424 
425  integer, intent(in) :: ixI^L, ixO^L
426  double precision, intent(in) :: x(ixi^s,1:ndim), w(ixi^s,1:nw)
427  !! qd store the heat conduction energy changing rate
428  double precision, intent(out) :: qd(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
429  double precision, dimension(ixI^S,1:ndim) :: qvec
430 
431  double precision, dimension(ixI^S,1:ndir) :: mf,Bc,Bcf
432  double precision, dimension(ixI^S,1:ndim) :: gradT
433  double precision, dimension(ixI^S) :: Te,ka,kaf,ke,kef,qdd,qe,Binv,minq,maxq,Bnorm
434  double precision :: alpha,dxinv(ndim)
435  integer, dimension(ndim) :: lowindex
436  integer :: idims,idir,ix^D,ix^L,ixC^L,ixA^L,ixB^L
437 
438  ! coefficient of limiting on normal component
439  if(ndim<3) then
440  alpha=0.75d0
441  else
442  alpha=0.85d0
443  end if
444  ix^l=ixo^l^ladd1;
445 
446  dxinv=1.d0/dxlevel
447 
448  ! tmp1 store internal energy
449  if(solve_internal_e) then
450  tmp1(ixi^s)=w(ixi^s,e_)
451  else
452  ! tmp2 store kinetic+magnetic energy before addition of heat conduction source
453  tmp2(ixi^s) = 0.5d0 * (sum(w(ixi^s,iw_mom(:))**2,dim=ndim+1)/w(ixi^s,rho_) + &
454  sum(w(ixi^s,iw_mag(:))**2,dim=ndim+1))
455  tmp1(ixi^s)=w(ixi^s,e_)-tmp2(ixi^s)
456  end if
457 
458  ! Clip off negative pressure if small_pressure is set
459  if(small_values_method=='error') then
460  if (any(tmp1(ixi^s)<small_e) .and. .not.crash) then
461  lowindex=minloc(tmp1(ixi^s))
462  ^d&lowindex(^d)=lowindex(^d)+iximin^d-1;
463  write(*,*)'too low internal energy = ',minval(tmp1(ixi^s)),' at x=',&
464  x(^d&lowindex(^d),1:ndim),lowindex,' with limit=',small_e,' on time=',global_time, ' it=',it
465  write(*,*) 'w',w(^d&lowindex(^d),:)
466  crash=.true.
467  end if
468  else
469  {do ix^db=iximin^db,iximax^db\}
470  if(tmp1(ix^d)<small_e) then
471  tmp1(ix^d)=small_e
472  end if
473  {end do\}
474  end if
475  ! compute the temperature
476  te(ixi^s)=tmp1(ixi^s)*(tc_gamma-one)/w(ixi^s,rho_)
477  ! B vector
478  if(b0field) then
479  mf(ixi^s,:)=w(ixi^s,iw_mag(:))+block%B0(ixi^s,:,0)
480  else
481  mf(ixi^s,:)=w(ixi^s,iw_mag(:));
482  end if
483  ! |B|
484  binv(ix^s)=dsqrt(sum(mf(ix^s,:)**2,dim=ndim+1))
485  where(binv(ix^s)/=0.d0)
486  binv(ix^s)=1.d0/binv(ix^s)
487  elsewhere
488  binv(ix^s)=bigdouble
489  end where
490  ! b unit vector: magnetic field direction vector
491  do idims=1,ndim
492  mf(ix^s,idims)=mf(ix^s,idims)*binv(ix^s)
493  end do
494  ! ixC is cell-corner index
495  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
496  ! b unit vector at cell corner
497  bc=0.d0
498  {do ix^db=0,1\}
499  ixamin^d=ixcmin^d+ix^d;
500  ixamax^d=ixcmax^d+ix^d;
501  bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)+mf(ixa^s,1:ndim)
502  {end do\}
503  bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)*0.5d0**ndim
504  if(tc_constant) then
505  if(tc_perpendicular) then
506  ka(ixc^s)=tc_k_para-tc_k_perp
507  ke(ixc^s)=tc_k_perp
508  else
509  ka(ixc^s)=tc_k_para
510  end if
511  else
512  ! conductivity at cell center
513  minq(ix^s)=tc_k_para*te(ix^s)**2.5d0
514  ka=0.d0
515  {do ix^db=0,1\}
516  ixbmin^d=ixcmin^d+ix^d;
517  ixbmax^d=ixcmax^d+ix^d;
518  ka(ixc^s)=ka(ixc^s)+minq(ixb^s)
519  {end do\}
520  ! cell corner conductivity
521  ka(ixc^s)=0.5d0**ndim*ka(ixc^s)
522  ! compensate with perpendicular conductivity
523  if(tc_perpendicular) then
524  minq(ix^s)=tc_k_perp*w(ix^s,rho_)**2*binv(ix^s)**2/dsqrt(te(ix^s))
525  ke=0.d0
526  {do ix^db=0,1\}
527  ixbmin^d=ixcmin^d+ix^d;
528  ixbmax^d=ixcmax^d+ix^d;
529  ke(ixc^s)=ke(ixc^s)+minq(ixb^s)
530  {end do\}
531  ! cell corner conductivity: k_parallel-k_perpendicular
532  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
533  where(ke(ixc^s)<ka(ixc^s))
534  ka(ixc^s)=ka(ixc^s)-ke(ixc^s)
535  elsewhere
536  ka(ixc^s)=0.d0
537  ke(ixc^s)=ka(ixc^s)
538  end where
539  end if
540  end if
541  ! T gradient at cell faces
542  gradt=0.d0
543  do idims=1,ndim
544  ixbmin^d=ixmin^d;
545  ixbmax^d=ixmax^d-kr(idims,^d);
546  call gradientc(te,ixi^l,ixb^l,idims,minq)
547  gradt(ixb^s,idims)=minq(ixb^s)
548  end do
549  if(tc_slope_limiter=='no') then
550  ! calculate thermal conduction flux with symmetric scheme
551  do idims=1,ndim
552  qd=0.d0
553  {do ix^db=0,1 \}
554  if({ ix^d==0 .and. ^d==idims | .or.}) then
555  ixbmin^d=ixcmin^d+ix^d;
556  ixbmax^d=ixcmax^d+ix^d;
557  qd(ixc^s)=qd(ixc^s)+gradt(ixb^s,idims)
558  end if
559  {end do\}
560  ! temperature gradient at cell corner
561  qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
562  end do
563  ! b grad T at cell corner
564  qd(ixc^s)=sum(qvec(ixc^s,1:ndim)*bc(ixc^s,1:ndim),dim=ndim+1)
565  do idims=1,ndim
566  ! TC flux at cell corner
567  gradt(ixc^s,idims)=ka(ixc^s)*bc(ixc^s,idims)*qd(ixc^s)
568  if(tc_perpendicular) gradt(ixc^s,idims)=gradt(ixc^s,idims)+ke(ixc^s)*qvec(ixc^s,idims)
569  end do
570  ! TC flux at cell face
571  qvec=0.d0
572  do idims=1,ndim
573  ixb^l=ixo^l-kr(idims,^d);
574  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
575  {do ix^db=0,1 \}
576  if({ ix^d==0 .and. ^d==idims | .or.}) then
577  ixbmin^d=ixamin^d-ix^d;
578  ixbmax^d=ixamax^d-ix^d;
579  qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
580  end if
581  {end do\}
582  qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
583  if(tc_saturate) then
584  ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
585  ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
586  bcf=0.d0
587  {do ix^db=0,1 \}
588  if({ ix^d==0 .and. ^d==idims | .or.}) then
589  ixbmin^d=ixamin^d-ix^d;
590  ixbmax^d=ixamax^d-ix^d;
591  bcf(ixa^s,idims)=bcf(ixa^s,idims)+bc(ixb^s,idims)
592  end if
593  {end do\}
594  ! averaged b at face centers
595  bcf(ixa^s,idims)=bcf(ixa^s,idims)*0.5d0**(ndim-1)
596  ixb^l=ixa^l+kr(idims,^d);
597  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))
598  {do ix^db=ixamin^db,ixamax^db\}
599  if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
600  qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
601  end if
602  {end do\}
603  end if
604  end do
605  else
606  ! calculate thermal conduction flux with slope-limited symmetric scheme
607  qvec=0.d0
608  do idims=1,ndim
609  ixb^l=ixo^l-kr(idims,^d);
610  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
611  ! calculate normal of magnetic field
612  ixb^l=ixa^l+kr(idims,^d);
613  bnorm(ixa^s)=0.5d0*(mf(ixa^s,idims)+mf(ixb^s,idims))
614  bcf=0.d0
615  kaf=0.d0
616  kef=0.d0
617  {do ix^db=0,1 \}
618  if({ ix^d==0 .and. ^d==idims | .or.}) then
619  ixbmin^d=ixamin^d-ix^d;
620  ixbmax^d=ixamax^d-ix^d;
621  bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)+bc(ixb^s,1:ndim)
622  kaf(ixa^s)=kaf(ixa^s)+ka(ixb^s)
623  if(tc_perpendicular) kef(ixa^s)=kef(ixa^s)+ke(ixb^s)
624  end if
625  {end do\}
626  ! averaged b at face centers
627  bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)*0.5d0**(ndim-1)
628  ! averaged thermal conductivity at face centers
629  kaf(ixa^s)=kaf(ixa^s)*0.5d0**(ndim-1)
630  if(tc_perpendicular) kef(ixa^s)=kef(ixa^s)*0.5d0**(ndim-1)
631  ! limited normal component
632  minq(ixa^s)=min(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
633  maxq(ixa^s)=max(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
634  ! eq (19)
635  qdd=0.d0
636  {do ix^db=0,1 \}
637  if({ ix^d==0 .and. ^d==idims | .or.}) then
638  ixbmin^d=ixcmin^d+ix^d;
639  ixbmax^d=ixcmax^d+ix^d;
640  qdd(ixc^s)=qdd(ixc^s)+gradt(ixb^s,idims)
641  end if
642  {end do\}
643  ! temperature gradient at cell corner
644  qdd(ixc^s)=qdd(ixc^s)*0.5d0**(ndim-1)
645  ! eq (21)
646  qe=0.d0
647  {do ix^db=0,1 \}
648  qd(ixc^s)=qdd(ixc^s)
649  if({ ix^d==0 .and. ^d==idims | .or.}) then
650  ixbmin^d=ixamin^d-ix^d;
651  ixbmax^d=ixamax^d-ix^d;
652  where(qd(ixb^s)<=minq(ixa^s))
653  qd(ixb^s)=minq(ixa^s)
654  elsewhere(qd(ixb^s)>=maxq(ixa^s))
655  qd(ixb^s)=maxq(ixa^s)
656  end where
657  qvec(ixa^s,idims)=qvec(ixa^s,idims)+bc(ixb^s,idims)**2*qd(ixb^s)
658  if(tc_perpendicular) qe(ixa^s)=qe(ixa^s)+qd(ixb^s)
659  end if
660  {end do\}
661  qvec(ixa^s,idims)=kaf(ixa^s)*qvec(ixa^s,idims)*0.5d0**(ndim-1)
662  ! add normal flux from perpendicular conduction
663  if(tc_perpendicular) qvec(ixa^s,idims)=qvec(ixa^s,idims)+kef(ixa^s)*qe(ixa^s)*0.5d0**(ndim-1)
664  ! limited transverse component, eq (17)
665  ixbmin^d=ixamin^d;
666  ixbmax^d=ixamax^d+kr(idims,^d);
667  do idir=1,ndim
668  if(idir==idims) cycle
669  qd(ixi^s)=slope_limiter(gradt(ixi^s,idir),ixi^l,ixb^l,idir,-1)
670  qd(ixi^s)=slope_limiter(qd,ixi^l,ixa^l,idims,1)
671  qvec(ixa^s,idims)=qvec(ixa^s,idims)+kaf(ixa^s)*bnorm(ixa^s)*bcf(ixa^s,idir)*qd(ixa^s)
672  end do
673 
674  ! consider magnetic null point
675  !where(Binv(ixA^S)==0.d0)
676  ! qvec(ixA^S,idims)=tc_k_para*(0.5d0*(Te(ixA^S)+Te(ixB^S)))**2.5d0*gradT(ixA^S,idims)
677  !end where
678 
679  if(tc_saturate) then
680  ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
681  ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
682  ixb^l=ixa^l+kr(idims,^d);
683  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))
684  {do ix^db=ixamin^db,ixamax^db\}
685  if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
686  ! write(*,*) 'it',it,qvec(ix^D,idims),qd(ix^D),' TC saturated at ',&
687  ! x(ix^D,:),' rho',w(ix^D,rho_),' Te',Te(ix^D)
688  qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
689  end if
690  {end do\}
691  end if
692  end do
693  end if
694 
695  qd=0.d0
696  if(slab_uniform) then
697  do idims=1,ndim
698  qvec(ix^s,idims)=dxinv(idims)*qvec(ix^s,idims)
699  ixb^l=ixo^l-kr(idims,^d);
700  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
701  end do
702  else
703  do idims=1,ndim
704  qvec(ix^s,idims)=qvec(ix^s,idims)*block%surfaceC(ix^s,idims)
705  ixb^l=ixo^l-kr(idims,^d);
706  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
707  end do
708  qd(ixo^s)=qd(ixo^s)/block%dvolume(ixo^s)
709  end if
710 
711  end subroutine mhd_get_heatconduct
712 
713  !> Calculate gradient of a scalar q at cell interfaces in direction idir
714  subroutine gradientc(q,ixI^L,ixO^L,idir,gradq)
716  use mod_geometry
717 
718  integer, intent(in) :: ixI^L, ixO^L, idir
719  double precision, intent(in) :: q(ixi^s)
720  double precision, intent(inout) :: gradq(ixi^s)
721  integer :: jxO^L
722 
723  associate(x=>block%x)
724 
725  jxo^l=ixo^l+kr(idir,^d);
726  select case(coordinate)
727  case(cartesian)
728  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/dxlevel(idir)
729  case(cartesian_stretched)
730  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,idir)-x(ixo^s,idir))
731  case(spherical)
732  select case(idir)
733  case(1)
734  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,1)-x(ixo^s,1))
735  {^nooned
736  case(2)
737  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/( (x(jxo^s,2)-x(ixo^s,2))*x(ixo^s,1) )
738  }
739  {^ifthreed
740  case(3)
741  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)) )
742  }
743  end select
744  case(cylindrical)
745  if(idir==phi_) then
746  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/((x(jxo^s,phi_)-x(ixo^s,phi_))*x(ixo^s,r_))
747  else
748  gradq(ixo^s)=(q(jxo^s)-q(ixo^s))/(x(jxo^s,idir)-x(ixo^s,idir))
749  end if
750  case default
751  call mpistop('Unknown geometry')
752  end select
753 
754  end associate
755  end subroutine gradientc
756 
757  function slope_limiter(f,ixI^L,ixO^L,idims,pm) result(lf)
759  integer, intent(in) :: ixI^L, ixO^L, idims, pm
760  double precision, intent(in) :: f(ixi^s)
761  double precision :: lf(ixi^s)
762 
763  double precision :: signf(ixi^s)
764  integer :: ixB^L
765 
766  ixb^l=ixo^l+pm*kr(idims,^d);
767  signf(ixo^s)=sign(1.d0,f(ixo^s))
768  select case(tc_slope_limiter)
769  case('minmod')
770  ! minmod limiter
771  lf(ixo^s)=signf(ixo^s)*max(0.d0,min(abs(f(ixo^s)),signf(ixo^s)*f(ixb^s)))
772  case ('MC')
773  ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
774  lf(ixo^s)=two*signf(ixo^s)* &
775  max(zero,min(dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s),&
776  signf(ixo^s)*quarter*(f(ixb^s)+f(ixo^s))))
777  case ('superbee')
778  ! Roes superbee limiter (eq.3.51i)
779  lf(ixo^s)=signf(ixo^s)* &
780  max(zero,min(two*dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s)),&
781  min(dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s)))
782  case ('koren')
783  ! Barry Koren Right variant
784  lf(ixo^s)=signf(ixo^s)* &
785  max(zero,min(two*dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s),&
786  (two*f(ixb^s)*signf(ixo^s)+dabs(f(ixo^s)))*third))
787  case default
788  call mpistop("Unknown slope limiter for thermal conduction")
789  end select
790 
791  end function slope_limiter
792 
793  subroutine mhd_getdt_heatconduct(w,ixI^L,ixO^L,dtnew,dx^D,x)
794  !Check diffusion time limit dt < tc_dtpar*dx_i**2/((gamma-1)*tc_k_para_i/rho)
795  !where tc_k_para_i=tc_k_para*B_i**2/B**2
796  !and T=p/rho
798  use mod_physics
799 
800  integer, intent(in) :: ixI^L, ixO^L
801  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
802  ! note that depending on small_values_method=='error' etc, w values may change
803  ! through call to getpthermal
804  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
805 
806  double precision :: dxinv(1:ndim),mf(ixi^s,1:ndir)
807  double precision :: tmp2(ixi^s),tmp(ixi^s),Te(ixi^s),B2(ixi^s)
808  double precision :: dtdiff_tcond, dtdiff_tsat
809  integer :: idim,ix^D
810 
811  ^d&dxinv(^d)=one/dx^d;
812 
813  call phys_get_pthermal(w,x,ixi^l,ixo^l,tmp)
814 
815  !temperature
816  te(ixo^s)=tmp(ixo^s)/w(ixo^s,rho_)
817  !tc_k_para_i
818  if(tc_constant) then
819  tmp(ixo^s)=tc_k_para
820  else
821  tmp(ixo^s)=tc_k_para*dsqrt(te(ixo^s)**5)/w(ixo^s,rho_)
822  end if
823 
824  ! B
825  if(b0field) then
826  mf(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
827  else
828  mf(ixo^s,:)=w(ixo^s,iw_mag(:))
829  end if
830  ! B^-2
831  b2(ixo^s)=sum(mf(ixo^s,:)**2,dim=ndim+1)
832  ! B_i**2/B**2
833  where(b2(ixo^s)/=0.d0)
834  ^d&mf(ixo^s,^d)=mf(ixo^s,^d)**2/b2(ixo^s);
835  elsewhere
836  ^d&mf(ixo^s,^d)=1.d0;
837  end where
838 
839  if(tc_saturate) b2(ixo^s)=22.d0*dsqrt(te(ixo^s))
840 
841  do idim=1,ndim
842  tmp2(ixo^s)=tmp(ixo^s)*mf(ixo^s,idim)
843  if(tc_saturate) then
844  where(tmp2(ixo^s)>b2(ixo^s))
845  tmp2(ixo^s)=b2(ixo^s)
846  end where
847  end if
848  ! dt< tc_dtpar * dx_idim**2/((gamma-1)*tc_k_para_i/rho*B_i**2/B**2)
849  dtdiff_tcond=tc_dtpar/(tc_gamma-1.d0)/maxval(tmp2(ixo^s)*dxinv(idim)**2)
850  ! limit the time step
851  dtnew=min(dtnew,dtdiff_tcond)
852  end do
853  dtnew=dtnew/dble(ndim)
854 
855  end subroutine mhd_getdt_heatconduct
856 
857  subroutine hd_get_heatconduct(qd,tmp1,tmp2,ixI^L,ixO^L,w,x,qvec)
860 
861  integer, intent(in) :: ixI^L, ixO^L
862  double precision, intent(in) :: x(ixi^s,1:ndim), w(ixi^s,1:nw)
863  !! tmp store the heat conduction energy changing rate
864  double precision, intent(out) :: qd(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
865  double precision :: qvec(ixi^s,1:ndim),gradT(ixi^s,1:ndim),Te(ixi^s),ke(ixi^s)
866  double precision :: dxinv(ndim)
867  integer, dimension(ndim) :: lowindex
868  integer :: idims,ix^D,ix^L,ixC^L,ixA^L,ixB^L
869 
870  ix^l=ixo^l^ladd1;
871  ! ixC is cell-corner index
872  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
873 
874  dxinv=1.d0/dxlevel
875 
876  ! tmp1 store internal energy
877  if(solve_internal_e) then
878  tmp1(ixi^s)=w(ixi^s,e_)
879  else
880  ! tmp2 store kinetic energy before addition of heat conduction source
881  tmp2(ixi^s) = 0.5d0*sum(w(ixi^s,iw_mom(:))**2,dim=ndim+1)/w(ixi^s,rho_)
882  tmp1(ixi^s)=w(ixi^s,e_)-tmp2(ixi^s)
883  end if
884 
885  ! Clip off negative pressure if small_pressure is set
886  if(small_values_method=='error') then
887  if (any(tmp1(ixi^s)<small_e) .and. .not.crash) then
888  lowindex=minloc(tmp1(ixi^s))
889  ^d&lowindex(^d)=lowindex(^d)+iximin^d-1;
890  write(*,*)'too low internal energy = ',minval(tmp1(ixi^s)),' at x=',&
891  x(^d&lowindex(^d),1:ndim),lowindex,' with limit=',small_e,' on time=',global_time, ' it=',it
892  write(*,*) 'w',w(^d&lowindex(^d),:)
893  crash=.true.
894  end if
895  else
896  {do ix^db=iximin^db,iximax^db\}
897  if(tmp1(ix^d)<small_e) then
898  tmp1(ix^d)=small_e
899  end if
900  {end do\}
901  end if
902 
903  ! compute temperature before source addition
904  te(ixi^s)=tmp1(ixi^s)/w(ixi^s,rho_)*(tc_gamma-1.d0)
905 
906  ! T gradient at cell faces
907  gradt=0.d0
908  do idims=1,ndim
909  ixbmin^d=ixmin^d;
910  ixbmax^d=ixmax^d-kr(idims,^d);
911  call gradientc(te,ixi^l,ixb^l,idims,ke)
912  gradt(ixb^s,idims)=ke(ixb^s)
913  end do
914  ! calculate thermal conduction flux with symmetric scheme
915  do idims=1,ndim
916  qd=0.d0
917  {do ix^db=0,1 \}
918  if({ ix^d==0 .and. ^d==idims | .or.}) then
919  ixbmin^d=ixcmin^d+ix^d;
920  ixbmax^d=ixcmax^d+ix^d;
921  qd(ixc^s)=qd(ixc^s)+gradt(ixb^s,idims)
922  end if
923  {end do\}
924  ! temperature gradient at cell corner
925  qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
926  end do
927  ! conductivity at cell center
928  qd(ix^s)=tc_k_para*dsqrt(te(ix^s))**5
929  ke=0.d0
930  {do ix^db=0,1\}
931  ixbmin^d=ixcmin^d+ix^d;
932  ixbmax^d=ixcmax^d+ix^d;
933  ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
934  {end do\}
935  ! cell corner conductivity
936  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
937  ! cell corner conduction flux
938  do idims=1,ndim
939  gradt(ixc^s,idims)=ke(ixc^s)*qvec(ixc^s,idims)
940  end do
941 
942  if(tc_saturate) then
943  ! consider saturation with unsigned saturated TC flux = 5 phi rho c**3
944  ! saturation flux at cell center
945  qd(ix^s)=5.d0*w(ix^s,rho_)*dsqrt(te(ix^s)**3)
946  ke=0.d0
947  {do ix^db=0,1\}
948  ixbmin^d=ixcmin^d+ix^d;
949  ixbmax^d=ixcmax^d+ix^d;
950  ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
951  {end do\}
952  ! cell corner saturation flux
953  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
954  ! magnitude of cell corner conduction flux
955  qd(ixc^s)=norm2(gradt(ixc^s,:),dim=ndim+1)
956  {do ix^db=ixcmin^db,ixcmax^db\}
957  if(qd(ix^d)>ke(ix^d)) then
958  ke(ix^d)=ke(ix^d)/qd(ix^d)
959  do idims=1,ndim
960  gradt(ix^d,idims)=ke(ix^d)*gradt(ix^d,idims)
961  end do
962  end if
963  {end do\}
964  end if
965 
966  ! conductionflux at cell face
967  qvec=0.d0
968  do idims=1,ndim
969  ixb^l=ixo^l-kr(idims,^d);
970  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
971  {do ix^db=0,1 \}
972  if({ ix^d==0 .and. ^d==idims | .or.}) then
973  ixbmin^d=ixamin^d-ix^d;
974  ixbmax^d=ixamax^d-ix^d;
975  qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
976  end if
977  {end do\}
978  qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
979  end do
980 
981  qd=0.d0
982  if(slab_uniform) then
983  do idims=1,ndim
984  qvec(ix^s,idims)=dxinv(idims)*qvec(ix^s,idims)
985  ixb^l=ixo^l-kr(idims,^d);
986  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
987  end do
988  else
989  do idims=1,ndim
990  qvec(ix^s,idims)=qvec(ix^s,idims)*block%surfaceC(ix^s,idims)
991  ixb^l=ixo^l-kr(idims,^d);
992  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
993  end do
994  qd(ixo^s)=qd(ixo^s)/block%dvolume(ixo^s)
995  end if
996 
997  end subroutine hd_get_heatconduct
998 
999  subroutine hd_getdt_heatconduct(w,ixI^L,ixO^L,dtnew,dx^D,x)
1000  ! Check diffusion time limit dt < tc_dtpar * dx_i**2 / ((gamma-1)*tc_k_para_i/rho)
1002  use mod_physics
1003 
1004  integer, intent(in) :: ixI^L, ixO^L
1005  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
1006  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
1007 
1008  double precision :: dxinv(1:ndim), tmp(ixi^s), Te(ixi^s)
1009  double precision :: dtdiff_tcond,dtdiff_tsat
1010  integer :: idim,ix^D
1011 
1012  ^d&dxinv(^d)=one/dx^d;
1013 
1014  call phys_get_pthermal(w,x,ixi^l,ixo^l,tmp)
1015 
1016  te(ixo^s)=tmp(ixo^s)/w(ixo^s,rho_)
1017  tmp(ixo^s)=(tc_gamma-one)*tc_k_para*dsqrt((te(ixo^s))**5)/w(ixo^s,rho_)
1018 
1019  do idim=1,ndim
1020  ! dt< tc_dtpar * dx_idim**2/((gamma-1)*tc_k_para_idim/rho)
1021  dtdiff_tcond=tc_dtpar/maxval(tmp(ixo^s)*dxinv(idim)**2)
1022  if(tc_saturate) then
1023  ! dt< tc_dtpar* dx_idim**2/((gamma-1)*sqrt(Te)*5*phi)
1024  dtdiff_tsat=tc_dtpar/maxval((tc_gamma-1.d0)*dsqrt(te(ixo^s))*&
1025  5.d0*dxinv(idim)**2)
1026  ! choose the slower flux (bigger time scale) between classic and saturated
1027  dtdiff_tcond=max(dtdiff_tcond,dtdiff_tsat)
1028  end if
1029  ! limit the time step
1030  dtnew=min(dtnew,dtdiff_tcond)
1031  end do
1032  dtnew=dtnew/dble(ndim)
1033 
1034  end subroutine hd_getdt_heatconduct
1035 
1036 end module mod_thermal_conduction
integer, parameter cylindrical
Definition: mod_geometry.t:9
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_f
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p1
integer, parameter cartesian_stretched
Definition: mod_geometry.t:8
subroutine evolve_stepj(igrid, qcmu, qcmut, qcnu, qcnut, qdt, ixIL, ixOL, w1, w2, w, x, wold)
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_f
update ghost cells of all blocks including physical boundaries
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
integer, dimension(:^d &,:^d &), pointer type_recv_srl
double precision function, dimension(ixi^s) slope_limiter(f, ixIL, ixOL, idims, pm)
subroutine create_bc_mpi_datatype(nwstart, nwbc)
double precision unit_length
Physical scaling factor for length.
integer, dimension(:^d &,:^d &), pointer type_send_r
double precision unit_density
Physical scaling factor for density.
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:14
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, dimension(:), allocatable type_gradient_limiter
Type of slope limiter used for computing gradients or divergences, when typegrad or typediv are set t...
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
logical fix_conserve_at_step
Whether to conserve fluxes at the current partial step.
integer, public s
Number of sub-steps of supertime stepping.
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
integer, parameter plevel_
integer ndir
Number of spatial dimensions (components) for vector variables.
integer coordinate
Definition: mod_geometry.t:6
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
procedure(get_heatconduct), pointer phys_get_heatconduct
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_f
Module for flux conservation near refinement boundaries.
integer, dimension(:^d &,:^d &), pointer type_recv_r
integer, dimension(:^d &,:^d &), pointer type_recv_p
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_f
logical b0field
split magnetic field as background B0 field
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
double precision dt_tc
Time step of thermal conduction.
subroutine, public sendflux(idimLIM)
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p1
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p1
subroutine gradientc(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q at cell interfaces in direction idir.
subroutine hd_get_heatconduct(qd, tmp1, tmp2, ixIL, ixOL, w, x, qvec)
subroutine mhd_getdt_heatconduct(w, ixIL, ixOL, dtnew, dxD, x)
double precision small_pressure
procedure(getdt_heatconduct), pointer phys_getdt_heatconduct
integer tc_ncycles
Maximal substeps of TC within one fluid time step to limit fluid time step.
logical stagger_grid
True for using stagger grid.
integer, parameter cartesian
Definition: mod_geometry.t:7
Thermal conduction for HD and MHD.
double precision unit_velocity
Physical scaling factor for velocity.
logical solve_internal_e
Solve polytropic process instead of solving total energy.
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, parameter unitpar
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
double precision, public tc_k_para
Coefficient of thermal conductivity (parallel to magnetic field)
double precision global_time
The global simulation time.
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, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_f
subroutine tc_params_read(files)
Read this module"s parameters from a file.
double precision, dimension(ndim) dxlevel
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
logical crash
Save a snapshot before crash a run met unphysical values.
integer, dimension(:^d &,:^d &), pointer type_send_p
subroutine evolve_step1(igrid, qcmut, qdt, ixIL, ixOL, w1, w, x, wold)
subroutine hd_getdt_heatconduct(w, ixIL, ixOL, dtnew, dxD, x)
subroutine, public recvflux(idimLIM)
procedure(thermal_conduction), pointer phys_thermal_conduction
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
integer, dimension(:,:), allocatable node
subroutine thermal_conduction_init(phys_gamma)
Initialize the module.
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
integer, parameter spherical
Definition: mod_geometry.t:10
integer, dimension(:^d &,:^d &), pointer type_send_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p1
integer it
Number of time steps taken.
character(len=20), public small_values_method
How to handle small values.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_numberdensity
Physical scaling factor for number density.
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition: mod_physics.t:59
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p1
subroutine mhd_get_heatconduct(qd, tmp1, tmp2, ixIL, ixOL, w, x, qvec)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p1