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 ! Number of variables need reconstruction in w
285 nw_recon=nwflux
286
287 ! set temperature as an auxiliary variable to get ionization degree
289 te_ = var_set_auxvar('Te','Te')
290 else
291 te_ = -1
292 end if
293
294 ! set cutoff temperature when using the TRAC method, as well as an auxiliary weight
295 tweight_ = -1
296 if(ffhd_trac) then
297 tcoff_ = var_set_wextra()
298 iw_tcoff=tcoff_
299 if(ffhd_trac_type .ge. 3) then
300 tweight_ = var_set_wextra()
301 iw_tweight=tweight_
302 end if
303 else
304 tcoff_ = -1
305 end if
306
307 nvector = 0 ! No. vector vars
308
309 ! Check whether custom flux types have been defined
310 if(.not. allocated(flux_type)) then
311 allocate(flux_type(ndir, nwflux))
312 flux_type = flux_default
313 else if(any(shape(flux_type) /= [ndir, nwflux])) then
314 call mpistop("phys_check error: flux_type has wrong shape")
315 end if
316
317 phys_get_dt => ffhd_get_dt
318 phys_get_cmax => ffhd_get_cmax_origin
319 phys_get_a2max => ffhd_get_a2max
320 phys_get_cs2max => ffhd_get_cs2max
321 phys_get_tcutoff => ffhd_get_tcutoff
322 phys_get_cbounds => ffhd_get_cbounds
323 phys_to_primitive => ffhd_to_primitive_origin
324 ffhd_to_primitive => ffhd_to_primitive_origin
325 phys_to_conserved => ffhd_to_conserved_origin
326 ffhd_to_conserved => ffhd_to_conserved_origin
327 phys_get_flux => ffhd_get_flux
328 phys_get_v => ffhd_get_v_origin
329 ffhd_get_v => ffhd_get_v_origin
330 phys_get_rho => ffhd_get_rho
331 ffhd_kin_en => ffhd_kin_en_origin
332 !> need to check source geom here
333 !phys_add_source_geom => ffhd_add_source_geom
334 phys_add_source => ffhd_add_source
335 phys_check_params => ffhd_check_params
336 phys_write_info => ffhd_write_info
337 phys_handle_small_values => ffhd_handle_small_values_origin
338 ffhd_handle_small_values => ffhd_handle_small_values_origin
339 phys_check_w => ffhd_check_w_origin
340
341 if(.not.ffhd_energy) then
342 phys_get_pthermal => ffhd_get_pthermal_iso
343 ffhd_get_pthermal => ffhd_get_pthermal_iso
344 else
345 phys_get_pthermal => ffhd_get_pthermal_origin
346 ffhd_get_pthermal => ffhd_get_pthermal_origin
347 end if
348
349 ! choose Rfactor in ideal gas law
351 ffhd_get_rfactor=>rfactor_from_temperature_ionization
352 phys_update_temperature => ffhd_update_temperature
353 else if(associated(usr_rfactor)) then
354 ffhd_get_rfactor=>usr_rfactor
355 else
356 ffhd_get_rfactor=>rfactor_from_constant_ionization
357 end if
358
360 ffhd_get_temperature => ffhd_get_temperature_from_te
361 else
362 ffhd_get_temperature => ffhd_get_temperature_from_etot
363 end if
364
365 ! derive units from basic units
366 call ffhd_physical_units()
367
370 end if
371 if(.not. ffhd_energy .and. ffhd_thermal_conduction) then
372 call mpistop("thermal conduction needs ffhd_energy=T")
373 end if
375 call mpistop("hyperbolic thermal conduction needs ffhd_energy=T")
376 end if
377 if(.not. ffhd_energy .and. ffhd_radiative_cooling) then
378 call mpistop("radiative cooling needs ffhd_energy=T")
379 end if
380
381 ! initialize thermal conduction module
383 call sts_init()
385
386 allocate(tc_fl)
387 call tc_get_hd_params(tc_fl,tc_params_read_ffhd)
388 call add_sts_method(ffhd_get_tc_dt_ffhd,ffhd_sts_set_source_tc_ffhd,e_,1,e_,1,.false.)
389 tc_fl%get_temperature_from_conserved => ffhd_get_temperature_from_etot
390 tc_fl%get_temperature_from_eint => ffhd_get_temperature_from_eint
392 call set_error_handling_to_head(ffhd_tc_handle_small_e)
393 tc_fl%get_rho => ffhd_get_rho
394 tc_fl%e_ = e_
395 tc_fl%Tcoff_ = tcoff_
396 end if
397
398 ! Initialize radiative cooling module
401 allocate(rc_fl)
402 call radiative_cooling_init(rc_fl,rc_params_read)
403 rc_fl%get_rho => ffhd_get_rho
404 rc_fl%get_pthermal => ffhd_get_pthermal
405 rc_fl%get_var_Rfactor => ffhd_get_rfactor
406 rc_fl%e_ = e_
407 rc_fl%Tcoff_ = tcoff_
408 rc_fl%has_equi = .false.
409 end if
410 allocate(te_fl_ffhd)
411 te_fl_ffhd%get_rho=> ffhd_get_rho
412 te_fl_ffhd%get_pthermal=> ffhd_get_pthermal
413 te_fl_ffhd%get_var_Rfactor => ffhd_get_rfactor
414{^ifthreed
415 phys_te_images => ffhd_te_images
416}
417 ! Initialize viscosity module
418 if(ffhd_viscosity) call viscosity_init(phys_wider_stencil)
419
420 ! Initialize gravity module
421 if(ffhd_gravity) then
422 call gravity_init()
423 end if
424
425 ! initialize ionization degree table
427 end subroutine ffhd_phys_init
428
429{^ifthreed
430 subroutine ffhd_te_images
433
434 select case(convert_type)
435 case('EIvtiCCmpi','EIvtuCCmpi')
437 case('ESvtiCCmpi','ESvtuCCmpi')
439 case('SIvtiCCmpi','SIvtuCCmpi')
441 case('WIvtiCCmpi','WIvtuCCmpi')
443 case default
444 call mpistop("Error in synthesize emission: Unknown convert_type")
445 end select
446 end subroutine ffhd_te_images
447}
448
449 subroutine ffhd_sts_set_source_tc_ffhd(ixI^L,ixO^L,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux)
453 integer, intent(in) :: ixi^l, ixo^l, igrid, nflux
454 double precision, intent(in) :: x(ixi^s,1:ndim)
455 double precision, intent(inout) :: wres(ixi^s,1:nw), w(ixi^s,1:nw)
456 double precision, intent(in) :: my_dt
457 logical, intent(in) :: fix_conserve_at_step
458 call sts_set_source_tc_mhd(ixi^l,ixo^l,w,x,wres,fix_conserve_at_step,my_dt,igrid,nflux,tc_fl)
459 end subroutine ffhd_sts_set_source_tc_ffhd
460
461 function ffhd_get_tc_dt_ffhd(w,ixI^L,ixO^L,dx^D,x) result(dtnew)
462 !Check diffusion time limit dt < dx_i**2/((gamma-1)*tc_k_para_i/rho)
463 !where tc_k_para_i=tc_k_para*B_i**2/B**2
464 !and T=p/rho
467
468 integer, intent(in) :: ixi^l, ixo^l
469 double precision, intent(in) :: dx^d, x(ixi^s,1:ndim)
470 double precision, intent(in) :: w(ixi^s,1:nw)
471 double precision :: dtnew
472
473 dtnew=get_tc_dt_mhd(w,ixi^l,ixo^l,dx^d,x,tc_fl)
474 end function ffhd_get_tc_dt_ffhd
475
476 subroutine ffhd_tc_handle_small_e(w, x, ixI^L, ixO^L, step)
478
479 integer, intent(in) :: ixi^l,ixo^l
480 double precision, intent(inout) :: w(ixi^s,1:nw)
481 double precision, intent(in) :: x(ixi^s,1:ndim)
482 integer, intent(in) :: step
483 character(len=140) :: error_msg
484
485 write(error_msg,"(a,i3)") "Thermal conduction step ", step
486 call ffhd_handle_small_ei(w,x,ixi^l,ixo^l,e_,error_msg)
487 end subroutine ffhd_tc_handle_small_e
488
489 subroutine tc_params_read_ffhd(fl)
491 type(tc_fluid), intent(inout) :: fl
492 integer :: n
493 ! list parameters
494 logical :: tc_saturate=.false.
495 double precision :: tc_k_para=0d0
496 character(len=std_len) :: tc_slope_limiter="MC"
497
498 namelist /tc_list/ tc_saturate, tc_slope_limiter, tc_k_para
499
500 do n = 1, size(par_files)
501 open(unitpar, file=trim(par_files(n)), status="old")
502 read(unitpar, tc_list, end=111)
503111 close(unitpar)
504 end do
505
506 fl%tc_saturate = tc_saturate
507 fl%tc_k_para = tc_k_para
508 select case(tc_slope_limiter)
509 case ('no','none')
510 fl%tc_slope_limiter = 0
511 case ('MC')
512 ! montonized central limiter Woodward and Collela limiter (eq.3.51h), a factor of 2 is pulled out
513 fl%tc_slope_limiter = 1
514 case('minmod')
515 ! minmod limiter
516 fl%tc_slope_limiter = 2
517 case ('superbee')
518 ! Roes superbee limiter (eq.3.51i)
519 fl%tc_slope_limiter = 3
520 case ('koren')
521 ! Barry Koren Right variant
522 fl%tc_slope_limiter = 4
523 case default
524 call mpistop("Unknown tc_slope_limiter, choose MC, minmod")
525 end select
526 end subroutine tc_params_read_ffhd
527
528 subroutine rc_params_read(fl)
530 use mod_constants, only: bigdouble
531 type(rc_fluid), intent(inout) :: fl
532 integer :: n
533 integer :: ncool = 4000
534 double precision :: cfrac=0.1d0
535
536 !> Name of cooling curve
537 character(len=std_len) :: coolcurve='JCcorona'
538
539 !> Name of cooling method
540 character(len=std_len) :: coolmethod='exact'
541
542 !> Fixed temperature not lower than tlow
543 logical :: tfix=.false.
544
545 !> Lower limit of temperature
546 double precision :: tlow=bigdouble
547
548 !> Add cooling source in a split way (.true.) or un-split way (.false.)
549 logical :: rc_split=.false.
550 logical :: rad_cut=.false.
551 double precision :: rad_cut_hgt=0.5d0
552 double precision :: rad_cut_dey=0.15d0
553
554 namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix, rc_split, rad_cut, rad_cut_hgt, rad_cut_dey
555
556 do n = 1, size(par_files)
557 open(unitpar, file=trim(par_files(n)), status="old")
558 read(unitpar, rc_list, end=111)
559111 close(unitpar)
560 end do
561
562 fl%ncool=ncool
563 fl%coolcurve=coolcurve
564 fl%coolmethod=coolmethod
565 fl%tlow=tlow
566 fl%Tfix=tfix
567 fl%rc_split=rc_split
568 fl%cfrac=cfrac
569 fl%rad_cut=rad_cut
570 fl%rad_cut_hgt=rad_cut_hgt
571 fl%rad_cut_dey=rad_cut_dey
572 end subroutine rc_params_read
573
574 subroutine ffhd_check_params
578
579 gamma_1=ffhd_gamma-1.d0
580 if (.not. ffhd_energy) then
581 if (ffhd_gamma <= 0.0d0) call mpistop ("Error: ffhd_gamma <= 0")
582 if (ffhd_adiab < 0.0d0) call mpistop ("Error: ffhd_adiab < 0")
584 else
585 if (ffhd_gamma <= 0.0d0 .or. ffhd_gamma == 1.0d0) &
586 call mpistop ("Error: ffhd_gamma <= 0 or ffhd_gamma == 1")
587 inv_gamma_1=1.d0/gamma_1
588 small_e = small_pressure * inv_gamma_1
589 end if
590
591 if (number_equi_vars > 0 .and. .not. associated(usr_set_equi_vars)) then
592 call mpistop("usr_set_equi_vars has to be implemented in the user file")
593 end if
594 end subroutine ffhd_check_params
595
596 subroutine ffhd_physical_units()
598 double precision :: mp,kb
599 double precision :: a,b
600
601 if(si_unit) then
602 mp=mp_si
603 kb=kb_si
604 else
605 mp=mp_cgs
606 kb=kb_cgs
607 end if
608 if(eq_state_units) then
609 a = 1d0 + 4d0 * he_abundance
611 b = 2+.3d0
612 else
613 b = 1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0)
614 end if
615 rr = 1d0
616 else
617 a = 1d0
618 b = 1d0
619 rr = (1d0 + h_ion_fr + he_abundance*(he_ion_fr*(he_ion_fr2 + 1d0)+1d0))/(1d0 + 4d0 * he_abundance)
620 end if
621 if(unit_density/=1.d0) then
623 else
625 end if
626 if(unit_velocity/=1.d0) then
629 else if(unit_pressure/=1.d0) then
632 else
635 end if
636 if(unit_time/=1.d0) then
638 else
640 end if
642 end subroutine ffhd_physical_units
643
644 subroutine ffhd_check_w_origin(primitive,ixI^L,ixO^L,w,flag)
646 logical, intent(in) :: primitive
647 integer, intent(in) :: ixi^l, ixo^l
648 double precision, intent(in) :: w(ixi^s,nw)
649 double precision :: tmp(ixi^s)
650 logical, intent(inout) :: flag(ixi^s,1:nw)
651
652 flag=.false.
653 where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
654
655 if(ffhd_energy) then
656 if(primitive) then
657 where(w(ixo^s,e_) < small_pressure) flag(ixo^s,e_) = .true.
658 else
659 tmp(ixo^s)=w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l)
660 where(tmp(ixo^s) < small_e) flag(ixo^s,e_) = .true.
661 end if
662 end if
663 end subroutine ffhd_check_w_origin
664
665 subroutine ffhd_to_conserved_origin(ixI^L,ixO^L,w,x)
667 integer, intent(in) :: ixi^l, ixo^l
668 double precision, intent(inout) :: w(ixi^s, nw)
669 double precision, intent(in) :: x(ixi^s, 1:ndim)
670 double precision :: inv_gamma2(ixo^s)
671 integer :: idir
672
673 if(ffhd_energy) then
674 w(ixo^s,e_)=w(ixo^s,p_)*inv_gamma_1+half*w(ixo^s,mom(1))**2*w(ixo^s,rho_)
675 end if
676 w(ixo^s,mom(1))=w(ixo^s,rho_)*w(ixo^s,mom(1))
677 end subroutine ffhd_to_conserved_origin
678
679 subroutine ffhd_to_primitive_origin(ixI^L,ixO^L,w,x)
681 integer, intent(in) :: ixi^l, ixo^l
682 double precision, intent(inout) :: w(ixi^s, nw)
683 double precision, intent(in) :: x(ixi^s, 1:ndim)
684 double precision :: inv_rho(ixo^s), gamma2(ixo^s)
685
686 if(fix_small_values) then
687 !> fix small values preventing NaN numbers in the following converting
688 call ffhd_handle_small_values(.false., w, x, ixi^l, ixo^l, 'ffhd_to_primitive_origin')
689 end if
690
691 w(ixo^s,mom(1)) = w(ixo^s,mom(1))/w(ixo^s,rho_)
692 if(ffhd_energy) then
693 w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-half*w(ixo^s,rho_)*w(ixo^s,mom(1))**2)
694 end if
695 end subroutine ffhd_to_primitive_origin
696
697 subroutine ffhd_ei_to_e(ixI^L,ixO^L,w,x)
699 integer, intent(in) :: ixi^l, ixo^l
700 double precision, intent(inout) :: w(ixi^s, nw)
701 double precision, intent(in) :: x(ixi^s, 1:ndim)
702
703 w(ixi^s,e_)=w(ixi^s,e_)+ffhd_kin_en(w,ixi^l,ixi^l)
704 end subroutine ffhd_ei_to_e
705
706 subroutine ffhd_e_to_ei(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
712 w(ixi^s,e_)=w(ixi^s,e_)-ffhd_kin_en(w,ixi^l,ixi^l)
713 if(fix_small_values) then
714 call ffhd_handle_small_ei(w,x,ixi^l,ixi^l,e_,'ffhd_e_to_ei')
715 end if
716 end subroutine ffhd_e_to_ei
717
718 subroutine ffhd_handle_small_values_origin(primitive, w, x, ixI^L, ixO^L, subname)
721 logical, intent(in) :: primitive
722 integer, intent(in) :: ixi^l,ixo^l
723 double precision, intent(inout) :: w(ixi^s,1:nw)
724 double precision, intent(in) :: x(ixi^s,1:ndim)
725 character(len=*), intent(in) :: subname
726
727 logical :: flag(ixi^s,1:nw)
728 double precision :: tmp2(ixi^s)
729
730 call phys_check_w(primitive, ixi^l, ixi^l, w, flag)
731
732 if(any(flag)) then
733 select case (small_values_method)
734 case ("replace")
735 where(flag(ixo^s,rho_)) w(ixo^s,rho_) = small_density
736 if(small_values_fix_iw(mom(1))) then
737 where(flag(ixo^s,rho_)) w(ixo^s, mom(1)) = 0.0d0
738 end if
739 if(ffhd_energy) then
740 if(primitive) then
741 where(flag(ixo^s,e_)) w(ixo^s,p_) = small_pressure
742 else
743 where(flag(ixo^s,e_))
744 w(ixo^s,e_) = small_e+ffhd_kin_en(w,ixi^l,ixo^l)
745 end where
746 end if
747 end if
748 case ("average")
749 call small_values_average(ixi^l, ixo^l, w, x, flag, rho_)
750 if(ffhd_energy) then
751 if(primitive) then
752 call small_values_average(ixi^l, ixo^l, w, x, flag, p_)
753 else
754 w(ixi^s,e_)=w(ixi^s,e_)-ffhd_kin_en(w,ixi^l,ixi^l)
755 call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
756 w(ixi^s,e_)=w(ixi^s,e_)+ffhd_kin_en(w,ixi^l,ixi^l)
757 end if
758 end if
759 case default
760 if(.not.primitive) then
761 if(ffhd_energy) then
762 w(ixo^s,p_)=gamma_1*(w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l))
763 end if
764 w(ixo^s,mom(1))=w(ixo^s,mom(1))/w(ixo^s,rho_)
765 end if
766 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
767 end select
768 end if
769 end subroutine ffhd_handle_small_values_origin
770
771 subroutine ffhd_get_v_origin(w,x,ixI^L,ixO^L,v)
773 integer, intent(in) :: ixi^l, ixo^l
774 double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
775 double precision, intent(out) :: v(ixi^s,ndir)
776 double precision :: rho(ixi^s)
777 integer :: idir
778
779 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
780 rho(ixo^s)=1.d0/rho(ixo^s)
781 do idir=1,ndir
782 v(ixo^s,ndir) = w(ixo^s,mom(1))*block%B0(ixo^s,idir,0)*rho(ixo^s)
783 end do
784 end subroutine ffhd_get_v_origin
785
786 subroutine ffhd_get_v_idim(w,x,ixI^L,ixO^L,idim,v)
788 integer, intent(in) :: ixi^l, ixo^l, idim
789 double precision, intent(in) :: w(ixi^s,nw), x(ixi^s,1:ndim)
790 double precision, intent(out) :: v(ixi^s)
791 double precision :: rho(ixi^s)
792
793 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
794 v(ixo^s) = (w(ixo^s, mom(1))*block%B0(ixo^s,idim,0)) / rho(ixo^s)
795 end subroutine ffhd_get_v_idim
796
797 subroutine ffhd_get_cmax_origin(w,x,ixI^L,ixO^L,idim,cmax)
799 integer, intent(in) :: ixi^l, ixo^l, idim
800 ! w in primitive form
801 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
802 double precision, intent(inout) :: cmax(ixi^s)
803
804 if(ffhd_energy) then
805 cmax(ixo^s)=sqrt(ffhd_gamma*w(ixo^s,p_)/w(ixo^s,rho_))*abs(block%B0(ixo^s,idim,idim))
806 else
807 cmax(ixo^s)=sqrt(ffhd_gamma*ffhd_adiab*w(ixo^s,rho_)**gamma_1)*abs(block%B0(ixo^s,idim,idim))
808 end if
809 cmax(ixo^s)=abs(w(ixo^s,mom(1))*block%B0(ixo^s,idim,0))+cmax(ixo^s)
810
811 end subroutine ffhd_get_cmax_origin
812
813 subroutine ffhd_get_cs2max(w,x,ixI^L,ixO^L,cs2max)
815 integer, intent(in) :: ixi^l, ixo^l
816 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
817 double precision, intent(inout) :: cs2max
818 double precision :: cs2(ixi^s)
819
820 call ffhd_get_csound2(w,x,ixi^l,ixo^l,cs2)
821 cs2max=maxval(cs2(ixo^s))
822 end subroutine ffhd_get_cs2max
823
824 subroutine ffhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
826 integer, intent(in) :: ixi^l, ixo^l
827 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
828 double precision, intent(inout) :: a2max(ndim)
829 double precision :: a2(ixi^s,ndim,nw)
830 integer :: gxo^l,hxo^l,jxo^l,kxo^l,i,j
831
832 a2=zero
833 do i = 1,ndim
834 !> 4th order
835 hxo^l=ixo^l-kr(i,^d);
836 gxo^l=hxo^l-kr(i,^d);
837 jxo^l=ixo^l+kr(i,^d);
838 kxo^l=jxo^l+kr(i,^d);
839 a2(ixo^s,i,1:nw)=abs(-w(kxo^s,1:nw)+16.d0*w(jxo^s,1:nw)&
840 -30.d0*w(ixo^s,1:nw)+16.d0*w(hxo^s,1:nw)-w(gxo^s,1:nw))
841 a2max(i)=maxval(a2(ixo^s,i,1:nw))/12.d0/dxlevel(i)**2
842 end do
843 end subroutine ffhd_get_a2max
844
845 subroutine ffhd_get_tcutoff(ixI^L,ixO^L,w,x,Tco_local,Tmax_local)
847 use mod_geometry
848 integer, intent(in) :: ixi^l,ixo^l
849 double precision, intent(in) :: x(ixi^s,1:ndim)
850 double precision, intent(inout) :: w(ixi^s,1:nw)
851 double precision, intent(out) :: tco_local,tmax_local
852 double precision, parameter :: trac_delta=0.25d0
853 double precision :: tmp1(ixi^s),te(ixi^s),lts(ixi^s)
854 double precision, dimension(ixI^S,1:ndir) :: bunitvec
855 double precision, dimension(ixI^S,1:ndim) :: gradt
856 double precision :: bdir(ndim)
857 double precision :: ltrc,ltrp,altr(ixi^s)
858 integer :: idims,jxo^l,hxo^l,ixa^d,ixb^d
859 integer :: jxp^l,hxp^l,ixp^l,ixq^l
860 logical :: lrlt(ixi^s)
861
862 call ffhd_get_temperature(w,x,ixi^l,ixi^l,te)
863 tco_local=zero
864 tmax_local=maxval(te(ixo^s))
865
866 {^ifoned
867 select case(ffhd_trac_type)
868 case(0)
869 !> test case, fixed cutoff temperature
870 block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
871 case(1)
872 hxo^l=ixo^l-1;
873 jxo^l=ixo^l+1;
874 lts(ixo^s)=0.5d0*abs(te(jxo^s)-te(hxo^s))/te(ixo^s)
875 lrlt=.false.
876 where(lts(ixo^s) > trac_delta)
877 lrlt(ixo^s)=.true.
878 end where
879 if(any(lrlt(ixo^s))) then
880 tco_local=maxval(te(ixo^s), mask=lrlt(ixo^s))
881 end if
882 case(2)
883 !> iijima et al. 2021, LTRAC method
884 ltrc=1.5d0
885 ltrp=4.d0
886 ixp^l=ixo^l^ladd1;
887 hxo^l=ixo^l-1;
888 jxo^l=ixo^l+1;
889 hxp^l=ixp^l-1;
890 jxp^l=ixp^l+1;
891 lts(ixp^s)=0.5d0*abs(te(jxp^s)-te(hxp^s))/te(ixp^s)
892 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
893 lts(ixo^s)=0.25d0*(lts(jxo^s)+two*lts(ixo^s)+lts(hxo^s))
894 block%wextra(ixo^s,tcoff_)=te(ixo^s)*lts(ixo^s)**0.4d0
895 case default
896 call mpistop("ffhd_trac_type not allowed for 1D simulation")
897 end select
898 }
899 {^nooned
900 select case(ffhd_trac_type)
901 case(0)
902 !> test case, fixed cutoff temperature
903 if(slab_uniform) then
904 !> assume cgs units
905 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)
906 else
907 block%wextra(ixi^s,tcoff_)=2.5d5/unit_temperature
908 end if
909 case(1,4,6)
910 do idims=1,ndim
911 call gradient(te,ixi^l,ixo^l,idims,tmp1)
912 gradt(ixo^s,idims)=tmp1(ixo^s)
913 end do
914 bunitvec(ixo^s,:)=block%B0(ixo^s,:,0)
915 if(ffhd_trac_type .gt. 1) then
916 ! B direction at cell center
917 bdir=zero
918 {do ixa^d=0,1\}
919 ixb^d=(ixomin^d+ixomax^d-1)/2+ixa^d;
920 bdir(1:ndim)=bdir(1:ndim)+bunitvec(ixb^d,1:ndim)
921 {end do\}
922 if(sum(bdir(:)**2) .gt. zero) then
923 bdir(1:ndim)=bdir(1:ndim)/dsqrt(sum(bdir(:)**2))
924 end if
925 block%special_values(3:ndim+2)=bdir(1:ndim)
926 end if
927 tmp1(ixo^s)=dsqrt(sum(bunitvec(ixo^s,:)**2,dim=ndim+1))
928 where(tmp1(ixo^s)/=0.d0)
929 tmp1(ixo^s)=1.d0/tmp1(ixo^s)
930 else where
931 tmp1(ixo^s)=bigdouble
932 end where
933 ! b unit vector: magnetic field direction vector
934 do idims=1,ndim
935 bunitvec(ixo^s,idims)=bunitvec(ixo^s,idims)*tmp1(ixo^s)
936 end do
937 ! temperature length scale inversed
938 lts(ixo^s)=abs(sum(gradt(ixo^s,1:ndim)*bunitvec(ixo^s,1:ndim),dim=ndim+1))/te(ixo^s)
939 ! fraction of cells size to temperature length scale
940 if(slab_uniform) then
941 lts(ixo^s)=minval(dxlevel)*lts(ixo^s)
942 else
943 lts(ixo^s)=minval(block%ds(ixo^s,:),dim=ndim+1)*lts(ixo^s)
944 end if
945 lrlt=.false.
946 where(lts(ixo^s) > trac_delta)
947 lrlt(ixo^s)=.true.
948 end where
949 if(any(lrlt(ixo^s))) then
950 block%special_values(1)=maxval(te(ixo^s), mask=lrlt(ixo^s))
951 else
952 block%special_values(1)=zero
953 end if
954 block%special_values(2)=tmax_local
955 case(2)
956 !> iijima et al. 2021, LTRAC method
957 ltrc=1.5d0
958 ltrp=4.d0
959 ixp^l=ixo^l^ladd2;
960 do idims=1,ndim
961 ixq^l=ixp^l;
962 hxp^l=ixp^l;
963 jxp^l=ixp^l;
964 select case(idims)
965 {case(^d)
966 ixqmin^d=ixqmin^d+1
967 ixqmax^d=ixqmax^d-1
968 hxpmax^d=ixpmin^d
969 jxpmin^d=ixpmax^d
970 \}
971 end select
972 call gradient(te,ixi^l,ixq^l,idims,gradt(ixi^s,idims))
973 call gradientf(te,x,ixi^l,hxp^l,idims,gradt(ixi^s,idims),nghostcells,.true.)
974 call gradientf(te,x,ixi^l,jxp^l,idims,gradt(ixi^s,idims),nghostcells,.false.)
975 end do
976 bunitvec(ixp^s,:)=block%B0(ixp^s,:,0)
977 lts(ixp^s)=abs(sum(gradt(ixp^s,1:ndim)*bunitvec(ixp^s,1:ndim),dim=ndim+1))/te(ixp^s)
978 if(slab_uniform) then
979 lts(ixp^s)=minval(dxlevel)*lts(ixp^s)
980 else
981 lts(ixp^s)=minval(block%ds(ixp^s,:),dim=ndim+1)*lts(ixp^s)
982 end if
983 lts(ixp^s)=max(one, (exp(lts(ixp^s))/ltrc)**ltrp)
984
985 altr=zero
986 ixp^l=ixo^l^ladd1;
987 do idims=1,ndim
988 hxo^l=ixp^l-kr(idims,^d);
989 jxo^l=ixp^l+kr(idims,^d);
990 altr(ixp^s)=altr(ixp^s)+0.25d0*(lts(hxo^s)+two*lts(ixp^s)+lts(jxo^s))*bunitvec(ixp^s,idims)**2
991 end do
992 block%wextra(ixp^s,tcoff_)=te(ixp^s)*altr(ixp^s)**0.4d0
993 case(3,5)
994 !> do nothing here
995 case default
996 call mpistop("unknown ffhd_trac_type")
997 end select
998 }
999 end subroutine ffhd_get_tcutoff
1000
1001 subroutine ffhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
1003 integer, intent(in) :: ixi^l, ixo^l, idim
1004 double precision, intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
1005 double precision, intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
1006 double precision, intent(in) :: x(ixi^s,1:ndim)
1007 double precision, intent(inout) :: cmax(ixi^s,1:number_species)
1008 double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
1009 double precision, intent(in) :: hspeed(ixi^s,1:number_species)
1010 double precision :: wmean(ixi^s,nw)
1011 double precision, dimension(ixI^S) :: umean, dmean, csoundl, csoundr, tmp1,tmp2,tmp3
1012 integer :: ix^d
1013
1014 select case (boundspeed)
1015 case (1)
1016 ! This implements formula (10.52) from "Riemann Solvers and Numerical
1017 ! Methods for Fluid Dynamics" by Toro.
1018 tmp1(ixo^s)=sqrt(wlp(ixo^s,rho_))
1019 tmp2(ixo^s)=sqrt(wrp(ixo^s,rho_))
1020 tmp3(ixo^s)=1.d0/(tmp1(ixo^s)+tmp2(ixo^s))
1021 umean(ixo^s)=(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)*tmp1(ixo^s)&
1022 +wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)*tmp2(ixo^s))*tmp3(ixo^s)
1023 call ffhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1024 call ffhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1025 dmean(ixo^s)=(tmp1(ixo^s)*csoundl(ixo^s)+tmp2(ixo^s)*csoundr(ixo^s)) * &
1026 tmp3(ixo^s) + 0.5d0*tmp1(ixo^s)*tmp2(ixo^s)*tmp3(ixo^s)**2 * &
1027 (wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)-wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))**2
1028 dmean(ixo^s)=dsqrt(dmean(ixo^s))
1029 if(present(cmin)) then
1030 cmin(ixo^s,1)=umean(ixo^s)-dmean(ixo^s)
1031 cmax(ixo^s,1)=umean(ixo^s)+dmean(ixo^s)
1032 else
1033 cmax(ixo^s,1)=abs(umean(ixo^s))+dmean(ixo^s)
1034 end if
1035 case (2)
1036 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
1037 tmp1(ixo^s)=wmean(ixo^s,mom(1))*block%B0(ixo^s,idim,idim)/wmean(ixo^s,rho_)
1038 call ffhd_get_csound(wmean,x,ixi^l,ixo^l,idim,csoundr)
1039 if(present(cmin)) then
1040 cmax(ixo^s,1)=max(tmp1(ixo^s)+csoundr(ixo^s),zero)
1041 cmin(ixo^s,1)=min(tmp1(ixo^s)-csoundr(ixo^s),zero)
1042 else
1043 cmax(ixo^s,1)=abs(tmp1(ixo^s))+csoundr(ixo^s)
1044 end if
1045 case (3)
1046 ! Miyoshi 2005 JCP 208, 315 equation (67)
1047 call ffhd_get_csound_prim(wlp,x,ixi^l,ixo^l,idim,csoundl)
1048 call ffhd_get_csound_prim(wrp,x,ixi^l,ixo^l,idim,csoundr)
1049 csoundl(ixo^s)=max(csoundl(ixo^s),csoundr(ixo^s))
1050 if(present(cmin)) then
1051 cmin(ixo^s,1)=min(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
1052 wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))-csoundl(ixo^s)
1053 cmax(ixo^s,1)=max(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
1054 wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))+csoundl(ixo^s)
1055 else
1056 cmax(ixo^s,1)=max(wlp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim),&
1057 wrp(ixo^s,mom(1))*block%B0(ixo^s,idim,idim))+csoundl(ixo^s)
1058 end if
1059 end select
1060 end subroutine ffhd_get_cbounds
1061
1062 subroutine ffhd_get_csound(w,x,ixI^L,ixO^L,idim,csound)
1064 integer, intent(in) :: ixi^l, ixo^l, idim
1065 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
1066 double precision, intent(out):: csound(ixi^s)
1067
1068 call ffhd_get_csound2(w,x,ixi^l,ixo^l,csound)
1069 csound(ixo^s) = dsqrt(csound(ixo^s))*abs(block%B0(ixo^s,idim,idim))
1070 end subroutine ffhd_get_csound
1071
1072 !> Calculate fast magnetosonic wave speed
1073 subroutine ffhd_get_csound_prim(w,x,ixI^L,ixO^L,idim,csound)
1075
1076 integer, intent(in) :: ixi^l, ixo^l, idim
1077 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
1078 double precision, intent(out):: csound(ixi^s)
1079
1080 if(ffhd_energy) then
1081 csound(ixo^s)=ffhd_gamma*w(ixo^s,e_)/w(ixo^s,rho_)
1082 else
1083 csound(ixo^s)=ffhd_gamma*ffhd_adiab*w(ixo^s,rho_)**gamma_1
1084 end if
1085 csound(ixo^s) = dsqrt(csound(ixo^s))
1086 end subroutine ffhd_get_csound_prim
1087
1088 subroutine ffhd_get_pthermal_iso(w,x,ixI^L,ixO^L,pth)
1090
1091 integer, intent(in) :: ixi^l, ixo^l
1092 double precision, intent(in) :: w(ixi^s,nw)
1093 double precision, intent(in) :: x(ixi^s,1:ndim)
1094 double precision, intent(out):: pth(ixi^s)
1095
1096 call ffhd_get_rho(w,x,ixi^l,ixo^l,pth)
1097 pth(ixo^s)=ffhd_adiab*pth(ixo^s)**ffhd_gamma
1098 end subroutine ffhd_get_pthermal_iso
1099
1100 subroutine ffhd_get_pthermal_origin(w,x,ixI^L,ixO^L,pth)
1103 integer, intent(in) :: ixi^l, ixo^l
1104 double precision, intent(in) :: w(ixi^s,nw)
1105 double precision, intent(in) :: x(ixi^s,1:ndim)
1106 double precision, intent(out):: pth(ixi^s)
1107 integer :: iw, ix^d
1108
1109 pth(ixo^s)=gamma_1*(w(ixo^s,e_)-ffhd_kin_en(w,ixi^l,ixo^l))
1110 if (fix_small_values) then
1111 {do ix^db= ixo^lim^db\}
1112 if(pth(ix^d)<small_pressure) then
1113 pth(ix^d)=small_pressure
1114 end if
1115 {end do^D&\}
1116 elseif(check_small_values) then
1117 {do ix^db= ixo^lim^db\}
1118 if(pth(ix^d)<small_pressure) then
1119 write(*,*) "Error: small value of gas pressure",pth(ix^d),&
1120 " encountered when call ffhd_get_pthermal"
1121 write(*,*) "Iteration: ", it, " Time: ", global_time
1122 write(*,*) "Location: ", x(ix^d,:)
1123 write(*,*) "Cell number: ", ix^d
1124 do iw=1,nw
1125 write(*,*) trim(cons_wnames(iw)),": ",w(ix^d,iw)
1126 end do
1127 ! use erroneous arithmetic operation to crash the run
1128 if(trace_small_values) write(*,*) sqrt(pth(ix^d)-bigdouble)
1129 write(*,*) "Saving status at the previous time step"
1130 crash=.true.
1131 end if
1132 {end do^D&\}
1133 end if
1134 end subroutine ffhd_get_pthermal_origin
1135
1136 subroutine ffhd_get_temperature_from_te(w, x, ixI^L, ixO^L, res)
1138 integer, intent(in) :: ixi^l, ixo^l
1139 double precision, intent(in) :: w(ixi^s, 1:nw)
1140 double precision, intent(in) :: x(ixi^s, 1:ndim)
1141 double precision, intent(out):: res(ixi^s)
1142
1143 res(ixo^s) = w(ixo^s, te_)
1144 end subroutine ffhd_get_temperature_from_te
1145
1146 subroutine ffhd_get_temperature_from_eint(w, x, ixI^L, ixO^L, res)
1148 integer, intent(in) :: ixi^l, ixo^l
1149 double precision, intent(in) :: w(ixi^s, 1:nw)
1150 double precision, intent(in) :: x(ixi^s, 1:ndim)
1151 double precision, intent(out):: res(ixi^s)
1152 double precision :: r(ixi^s)
1153
1154 call ffhd_get_rfactor(w,x,ixi^l,ixo^l,r)
1155 res(ixo^s) = gamma_1 * w(ixo^s, e_)/(w(ixo^s,rho_)*r(ixo^s))
1156 end subroutine ffhd_get_temperature_from_eint
1157
1158 subroutine ffhd_get_temperature_from_etot(w, x, ixI^L, ixO^L, res)
1160 integer, intent(in) :: ixi^l, ixo^l
1161 double precision, intent(in) :: w(ixi^s, 1:nw)
1162 double precision, intent(in) :: x(ixi^s, 1:ndim)
1163 double precision, intent(out):: res(ixi^s)
1164
1165 double precision :: r(ixi^s)
1166
1167 call ffhd_get_rfactor(w,x,ixi^l,ixo^l,r)
1168 call ffhd_get_pthermal(w,x,ixi^l,ixo^l,res)
1169 res(ixo^s)=res(ixo^s)/(r(ixo^s)*w(ixo^s,rho_))
1170 end subroutine ffhd_get_temperature_from_etot
1171
1172 subroutine ffhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
1174 integer, intent(in) :: ixi^l, ixo^l
1175 double precision, intent(in) :: w(ixi^s,nw)
1176 double precision, intent(in) :: x(ixi^s,1:ndim)
1177 double precision, intent(out) :: csound2(ixi^s)
1178 double precision :: rho(ixi^s)
1179
1180 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
1181 if(ffhd_energy) then
1182 call ffhd_get_pthermal(w,x,ixi^l,ixo^l,csound2)
1183 csound2(ixo^s)=ffhd_gamma*csound2(ixo^s)/rho(ixo^s)
1184 else
1185 csound2(ixo^s)=ffhd_gamma*ffhd_adiab*rho(ixo^s)**gamma_1
1186 end if
1187 end subroutine ffhd_get_csound2
1188
1189 subroutine ffhd_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
1191 use mod_geometry
1192 integer, intent(in) :: ixi^l, ixo^l, idim
1193 ! conservative w
1194 double precision, intent(in) :: wc(ixi^s,nw)
1195 ! primitive w
1196 double precision, intent(in) :: w(ixi^s,nw)
1197 double precision, intent(in) :: x(ixi^s,1:ndim)
1198 double precision,intent(out) :: f(ixi^s,nwflux)
1199 double precision :: ptotal(ixo^s)
1200 double precision :: tmp(ixi^s)
1201 integer :: idirmin, iw, idir, jdir, kdir
1202 double precision, dimension(ixI^S) :: te,tau,sigt
1203
1204 f(ixo^s,rho_)=w(ixo^s,mom(1))*w(ixo^s,rho_)*block%B0(ixo^s,idim,idim)
1205
1206 if(ffhd_energy) then
1207 ptotal(ixo^s)=w(ixo^s,p_)
1208 else
1209 ptotal(ixo^s)=ffhd_adiab*w(ixo^s,rho_)**ffhd_gamma
1210 end if
1211
1212 ! Get flux of momentum
1213 f(ixo^s,mom(1))=(wc(ixo^s,mom(1))*w(ixo^s,mom(1))+ptotal(ixo^s))*block%B0(ixo^s,idim,idim)
1214
1215 ! Get flux of energy
1216 if(ffhd_energy) then
1217 f(ixo^s,e_)=w(ixo^s,mom(1))*(wc(ixo^s,e_)+ptotal(ixo^s))*block%B0(ixo^s,idim,idim)
1219 f(ixo^s,e_)=f(ixo^s,e_)+w(ixo^s,q_)*block%B0(ixo^s,idim,idim)
1220 f(ixo^s,q_)=zero
1221 end if
1222 end if
1223 end subroutine ffhd_get_flux
1224
1225 subroutine ffhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
1230 integer, intent(in) :: ixi^l, ixo^l
1231 double precision, intent(in) :: qdt,dtfactor
1232 double precision, intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:ndim)
1233 double precision, intent(inout) :: w(ixi^s,1:nw)
1234 logical, intent(in) :: qsourcesplit
1235 logical, intent(inout) :: active
1236
1237 if (.not. qsourcesplit) then
1238 active = .true.
1239 call add_punitb(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
1241 call add_hypertc_source(qdt,ixi^l,ixo^l,wct,w,x,wctprim)
1242 end if
1243 end if
1244
1245 if(ffhd_radiative_cooling) then
1246 call radiative_cooling_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
1247 w,x,qsourcesplit,active, rc_fl)
1248 end if
1249
1250 if(ffhd_viscosity) then
1251 call viscosity_add_source(qdt,ixi^l,ixo^l,wct,&
1252 w,x,ffhd_energy,qsourcesplit,active)
1253 end if
1254
1255 if(ffhd_gravity) then
1256 call gravity_add_source(qdt,ixi^l,ixo^l,wct,wctprim,&
1257 w,x,ffhd_energy,.false.,qsourcesplit,active)
1258 end if
1259
1260 ! update temperature from new pressure, density, and old ionization degree
1262 if(.not.qsourcesplit) then
1263 active = .true.
1264 call ffhd_update_temperature(ixi^l,ixo^l,wct,w,x)
1265 end if
1266 end if
1267 end subroutine ffhd_add_source
1268
1269 subroutine add_punitb(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
1271 use mod_geometry
1272 integer, intent(in) :: ixi^l,ixo^l
1273 double precision, intent(in) :: qdt
1274 double precision, intent(in) :: wct(ixi^s,1:nw),x(ixi^s,1:ndim)
1275 double precision, intent(in) :: wctprim(ixi^s,1:nw)
1276 double precision, intent(inout) :: w(ixi^s,1:nw)
1277
1278 integer :: idims,hxo^l
1279 double precision :: divb(ixi^s)
1280
1281 divb=zero
1282 if(slab_uniform) then
1283 do idims=1,ndim
1284 hxo^l=ixo^l-kr(idims,^d);
1285 divb(ixo^s)=divb(ixo^s)+(block%B0(ixo^s,idims,idims)-block%B0(hxo^s,idims,idims))/dxlevel(idims)
1286 end do
1287 else
1288 call divvector(block%B0(ixi^s,1:ndir,0),ixi^l,ixo^l,divb)
1289 end if
1290 w(ixo^s,mom(1))=w(ixo^s,mom(1))+qdt*wctprim(ixo^s,p_)*divb(ixo^s)
1291 end subroutine add_punitb
1292
1293 subroutine ffhd_get_rho(w,x,ixI^L,ixO^L,rho)
1295 integer, intent(in) :: ixi^l, ixo^l
1296 double precision, intent(in) :: w(ixi^s,1:nw),x(ixi^s,1:ndim)
1297 double precision, intent(out) :: rho(ixi^s)
1298
1299 rho(ixo^s) = w(ixo^s,rho_)
1300 end subroutine ffhd_get_rho
1301
1302 subroutine ffhd_handle_small_ei(w, x, ixI^L, ixO^L, ie, subname)
1305 integer, intent(in) :: ixi^l,ixo^l, ie
1306 double precision, intent(inout) :: w(ixi^s,1:nw)
1307 double precision, intent(in) :: x(ixi^s,1:ndim)
1308 character(len=*), intent(in) :: subname
1309 integer :: idir
1310 logical :: flag(ixi^s,1:nw)
1311 double precision :: rho(ixi^s)
1312
1313 flag=.false.
1314 where(w(ixo^s,ie)<small_e) flag(ixo^s,ie)=.true.
1315 if(any(flag(ixo^s,ie))) then
1316 select case (small_values_method)
1317 case ("replace")
1318 where(flag(ixo^s,ie)) w(ixo^s,ie)=small_e
1319 case ("average")
1320 call small_values_average(ixi^l, ixo^l, w, x, flag, ie)
1321 case default
1322 w(ixo^s,e_)=w(ixo^s,e_)*gamma_1
1323 call ffhd_get_rho(w,x,ixi^l,ixo^l,rho)
1324 w(ixo^s,mom(1)) = w(ixo^s,mom(1))/rho(ixo^s)
1325 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
1326 end select
1327 end if
1328 end subroutine ffhd_handle_small_ei
1329
1330 subroutine ffhd_update_temperature(ixI^L,ixO^L,wCT,w,x)
1333 integer, intent(in) :: ixi^l, ixo^l
1334 double precision, intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
1335 double precision, intent(inout) :: w(ixi^s,1:nw)
1336 double precision :: iz_h(ixo^s),iz_he(ixo^s), pth(ixi^s)
1337
1338 call ionization_degree_from_temperature(ixi^l,ixo^l,wct(ixi^s,te_),iz_h,iz_he)
1339 call ffhd_get_pthermal(w,x,ixi^l,ixo^l,pth)
1340 w(ixo^s,te_)=(2.d0+3.d0*he_abundance)*pth(ixo^s)/(w(ixo^s,rho_)*(1.d0+iz_h(ixo^s)+&
1341 he_abundance*(iz_he(ixo^s)*(iz_he(ixo^s)+1.d0)+1.d0)))
1342 end subroutine ffhd_update_temperature
1343
1344 subroutine ffhd_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1346 use mod_usr_methods
1349 use mod_gravity, only: gravity_get_dt
1350 integer, intent(in) :: ixi^l, ixo^l
1351 double precision, intent(inout) :: dtnew
1352 double precision, intent(in) :: dx^d
1353 double precision, intent(in) :: w(ixi^s,1:nw)
1354 double precision, intent(in) :: x(ixi^s,1:ndim)
1355 integer :: idirmin,idim
1356 double precision :: dxarr(ndim)
1357 double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1358
1359 dtnew = bigdouble
1360
1361 if(ffhd_radiative_cooling) then
1362 call cooling_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x,rc_fl)
1363 end if
1364
1365 if(ffhd_viscosity) then
1366 call viscosity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1367 end if
1368
1369 if(ffhd_gravity) then
1370 call gravity_get_dt(w,ixi^l,ixo^l,dtnew,dx^d,x)
1371 end if
1372 end subroutine ffhd_get_dt
1373
1374! subroutine ffhd_add_source_geom(qdt,dtfactor,ixI^L,ixO^L,wCT,w,x)
1375! use mod_global_parameters
1376! use mod_geometry
1377!
1378! integer, intent(in) :: ixI^L, ixO^L
1379! double precision, intent(in) :: qdt, dtfactor,x(ixI^S,1:ndim)
1380! double precision, intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1381!
1382! integer :: iw,idir, h1x^L{^NOONED, h2x^L}
1383! double precision :: tmp(ixI^S),tmp1(ixI^S),tmp2(ixI^S),invrho(ixO^S),invr(ixO^S)
1384!
1385! integer :: mr_,mphi_ ! Polar var. names
1386! integer :: br_,bphi_
1387!
1388! mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1389! br_=mag(1); bphi_=mag(1)-1+phi_
1390!
1391! ! 1/rho
1392! invrho(ixO^S)=1.d0/wCT(ixO^S,rho_)
1393! ! include dt in invr, invr is always used with qdt
1394! if(local_timestep) then
1395! invr(ixO^S) = block%dt(ixO^S) * dtfactor/x(ixO^S,1)
1396! else
1397! invr(ixO^S) = qdt/x(ixO^S,1)
1398! end if
1399!
1400!
1401! select case (coordinate)
1402! case (cylindrical)
1403! call ffhd_get_p_total(wCT,x,ixI^L,ixO^L,tmp)
1404! if(phi_>0) then
1405! w(ixO^S,mr_)=w(ixO^S,mr_)+invr(ixO^S)*(tmp(ixO^S)-&
1406! wCT(ixO^S,bphi_)**2+wCT(ixO^S,mphi_)**2*invrho(ixO^S))
1407! w(ixO^S,mphi_)=w(ixO^S,mphi_)+invr(ixO^S)*(&
1408! -wCT(ixO^S,mphi_)*wCT(ixO^S,mr_)*invrho(ixO^S) &
1409! +wCT(ixO^S,bphi_)*wCT(ixO^S,br_))
1410! if(.not.stagger_grid) then
1411! w(ixO^S,bphi_)=w(ixO^S,bphi_)+invr(ixO^S)*&
1412! (wCT(ixO^S,bphi_)*wCT(ixO^S,mr_) &
1413! -wCT(ixO^S,br_)*wCT(ixO^S,mphi_)) &
1414! *invrho(ixO^S)
1415! end if
1416! else
1417! w(ixO^S,mr_)=w(ixO^S,mr_)+invr(ixO^S)*tmp(ixO^S)
1418! end if
1419! if(ffhd_glm) w(ixO^S,br_)=w(ixO^S,br_)+wCT(ixO^S,psi_)*invr(ixO^S)
1420! case (spherical)
1421! h1x^L=ixO^L-kr(1,^D); {^NOONED h2x^L=ixO^L-kr(2,^D);}
1422! call ffhd_get_p_total(wCT,x,ixI^L,ixO^L,tmp1)
1423! ! m1
1424! tmp(ixO^S)=tmp1(ixO^S)*x(ixO^S,1) &
1425! *(block%surfaceC(ixO^S,1)-block%surfaceC(h1x^S,1))/block%dvolume(ixO^S)
1426! do idir=2,ndir
1427! tmp(ixO^S)=tmp(ixO^S)+wCT(ixO^S,mom(idir))**2*invrho(ixO^S)-wCT(ixO^S,mag(idir))**2
1428! end do
1429! w(ixO^S,mom(1))=w(ixO^S,mom(1))+tmp(ixO^S)*invr(ixO^S)
1430! ! b1
1431! if(ffhd_glm) then
1432! w(ixO^S,mag(1))=w(ixO^S,mag(1))+invr(ixO^S)*2.0d0*wCT(ixO^S,psi_)
1433! end if
1434!
1435! {^NOONED
1436! ! m2
1437! ! This will make hydrostatic p=const an exact solution
1438! if(local_timestep) then
1439! tmp(ixO^S) = block%dt(ixO^S) * tmp1(ixO^S)
1440! else
1441! tmp(ixO^S) = qdt * tmp1(ixO^S)
1442! end if
1443! w(ixO^S,mom(2))=w(ixO^S,mom(2))+tmp(ixO^S) &
1444! *(block%surfaceC(ixO^S,2)-block%surfaceC(h2x^S,2)) &
1445! /block%dvolume(ixO^S)
1446! tmp(ixO^S)=-(wCT(ixO^S,mom(1))*wCT(ixO^S,mom(2))*invrho(ixO^S) &
1447! -wCT(ixO^S,mag(1))*wCT(ixO^S,mag(2)))
1448! if(ndir==3) then
1449! tmp(ixO^S)=tmp(ixO^S)+(wCT(ixO^S,mom(3))**2*invrho(ixO^S) &
1450! -wCT(ixO^S,mag(3))**2)*dcos(x(ixO^S,2))/dsin(x(ixO^S,2))
1451! end if
1452! w(ixO^S,mom(2))=w(ixO^S,mom(2))+tmp(ixO^S)*invr(ixO^S)
1453! ! b2
1454! if(.not.stagger_grid) then
1455! tmp(ixO^S)=(wCT(ixO^S,mom(1))*wCT(ixO^S,mag(2)) &
1456! -wCT(ixO^S,mom(2))*wCT(ixO^S,mag(1)))*invrho(ixO^S)
1457! if(ffhd_glm) then
1458! tmp(ixO^S)=tmp(ixO^S) &
1459! + dcos(x(ixO^S,2))/dsin(x(ixO^S,2))*wCT(ixO^S,psi_)
1460! end if
1461! w(ixO^S,mag(2))=w(ixO^S,mag(2))+tmp(ixO^S)*invr(ixO^S)
1462! end if
1463! }
1464!
1465! if(ndir==3) then
1466! ! m3
1467! tmp(ixO^S)=-(wCT(ixO^S,mom(3))*wCT(ixO^S,mom(1))*invrho(ixO^S) &
1468! -wCT(ixO^S,mag(3))*wCT(ixO^S,mag(1))) {^NOONED &
1469! -(wCT(ixO^S,mom(2))*wCT(ixO^S,mom(3))*invrho(ixO^S) &
1470! -wCT(ixO^S,mag(2))*wCT(ixO^S,mag(3))) &
1471! *dcos(x(ixO^S,2))/dsin(x(ixO^S,2)) }
1472! w(ixO^S,mom(3))=w(ixO^S,mom(3))+tmp(ixO^S)*invr(ixO^S)
1473! ! b3
1474! if(.not.stagger_grid) then
1475! tmp(ixO^S)=(wCT(ixO^S,mom(1))*wCT(ixO^S,mag(3)) &
1476! -wCT(ixO^S,mom(3))*wCT(ixO^S,mag(1)))*invrho(ixO^S) {^NOONED &
1477! -(wCT(ixO^S,mom(3))*wCT(ixO^S,mag(2)) &
1478! -wCT(ixO^S,mom(2))*wCT(ixO^S,mag(3)))*dcos(x(ixO^S,2)) &
1479! /(wCT(ixO^S,rho_)*dsin(x(ixO^S,2))) }
1480! w(ixO^S,mag(3))=w(ixO^S,mag(3))+tmp(ixO^S)*invr(ixO^S)
1481! end if
1482! end if
1483! end select
1484! end subroutine ffhd_add_source_geom
1485
1486 function ffhd_kin_en_origin(w, ixI^L, ixO^L, inv_rho) result(ke)
1487 use mod_global_parameters, only: nw, ndim,block
1488 integer, intent(in) :: ixi^l, ixo^l
1489 double precision, intent(in) :: w(ixi^s, nw)
1490 double precision :: ke(ixo^s)
1491 double precision, intent(in), optional :: inv_rho(ixo^s)
1492
1493 if(present(inv_rho)) then
1494 ke(ixo^s)=0.5d0*w(ixo^s,mom(1))**2*inv_rho(ixo^s)
1495 else
1496 ke(ixo^s)=0.5d0*w(ixo^s,mom(1))**2/w(ixo^s,rho_)
1497 end if
1498 end function ffhd_kin_en_origin
1499
1500 subroutine rfactor_from_temperature_ionization(w,x,ixI^L,ixO^L,Rfactor)
1503 integer, intent(in) :: ixi^l, ixo^l
1504 double precision, intent(in) :: w(ixi^s,1:nw)
1505 double precision, intent(in) :: x(ixi^s,1:ndim)
1506 double precision, intent(out):: rfactor(ixi^s)
1507 double precision :: iz_h(ixo^s),iz_he(ixo^s)
1508
1509 call ionization_degree_from_temperature(ixi^l,ixo^l,w(ixi^s,te_),iz_h,iz_he)
1510 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)
1511 end subroutine rfactor_from_temperature_ionization
1512
1513 subroutine rfactor_from_constant_ionization(w,x,ixI^L,ixO^L,Rfactor)
1515 integer, intent(in) :: ixi^l, ixo^l
1516 double precision, intent(in) :: w(ixi^s,1:nw)
1517 double precision, intent(in) :: x(ixi^s,1:ndim)
1518 double precision, intent(out):: rfactor(ixi^s)
1519
1520 rfactor(ixo^s)=rr
1521 end subroutine rfactor_from_constant_ionization
1522
1523 subroutine get_tau(ixI^L,ixO^L,w,Te,tau,sigT5)
1525 integer, intent(in) :: ixi^l, ixo^l
1526 double precision, dimension(ixI^S,1:nw), intent(in) :: w
1527 double precision, dimension(ixI^S), intent(in) :: te
1528 double precision, dimension(ixI^S), intent(out) :: tau,sigt5
1529 integer :: ix^d
1530 double precision :: dxmin,taumin
1531 double precision, dimension(ixI^S) :: sigt7,eint
1532
1533 taumin=4.d0
1534 !> w supposed to be wCTprim here
1535 if(ffhd_trac) then
1536 where(te(ixo^s) .lt. block%wextra(ixo^s,tcoff_))
1537 sigt5(ixo^s)=hypertc_kappa*sqrt(block%wextra(ixo^s,tcoff_)**5)
1538 sigt7(ixo^s)=sigt5(ixo^s)*block%wextra(ixo^s,tcoff_)
1539 else where
1540 sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
1541 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
1542 end where
1543 else
1544 sigt5(ixo^s)=hypertc_kappa*sqrt(te(ixo^s)**5)
1545 sigt7(ixo^s)=sigt5(ixo^s)*te(ixo^s)
1546 end if
1547 eint(ixo^s)=w(ixo^s,p_)/(ffhd_gamma-one)
1548 tau(ixo^s)=max(taumin*dt,sigt7(ixo^s)/eint(ixo^s)/cs2max_global)
1549 end subroutine get_tau
1550
1551 subroutine add_hypertc_source(qdt,ixI^L,ixO^L,wCT,w,x,wCTprim)
1553 integer, intent(in) :: ixi^l,ixo^l
1554 double precision, intent(in) :: qdt
1555 double precision, dimension(ixI^S,1:ndim), intent(in) :: x
1556 double precision, dimension(ixI^S,1:nw), intent(in) :: wct,wctprim
1557 double precision, dimension(ixI^S,1:nw), intent(inout) :: w
1558 integer :: idims
1559 integer :: hxc^l,hxo^l,ixc^l,jxc^l,jxo^l,kxc^l
1560 double precision :: invdx
1561 double precision, dimension(ixI^S) :: te,tau,sigt,htc_qsrc,tface
1562 double precision, dimension(ixI^S) :: htc_esrc
1563
1564 te(ixi^s)=wctprim(ixi^s,p_)/wct(ixi^s,rho_)
1565 call get_tau(ixi^l,ixo^l,wctprim,te,tau,sigt)
1566 htc_qsrc=zero
1567 do idims=1,ndim
1568 invdx=1.d0/dxlevel(idims)
1569 ixc^l=ixo^l;
1570 ixcmin^d=ixomin^d-kr(idims,^d);ixcmax^d=ixomax^d;
1571 jxc^l=ixc^l+kr(idims,^d);
1572 kxc^l=jxc^l+kr(idims,^d);
1573 hxc^l=ixc^l-kr(idims,^d);
1574 hxo^l=ixo^l-kr(idims,^d);
1575 tface(ixc^s)=(7.d0*(te(ixc^s)+te(jxc^s))-(te(hxc^s)+te(kxc^s)))/12.d0
1576 htc_qsrc(ixo^s)=htc_qsrc(ixo^s)+sigt(ixo^s)*block%B0(ixo^s,idims,0)*(tface(ixo^s)-tface(hxo^s))*invdx
1577 end do
1578 htc_qsrc(ixo^s)=(htc_qsrc(ixo^s)+wct(ixo^s,q_))/tau(ixo^s)
1579 w(ixo^s,q_)=w(ixo^s,q_)-qdt*htc_qsrc(ixo^s)
1580 end subroutine add_hypertc_source
1581end 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)