MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_ffhd_phys.t
Go to the documentation of this file.
1!> Frozen-field hydrodynamics module
3
4#include "amrvac.h"
5
6 use mod_global_parameters, only: std_len, const_c
10 use mod_physics
11 use mod_comm_lib, only: mpistop
12
13 implicit none
14 private
15
16 !> Whether an energy equation is used
17 logical, public, protected :: ffhd_energy = .true.
18
19 !> Whether thermal conduction is used
20 logical, public, protected :: ffhd_thermal_conduction = .false.
21 !> Whether hyperbolic type thermal conduction is used
22 logical, public, protected :: ffhd_hyperbolic_thermal_conduction = .false.
23 !> type of fluid for thermal conduction
24 type(tc_fluid), public, allocatable :: tc_fl
25 !> type of fluid for thermal emission synthesis
26 type(te_fluid), public, allocatable :: te_fl_ffhd
27
28 !> Whether radiative cooling is added
29 logical, public, protected :: ffhd_radiative_cooling = .false.
30 !> type of fluid for radiative cooling
31 type(rc_fluid), public, allocatable :: rc_fl
32
33 !> Whether viscosity is added
34 logical, public, protected :: ffhd_viscosity = .false.
35
36 !> Whether gravity is added
37 logical, public, protected :: ffhd_gravity = .false.
38
39 !> Whether TRAC method is used
40 logical, public, protected :: ffhd_trac = .false.
41
42 !> Which TRAC method is used
43 integer, public, protected :: ffhd_trac_type=1
44
45 !> Height of the mask used in the TRAC method
46 double precision, public, protected :: ffhd_trac_mask = 0.d0
47
48 !> Distance between two adjacent traced magnetic field lines (in finest cell size)
49 integer, public, protected :: ffhd_trac_finegrid=4
50
51 !> Whether plasma is partially ionized
52 logical, public, protected :: ffhd_partial_ionization = .false.
53
54 !> Index of the density (in the w array)
55 integer, public, protected :: rho_
56
57 !> Indices of the momentum density
58 integer, allocatable, public, protected :: mom(:)
59
60 !> Index of the energy density (-1 if not present)
61 integer, public, protected :: e_
62
63 !> Index of the gas pressure (-1 if not present) should equal e_
64 integer, public, protected :: p_
65
66 !> Indices of temperature
67 integer, public, protected :: te_
68
69 !> Index of the cutoff temperature for the TRAC method
70 integer, public, protected :: tcoff_
71 integer, public, protected :: tweight_
72 integer, public, protected :: q_
73
74 !> The adiabatic index
75 double precision, public :: ffhd_gamma = 5.d0/3.0d0
76
77 !> The adiabatic constant
78 double precision, public :: ffhd_adiab = 1.0d0
79
80 !> The small_est allowed energy
81 double precision, protected :: small_e
82
83 !> The thermal conductivity kappa in hyperbolic thermal conduction
84 double precision, public :: hypertc_kappa
85
86 !> Helium abundance over Hydrogen
87 double precision, public, protected :: he_abundance=0.1d0
88 !> Ionization fraction of H
89 !> H_ion_fr = H+/(H+ + H)
90 double precision, public, protected :: h_ion_fr=1d0
91 !> Ionization fraction of He
92 !> He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
93 double precision, public, protected :: he_ion_fr=1d0
94 !> Ratio of number He2+ / number He+ + He2+
95 !> He_ion_fr2 = He2+/(He2+ + He+)
96 double precision, public, protected :: he_ion_fr2=1d0
97 ! used for eq of state when it is not defined by units,
98 ! the units do not contain terms related to ionization fraction
99 ! and it is p = RR * rho * T
100 double precision, public, protected :: rr=1d0
101 ! remove the below flag and assume default value = .false.
102 ! when eq state properly implemented everywhere
103 ! and not anymore through units
104 logical, public, protected :: eq_state_units = .true.
105
106 !> gamma minus one and its inverse
107 double precision :: gamma_1, inv_gamma_1
108
109 !define the subroutine interface for the ambipolar mask
110 abstract interface
111
112 function fun_kin_en(w, ixI^L, ixO^L, inv_rho) result(ke)
113 use mod_global_parameters, only: nw, ndim,block
114 integer, intent(in) :: ixi^l, ixo^l
115 double precision, intent(in) :: w(ixi^s, nw)
116 double precision :: ke(ixo^s)
117 double precision, intent(in), optional :: inv_rho(ixo^s)
118 end function fun_kin_en
119
120 end interface
121
122 procedure(sub_convert), pointer :: ffhd_to_primitive => null()
123 procedure(sub_convert), pointer :: ffhd_to_conserved => null()
124 procedure(sub_small_values), pointer :: ffhd_handle_small_values => null()
125 procedure(sub_get_pthermal), pointer :: ffhd_get_pthermal => null()
126 procedure(sub_get_pthermal), pointer :: ffhd_get_rfactor => null()
127 procedure(sub_get_pthermal), pointer :: ffhd_get_temperature => null()
128 procedure(sub_get_v), pointer :: ffhd_get_v => null()
129 procedure(fun_kin_en), pointer :: ffhd_kin_en => null()
130 ! Public methods
131 public :: ffhd_phys_init
132 public :: ffhd_kin_en
133 public :: ffhd_get_pthermal
134 public :: ffhd_get_temperature
135 public :: ffhd_get_v
136 public :: ffhd_get_rho
137 public :: ffhd_get_v_idim
138 public :: ffhd_to_conserved
139 public :: ffhd_to_primitive
140 public :: ffhd_get_csound2
141 public :: ffhd_e_to_ei
142 public :: ffhd_ei_to_e
143
144contains
145
146 subroutine ffhd_read_params(files)
148 use mod_particles, only: particles_eta, particles_etah
149 character(len=*), intent(in) :: files(:)
150 integer :: n
151
152 namelist /ffhd_list/ ffhd_energy, ffhd_gamma, ffhd_adiab,&
156
157 do n = 1, size(files)
158 open(unitpar, file=trim(files(n)), status="old")
159 read(unitpar, ffhd_list, end=111)
160111 close(unitpar)
161 end do
162 end subroutine ffhd_read_params
163
164 !> Write this module's parameters to a snapsoht
165 subroutine ffhd_write_info(fh)
167 integer, intent(in) :: fh
168 integer, parameter :: n_par = 1
169 double precision :: values(n_par)
170 character(len=name_len) :: names(n_par)
171 integer, dimension(MPI_STATUS_SIZE) :: st
172 integer :: er
173
174 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
175
176 names(1) = "gamma"
177 values(1) = ffhd_gamma
178 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
179 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
180 end subroutine ffhd_write_info
181
182 subroutine ffhd_phys_init()
187 use mod_gravity, only: gravity_init
192 integer :: itr, idir
193
194 call ffhd_read_params(par_files)
195
196 if(.not. ffhd_energy) then
199 if(mype==0) write(*,*) 'WARNING: set ffhd_thermal_conduction=F when ffhd_energy=F'
200 end if
203 if(mype==0) write(*,*) 'WARNING: set ffhd_hyperbolic_thermal_conduction=F when ffhd_energy=F'
204 end if
207 if(mype==0) write(*,*) 'WARNING: set ffhd_radiative_cooling=F when ffhd_energy=F'
208 end if
209 if(ffhd_trac) then
210 ffhd_trac=.false.
211 if(mype==0) write(*,*) 'WARNING: set ffhd_trac=F when ffhd_energy=F'
212 end if
215 if(mype==0) write(*,*) 'WARNING: set ffhd_partial_ionization=F when ffhd_energy=F'
216 end if
217 end if
218 if(.not.eq_state_units) then
221 if(mype==0) write(*,*) 'WARNING: set ffhd_partial_ionization=F when eq_state_units=F'
222 end if
223 end if
224
227 if(mype==0) write(*,*) 'WARNING: turn off parabolic TC when using hyperbolic TC'
228 end if
229
230 physics_type = "ffhd"
231 phys_energy=ffhd_energy
232 phys_internal_e=.false.
235 phys_partial_ionization=ffhd_partial_ionization
236
237 phys_gamma = ffhd_gamma
238 phys_total_energy=ffhd_energy
240
241 {^ifoned
242 if(ffhd_trac .and. ffhd_trac_type .gt. 2) then
244 if(mype==0) write(*,*) 'WARNING: reset ffhd_trac_type=1 for 1D simulation'
245 end if
246 }
247 if(ffhd_trac .and. ffhd_trac_type .le. 4) then
248 ffhd_trac_mask=bigdouble
249 if(mype==0) write(*,*) 'WARNING: set ffhd_trac_mask==bigdouble for global TRAC method'
250 end if
252
253 allocate(start_indices(number_species),stop_indices(number_species))
254 ! set the index of the first flux variable for species 1
255 start_indices(1)=1
256 ! Determine flux variables
257 rho_ = var_set_rho()
258
259 allocate(mom(1))
260 mom(:) = var_set_momentum(1)
261
262 ! Set index of energy variable
263 if(ffhd_energy) then
264 e_ = var_set_energy() ! energy density
265 p_ = e_ ! gas pressure
266 else
267 e_ = -1
268 p_ = -1
269 end if
270
272 q_ = var_set_q()
273 need_global_cs2max=.true.
274 else
275 q_=-1
276 end if
277
278 ! set number of variables which need update ghostcells
279 nwgc=nwflux
280
281 ! set the index of the last flux variable for species 1
282 stop_indices(1)=nwflux
283
284 ! set temperature as an auxiliary variable to get ionization degree
286 te_ = var_set_auxvar('Te','Te')
287 else
288 te_ = -1
289 end if
290
291 ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
292 tweight_ = -1
293 if(ffhd_trac) then
294 tcoff_ = var_set_wextra()
295 iw_tcoff=tcoff_
296 if(ffhd_trac_type .ge. 3) then
297 tweight_ = var_set_wextra()
298 iw_tweight=tweight_
299 end if
300 else
301 tcoff_ = -1
302 end if
303
304 nvector = 0 ! No. vector vars
305
306 ! Check whether custom flux types have been defined
307 if(.not. allocated(flux_type)) then
308 allocate(flux_type(ndir, nwflux))
309 flux_type = flux_default
310 else if(any(shape(flux_type) /= [ndir, nwflux])) then
311 call mpistop("phys_check error: flux_type has wrong shape")
312 end if
313
314 phys_get_dt => ffhd_get_dt
315 phys_get_cmax => ffhd_get_cmax_origin
316 phys_get_a2max => ffhd_get_a2max
317 phys_get_cs2max => ffhd_get_cs2max
318 phys_get_tcutoff => ffhd_get_tcutoff
319 phys_get_cbounds => ffhd_get_cbounds
320 phys_to_primitive => ffhd_to_primitive_origin
321 ffhd_to_primitive => ffhd_to_primitive_origin
322 phys_to_conserved => ffhd_to_conserved_origin
323 ffhd_to_conserved => ffhd_to_conserved_origin
324 phys_get_flux => ffhd_get_flux
325 phys_get_v => ffhd_get_v_origin
326 ffhd_get_v => ffhd_get_v_origin
327 phys_get_rho => ffhd_get_rho
328 ffhd_kin_en => ffhd_kin_en_origin
329 !> need to check source geom here
330 !phys_add_source_geom => ffhd_add_source_geom
331 phys_add_source => ffhd_add_source
332 phys_check_params => ffhd_check_params
333 phys_write_info => ffhd_write_info
334 phys_handle_small_values => ffhd_handle_small_values_origin
335 ffhd_handle_small_values => ffhd_handle_small_values_origin
336 phys_check_w => ffhd_check_w_origin
337
338 if(.not.ffhd_energy) then
339 phys_get_pthermal => ffhd_get_pthermal_iso
340 ffhd_get_pthermal => ffhd_get_pthermal_iso
341 else
342 phys_get_pthermal => ffhd_get_pthermal_origin
343 ffhd_get_pthermal => ffhd_get_pthermal_origin
344 end if
345
346 ! choose Rfactor in ideal gas law
348 ffhd_get_rfactor=>rfactor_from_temperature_ionization
349 phys_update_temperature => ffhd_update_temperature
350 else if(associated(usr_rfactor)) then
351 ffhd_get_rfactor=>usr_rfactor
352 else
353 ffhd_get_rfactor=>rfactor_from_constant_ionization
354 end if
355
357 ffhd_get_temperature => ffhd_get_temperature_from_te
358 else
359 ffhd_get_temperature => ffhd_get_temperature_from_etot
360 end if
361
362 ! derive units from basic units
363 call ffhd_physical_units()
364
367 end if
368 if(.not. ffhd_energy .and. ffhd_thermal_conduction) then
369 call mpistop("thermal conduction needs ffhd_energy=T")
370 end if
372 call mpistop("hyperbolic thermal conduction needs ffhd_energy=T")
373 end if
374 if(.not. ffhd_energy .and. ffhd_radiative_cooling) then
375 call mpistop("radiative cooling needs ffhd_energy=T")
376 end if
377
378 ! initialize thermal conduction module
380 call sts_init()
382
383 allocate(tc_fl)
384 call tc_get_hd_params(tc_fl,tc_params_read_ffhd)
385 call add_sts_method(ffhd_get_tc_dt_ffhd,ffhd_sts_set_source_tc_ffhd,e_,1,e_,1,.false.)
386 tc_fl%get_temperature_from_conserved => ffhd_get_temperature_from_etot
387 tc_fl%get_temperature_from_eint => ffhd_get_temperature_from_eint
389 call set_error_handling_to_head(ffhd_tc_handle_small_e)
390 tc_fl%get_rho => ffhd_get_rho
391 tc_fl%e_ = e_
392 tc_fl%Tcoff_ = tcoff_
393 end if
394
395 ! Initialize radiative cooling module
398 allocate(rc_fl)
399 call radiative_cooling_init(rc_fl,rc_params_read)
400 rc_fl%get_rho => ffhd_get_rho
401 rc_fl%get_pthermal => ffhd_get_pthermal
402 rc_fl%get_var_Rfactor => ffhd_get_rfactor
403 rc_fl%e_ = e_
404 rc_fl%Tcoff_ = tcoff_
405 rc_fl%has_equi = .false.
406 end if
407 allocate(te_fl_ffhd)
408 te_fl_ffhd%get_rho=> ffhd_get_rho
409 te_fl_ffhd%get_pthermal=> ffhd_get_pthermal
410 te_fl_ffhd%get_var_Rfactor => ffhd_get_rfactor
411{^ifthreed
412 phys_te_images => ffhd_te_images
413}
414 ! Initialize viscosity module
415 if(ffhd_viscosity) call viscosity_init(phys_wider_stencil)
416
417 ! Initialize gravity module
418 if(ffhd_gravity) then
419 call gravity_init()
420 end if
421
422 ! initialize ionization degree table
424 end subroutine ffhd_phys_init
425
426{^ifthreed
427 subroutine ffhd_te_images
430
431 select case(convert_type)
432 case('EIvtiCCmpi','EIvtuCCmpi')
434 case('ESvtiCCmpi','ESvtuCCmpi')
436 case('SIvtiCCmpi','SIvtuCCmpi')
438 case('WIvtiCCmpi','WIvtuCCmpi')
440 case default
441 call mpistop("Error in synthesize emission: Unknown convert_type")
442 end select
443 end subroutine ffhd_te_images
444}
445
446 subroutine ffhd_sts_set_source_tc_ffhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
450 integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
451 double precision, intent(in) :: x(ixi^s,1:ndim)
452 double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
453 double precision, intent(in) :: my_dt
454 logical, intent(in) :: fix_conserve_at_step
455 call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
456 end subroutine ffhd_sts_set_source_tc_ffhd
457
458 function ffhd_get_tc_dt_ffhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
459 !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
460 !where tc_k_para_i=tc_k_para*B_i**2/B**2
461 !and T=p/rho
464
465 integer, intent(in) :: ixi^l, ixo^l
466 double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
467 double precision, intent(in) :: w(ixi^s,1:nw)
468 double precision :: dtnew
469
470 dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
471 end function ffhd_get_tc_dt_ffhd
472
473 subroutine ffhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
475
476 integer, intent(in) :: ixi^l,ixo^l
477 double precision, intent(inout) :: w(ixi^s,1:nw)
478 double precision, intent(in) :: x(ixi^s,1:ndim)
479 integer, intent(in) :: step
480 character(len=140) :: error_msg
481
482 write(error_msg,"(a,i3)") "Thermal conduction step ", step
483 call ffhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,error_msg)
484 end subroutine ffhd_tc_handle_small_e
485
486 subroutine tc_params_read_ffhd(fl)
488 type(tc_fluid), intent(inout) :: fl
489 integer :: n
490 ! list parameters
491 logical :: tc_saturate=.false.
492 double precision :: tc_k_para=0d0
493 character(len=std_len) :: tc_slope_limiter="MC"
494
495 namelist /tc_list/ tc_saturate, tc_slope_limiter, tc_k_para
496
497 do n = 1, size(par_files)
498 open(unitpar, file=trim(par_files(n)), status="old")
499 read(unitpar, tc_list, end=111)
500111 close(unitpar)
501 end do
502
503 fl%tc_saturate = tc_saturate
504 fl%tc_k_para = tc_k_para
505 select case(tc_slope_limiter)
506 case ('no','none')
507 fl%tc_slope_limiter = 0
508 case ('MC')
509 ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
510 fl%tc_slope_limiter = 1
511 case('minmod')
512 ! minmod limiter
513 fl%tc_slope_limiter = 2
514 case ('superbee')
515 ! Roes superbee limiter (eq.3.51i)
516 fl%tc_slope_limiter = 3
517 case ('koren')
518 ! Barry Koren Right variant
519 fl%tc_slope_limiter = 4
520 case default
521 call mpistop("Unknown tc_slope_limiter, choose MC, minmod")
522 end select
523 end subroutine tc_params_read_ffhd
524
525 subroutine rc_params_read(fl)
527 use mod_constants, only: bigdouble
528 type(rc_fluid), intent(inout) :: fl
529 integer :: n
530 integer :: ncool = 4000
531 double precision :: cfrac=0.1d0
532
533 !> Name of cooling curve
534 character(len=std_len) :: coolcurve='JCcorona'
535
536 !> Name of cooling method
537 character(len=std_len) :: coolmethod='exact'
538
539 !> Fixed temperature not lower than tlow
540 logical :: tfix=.false.
541
542 !> Lower limit of temperature
543 double precision :: tlow=bigdouble
544
545 !> Add cooling source in a split way (.true.) or un-split way (.false.)
546 logical :: rc_split=.false.
547 logical :: rad_cut=.false.
548 double precision :: rad_cut_hgt=0.5d0
549 double precision :: rad_cut_dey=0.15d0
550
551 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split, rad_cut, rad_cut_hgt, rad_cut_dey
552
553 do n = 1, size(par_files)
554 open(unitpar, file=trim(par_files(n)), status="old")
555 read(unitpar, rc_list, end=111)
556111 close(unitpar)
557 end do
558
559 fl%ncool=ncool
560 fl%coolcurve=coolcurve
561 fl%coolmethod=coolmethod
562 fl%tlow=tlow
563 fl%Tfix=tfix
564 fl%rc_split=rc_split
565 fl%cfrac=cfrac
566 fl%rad_cut=rad_cut
567 fl%rad_cut_hgt=rad_cut_hgt
568 fl%rad_cut_dey=rad_cut_dey
569 end subroutine rc_params_read
570
571 subroutine ffhd_check_params
575
576 gamma_1=ffhd_gamma-1.d0
577 if (.not. ffhd_energy) then
578 if (ffhd_gamma <= 0.0d0) call mpistop ("Error: ffhd_gamma <= 0")
579 if (ffhd_adiab < 0.0d0) call mpistop ("Error: ffhd_adiab < 0")
581 else
582 if (ffhd_gamma <= 0.0d0 .or. ffhd_gamma == 1.0d0) &
583 call mpistop ("Error: ffhd_gamma <= 0 or ffhd_gamma == 1")
584 inv_gamma_1=1.d0/gamma_1
585 small_e = small_pressure * inv_gamma_1
586 end if
587
588 if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
589 call mpistop("usr_set_equi_vars has to be implemented in the user file")
590 end if
591 end subroutine ffhd_check_params
592
593 subroutine ffhd_physical_units()
595 double precision :: mp,kb
596 double precision :: a,b
597
598 if(si_unit) then
599 mp=mp_si
600 kb=kb_si
601 else
602 mp=mp_cgs
603 kb=kb_cgs
604 end if
605 if(eq_state_units) then
606 a=1d0+4d0*he_abundance
608 b=1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0)
609 else
610 b=2d0+3d0*he_abundance
611 end if
612 rr=1d0
613 else
614 a=1d0
615 b=1d0
616 rr=(1d0+h_ion_fr+he_abundance*(he_ion_fr*(he_ion_fr2+1d0)+1d0))/(1d0+4d0*he_abundance)
617 end if
618 if(unit_density/=1.d0 .or. unit_numberdensity/=1.d0) then
619 if(unit_density/=1.d0) then
621 else if(unit_numberdensity/=1.d0) then
623 end if
624 if(unit_temperature/=1.d0) then
627 if(unit_length/=1.d0) then
629 else if(unit_time/=1.d0) then
631 end if
632 else if(unit_pressure/=1.d0) then
635 if(unit_length/=1.d0) then
637 else if(unit_time/=1.d0) then
639 end if
640 else if(unit_velocity/=1.d0) then
643 if(unit_length/=1.d0) then
645 else if(unit_time/=1.d0) then
647 end if
648 else if(unit_time/=1.d0) then
652 end if
653 else if(unit_temperature/=1.d0) then
654 ! units of temperature and velocity are dependent
655 if(unit_pressure/=1.d0) then
659 if(unit_length/=1.d0) then
661 else if(unit_time/=1.d0) then
663 end if
664 end if
665 else if(unit_pressure/=1.d0) then
666 if(unit_velocity/=1.d0) then
670 if(unit_length/=1.d0) then
672 else if(unit_time/=1.d0) then
674 end if
675 else if(unit_time/=0.d0) then
680 end if
681 end if
683 end subroutine ffhd_physical_units
684
685 subroutine ffhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
687 logical, intent(in) :: primitive
688 integer, intent(in) :: ixi^l, ixo^l
689 double precision, intent(in) :: w(ixi^s,nw)
690 double precision :: tmp(ixi^s)
691 logical, intent(inout) :: flag(ixi^s,1:nw)
692
693 flag=.false.
694 where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
695
696 if(ffhd_energy) then
697 if(primitive) then
698 where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
699 else
700 tmp(ixo^s)=w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l)
701 where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
702 end if
703 end if
704 end subroutine ffhd_check_w_origin
705
706 subroutine ffhd_to_conserved_origin(ixI^L,ixO^L,w,x)
708 integer, intent(in) :: ixi^l, ixo^l
709 double precision, intent(inout) :: w(ixi^s, nw)
710 double precision, intent(in) :: x(ixi^s, 1:ndim)
711 double precision :: inv_gamma2(ixo^s)
712 integer :: idir
713
714 if(ffhd_energy) then
715 w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1+half*w(ixo^s,mom(1))**2*w(ixo^s,rho_)
716 end if
717 w(ixo^s,mom(1))=w(ixo^s,rho_)*w(ixo^s,mom(1))
718 end subroutine ffhd_to_conserved_origin
719
720 subroutine ffhd_to_primitive_origin(ixI^L,ixO^L,w,x)
722 integer, intent(in) :: ixi^l, ixo^l
723 double precision, intent(inout) :: w(ixi^s, nw)
724 double precision, intent(in) :: x(ixi^s, 1:ndim)
725 double precision :: inv_rho(ixo^s), gamma2(ixo^s)
726
727 if(fix_small_values) then
728 !> fix small values preventing NaN numbers in the following converting
729 call ffhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'ffhd_to_primitive_origin')
730 end if
731
732 w(ixo^s,mom(1)) = w(ixo^s,mom(1))/w(ixo^s,rho_)
733 if(ffhd_energy) then
734 w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-half*w(ixo^s,rho_)*w(ixo^s,mom(1))**2)
735 end if
736 end subroutine ffhd_to_primitive_origin
737
738 subroutine ffhd_ei_to_e(ixI^L,ixO^L,w,x)
740 integer, intent(in) :: ixi^l, ixo^l
741 double precision, intent(inout) :: w(ixi^s, nw)
742 double precision, intent(in) :: x(ixi^s, 1:ndim)
743
744 w(ixi^s,e_)=w(ixi^s,e_)+ffhd_kin_en(w,ixi^l,ixi^l)
745 end subroutine ffhd_ei_to_e
746
747 subroutine ffhd_e_to_ei(ixI^L,ixO^L,w,x)
749 integer, intent(in) :: ixi^l, ixo^l
750 double precision, intent(inout) :: w(ixi^s, nw)
751 double precision, intent(in) :: x(ixi^s, 1:ndim)
752
753 w(ixi^s,e_)=w(ixi^s,e_)-ffhd_kin_en(w,ixi^l,ixi^l)
754 if(fix_small_values) then
755 call ffhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'ffhd_e_to_ei')
756 end if
757 end subroutine ffhd_e_to_ei
758
759 subroutine ffhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
762 logical, intent(in) :: primitive
763 integer, intent(in) :: ixi^l,ixo^l
764 double precision, intent(inout) :: w(ixi^s,1:nw)
765 double precision, intent(in) :: x(ixi^s,1:ndim)
766 character(len=*), intent(in) :: subname
767
768 logical :: flag(ixi^s,1:nw)
769 double precision :: tmp2(ixi^s)
770
771 call phys_check_w(primitive, ixi^l, ixi^l, w, flag)
772
773 if(any(flag)) then
774 select case (small_values_method)
775 case ("replace")
776 where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
777 if(small_values_fix_iw(mom(1))) then
778 where(flag(ixo^s,rho_)) w(ixo^s, mom(1)) = 0.0d0
779 end if
780 if(ffhd_energy) then
781 if(primitive) then
782 where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
783 else
784 where(flag(ixo^s,e_))
785 w(ixo^s,e_) = small_e+ffhd_kin_en(w,ixi^l,ixo^l)
786 end where
787 end if
788 end if
789 case ("average")
790 call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
791 if(ffhd_energy) then
792 if(primitive) then
793 call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
794 else
795 w(ixi^s,e_)=w(ixi^s,e_)-ffhd_kin_en(w,ixi^l,ixi^l)
796 call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
797 w(ixi^s,e_)=w(ixi^s,e_)+ffhd_kin_en(w,ixi^l,ixi^l)
798 end if
799 end if
800 case default
801 if(.not.primitive) then
802 if(ffhd_energy) then
803 w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l))
804 end if
805 w(ixo^s,mom(1))=w(ixo^s,mom(1))/w(ixo^s,rho_)
806 end if
807 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
808 end select
809 end if
810 end subroutine ffhd_handle_small_values_origin
811
812 subroutine ffhd_get_v_origin(w,x,ixI^L,ixO^L,v)
814 integer, intent(in) :: ixi^l, ixo^l
815 double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
816 double precision, intent(out) :: v(ixi^s,ndir)
817 double precision :: rho(ixi^s)
818 integer :: idir
819
820 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
821 rho(ixo^s)=1.d0/rho(ixo^s)
822 do idir=1,ndir
823 v(ixo^s,ndir) = w(ixo^s,mom(1))*block%B0(ixo^s,idir,0)*rho(ixo^s)
824 end do
825 end subroutine ffhd_get_v_origin
826
827 subroutine ffhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
829 integer, intent(in) :: ixi^l, ixo^l, idim
830 double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
831 double precision, intent(out) :: v(ixi^s)
832 double precision :: rho(ixi^s)
833
834 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
835 v(ixo^s) = (w(ixo^s, mom(1))*block%B0(ixo^s,idim,0)) / rho(ixo^s)
836 end subroutine ffhd_get_v_idim
837
838 subroutine ffhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
840 integer, intent(in) :: ixi^l, ixo^l, idim
841 ! w in primitive form
842 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
843 double precision, intent(inout) :: cmax(ixi^s)
844
845 if(ffhd_energy) then
846 cmax(ixo^s)=sqrt(ffhd_gamma*w(ixo^s,p_)/w(ixo^s,rho_))*abs(block%B0(ixo^s,idim,idim))
847 else
848 cmax(ixo^s)=sqrt(ffhd_gamma*ffhd_adiab*w(ixo^s,rho_)**gamma_1)*abs(block%B0(ixo^s,idim,idim))
849 end if
850 cmax(ixo^s)=abs(w(ixo^s,mom(1))*block%B0(ixo^s,idim,0))+cmax(ixo^s)
851
852 end subroutine ffhd_get_cmax_origin
853
854 subroutine ffhd_get_cs2max(w,x,ixI^L,ixO^L,cs2max)
856 integer, intent(in) :: ixi^l, ixo^l
857 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
858 double precision, intent(inout) :: cs2max
859 double precision :: cs2(ixi^s)
860
861 call ffhd_get_csound2(w,x,ixi^l,ixo^l,cs2)
862 cs2max=maxval(cs2(ixo^s))
863 end subroutine ffhd_get_cs2max
864
865 subroutine ffhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
867 integer, intent(in) :: ixi^l, ixo^l
868 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
869 double precision, intent(inout) :: a2max(ndim)
870 double precision :: a2(ixi^s,ndim,nw)
871 integer :: gxo^l,hxo^l,jxo^l,kxo^l,i,j
872
873 a2=zero
874 do i = 1,ndim
875 !> 4th order
876 hxo^l=ixo^l-kr(i,^d);
877 gxo^l=hxo^l-kr(i,^d);
878 jxo^l=ixo^l+kr(i,^d);
879 kxo^l=jxo^l+kr(i,^d);
880 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
881 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
882 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
883 end do
884 end subroutine ffhd_get_a2max
885
886 subroutine ffhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
888 use mod_geometry
889 integer, intent(in) :: ixi^l,ixo^l
890 double precision, intent(in) :: x(ixi^s,1:ndim)
891 double precision, intent(inout) :: w(ixi^s,1:nw)
892 double precision, intent(out) :: tco_local,tmax_local
893 double precision, parameter :: trac_delta=0.25d0
894 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
895 double precision, dimension(ixI^S,1:ndir) :: bunitvec
896 double precision, dimension(ixI^S,1:ndim) :: gradt
897 double precision :: bdir(ndim)
898 double precision :: ltrc,ltrp,altr(ixi^s)
899 integer :: idims,jxo^l,hxo^l,ixa^d,ixb^d
900 integer :: jxp^l,hxp^l,ixp^l,ixq^l
901 logical :: lrlt(ixi^s)
902
903 call ffhd_get_temperature(w,x,ixi^l,ixi^l,te)
904 tco_local=zero
905 tmax_local=maxval(te(ixo^s))
906
907 {^ifoned
908 select case(ffhd_trac_type)
909 case(0)
910 !> test case, fixed cutoff temperature
911 block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
912 case(1)
913 hxo^l=ixo^l-1;
914 jxo^l=ixo^l+1;
915 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
916 lrlt=.false.
917 where(lts(ixo^s) > trac_delta)
918 lrlt(ixo^s)=.true.
919 end where
920 if(any(lrlt(ixo^s))) then
921 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
922 end if
923 case(2)
924 !> iijima et al. 2021, LTRAC method
925 ltrc=1.5d0
926 ltrp=4.d0
927 ixp^l=ixo^l^ladd1;
928 hxo^l=ixo^l-1;
929 jxo^l=ixo^l+1;
930 hxp^l=ixp^l-1;
931 jxp^l=ixp^l+1;
932 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
933 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
934 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
935 block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
936 case default
937 call mpistop("ffhd_trac_type not allowed for 1D simulation")
938 end select
939 }
940 {^nooned
941 select case(ffhd_trac_type)
942 case(0)
943 !> test case, fixed cutoff temperature
944 if(slab_uniform) then
945 !> assume cgs units
946 block%wextra(ixi^s,tcoff_)=max(min(3.d5/unit_temperature,6.d5/unit_temperature-3.d-4/unit_temperature*unit_length*x(ixi^s,ndim)),zero)
947 else
948 block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
949 end if
950 case(1,4,6)
951 do idims=1,ndim
952 call gradient(te,ixi^l,ixo^l,idims,tmp1)
953 gradt(ixo^s,idims)=tmp1(ixo^s)
954 end do
955 bunitvec(ixo^s,:)=block%B0(ixo^s,:,0)
956 if(ffhd_trac_type .gt. 1) then
957 ! B direction at cell center
958 bdir=zero
959 {do ixa^d=0,1\}
960 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
961 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
962 {end do\}
963 if(sum(bdir(:)**2) .gt. zero) then
964 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
965 end if
966 block%special_values(3:ndim+2)=bdir(1:ndim)
967 end if
968 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
969 where(tmp1(ixo^s)/=0.d0)
970 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
971 else where
972 tmp1(ixo^s)=bigdouble
973 end where
974 ! b unit vector: magnetic field direction vector
975 do idims=1,ndim
976 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
977 end do
978 ! temperature length scale inversed
979 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
980 ! fraction of cells size to temperature length scale
981 if(slab_uniform) then
982 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
983 else
984 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
985 end if
986 lrlt=.false.
987 where(lts(ixo^s) > trac_delta)
988 lrlt(ixo^s)=.true.
989 end where
990 if(any(lrlt(ixo^s))) then
991 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
992 else
993 block%special_values(1)=zero
994 end if
995 block%special_values(2)=tmax_local
996 case(2)
997 !> iijima et al. 2021, LTRAC method
998 ltrc=1.5d0
999 ltrp=4.d0
1000 ixp^l=ixo^l^ladd2;
1001 do idims=1,ndim
1002 ixq^l=ixp^l;
1003 hxp^l=ixp^l;
1004 jxp^l=ixp^l;
1005 select case(idims)
1006 {case(^d)
1007 ixqmin^d=ixqmin^d+1
1008 ixqmax^d=ixqmax^d-1
1009 hxpmax^d=ixpmin^d
1010 jxpmin^d=ixpmax^d
1011 \}
1012 end select
1013 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
1014 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
1015 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
1016 end do
1017 bunitvec(ixp^s,:)=block%B0(ixp^s,:,0)
1018 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
1019 if(slab_uniform) then
1020 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
1021 else
1022 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
1023 end if
1024 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
1025
1026 altr=zero
1027 ixp^l=ixo^l^ladd1;
1028 do idims=1,ndim
1029 hxo^l=ixp^l-kr(idims,^d);
1030 jxo^l=ixp^l+kr(idims,^d);
1031 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
1032 end do
1033 block%wextra(ixp^s,tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
1034 case(3,5)
1035 !> do nothing here
1036 case default
1037 call mpistop("unknown ffhd_trac_type")
1038 end select
1039 }
1040 end subroutine ffhd_get_tcutoff
1041
1042 subroutine ffhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
1044 integer, intent(in) :: ixi^l, ixo^l, idim
1045 double precision, intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
1046 double precision, intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
1047 double precision, intent(in) :: x(ixi^s,1:ndim)
1048 double precision, intent(inout) :: cmax(ixi^s,1:number_species)
1049 double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
1050 double precision, intent(in) :: hspeed(ixi^s,1:number_species)
1051 double precision :: wmean(ixi^s,nw)
1052 double precision, dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
1053 integer :: ix^d
1054
1055 select case (boundspeed)
1056 case (1)
1057 ! This implements formula (10.52) from "Riemann Solvers and Numerical
1058 ! Methods for Fluid Dynamics" by Toro.
1059 tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
1060 tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
1061 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
1062 umean(ixo^s)=(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)*tmp1(ixo^s)&
1063 +wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)*tmp2(ixo^s))*tmp3(ixo^s)
1064 call ffhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1065 call ffhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1066 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
1067 tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
1068 (wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)-wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))**2
1069 dmean(ixo^s)=dsqrt(dmean(ixo^s))
1070 if(present(cmin)) then
1071 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
1072 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
1073 else
1074 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
1075 end if
1076 case (2)
1077 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1078 tmp1(ixo^s)=wmean(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)/wmean(ixo^s,rho_)
1079 call ffhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
1080 if(present(cmin)) then
1081 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
1082 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
1083 else
1084 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
1085 end if
1086 case (3)
1087 ! Miyoshi 2005 JCP 208, 315 equation (67)
1088 call ffhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1089 call ffhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1090 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
1091 if(present(cmin)) then
1092 cmin(ixo^s,1)=min(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
1093 wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))-csoundl(ixo^s)
1094 cmax(ixo^s,1)=max(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
1095 wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))+csoundl(ixo^s)
1096 else
1097 cmax(ixo^s,1)=max(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
1098 wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))+csoundl(ixo^s)
1099 end if
1100 end select
1101 end subroutine ffhd_get_cbounds
1102
1103 subroutine ffhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
1105 integer, intent(in) :: ixi^l, ixo^l, idim
1106 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
1107 double precision, intent(out):: csound(ixi^s)
1108
1109 call ffhd_get_csound2(w,x,ixi^l,ixo^l,csound)
1110 csound(ixo^s) = dsqrt(csound(ixo^s))*abs(block%B0(ixo^s,idim,idim))
1111 end subroutine ffhd_get_csound
1112
1113 !> Calculate fast magnetosonic wave speed
1114 subroutine ffhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
1116
1117 integer, intent(in) :: ixi^l, ixo^l, idim
1118 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
1119 double precision, intent(out):: csound(ixi^s)
1120
1121 if(ffhd_energy) then
1122 csound(ixo^s)=ffhd_gamma*w(ixo^s,e_)/w(ixo^s,rho_)
1123 else
1124 csound(ixo^s)=ffhd_gamma*ffhd_adiab*w(ixo^s,rho_)**gamma_1
1125 end if
1126 csound(ixo^s) = dsqrt(csound(ixo^s))
1127 end subroutine ffhd_get_csound_prim
1128
1129 subroutine ffhd_get_pthermal_iso(w,x,ixI^L,ixO^L,pth)
1131
1132 integer, intent(in) :: ixi^l, ixo^l
1133 double precision, intent(in) :: w(ixi^s,nw)
1134 double precision, intent(in) :: x(ixi^s,1:ndim)
1135 double precision, intent(out):: pth(ixi^s)
1136
1137 call ffhd_get_rho(w,x,ixi^l,ixo^l,pth)
1138 pth(ixo^s)=ffhd_adiab*pth(ixo^s)**ffhd_gamma
1139 end subroutine ffhd_get_pthermal_iso
1140
1141 subroutine ffhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
1144 integer, intent(in) :: ixi^l, ixo^l
1145 double precision, intent(in) :: w(ixi^s,nw)
1146 double precision, intent(in) :: x(ixi^s,1:ndim)
1147 double precision, intent(out):: pth(ixi^s)
1148 integer :: iw, ix^d
1149
1150 pth(ixo^s)=gamma_1*(w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l))
1151 if (fix_small_values) then
1152 {do ix^db= ixo^lim^db\}
1153 if(pth(ix^d)<small_pressure) then
1154 pth(ix^d)=small_pressure
1155 end if
1156 {end do^D&\}
1157 elseif(check_small_values) then
1158 {do ix^db= ixo^lim^db\}
1159 if(pth(ix^d)<small_pressure) then
1160 write(*,*) "Error: small value of gas pressure",pth(ix^d),&
1161 " encountered when call ffhd_get_pthermal"
1162 write(*,*) "Iteration: ", it, " Time: ", global_time
1163 write(*,*) "Location: ", x(ix^d,:)
1164 write(*,*) "Cell number: ", ix^d
1165 do iw=1,nw
1166 write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
1167 end do
1168 ! use erroneous arithmetic operation to crash the run
1169 if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
1170 write(*,*) "Saving status at the previous time step"
1171 crash=.true.
1172 end if
1173 {end do^D&\}
1174 end if
1175 end subroutine ffhd_get_pthermal_origin
1176
1177 subroutine ffhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
1179 integer, intent(in) :: ixi^l, ixo^l
1180 double precision, intent(in) :: w(ixi^s, 1:nw)
1181 double precision, intent(in) :: x(ixi^s, 1:ndim)
1182 double precision, intent(out):: res(ixi^s)
1183
1184 res(ixo^s) = w(ixo^s, te_)
1185 end subroutine ffhd_get_temperature_from_te
1186
1187 subroutine ffhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1189 integer, intent(in) :: ixi^l, ixo^l
1190 double precision, intent(in) :: w(ixi^s, 1:nw)
1191 double precision, intent(in) :: x(ixi^s, 1:ndim)
1192 double precision, intent(out):: res(ixi^s)
1193 double precision :: r(ixi^s)
1194
1195 call ffhd_get_rfactor(w,x,ixi^l,ixo^l,r)
1196 res(ixo^s) = gamma_1 * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
1197 end subroutine ffhd_get_temperature_from_eint
1198
1199 subroutine ffhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1201 integer, intent(in) :: ixi^l, ixo^l
1202 double precision, intent(in) :: w(ixi^s, 1:nw)
1203 double precision, intent(in) :: x(ixi^s, 1:ndim)
1204 double precision, intent(out):: res(ixi^s)
1205
1206 double precision :: r(ixi^s)
1207
1208 call ffhd_get_rfactor(w,x,ixi^l,ixo^l,r)
1209 call ffhd_get_pthermal(w,x,ixi^l,ixo^l,res)
1210 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
1211 end subroutine ffhd_get_temperature_from_etot
1212
1213 subroutine ffhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
1215 integer, intent(in) :: ixi^l, ixo^l
1216 double precision, intent(in) :: w(ixi^s,nw)
1217 double precision, intent(in) :: x(ixi^s,1:ndim)
1218 double precision, intent(out) :: csound2(ixi^s)
1219 double precision :: rho(ixi^s)
1220
1221 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
1222 if(ffhd_energy) then
1223 call ffhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
1224 csound2(ixo^s)=ffhd_gamma*csound2(ixo^s)/rho(ixo^s)
1225 else
1226 csound2(ixo^s)=ffhd_gamma*ffhd_adiab*rho(ixo^s)**gamma_1
1227 end if
1228 end subroutine ffhd_get_csound2
1229
1230 subroutine ffhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
1232 use mod_geometry
1233 integer, intent(in) :: ixi^l, ixo^l, idim
1234 ! conservative w
1235 double precision, intent(in) :: wc(ixi^s,nw)
1236 ! primitive w
1237 double precision, intent(in) :: w(ixi^s,nw)
1238 double precision, intent(in) :: x(ixi^s,1:ndim)
1239 double precision,intent(out) :: f(ixi^s,nwflux)
1240 double precision :: ptotal(ixo^s)
1241 double precision :: tmp(ixi^s)
1242 integer :: idirmin, iw, idir, jdir, kdir
1243 double precision, dimension(ixI^S) :: te,tau,sigt
1244
1245 f(ixo^s,rho_)=w(ixo^s,mom(1))*w(ixo^s,rho_)*block%B0(ixo^s,idim,idim)
1246
1247 if(ffhd_energy) then
1248 ptotal(ixo^s)=w(ixo^s,p_)
1249 else
1250 ptotal(ixo^s)=ffhd_adiab*w(ixo^s,rho_)**ffhd_gamma
1251 end if
1252
1253 ! Get flux of momentum
1254 f(ixo^s,mom(1))=(wc(ixo^s,mom(1))*w(ixo^s,mom(1))+ptotal(ixo^s))*block%B0(ixo^s,idim,idim)
1255
1256 ! Get flux of energy
1257 if(ffhd_energy) then
1258 f(ixo^s,e_)=w(ixo^s,mom(1))*(wc(ixo^s,e_)+ptotal(ixo^s))*block%B0(ixo^s,idim,idim)
1260 f(ixo^s,e_)=f(ixo^s,e_)+w(ixo^s,q_)*block%B0(ixo^s,idim,idim)
1261 f(ixo^s,q_)=zero
1262 end if
1263 end if
1264 end subroutine ffhd_get_flux
1265
1266 subroutine ffhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
1271 integer, intent(in) :: ixi^l, ixo^l
1272 double precision, intent(in) :: qdt,dtfactor
1273 double precision, intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:ndim)
1274 double precision, intent(inout) :: w(ixi^s,1:nw)
1275 logical, intent(in) :: qsourcesplit
1276 logical, intent(inout) :: active
1277
1278 if (.not. qsourcesplit) then
1279 active = .true.
1280 call add_punitb(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
1282 call add_hypertc_source(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
1283 end if
1284 end if
1285
1286 if(ffhd_radiative_cooling) then
1287 call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
1288 w,x,qsourcesplit,active, rc_fl)
1289 end if
1290
1291 if(ffhd_viscosity) then
1292 call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
1293 w,x,ffhd_energy,qsourcesplit,active)
1294 end if
1295
1296 if(ffhd_gravity) then
1297 call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
1298 w,x,ffhd_energy,.false.,qsourcesplit,active)
1299 end if
1300
1301 ! update temperature from new pressure, density, and old ionization degree
1303 if(.not.qsourcesplit) then
1304 active = .true.
1305 call ffhd_update_temperature(ixi^l,ixo^l,wct,w,x)
1306 end if
1307 end if
1308 end subroutine ffhd_add_source
1309
1310 subroutine add_punitb(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
1312 use mod_geometry
1313 integer, intent(in) :: ixi^l,ixo^l
1314 double precision, intent(in) :: qdt
1315 double precision, intent(in) :: wct(ixi^s,1:nw),x(ixi^s,1:ndim)
1316 double precision, intent(in) :: wctprim(ixi^s,1:nw)
1317 double precision, intent(inout) :: w(ixi^s,1:nw)
1318
1319 integer :: idims,hxo^l
1320 double precision :: divb(ixi^s)
1321
1322 divb=zero
1323 if(slab_uniform) then
1324 do idims=1,ndim
1325 hxo^l=ixo^l-kr(idims,^d);
1326 divb(ixo^s)=divb(ixo^s)+(block%B0(ixo^s,idims,idims)-block%B0(hxo^s,idims,idims))/dxlevel(idims)
1327 end do
1328 else
1329 call divvector(block%B0(ixi^s,1:ndir,0),ixi^l,ixo^l,divb)
1330 end if
1331 w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*wctprim(ixo^s,p_)*divb(ixo^s)
1332 end subroutine add_punitb
1333
1334 subroutine ffhd_get_rho(w,x,ixI^L,ixO^L,rho)
1336 integer, intent(in) :: ixi^l, ixo^l
1337 double precision, intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
1338 double precision, intent(out) :: rho(ixi^s)
1339
1340 rho(ixo^s) = w(ixo^s,rho_)
1341 end subroutine ffhd_get_rho
1342
1343 subroutine ffhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
1346 integer, intent(in) :: ixi^l,ixo^l, ie
1347 double precision, intent(inout) :: w(ixi^s,1:nw)
1348 double precision, intent(in) :: x(ixi^s,1:ndim)
1349 character(len=*), intent(in) :: subname
1350 integer :: idir
1351 logical :: flag(ixi^s,1:nw)
1352 double precision :: rho(ixi^s)
1353
1354 flag=.false.
1355 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
1356 if(any(flag(ixo^s,ie))) then
1357 select case (small_values_method)
1358 case ("replace")
1359 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
1360 case ("average")
1361 call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
1362 case default
1363 w(ixo^s,e_)=w(ixo^s,e_)*gamma_1
1364 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
1365 w(ixo^s,mom(1)) = w(ixo^s,mom(1))/rho(ixo^s)
1366 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
1367 end select
1368 end if
1369 end subroutine ffhd_handle_small_ei
1370
1371 subroutine ffhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
1374 integer, intent(in) :: ixi^l, ixo^l
1375 double precision, intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
1376 double precision, intent(inout) :: w(ixi^s,1:nw)
1377 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
1378
1379 call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
1380 call ffhd_get_pthermal(w,x,ixi^l,ixo^l,pth)
1381 w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
1382 he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
1383 end subroutine ffhd_update_temperature
1384
1385 subroutine ffhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1387 use mod_usr_methods
1390 use mod_gravity, only: gravity_get_dt
1391 integer, intent(in) :: ixi^l, ixo^l
1392 double precision, intent(inout) :: dtnew
1393 double precision, intent(in) :: dx^d
1394 double precision, intent(in) :: w(ixi^s,1:nw)
1395 double precision, intent(in) :: x(ixi^s,1:ndim)
1396 integer :: idirmin,idim
1397 double precision :: dxarr(ndim)
1398 double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1399
1400 dtnew = bigdouble
1401
1402 if(ffhd_radiative_cooling) then
1403 call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl)
1404 end if
1405
1406 if(ffhd_viscosity) then
1407 call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1408 end if
1409
1410 if(ffhd_gravity) then
1411 call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1412 end if
1413 end subroutine ffhd_get_dt
1414
1415! subroutine ffhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
1416! use mod_global_parameters
1417! use mod_geometry
1418!
1419! integer, intent(in) :: ixI^L, ixO^L
1420! double precision, intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
1421! double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1422!
1423! integer :: iw,idir, h1x^L{^NOONED, h2x^L}
1424! double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
1425!
1426! integer :: mr_,mphi_ ! Polar var. names
1427! integer :: br_,bphi_
1428!
1429! mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1430! br_=mag(1); bphi_=mag(1)-1+phi_
1431!
1432! ! 1/rho
1433! invrho(ixO^S)=1.d0/wCT(ixO^S,rho_)
1434! ! include dt in invr, invr is always used with qdt
1435! if(local_timestep) then
1436! invr(ixO^S) = block%dt(ixO^S) * dtfactor/x(ixO^S,1)
1437! else
1438! invr(ixO^S) = qdt/x(ixO^S,1)
1439! end if
1440!
1441!
1442! select case (coordinate)
1443! case (cylindrical)
1444! call ffhd_get_p_total(wCT,x,ixI^L,ixO^L,tmp)
1445! if(phi_>0) then
1446! w(ixO^S,mr_)=w(ixO^S,mr_)+invr(ixO^S)*(tmp(ixO^S)-&
1447! wCT(ixO^S,bphi_)**2+wCT(ixO^S,mphi_)**2*invrho(ixO^S))
1448! w(ixO^S,mphi_)=w(ixO^S,mphi_)+invr(ixO^S)*(&
1449! -wCT(ixO^S,mphi_)*wCT(ixO^S,mr_)*invrho(ixO^S) &
1450! +wCT(ixO^S,bphi_)*wCT(ixO^S,br_))
1451! if(.not.stagger_grid) then
1452! w(ixO^S,bphi_)=w(ixO^S,bphi_)+invr(ixO^S)*&
1453! (wCT(ixO^S,bphi_)*wCT(ixO^S,mr_) &
1454! -wCT(ixO^S,br_)*wCT(ixO^S,mphi_)) &
1455! *invrho(ixO^S)
1456! end if
1457! else
1458! w(ixO^S,mr_)=w(ixO^S,mr_)+invr(ixO^S)*tmp(ixO^S)
1459! end if
1460! if(ffhd_glm) w(ixO^S,br_)=w(ixO^S,br_)+wCT(ixO^S,psi_)*invr(ixO^S)
1461! case (spherical)
1462! h1x^L=ixO^L-kr(1,^D); {^NOONED h2x^L=ixO^L-kr(2,^D);}
1463! call ffhd_get_p_total(wCT,x,ixI^L,ixO^L,tmp1)
1464! ! m1
1465! tmp(ixO^S)=tmp1(ixO^S)*x(ixO^S,1) &
1466! *(block%surfaceC(ixO^S,1)-block%surfaceC(h1x^S,1))/block%dvolume(ixO^S)
1467! do idir=2,ndir
1468! tmp(ixO^S)=tmp(ixO^S)+wCT(ixO^S,mom(idir))**2*invrho(ixO^S)-wCT(ixO^S,mag(idir))**2
1469! end do
1470! w(ixO^S,mom(1))=w(ixO^S,mom(1))+tmp(ixO^S)*invr(ixO^S)
1471! ! b1
1472! if(ffhd_glm) then
1473! w(ixO^S,mag(1))=w(ixO^S,mag(1))+invr(ixO^S)*2.0d0*wCT(ixO^S,psi_)
1474! end if
1475!
1476! {^NOONED
1477! ! m2
1478! ! This will make hydrostatic p=const an exact solution
1479! if(local_timestep) then
1480! tmp(ixO^S) = block%dt(ixO^S) * tmp1(ixO^S)
1481! else
1482! tmp(ixO^S) = qdt * tmp1(ixO^S)
1483! end if
1484! w(ixO^S,mom(2))=w(ixO^S,mom(2))+tmp(ixO^S) &
1485! *(block%surfaceC(ixO^S,2)-block%surfaceC(h2x^S,2)) &
1486! /block%dvolume(ixO^S)
1487! tmp(ixO^S)=-(wCT(ixO^S,mom(1))*wCT(ixO^S,mom(2))*invrho(ixO^S) &
1488! -wCT(ixO^S,mag(1))*wCT(ixO^S,mag(2)))
1489! if(ndir==3) then
1490! tmp(ixO^S)=tmp(ixO^S)+(wCT(ixO^S,mom(3))**2*invrho(ixO^S) &
1491! -wCT(ixO^S,mag(3))**2)*dcos(x(ixO^S,2))/dsin(x(ixO^S,2))
1492! end if
1493! w(ixO^S,mom(2))=w(ixO^S,mom(2))+tmp(ixO^S)*invr(ixO^S)
1494! ! b2
1495! if(.not.stagger_grid) then
1496! tmp(ixO^S)=(wCT(ixO^S,mom(1))*wCT(ixO^S,mag(2)) &
1497! -wCT(ixO^S,mom(2))*wCT(ixO^S,mag(1)))*invrho(ixO^S)
1498! if(ffhd_glm) then
1499! tmp(ixO^S)=tmp(ixO^S) &
1500! + dcos(x(ixO^S,2))/dsin(x(ixO^S,2))*wCT(ixO^S,psi_)
1501! end if
1502! w(ixO^S,mag(2))=w(ixO^S,mag(2))+tmp(ixO^S)*invr(ixO^S)
1503! end if
1504! }
1505!
1506! if(ndir==3) then
1507! ! m3
1508! tmp(ixO^S)=-(wCT(ixO^S,mom(3))*wCT(ixO^S,mom(1))*invrho(ixO^S) &
1509! -wCT(ixO^S,mag(3))*wCT(ixO^S,mag(1))) {^NOONED &
1510! -(wCT(ixO^S,mom(2))*wCT(ixO^S,mom(3))*invrho(ixO^S) &
1511! -wCT(ixO^S,mag(2))*wCT(ixO^S,mag(3))) &
1512! *dcos(x(ixO^S,2))/dsin(x(ixO^S,2)) }
1513! w(ixO^S,mom(3))=w(ixO^S,mom(3))+tmp(ixO^S)*invr(ixO^S)
1514! ! b3
1515! if(.not.stagger_grid) then
1516! tmp(ixO^S)=(wCT(ixO^S,mom(1))*wCT(ixO^S,mag(3)) &
1517! -wCT(ixO^S,mom(3))*wCT(ixO^S,mag(1)))*invrho(ixO^S) {^NOONED &
1518! -(wCT(ixO^S,mom(3))*wCT(ixO^S,mag(2)) &
1519! -wCT(ixO^S,mom(2))*wCT(ixO^S,mag(3)))*dcos(x(ixO^S,2)) &
1520! /(wCT(ixO^S,rho_)*dsin(x(ixO^S,2))) }
1521! w(ixO^S,mag(3))=w(ixO^S,mag(3))+tmp(ixO^S)*invr(ixO^S)
1522! end if
1523! end if
1524! end select
1525! end subroutine ffhd_add_source_geom
1526
1527 function ffhd_kin_en_origin(w, ixI^L, ixO^L, inv_rho) result(ke)
1528 use mod_global_parameters, only: nw, ndim,block
1529 integer, intent(in) :: ixi^l, ixo^l
1530 double precision, intent(in) :: w(ixi^s, nw)
1531 double precision :: ke(ixo^s)
1532 double precision, intent(in), optional :: inv_rho(ixo^s)
1533
1534 if(present(inv_rho)) then
1535 ke(ixo^s)=0.5d0*w(ixo^s,mom(1))**2*inv_rho(ixo^s)
1536 else
1537 ke(ixo^s)=0.5d0*w(ixo^s,mom(1))**2/w(ixo^s,rho_)
1538 end if
1539 end function ffhd_kin_en_origin
1540
1541 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
1544 integer, intent(in) :: ixi^l, ixo^l
1545 double precision, intent(in) :: w(ixi^s,1:nw)
1546 double precision, intent(in) :: x(ixi^s,1:ndim)
1547 double precision, intent(out):: rfactor(ixi^s)
1548 double precision :: iz_h(ixo^s),iz_he(ixo^s)
1549
1550 call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,te_),iz_h,iz_he)
1551 rfactor(ixo^s)=(1.d0+iz_h(ixo^s)+0.1d0*(1.d0+iz_he(ixo^s)*(1.d0+iz_he(ixo^s))))/(2.d0+3.d0*he_abundance)
1552 end subroutine rfactor_from_temperature_ionization
1553
1554 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
1556 integer, intent(in) :: ixi^l, ixo^l
1557 double precision, intent(in) :: w(ixi^s,1:nw)
1558 double precision, intent(in) :: x(ixi^s,1:ndim)
1559 double precision, intent(out):: rfactor(ixi^s)
1560
1561 rfactor(ixo^s)=rr
1562 end subroutine rfactor_from_constant_ionization
1563
1564 subroutine get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
1566 integer, intent(in) :: ixi^l, ixo^l
1567 double precision, dimension(ixI^S,1:nw), intent(in) :: w
1568 double precision, dimension(ixI^S), intent(in) :: te
1569 double precision, dimension(ixI^S), intent(out) :: tau,sigt5
1570 integer :: ix^d
1571 double precision :: dxmin,taumin
1572 double precision, dimension(ixI^S) :: sigt7,eint
1573
1574 taumin=4.d0
1575 !> w supposed to be wCTprim here
1576 if(ffhd_trac) then
1577 where(te(ixo^s) .lt. block%wextra(ixo^s,tcoff_))
1578 sigt5(ixo^s)=hypertc_kappa*sqrt(block%wextra(ixo^s,tcoff_)**5)
1579 sigt7(ixo^s)=sigt5(ixo^s)*block%wextra(ixo^s,tcoff_)
1580 else where
1581 sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
1582 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
1583 end where
1584 else
1585 sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
1586 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
1587 end if
1588 eint(ixo^s)=w(ixo^s,p_)/(ffhd_gamma-one)
1589 tau(ixo^s)=max(taumin*dt,sigt7(ixo^s)/eint(ixo^s)/cs2max_global)
1590 end subroutine get_tau
1591
1592 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
1594 integer, intent(in) :: ixi^l,ixo^l
1595 double precision, intent(in) :: qdt
1596 double precision, dimension(ixI^S,1:ndim), intent(in) :: x
1597 double precision, dimension(ixI^S,1:nw), intent(in) :: wct,wctprim
1598 double precision, dimension(ixI^S,1:nw), intent(inout) :: w
1599 integer :: idims
1600 integer :: hxc^l,hxo^l,ixc^l,jxc^l,jxo^l,kxc^l
1601 double precision :: invdx
1602 double precision, dimension(ixI^S) :: te,tau,sigt,htc_qsrc,tface
1603 double precision, dimension(ixI^S) :: htc_esrc
1604
1605 te(ixi^s)=wctprim(ixi^s,p_)/wct(ixi^s,rho_)
1606 call get_tau(ixi^l,ixo^l,wctprim,te,tau,sigt)
1607 htc_qsrc=zero
1608 do idims=1,ndim
1609 invdx=1.d0/dxlevel(idims)
1610 ixc^l=ixo^l;
1611 ixcmin^d=ixomin^d-kr(idims,^d);ixcmax^d=ixomax^d;
1612 jxc^l=ixc^l+kr(idims,^d);
1613 kxc^l=jxc^l+kr(idims,^d);
1614 hxc^l=ixc^l-kr(idims,^d);
1615 hxo^l=ixo^l-kr(idims,^d);
1616 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
1617 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*block%B0(ixo^s,idims,0)*(tface(ixo^s)-tface(hxo^s))*invdx
1618 end do
1619 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,q_))/tau(ixo^s)
1620 w(ixo^s,q_)=w(ixo^s,q_)-qdt*htc_qsrc(ixo^s)
1621 end subroutine add_hypertc_source
1622end module mod_ffhd_phys
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter bigdouble
A very large real number.
subroutine add_convert_method(phys_convert_vars, nwc, dataset_names, file_suffix)
Definition mod_convert.t:59
Frozen-field hydrodynamics module.
integer, public, protected te_
Indices of temperature.
integer, public, protected ffhd_trac_type
Which TRAC method is used.
double precision, public hypertc_kappa
The thermal conductivity kappa in hyperbolic thermal conduction.
integer, public, protected e_
Index of the energy density (-1 if not present)
double precision, public ffhd_gamma
The adiabatic index.
logical, public, protected eq_state_units
double precision, public ffhd_adiab
The adiabatic constant.
procedure(sub_get_pthermal), pointer, public ffhd_get_temperature
logical, public, protected ffhd_hyperbolic_thermal_conduction
Whether hyperbolic type thermal conduction is used.
type(rc_fluid), allocatable, public rc_fl
type of fluid for radiative cooling
integer, public, protected rho_
Index of the density (in the w array)
integer, public, protected tcoff_
Index of the cutoff temperature for the TRAC method.
procedure(sub_get_pthermal), pointer, public ffhd_get_pthermal
subroutine, public ffhd_get_rho(w, x, ixil, ixol, rho)
logical, public, protected ffhd_partial_ionization
Whether plasma is partially ionized.
procedure(sub_convert), pointer, public ffhd_to_conserved
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
double precision, public, protected h_ion_fr
Ionization fraction of H H_ion_fr = H+/(H+ + H)
procedure(fun_kin_en), pointer, public ffhd_kin_en
subroutine, public ffhd_phys_init()
procedure(sub_get_v), pointer, public ffhd_get_v
subroutine, public ffhd_get_csound2(w, x, ixil, ixol, csound2)
logical, public, protected ffhd_energy
Whether an energy equation is used.
double precision, public, protected rr
type(tc_fluid), allocatable, public tc_fl
type of fluid for thermal conduction
type(te_fluid), allocatable, public te_fl_ffhd
type of fluid for thermal emission synthesis
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
logical, public, protected ffhd_viscosity
Whether viscosity is added.
logical, public, protected ffhd_radiative_cooling
Whether radiative cooling is added.
subroutine, public ffhd_get_v_idim(w, x, ixil, ixol, idim, v)
integer, public, protected q_
procedure(sub_convert), pointer, public ffhd_to_primitive
integer, public, protected tweight_
double precision, public, protected he_ion_fr2
Ratio of number He2+ / number He+ + He2+ He_ion_fr2 = He2+/(He2+ + He+)
logical, public, protected ffhd_trac
Whether TRAC method is used.
logical, public, protected ffhd_thermal_conduction
Whether thermal conduction is used.
subroutine, public ffhd_ei_to_e(ixil, ixol, w, x)
logical, public, protected ffhd_gravity
Whether gravity is added.
integer, public, protected p_
Index of the gas pressure (-1 if not present) should equal e_.
double precision, public, protected ffhd_trac_mask
Height of the mask used in the TRAC method.
double precision, public, protected he_ion_fr
Ionization fraction of He He_ion_fr = (He2+ + He+)/(He2+ + He+ + He)
integer, public, protected ffhd_trac_finegrid
Distance between two adjacent traced magnetic field lines (in finest cell size)
subroutine, public ffhd_e_to_ei(ixil, ixol, w, x)
Module for flux conservation near refinement boundaries.
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
subroutine divvector(qvec, ixil, ixol, divq, nth_in)
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision unit_numberdensity
Physical scaling factor for number density.
character(len=std_len) convert_type
Which format to use when converting.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
double precision dt
global time step
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision unit_velocity
Physical scaling factor for velocity.
logical b0field
split magnetic field as background B0 field
logical need_global_cs2max
global value for csound speed
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
logical phys_trac
Use TRAC for MHD or 1D HD.
logical fix_small_values
fix small values with average or replace methods
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
double precision cs2max_global
global largest cs2 for hyperbolic thermal conduction
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
integer, parameter unitconvert
integer number_equi_vars
number of equilibrium set variables, besides the mag field
Module for including gravity in (magneto)hydrodynamics simulations.
Definition mod_gravity.t:2
subroutine gravity_get_dt(w, ixil, ixol, dtnew, dxd, x)
Definition mod_gravity.t:87
subroutine gravity_init()
Initialize the module.
Definition mod_gravity.t:26
subroutine gravity_add_source(qdt, ixil, ixol, wct, wctprim, w, x, energy, rhov, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
Definition mod_gravity.t:43
module ionization degree - get ionization degree for given temperature
subroutine ionization_degree_from_temperature(ixil, ixol, te, iz_h, iz_he)
Module containing all the particle routines.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
module radiative cooling – add optically thin radiative cooling for HD and MHD
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
subroutine radiative_cooling_init(fl, read_params)
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Module for handling problematic values in simulations, such as negative pressures.
subroutine, public small_values_average(ixil, ixol, w, x, w_flag, windex)
logical, public trace_small_values
trace small values in the source file using traceback flag of compiler
subroutine, public small_values_error(wprim, x, ixil, ixol, w_flag, subname)
logical, dimension(:), allocatable, public small_values_fix_iw
Whether to apply small value fixes to certain variables.
character(len=20), public small_values_method
How to handle small values.
Generic supertimestepping method 1) in amrvac.par in sts_list set the following parameters which have...
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...
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...
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....
subroutine, public sts_init()
Initialize sts module.
Thermal conduction for HD and MHD or RHD and RMHD or twofl (plasma-neutral) module Adaptation of mod_...
subroutine, public tc_get_hd_params(fl, read_hd_params)
Init TC coefficients: HD case.
double precision function, public get_tc_dt_mhd(w, ixil, ixol, dxd, x, fl)
Get the explicut timestep for the TC (mhd implementation)
subroutine tc_init_params(phys_gamma)
subroutine, public sts_set_source_tc_mhd(ixil, ixol, w, x, wres, fix_conserve_at_step, my_dt, igrid, nflux, fl)
anisotropic thermal conduction with slope limited symmetric scheme Sharma 2007 Journal of Computation...
subroutine get_euv_image(qunit, fl)
subroutine get_sxr_image(qunit, fl)
subroutine get_euv_spectrum(qunit, fl)
subroutine get_whitelight_image(qunit, fl)
Module with all the methods that users can customize in AMRVAC.
procedure(rfactor), pointer usr_rfactor
procedure(set_equi_vars), pointer usr_set_equi_vars
The module add viscous source terms and check time step.
subroutine viscosity_add_source(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
subroutine viscosity_init(phys_wider_stencil)
Initialize the module.
subroutine viscosity_get_dt(w, ixil, ixol, dtnew, dxd, x)