MPI-AMRVAC  2.0
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
194 
195  ! point bc mpi datatype to partial type for thermalconduction
202  ! create bc mpi datatype for ghostcells update
203  if(first) then
204  call create_bc_mpi_datatype(e_-1,1)
205  first=.false.
206  end if
207 
208  call init_comm_fix_conserve(1,ndim,1)
210 
211  do iigrid=1,igridstail; igrid=igrids(iigrid);
212  if(.not. allocated(ps2(igrid)%w)) allocate(ps2(igrid)%w(ixg^t,1:nw))
213  if(.not. allocated(ps3(igrid)%w)) allocate(ps3(igrid)%w(ixg^t,1:nw))
214  ps1(igrid)%w=ps(igrid)%w
215  ps2(igrid)%w=ps(igrid)%w
216  ps3(igrid)%w=ps(igrid)%w
217  end do
218 
219  allocate(bj(0:s))
220  bj(0)=1.d0/3.d0
221  bj(1)=bj(0)
222  if(s>1) then
223  omega1=4.d0/dble(s**2+s-2)
224  cmut=omega1/3.d0
225  else
226  omega1=0.d0
227  cmut=1.d0
228  endif
229 
230  !$OMP PARALLEL DO PRIVATE(igrid)
231  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
232  block=>ps(igrid)
235  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
236  call evolve_step1(igrid,cmut,dt_tc,ixg^ll,ixm^ll,ps1(igrid)%w,ps(igrid)%w,&
237  ps(igrid)%x,ps3(igrid)%w)
238  end do
239  !$OMP END PARALLEL DO
240  ! fix conservation of AMR grid by replacing flux from finer neighbors
241  if (fix_conserve_at_step) then
242  call recvflux(1,ndim)
243  call sendflux(1,ndim)
244  call fix_conserve(ps1,1,ndim,e_,1)
245  end if
246  bcphys=.false.
247  call getbc(global_time,0.d0,ps1,e_-1,1)
248  if(s==1) then
249  do iigrid=1,igridstail; igrid=igrids(iigrid);
250  ps(igrid)%w(ixg^t,e_)=ps1(igrid)%w(ixg^t,e_)
251  end do
252  ! point bc mpi data type back to full type for (M)HD
259  bcphys=.true.
260  deallocate(bj)
261  return
262  endif
263  evenstep=.true.
264  do j=2,s
265  bj(j)=dble(j**2+j-2)/dble(2*j*(j+1))
266  cmu=dble(2*j-1)/dble(j)*bj(j)/bj(j-1)
267  cmut=omega1*cmu
268  cnu=dble(1-j)/dble(j)*bj(j)/bj(j-2)
269  cnut=(bj(j-1)-1.d0)*cmut
270  if(evenstep) then
271  !$OMP PARALLEL DO PRIVATE(igrid)
272  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
273  block=>ps(igrid)
276  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
277  call evolve_stepj(igrid,cmu,cmut,cnu,cnut,dt_tc,ixg^ll,ixm^ll,ps1(igrid)%w,&
278  ps2(igrid)%w,ps(igrid)%w,ps(igrid)%x,ps3(igrid)%w)
279  end do
280  !$OMP END PARALLEL DO
281  ! fix conservation of AMR grid by replacing flux from finer neighbors
282  if (fix_conserve_at_step) then
283  call recvflux(1,ndim)
284  call sendflux(1,ndim)
285  call fix_conserve(ps2,1,ndim,e_,1)
286  end if
287  call getbc(global_time,0.d0,ps2,e_-1,1)
288  evenstep=.false.
289  else
290  !$OMP PARALLEL DO PRIVATE(igrid)
291  do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
292  block=>ps(igrid)
295  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
296  call evolve_stepj(igrid,cmu,cmut,cnu,cnut,dt_tc,ixg^ll,ixm^ll,ps2(igrid)%w,&
297  ps1(igrid)%w,ps(igrid)%w,ps(igrid)%x,ps3(igrid)%w)
298  end do
299  !$OMP END PARALLEL DO
300  ! fix conservation of AMR grid by replacing flux from finer neighbors
301  if (fix_conserve_at_step) then
302  call recvflux(1,ndim)
303  call sendflux(1,ndim)
304  call fix_conserve(ps1,1,ndim,e_,1)
305  end if
306  call getbc(global_time,0.d0,ps1,e_-1,1)
307  evenstep=.true.
308  end if
309  end do
310  if(evenstep) then
311  do iigrid=1,igridstail; igrid=igrids(iigrid);
312  ps(igrid)%w(ixg^t,e_)=ps1(igrid)%w(ixg^t,e_)
313  end do
314  else
315  do iigrid=1,igridstail; igrid=igrids(iigrid);
316  ps(igrid)%w(ixg^t,e_)=ps2(igrid)%w(ixg^t,e_)
317  end do
318  end if
319  deallocate(bj)
320  ! point bc mpi data type back to full type for (M)HD
327  bcphys=.true.
328 
329  end subroutine do_thermal_conduction
330 
331  subroutine evolve_stepj(igrid,qcmu,qcmut,qcnu,qcnut,qdt,ixI^L,ixO^L,w1,w2,w,x,wold)
333  use mod_fix_conserve
334 
335  integer, intent(in) :: igrid,ixI^L,ixO^L
336  double precision, intent(in) :: qcmu,qcmut,qcnu,qcnut,qdt
337  double precision, intent(in) :: w1(ixi^s,1:nw),w(ixi^s,1:nw),wold(ixi^s,1:nw)
338  double precision, intent(in) :: x(ixi^s,1:ndim)
339  double precision, intent(inout) :: w2(ixi^s,1:nw)
340 
341  double precision :: fC(ixi^s,1:1,1:ndim)
342  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
343 
344  call phys_get_heatconduct(tmp,tmp1,tmp2,ixi^l,ixo^l,w1,x,fc)
345 
346  w2(ixo^s,e_)=qcmu*w1(ixo^s,e_)+qcnu*w2(ixo^s,e_)+(1.d0-qcmu-qcnu)*w(ixo^s,e_)&
347  +qcmut*qdt*tmp(ixo^s)+qcnut*wold(ixo^s,e_)
348 
349  if (fix_conserve_at_step) then
350  fc=qcmut*qdt*fc
351  call storeflux(igrid,fc,1,ndim,1)
352  end if
353 
354  end subroutine evolve_stepj
355 
356  subroutine evolve_step1(igrid,qcmut,qdt,ixI^L,ixO^L,w1,w,x,wold)
358  use mod_fix_conserve
360 
361  integer, intent(in) :: igrid,ixI^L,ixO^L
362  double precision, intent(in) :: qcmut, qdt, w(ixi^s,1:nw), x(ixi^s,1:ndim)
363  double precision, intent(out) ::w1(ixi^s,1:nw),wold(ixi^s,1:nw)
364 
365  double precision :: fC(ixi^s,1:1,1:ndim)
366  double precision :: tmp(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
367  integer :: lowindex(ndim), ix^D
368 
369  call phys_get_heatconduct(tmp,tmp1,tmp2,ixi^l,ixo^l,w,x,fc)
370 
371  wold(ixo^s,e_)=qdt*tmp(ixo^s)
372  ! update internal energy
373  tmp1(ixo^s) = tmp1(ixo^s) + qcmut*wold(ixo^s,e_)
374 
375  ! ensure you never trigger negative pressure
376  ! hence code up energy change with respect to kinetic and magnetic
377  ! part(nonthermal)
378  if(small_values_method=='error') then
379  if(any(tmp1(ixo^s)<small_e).and. .not.crash) then
380  lowindex=minloc(tmp1(ixo^s))
381  ^d&lowindex(^d)=lowindex(^d)+ixomin^d-1;
382  write(*,*)'too small internal energy = ',minval(tmp1(ixo^s)),'at x=',&
383  x(^d&lowindex(^d),1:ndim),lowindex,' with limit=',small_e,&
384  ' on time=',global_time,' step=',it, 'where w(1:nwflux)=',w(^d&lowindex(^d),1:nwflux)
385  crash=.true.
386  else
387  if(solve_internal_e) then
388  w1(ixo^s,e_) = tmp1(ixo^s)
389  else
390  w1(ixo^s,e_) = tmp2(ixo^s)+tmp1(ixo^s)
391  end if
392  end if
393  else
394  where(tmp1(ixo^s)<small_e)
395  tmp1(ixo^s)=small_e
396  endwhere
397  if(solve_internal_e) then
398  w1(ixo^s,e_) = tmp1(ixo^s)
399  else
400  w1(ixo^s,e_) = tmp2(ixo^s)+tmp1(ixo^s)
401  end if
402  end if
403 
404  if (fix_conserve_at_step) then
405  fc=qcmut*qdt*fc
406  call storeflux(igrid,fc,1,ndim,1)
407  end if
408 
409  end subroutine evolve_step1
410 
411  !> anisotropic thermal conduction with slope limited symmetric scheme
412  !> Sharma 2007 Journal of Computational Physics 227, 123
413  subroutine mhd_get_heatconduct(qd,tmp1,tmp2,ixI^L,ixO^L,w,x,qvec)
416 
417  integer, intent(in) :: ixI^L, ixO^L
418  double precision, intent(in) :: x(ixi^s,1:ndim), w(ixi^s,1:nw)
419  !! qd store the heat conduction energy changing rate
420  double precision, intent(out) :: qd(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
421  double precision, dimension(ixI^S,1:ndim) :: qvec
422 
423  double precision, dimension(ixI^S,1:ndir) :: mf,Bc,Bcf
424  double precision, dimension(ixI^S,1:ndim) :: gradT
425  double precision, dimension(ixI^S) :: Te,ka,kaf,ke,kef,qdd,qe,Binv,minq,maxq,Bnorm
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  ! tmp1 store internal energy
441  if(solve_internal_e) then
442  tmp1(ixi^s)=w(ixi^s,e_)
443  else
444  ! tmp2 store kinetic+magnetic energy before addition of heat conduction source
445  tmp2(ixi^s) = 0.5d0 * (sum(w(ixi^s,iw_mom(:))**2,dim=ndim+1)/w(ixi^s,rho_) + &
446  sum(w(ixi^s,iw_mag(:))**2,dim=ndim+1))
447  tmp1(ixi^s)=w(ixi^s,e_)-tmp2(ixi^s)
448  end if
449 
450  ! Clip off negative pressure if small_pressure is set
451  if(small_values_method=='error') then
452  if (any(tmp1(ixi^s)<small_e) .and. .not.crash) then
453  lowindex=minloc(tmp1(ixi^s))
454  ^d&lowindex(^d)=lowindex(^d)+iximin^d-1;
455  write(*,*)'too low internal energy = ',minval(tmp1(ixi^s)),' at x=',&
456  x(^d&lowindex(^d),1:ndim),lowindex,' with limit=',small_e,' on time=',global_time, ' it=',it
457  write(*,*) 'w',w(^d&lowindex(^d),:)
458  crash=.true.
459  end if
460  else
461  {do ix^db=iximin^db,iximax^db\}
462  if(tmp1(ix^d)<small_e) then
463  tmp1(ix^d)=small_e
464  end if
465  {end do\}
466  end if
467  ! compute the temperature
468  te(ixi^s)=tmp1(ixi^s)*(tc_gamma-one)/w(ixi^s,rho_)
469  ! B vector
470  if(b0field) then
471  mf(ixi^s,:)=w(ixi^s,iw_mag(:))+block%B0(ixi^s,:,0)
472  else
473  mf(ixi^s,:)=w(ixi^s,iw_mag(:));
474  end if
475  ! |B|
476  binv(ix^s)=dsqrt(sum(mf(ix^s,:)**2,dim=ndim+1))
477  where(binv(ix^s)/=0.d0)
478  binv(ix^s)=1.d0/binv(ix^s)
479  elsewhere
480  binv(ix^s)=bigdouble
481  end where
482  ! b unit vector: magnetic field direction vector
483  do idims=1,ndim
484  mf(ix^s,idims)=mf(ix^s,idims)*binv(ix^s)
485  end do
486  ! ixC is cell-corner index
487  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
488  ! b unit vector at cell corner
489  bc=0.d0
490  {do ix^db=0,1\}
491  ixamin^d=ixcmin^d+ix^d;
492  ixamax^d=ixcmax^d+ix^d;
493  bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)+mf(ixa^s,1:ndim)
494  {end do\}
495  bc(ixc^s,1:ndim)=bc(ixc^s,1:ndim)*0.5d0**ndim
496  if(tc_constant) then
497  if(tc_perpendicular) then
498  ka(ixc^s)=tc_k_para-tc_k_perp
499  ke(ixc^s)=tc_k_perp
500  else
501  ka(ixc^s)=tc_k_para
502  end if
503  else
504  ! conductivity at cell center
505  minq(ix^s)=tc_k_para*te(ix^s)**2.5d0
506  ka=0.d0
507  {do ix^db=0,1\}
508  ixbmin^d=ixcmin^d+ix^d;
509  ixbmax^d=ixcmax^d+ix^d;
510  ka(ixc^s)=ka(ixc^s)+minq(ixb^s)
511  {end do\}
512  ! cell corner conductivity
513  ka(ixc^s)=0.5d0**ndim*ka(ixc^s)
514  ! compensate with perpendicular conductivity
515  if(tc_perpendicular) then
516  minq(ix^s)=tc_k_perp*w(ix^s,rho_)**2*binv(ix^s)**2/dsqrt(te(ix^s))
517  ke=0.d0
518  {do ix^db=0,1\}
519  ixbmin^d=ixcmin^d+ix^d;
520  ixbmax^d=ixcmax^d+ix^d;
521  ke(ixc^s)=ke(ixc^s)+minq(ixb^s)
522  {end do\}
523  ! cell corner conductivity: k_parallel-k_perpendicular
524  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
525  where(ke(ixc^s)<ka(ixc^s))
526  ka(ixc^s)=ka(ixc^s)-ke(ixc^s)
527  elsewhere
528  ka(ixc^s)=0.d0
529  ke(ixc^s)=ka(ixc^s)
530  end where
531  end if
532  end if
533  ! T gradient at cell faces
534  gradt=0.d0
535  do idims=1,ndim
536  ixbmin^d=ixmin^d;
537  ixbmax^d=ixmax^d-kr(idims,^d);
538  ixa^l=ixb^l+kr(idims,^d);
539  gradt(ixb^s,idims)=(te(ixa^s)-te(ixb^s))*dxinv(idims)
540  end do
541  if(tc_slope_limiter=='no') then
542  ! calculate thermal conduction flux with symmetric scheme
543  do idims=1,ndim
544  qd=0.d0
545  {do ix^db=0,1 \}
546  if({ ix^d==0 .and. ^d==idims | .or.}) then
547  ixbmin^d=ixcmin^d+ix^d;
548  ixbmax^d=ixcmax^d+ix^d;
549  qd(ixc^s)=qd(ixc^s)+gradt(ixb^s,idims)
550  end if
551  {end do\}
552  ! temperature gradient at cell corner
553  qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
554  end do
555  ! b grad T at cell corner
556  qd(ixc^s)=sum(qvec(ixc^s,1:ndim)*bc(ixc^s,1:ndim),dim=ndim+1)
557  do idims=1,ndim
558  ! TC flux at cell corner
559  gradt(ixc^s,idims)=ka(ixc^s)*bc(ixc^s,idims)*qd(ixc^s)
560  if(tc_perpendicular) gradt(ixc^s,idims)=gradt(ixc^s,idims)+ke(ixc^s)*qvec(ixc^s,idims)
561  end do
562  ! TC flux at cell face
563  qvec=0.d0
564  do idims=1,ndim
565  ixb^l=ixo^l-kr(idims,^d);
566  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
567  {do ix^db=0,1 \}
568  if({ ix^d==0 .and. ^d==idims | .or.}) then
569  ixbmin^d=ixamin^d-ix^d;
570  ixbmax^d=ixamax^d-ix^d;
571  qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
572  end if
573  {end do\}
574  qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
575  if(tc_saturate) then
576  ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
577  ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
578  bcf=0.d0
579  {do ix^db=0,1 \}
580  if({ ix^d==0 .and. ^d==idims | .or.}) then
581  ixbmin^d=ixamin^d-ix^d;
582  ixbmax^d=ixamax^d-ix^d;
583  bcf(ixa^s,idims)=bcf(ixa^s,idims)+bc(ixb^s,idims)
584  end if
585  {end do\}
586  ! averaged b at face centers
587  bcf(ixa^s,idims)=bcf(ixa^s,idims)*0.5d0**(ndim-1)
588  ixb^l=ixa^l+kr(idims,^d);
589  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))
590  {do ix^db=ixamin^db,ixamax^db\}
591  if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
592  qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
593  end if
594  {end do\}
595  end if
596  end do
597  else
598  ! calculate thermal conduction flux with slope-limited symmetric scheme
599  qvec=0.d0
600  do idims=1,ndim
601  ixb^l=ixo^l-kr(idims,^d);
602  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
603  ! calculate normal of magnetic field
604  ixb^l=ixa^l+kr(idims,^d);
605  bnorm(ixa^s)=0.5d0*(mf(ixa^s,idims)+mf(ixb^s,idims))
606  bcf=0.d0
607  kaf=0.d0
608  kef=0.d0
609  {do ix^db=0,1 \}
610  if({ ix^d==0 .and. ^d==idims | .or.}) then
611  ixbmin^d=ixamin^d-ix^d;
612  ixbmax^d=ixamax^d-ix^d;
613  bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)+bc(ixb^s,1:ndim)
614  kaf(ixa^s)=kaf(ixa^s)+ka(ixb^s)
615  if(tc_perpendicular) kef(ixa^s)=kef(ixa^s)+ke(ixb^s)
616  end if
617  {end do\}
618  ! averaged b at face centers
619  bcf(ixa^s,1:ndim)=bcf(ixa^s,1:ndim)*0.5d0**(ndim-1)
620  ! averaged thermal conductivity at face centers
621  kaf(ixa^s)=kaf(ixa^s)*0.5d0**(ndim-1)
622  if(tc_perpendicular) kef(ixa^s)=kef(ixa^s)*0.5d0**(ndim-1)
623  ! limited normal component
624  minq(ixa^s)=min(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
625  maxq(ixa^s)=max(alpha*gradt(ixa^s,idims),gradt(ixa^s,idims)/alpha)
626  ! eq (19)
627  qdd=0.d0
628  {do ix^db=0,1 \}
629  if({ ix^d==0 .and. ^d==idims | .or.}) then
630  ixbmin^d=ixcmin^d+ix^d;
631  ixbmax^d=ixcmax^d+ix^d;
632  qdd(ixc^s)=qdd(ixc^s)+gradt(ixb^s,idims)
633  end if
634  {end do\}
635  ! temperature gradient at cell corner
636  qdd(ixc^s)=qdd(ixc^s)*0.5d0**(ndim-1)
637  ! eq (21)
638  qe=0.d0
639  {do ix^db=0,1 \}
640  qd(ixc^s)=qdd(ixc^s)
641  if({ ix^d==0 .and. ^d==idims | .or.}) then
642  ixbmin^d=ixamin^d-ix^d;
643  ixbmax^d=ixamax^d-ix^d;
644  where(qd(ixb^s)<=minq(ixa^s))
645  qd(ixb^s)=minq(ixa^s)
646  elsewhere(qd(ixb^s)>=maxq(ixa^s))
647  qd(ixb^s)=maxq(ixa^s)
648  end where
649  qvec(ixa^s,idims)=qvec(ixa^s,idims)+bc(ixb^s,idims)**2*qd(ixb^s)
650  if(tc_perpendicular) qe(ixa^s)=qe(ixa^s)+qd(ixb^s)
651  end if
652  {end do\}
653  qvec(ixa^s,idims)=kaf(ixa^s)*qvec(ixa^s,idims)*0.5d0**(ndim-1)
654  ! add normal flux from perpendicular conduction
655  if(tc_perpendicular) qvec(ixa^s,idims)=qvec(ixa^s,idims)+kef(ixa^s)*qe(ixa^s)*0.5d0**(ndim-1)
656  ! limited transverse component, eq (17)
657  ixbmin^d=ixamin^d;
658  ixbmax^d=ixamax^d+kr(idims,^d);
659  do idir=1,ndim
660  if(idir==idims) cycle
661  qd(ixi^s)=slope_limiter(gradt(ixi^s,idir),ixi^l,ixb^l,idir,-1)
662  qd(ixi^s)=slope_limiter(qd,ixi^l,ixa^l,idims,1)
663  qvec(ixa^s,idims)=qvec(ixa^s,idims)+kaf(ixa^s)*bnorm(ixa^s)*bcf(ixa^s,idir)*qd(ixa^s)
664  end do
665 
666  ! consider magnetic null point
667  !where(Binv(ixA^S)==0.d0)
668  ! qvec(ixA^S,idims)=tc_k_para*(0.5d0*(Te(ixA^S)+Te(ixB^S)))**2.5d0*gradT(ixA^S,idims)
669  !end where
670 
671  if(tc_saturate) then
672  ! consider saturation (Cowie and Mckee 1977 ApJ, 211, 135)
673  ! unsigned saturated TC flux = 5 phi rho c**3, c is isothermal sound speed
674  ixb^l=ixa^l+kr(idims,^d);
675  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))
676  {do ix^db=ixamin^db,ixamax^db\}
677  if(dabs(qvec(ix^d,idims))>qd(ix^d)) then
678  ! write(*,*) 'it',it,qvec(ix^D,idims),qd(ix^D),' TC saturated at ',&
679  ! x(ix^D,:),' rho',w(ix^D,rho_),' Te',Te(ix^D)
680  qvec(ix^d,idims)=sign(1.d0,qvec(ix^d,idims))*qd(ix^d)
681  end if
682  {end do\}
683  end if
684  end do
685  end if
686 
687  qd=0.d0
688  do idims=1,ndim
689  qvec(ix^s,idims)=dxinv(idims)*qvec(ix^s,idims)
690  ixb^l=ixo^l-kr(idims,^d);
691  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
692  end do
693 
694  end subroutine mhd_get_heatconduct
695 
696  function slope_limiter(f,ixI^L,ixO^L,idims,pm) result(lf)
698  integer, intent(in) :: ixI^L, ixO^L, idims, pm
699  double precision, intent(in) :: f(ixi^s)
700  double precision :: lf(ixi^s)
701 
702  double precision :: signf(ixi^s)
703  integer :: ixB^L
704 
705  ixb^l=ixo^l+pm*kr(idims,^d);
706  signf(ixo^s)=sign(1.d0,f(ixo^s))
707  select case(tc_slope_limiter)
708  case('minmod')
709  ! minmod limiter
710  lf(ixo^s)=signf(ixo^s)*max(0.d0,min(abs(f(ixo^s)),signf(ixo^s)*f(ixb^s)))
711  case ('MC')
712  ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
713  lf(ixo^s)=two*signf(ixo^s)* &
714  max(zero,min(dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s),&
715  signf(ixo^s)*quarter*(f(ixb^s)+f(ixo^s))))
716  case ('superbee')
717  ! Roes superbee limiter (eq.3.51i)
718  lf(ixo^s)=signf(ixo^s)* &
719  max(zero,min(two*dabs(f(ixo^s)),signf(ixo^s)*f(ixb^s)),&
720  min(dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s)))
721  case ('koren')
722  ! Barry Koren Right variant
723  lf(ixo^s)=signf(ixo^s)* &
724  max(zero,min(two*dabs(f(ixo^s)),two*signf(ixo^s)*f(ixb^s),&
725  (two*f(ixb^s)*signf(ixo^s)+dabs(f(ixo^s)))*third))
726  case default
727  call mpistop("Unknown slope limiter for thermal conduction")
728  end select
729 
730  end function slope_limiter
731 
732  subroutine mhd_getdt_heatconduct(w,ixI^L,ixO^L,dtnew,dx^D,x)
733  !Check diffusion time limit dt < tc_dtpar*dx_i**2/((gamma-1)*tc_k_para_i/rho)
734  !where tc_k_para_i=tc_k_para*B_i**2/B**2
735  !and T=p/rho
737  use mod_physics
738 
739  integer, intent(in) :: ixI^L, ixO^L
740  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
741  ! note that depending on small_values_method=='error' etc, w values may change
742  ! through call to getpthermal
743  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
744 
745  double precision :: dxinv(1:ndim),mf(ixi^s,1:ndir)
746  double precision :: tmp2(ixi^s),tmp(ixi^s),Te(ixi^s),B2(ixi^s)
747  double precision :: dtdiff_tcond, dtdiff_tsat
748  integer :: idim,ix^D
749 
750  ^d&dxinv(^d)=one/dx^d;
751 
752  call phys_get_pthermal(w,x,ixi^l,ixo^l,tmp)
753 
754  !temperature
755  te(ixo^s)=tmp(ixo^s)/w(ixo^s,rho_)
756  !tc_k_para_i
757  if(tc_constant) then
758  tmp(ixo^s)=tc_k_para
759  else
760  tmp(ixo^s)=tc_k_para*dsqrt(te(ixo^s)**5)/w(ixo^s,rho_)
761  end if
762 
763  ! B
764  if(b0field) then
765  mf(ixo^s,:)=w(ixo^s,iw_mag(:))+block%B0(ixo^s,:,0)
766  else
767  mf(ixo^s,:)=w(ixo^s,iw_mag(:))
768  end if
769  ! B^-2
770  b2(ixo^s)=sum(mf(ixo^s,:)**2,dim=ndim+1)
771  ! B_i**2/B**2
772  where(b2(ixo^s)/=0.d0)
773  ^d&mf(ixo^s,^d)=mf(ixo^s,^d)**2/b2(ixo^s);
774  elsewhere
775  ^d&mf(ixo^s,^d)=1.d0;
776  end where
777 
778  if(tc_saturate) b2(ixo^s)=22.d0*dsqrt(te(ixo^s))
779 
780  do idim=1,ndim
781  tmp2(ixo^s)=tmp(ixo^s)*mf(ixo^s,idim)
782  if(tc_saturate) then
783  where(tmp2(ixo^s)>b2(ixo^s))
784  tmp2(ixo^s)=b2(ixo^s)
785  end where
786  end if
787  ! dt< tc_dtpar * dx_idim**2/((gamma-1)*tc_k_para_i/rho*B_i**2/B**2)
788  dtdiff_tcond=tc_dtpar/(tc_gamma-1.d0)/maxval(tmp2(ixo^s)*dxinv(idim)**2)
789  ! limit the time step
790  dtnew=min(dtnew,dtdiff_tcond)
791  end do
792  dtnew=dtnew/dble(ndim)
793 
794  end subroutine mhd_getdt_heatconduct
795 
796  subroutine hd_get_heatconduct(qd,tmp1,tmp2,ixI^L,ixO^L,w,x,qvec)
799 
800  integer, intent(in) :: ixI^L, ixO^L
801  double precision, intent(in) :: x(ixi^s,1:ndim), w(ixi^s,1:nw)
802  !! tmp store the heat conduction energy changing rate
803  double precision, intent(out) :: qd(ixi^s),tmp1(ixi^s),tmp2(ixi^s)
804  double precision :: qvec(ixi^s,1:ndim),gradT(ixi^s,1:ndim),Te(ixi^s),ke(ixi^s)
805  double precision :: dxinv(ndim)
806  integer, dimension(ndim) :: lowindex
807  integer :: idims,ix^D,ix^L,ixC^L,ixA^L,ixB^L
808 
809  ix^l=ixo^l^ladd1;
810  ! ixC is cell-corner index
811  ixcmax^d=ixomax^d; ixcmin^d=ixomin^d-1;
812 
813  dxinv=1.d0/dxlevel
814 
815  ! tmp1 store internal energy
816  if(solve_internal_e) then
817  tmp1(ixi^s)=w(ixi^s,e_)
818  else
819  ! tmp2 store kinetic energy before addition of heat conduction source
820  tmp2(ixi^s) = 0.5d0*sum(w(ixi^s,iw_mom(:))**2,dim=ndim+1)/w(ixi^s,rho_)
821  tmp1(ixi^s)=w(ixi^s,e_)-tmp2(ixi^s)
822  end if
823 
824  ! Clip off negative pressure if small_pressure is set
825  if(small_values_method=='error') then
826  if (any(tmp1(ixi^s)<small_e) .and. .not.crash) then
827  lowindex=minloc(tmp1(ixi^s))
828  ^d&lowindex(^d)=lowindex(^d)+iximin^d-1;
829  write(*,*)'too low internal energy = ',minval(tmp1(ixi^s)),' at x=',&
830  x(^d&lowindex(^d),1:ndim),lowindex,' with limit=',small_e,' on time=',global_time, ' it=',it
831  write(*,*) 'w',w(^d&lowindex(^d),:)
832  crash=.true.
833  end if
834  else
835  {do ix^db=iximin^db,iximax^db\}
836  if(tmp1(ix^d)<small_e) then
837  tmp1(ix^d)=small_e
838  end if
839  {end do\}
840  end if
841 
842  ! compute temperature before source addition
843  te(ixi^s)=tmp1(ixi^s)/w(ixi^s,rho_)*(tc_gamma-1.d0)
844 
845  ! T gradient at cell faces
846  gradt=0.d0
847  do idims=1,ndim
848  ixbmin^d=ixmin^d;
849  ixbmax^d=ixmax^d-kr(idims,^d);
850  ixa^l=ixb^l+kr(idims,^d);
851  gradt(ixb^s,idims)=(te(ixa^s)-te(ixb^s))*dxinv(idims)
852  end do
853  ! calculate thermal conduction flux with symmetric scheme
854  do idims=1,ndim
855  qd=0.d0
856  {do ix^db=0,1 \}
857  if({ ix^d==0 .and. ^d==idims | .or.}) then
858  ixbmin^d=ixcmin^d+ix^d;
859  ixbmax^d=ixcmax^d+ix^d;
860  qd(ixc^s)=qd(ixc^s)+gradt(ixb^s,idims)
861  end if
862  {end do\}
863  ! temperature gradient at cell corner
864  qvec(ixc^s,idims)=qd(ixc^s)*0.5d0**(ndim-1)
865  end do
866  ! conductivity at cell center
867  qd(ix^s)=tc_k_para*dsqrt(te(ix^s))**5
868  ke=0.d0
869  {do ix^db=0,1\}
870  ixbmin^d=ixcmin^d+ix^d;
871  ixbmax^d=ixcmax^d+ix^d;
872  ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
873  {end do\}
874  ! cell corner conductivity
875  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
876  ! cell corner conduction flux
877  do idims=1,ndim
878  gradt(ixc^s,idims)=ke(ixc^s)*qvec(ixc^s,idims)
879  end do
880 
881  if(tc_saturate) then
882  ! consider saturation with unsigned saturated TC flux = 5 phi rho c**3
883  ! saturation flux at cell center
884  qd(ix^s)=5.d0*w(ix^s,rho_)*dsqrt(te(ix^s)**3)
885  ke=0.d0
886  {do ix^db=0,1\}
887  ixbmin^d=ixcmin^d+ix^d;
888  ixbmax^d=ixcmax^d+ix^d;
889  ke(ixc^s)=ke(ixc^s)+qd(ixb^s)
890  {end do\}
891  ! cell corner saturation flux
892  ke(ixc^s)=0.5d0**ndim*ke(ixc^s)
893  ! magnitude of cell corner conduction flux
894  qd(ixc^s)=norm2(gradt(ixc^s,:),dim=ndim+1)
895  {do ix^db=ixcmin^db,ixcmax^db\}
896  if(qd(ix^d)>ke(ix^d)) then
897  ke(ix^d)=ke(ix^d)/qd(ix^d)
898  do idims=1,ndim
899  gradt(ix^d,idims)=ke(ix^d)*gradt(ix^d,idims)
900  end do
901  end if
902  {end do\}
903  end if
904 
905  ! conductionflux at cell face
906  qvec=0.d0
907  do idims=1,ndim
908  ixb^l=ixo^l-kr(idims,^d);
909  ixamax^d=ixomax^d; ixamin^d=ixbmin^d;
910  {do ix^db=0,1 \}
911  if({ ix^d==0 .and. ^d==idims | .or.}) then
912  ixbmin^d=ixamin^d-ix^d;
913  ixbmax^d=ixamax^d-ix^d;
914  qvec(ixa^s,idims)=qvec(ixa^s,idims)+gradt(ixb^s,idims)
915  end if
916  {end do\}
917  qvec(ixa^s,idims)=qvec(ixa^s,idims)*0.5d0**(ndim-1)
918  end do
919 
920  qd=0.d0
921  do idims=1,ndim
922  qvec(ix^s,idims)=dxinv(idims)*qvec(ix^s,idims)
923  ixb^l=ixo^l-kr(idims,^d);
924  qd(ixo^s)=qd(ixo^s)+qvec(ixo^s,idims)-qvec(ixb^s,idims)
925  end do
926 
927  end subroutine hd_get_heatconduct
928 
929  subroutine hd_getdt_heatconduct(w,ixI^L,ixO^L,dtnew,dx^D,x)
930  ! Check diffusion time limit dt < tc_dtpar * dx_i**2 / ((gamma-1)*tc_k_para_i/rho)
932  use mod_physics
933 
934  integer, intent(in) :: ixI^L, ixO^L
935  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim)
936  double precision, intent(inout) :: w(ixi^s,1:nw), dtnew
937 
938  double precision :: dxinv(1:ndim), tmp(ixi^s), Te(ixi^s)
939  double precision :: dtdiff_tcond,dtdiff_tsat
940  integer :: idim,ix^D
941 
942  ^d&dxinv(^d)=one/dx^d;
943 
944  call phys_get_pthermal(w,x,ixi^l,ixo^l,tmp)
945 
946  te(ixo^s)=tmp(ixo^s)/w(ixo^s,rho_)
947  tmp(ixo^s)=(tc_gamma-one)*tc_k_para*dsqrt((te(ixo^s))**5)/w(ixo^s,rho_)
948 
949  do idim=1,ndim
950  ! dt< tc_dtpar * dx_idim**2/((gamma-1)*tc_k_para_idim/rho)
951  dtdiff_tcond=tc_dtpar/maxval(tmp(ixo^s)*dxinv(idim)**2)
952  if(tc_saturate) then
953  ! dt< tc_dtpar* dx_idim**2/((gamma-1)*sqrt(Te)*5*phi)
954  dtdiff_tsat=tc_dtpar/maxval((tc_gamma-1.d0)*dsqrt(te(ixo^s))*&
955  5.d0*dxinv(idim)**2)
956  ! choose the slower flux (bigger time scale) between classic and saturated
957  dtdiff_tcond=max(dtdiff_tcond,dtdiff_tsat)
958  end if
959  ! limit the time step
960  dtnew=min(dtnew,dtdiff_tcond)
961  end do
962  dtnew=dtnew/dble(ndim)
963 
964  end subroutine hd_getdt_heatconduct
965 
966 end module mod_thermal_conduction
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
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
integer, dimension(:), allocatable type_gradient_limiter
Type of slope limiter used for computing gradients or divergences, when typegrad or typediv are set t...
logical fix_conserve_at_step
Whether to conserve fluxes at the current partial step.
integer, public s
Number of sub-steps of supertime stepping.
integer, parameter plevel_
integer ndir
Number of spatial dimensions (components) for vector variables.
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
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 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.
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)
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.
subroutine, public storeflux(igrid, fC, idimLIM, nwfluxin)
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:160
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, 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:49
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