16 integer,
public,
protected ::
rho_
17 integer,
public,
protected ::
d_
20 integer,
allocatable,
public,
protected ::
mom(:)
23 integer,
allocatable,
public,
protected ::
tracer(:)
26 integer,
public,
protected ::
e_
28 integer,
public,
protected ::
p_
31 integer,
public,
protected ::
lfac_
34 integer,
public,
protected ::
xi_
54 double precision,
public ::
tolernr = 1.0d-9
55 double precision,
public ::
dmaxvel = 1.0d-7
74 subroutine srhd_read_params(files)
76 character(len=*),
intent(in) :: files(:)
84 open(
unitpar, file=trim(files(n)), status=
"old")
85 read(
unitpar, srhd_list,
end=111)
89 end subroutine srhd_read_params
92 subroutine srhd_write_info(fh)
94 integer,
intent(in) :: fh
95 integer,
parameter :: n_par = 1
96 double precision :: values(n_par)
97 character(len=name_len) :: names(n_par)
98 integer,
dimension(MPI_STATUS_SIZE) :: st
101 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
105 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
106 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
108 end subroutine srhd_write_info
118 physics_type =
"srhd"
120 phys_total_energy = .true.
124 phys_internal_e = .false.
125 phys_partial_ionization=.false.
131 allocate(start_indices(number_species),stop_indices(number_species))
140 mom(:) = var_set_momentum(
ndir)
143 e_ = var_set_energy()
147 call srhd_physical_units()
153 tracer(itr) = var_set_fluxvar(
"trc",
"trp", itr, need_bc=.false.)
158 xi_ = var_set_auxvar(
'xi',
'xi')
159 lfac_= var_set_auxvar(
'lfac',
'lfac')
165 stop_indices(1)=nwflux
171 if (.not.
allocated(flux_type))
then
172 allocate(flux_type(
ndir, nw))
173 flux_type = flux_default
174 else if (any(shape(flux_type) /= [
ndir, nw]))
then
175 call mpistop(
"phys_check error: flux_type has wrong shape")
179 allocate(iw_vector(nvector))
180 iw_vector(1) =
mom(1) - 1
183 phys_add_source => srhd_add_source
184 phys_get_dt => srhd_get_dt
186 phys_get_a2max => srhd_get_a2max
191 phys_get_cmax => srhd_get_cmax
192 phys_get_cbounds => srhd_get_cbounds
193 phys_get_flux => srhd_get_flux
194 phys_add_source_geom => srhd_add_source_geom
201 phys_get_v => srhd_get_v
202 phys_write_info => srhd_write_info
203 phys_handle_small_values => srhd_handle_small_values
216 call mpistop (
"Error: srhd_gamma <= 0 or srhd_gamma == 1")
224 call srhd_get_smallvalues_eos
227 write(*,*)
'------------------------------------------------------------'
228 write(*,*)
'Using EOS set via srhd_eos=',
srhd_eos
229 write(*,*)
'Maximal lorentz factor (via dmaxvel) is=',
lfacmax
233 write(*,*)
'------------------------------------------------------------'
238 subroutine srhd_physical_units
240 double precision :: mp,kb
253 call mpistop(
"Abort: must set positive values for unit length and numberdensity")
263 end subroutine srhd_physical_units
268 logical,
intent(in) :: primitive
269 integer,
intent(in) :: ixi^
l, ixo^
l
270 double precision,
intent(in) :: w(ixi^s, nw)
271 logical,
intent(inout) :: flag(ixi^s,1:nw)
281 where(w(ixo^s,
e_) <
small_e) flag(ixo^s,
e_) = .true.
287 subroutine srhd_check_w_aux(ixI^L, ixO^L, w, flag)
289 integer,
intent(in) :: ixi^
l, ixo^
l
290 double precision,
intent(in) :: w(ixi^s, nw)
291 logical,
intent(inout) :: flag(ixi^s,1:nw)
296 where(w(ixo^s,
lfac_) < one) flag(ixo^s,
lfac_) = .true.
298 if(any(flag(ixo^s,
xi_)))
then
299 write(*,*)
'auxiliary xi too low: abort'
300 call mpistop(
'auxiliary check failed')
302 if(any(flag(ixo^s,
lfac_)))
then
303 write(*,*)
'auxiliary lfac too low: abort'
304 call mpistop(
'auxiliary check failed')
307 end subroutine srhd_check_w_aux
313 integer,
intent(in) :: ixi^
l, ixo^
l
314 double precision,
intent(inout) :: w(ixi^s, nw)
315 double precision,
dimension(ixO^S) :: rho,rhoh,pth
318 rho(ixo^s) = sum(w(ixo^s,
mom(:))**2, dim=
ndim+1)
319 w(ixo^s,
lfac_) = dsqrt(1.0d0+rho(ixo^s))
321 rho(ixo^s)=w(ixo^s,
rho_)
322 pth(ixo^s)=w(ixo^s,
p_)
327 w(ixo^s,
xi_) = w(ixo^s,
lfac_)**2.0d0*rhoh(ixo^s)
337 integer,
intent(in) :: ixi^
l,ixo^
l
338 double precision,
intent(inout) :: w(ixi^s, nw)
339 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
341 integer :: ix^
d,ierror,idir
342 integer :: flag_error(ixo^s)
343 double precision :: ssqr
346 {
do ix^db=ixomin^db,ixomax^db\}
350 ssqr= ssqr+w(ix^
d,
mom(idir))**2
353 print *,
'entering con2prim with', w(ix^
d,
lfac_),w(ix^
d,
xi_), &
354 w(ix^
d,
d_),ssqr,w(ix^
d,
e_)
355 print *,
'in position',ix^
d,x(ix^
d,1:
ndim)
361 print *,
'entering con2prim with', w(ix^
d,
lfac_),w(ix^
d,
xi_), &
362 w(ix^
d,
d_),ssqr,w(ix^
d,
e_)
363 print *,
'in position',ix^
d,x(ix^
d,1:
ndim)
368 call con2prim_eos(w(ix^
d,
lfac_),w(ix^
d,
xi_), &
369 w(ix^
d,
d_),ssqr,w(ix^
d,
e_),ierror)
370 flag_error(ix^
d) = ierror
373 {
do ix^db=ixomin^db,ixomax^db\}
377 ssqr= ssqr+w(ix^
d,
mom(idir))**2
380 print *,
'entering con2prim with', w(ix^
d,
lfac_),w(ix^
d,
xi_), &
381 w(ix^
d,
d_),ssqr,w(ix^
d,
e_)
382 print *,
'in position',ix^
d,x(ix^
d,1:
ndim)
388 print *,
'entering con2prim with', w(ix^
d,
lfac_),w(ix^
d,
xi_), &
389 w(ix^
d,
d_),ssqr,w(ix^
d,
e_)
390 print *,
'in position',ix^
d,x(ix^
d,1:
ndim)
396 w(ix^
d,
d_),ssqr,w(ix^
d,
e_),ierror)
397 flag_error(ix^
d) = ierror
402 if(any(flag_error(ixo^s)/=0))
then
403 print *,flag_error(ixo^s)
404 call mpistop(
'Problem when getting auxiliaries')
414 integer,
intent(in) :: ixi^
l, ixo^
l
415 double precision,
intent(inout) :: w(ixi^s, nw)
416 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
418 double precision,
dimension(ixO^S) :: rhoh,rho,pth
422 rhoh(ixo^s) = sum(w(ixo^s,
mom(:))**2, dim=
ndim+1)
423 w(ixo^s,
lfac_) = dsqrt(1.0d0+rhoh(ixo^s))
425 rho(ixo^s)=w(ixo^s,
rho_)
426 pth(ixo^s)=w(ixo^s,
p_)
431 rhoh(ixo^s)= rhoh(ixo^s)*w(ixo^s,
lfac_)
433 w(ixo^s,
xi_) = w(ixo^s,
lfac_)*rhoh(ixo^s)
436 w(ixo^s,
d_)=w(ixo^s,
lfac_)*rho(ixo^s)
440 w(ixo^s,
mom(idir)) = rhoh(ixo^s)*w(ixo^s,
mom(idir))
444 w(ixo^s,
e_) = w(ixo^s,
xi_)-pth(ixo^s)-w(ixo^s,
d_)
455 integer,
intent(in) :: ixi^
l, ixo^
l
456 double precision,
intent(inout) :: w(ixi^s, nw)
457 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
459 double precision,
dimension(ixO^S) :: rho,rhoh,e
460 double precision,
dimension(ixI^S) :: pth
461 character(len=30) :: subname_loc
467 rho(ixo^s) = w(ixo^s,
d_)/w(ixo^s,
lfac_)
471 rhoh(ixo^s) = w(ixo^s,
xi_)/w(ixo^s,
lfac_)**2.0d0
472 call srhd_get_pressure_eos(ixi^
l,ixo^
l,rho,rhoh,pth,e)
474 w(ixo^s,
rho_)=rho(ixo^s)
477 w(ixo^s,
mom(idir)) = w(ixo^s,
lfac_)*w(ixo^s,
mom(idir))&
480 w(ixo^s,
p_)=pth(ixo^s)
484 /(rho(ixo^s)*w(ixo^s,
lfac_))
490 subroutine srhd_get_v(w,x,ixI^L,ixO^L,v)
492 integer,
intent(in) :: ixi^
l, ixo^
l
493 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
494 double precision,
intent(out) :: v(ixi^s,1:
ndir)
499 v(ixo^s,idir) = w(ixo^s,
mom(idir))/w(ixo^s,
xi_)
502 end subroutine srhd_get_v
507 subroutine srhd_get_csound2_rhoh(w,x,ixI^L,ixO^L,rhoh,csound2)
509 integer,
intent(in) :: ixi^
l, ixo^
l
510 double precision,
intent(in) :: w(ixi^s,nw),rhoh(ixo^s)
511 double precision,
intent(in) :: x(ixi^s,1:
ndim)
512 double precision,
intent(out) :: csound2(ixo^s)
514 double precision :: rho(ixo^s)
517 call srhd_get_csound2_eos(ixi^
l,ixo^
l,rho,rhoh,csound2)
519 end subroutine srhd_get_csound2_rhoh
526 integer,
intent(in) :: ixi^
l, ixo^
l
527 double precision,
intent(inout) :: w(ixi^s,nw)
528 double precision,
intent(in) :: x(ixi^s,1:
ndim)
529 double precision,
intent(out) :: csound2(ixo^s)
531 double precision :: rho(ixo^s),rhoh(ixo^s)
536 rho = w(ixo^s,
d_)/w(ixo^s,
lfac_)
537 rhoh = w(ixo^s,
xi_)/w(ixo^s,
lfac_)**2.0d0
538 call srhd_get_csound2_eos(ixi^
l,ixo^
l,rho,rhoh,csound2)
548 integer,
intent(in) :: ixi^
l, ixo^
l
549 double precision,
intent(in) :: w(ixi^s, nw)
550 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
551 double precision,
intent(out) :: pth(ixi^s)
554 double precision :: rho(ixo^s),rhoh(ixo^s),e(ixo^s)
557 rho = w(ixo^s,
d_)/w(ixo^s,
lfac_)
558 rhoh = w(ixo^s,
xi_)/w(ixo^s,
lfac_)**2.0d0
559 call srhd_get_pressure_eos(ixi^
l,ixo^
l,rho,rhoh,pth,e)
565 subroutine srhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
568 integer,
intent(in) :: ixi^
l, ixo^
l
569 double precision,
intent(in) :: qdt,dtfactor
570 double precision,
intent(in) :: wct(ixi^s, 1:nw),wctprim(ixi^s,1:nw), x(ixi^s, 1:
ndim)
571 double precision,
intent(inout) :: w(ixi^s, 1:nw)
572 logical,
intent(in) :: qsourcesplit
573 logical,
intent(inout) :: active
575 end subroutine srhd_add_source
578 subroutine srhd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
581 integer,
intent(in) :: ixi^
l, ixo^
l
582 double precision,
intent(in) ::
dx^
d, x(ixi^s, 1:^nd)
583 double precision,
intent(in) :: w(ixi^s, 1:nw)
584 double precision,
intent(inout) :: dtnew
588 end subroutine srhd_get_dt
592 subroutine srhd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
595 integer,
intent(in) :: ixi^
l, ixo^
l, idim
596 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s, 1:
ndim)
597 double precision,
intent(inout) :: cmax(ixi^s)
599 double precision :: wc(ixi^s,nw)
600 double precision,
dimension(ixO^S) :: csound2,tmp1,tmp2,v2
601 double precision,
dimension(ixI^S) :: vidim, cmin
603 logical :: flag(ixi^s,1:nw)
611 tmp1(ixo^s)=wc(ixo^s,
xi_)/wc(ixo^s,
lfac_)**2.0d0
612 v2(ixo^s)=1.0d0-1.0d0/wc(ixo^s,
lfac_)**2
613 call srhd_get_csound2_rhoh(wc,x,ixi^
l,ixo^
l,tmp1,csound2)
614 vidim(ixo^s) = wc(ixo^s,
mom(idim))/wc(ixo^s,
xi_)
615 tmp2(ixo^s)=vidim(ixo^s)**2.0d0
616 tmp1(ixo^s)=1.0d0-v2(ixo^s)*csound2(ixo^s) &
617 -tmp2(ixo^s)*(1.0d0-csound2(ixo^s))
618 tmp2(ixo^s)=dsqrt(csound2(ixo^s)*(one-v2(ixo^s))*tmp1(ixo^s))
619 tmp1(ixo^s)=vidim(ixo^s)*(one-csound2(ixo^s))
620 cmax(ixo^s)=(tmp1(ixo^s)+tmp2(ixo^s))/(one-v2(ixo^s)*csound2(ixo^s))
621 cmin(ixo^s)=(tmp1(ixo^s)-tmp2(ixo^s))/(one-v2(ixo^s)*csound2(ixo^s))
623 cmin(ixo^s) = max(cmin(ixo^s), - 1.0d0)
624 cmin(ixo^s) = min(cmin(ixo^s), 1.0d0)
625 cmax(ixo^s) = max(cmax(ixo^s), - 1.0d0)
626 cmax(ixo^s) = min(cmax(ixo^s), 1.0d0)
628 cmax(ixo^s) = max(dabs(cmax(ixo^s)),dabs(cmin(ixo^s)))
630 end subroutine srhd_get_cmax
632 subroutine srhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
635 integer,
intent(in) :: ixi^
l, ixo^
l
636 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
637 double precision,
intent(inout) :: a2max(
ndim)
638 double precision :: a2(ixi^s,
ndim,nw)
639 integer :: gxo^
l,hxo^
l,jxo^
l,kxo^
l,i
644 hxo^
l=ixo^
l-
kr(i,^
d);
645 gxo^
l=hxo^
l-
kr(i,^
d);
646 jxo^
l=ixo^
l+
kr(i,^
d);
647 kxo^
l=jxo^
l+
kr(i,^
d);
648 a2(ixo^s,i,1:nwflux)=dabs(-w(kxo^s,1:nwflux)+16.d0*w(jxo^s,1:nwflux)&
649 -30.d0*w(ixo^s,1:nwflux)+16.d0*w(hxo^s,1:nwflux)-w(gxo^s,1:nwflux))
650 a2max(i)=maxval(a2(ixo^s,i,1:nwflux))/12.d0/
dxlevel(i)**2
653 end subroutine srhd_get_a2max
656 subroutine srhd_get_cmax_loc(ixI^L,ixO^L,vidim,csound2,v2,cmax,cmin)
658 integer,
intent(in) :: ixi^
l, ixo^
l
659 double precision,
intent(in) :: vidim(ixi^s)
660 double precision,
intent(in),
dimension(ixO^S) :: csound2
661 double precision,
intent(in) :: v2(ixi^s)
662 double precision,
intent(out) :: cmax(ixi^s)
663 double precision,
intent(out) :: cmin(ixi^s)
665 double precision,
dimension(ixI^S):: tmp1,tmp2
667 tmp2(ixo^s)=vidim(ixo^s)**2.0d0
668 tmp1(ixo^s)=1.0d0-v2(ixo^s)*csound2(ixo^s) &
669 -tmp2(ixo^s)*(1.0d0-csound2(ixo^s))
670 tmp2(ixo^s)=dsqrt(csound2(ixo^s)*(one-v2(ixo^s))*tmp1(ixo^s))
671 tmp1(ixo^s)=vidim(ixo^s)*(one-csound2(ixo^s))
672 cmax(ixo^s)=(tmp1(ixo^s)+tmp2(ixo^s))/(one-v2(ixo^s)*csound2(ixo^s))
674 cmax(ixo^s) = max(cmax(ixo^s), - 1.0d0)
675 cmax(ixo^s) = min(cmax(ixo^s), 1.0d0)
676 cmin(ixo^s)=(tmp1(ixo^s)-tmp2(ixo^s))/(one-v2(ixo^s)*csound2(ixo^s))
678 cmin(ixo^s) = max(cmin(ixo^s), - 1.0d0)
679 cmin(ixo^s) = min(cmin(ixo^s), 1.0d0)
681 end subroutine srhd_get_cmax_loc
685 subroutine srhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
689 integer,
intent(in) :: ixi^
l, ixo^
l, idim
691 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
693 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
694 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
696 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
699 double precision :: wmean(ixi^s,
nw)
700 double precision,
dimension(ixO^S) :: csound2,tmp1,tmp2,tmp3
701 double precision,
dimension(ixI^S) :: vidim,cmaxl,cmaxr,cminl,cminr,v2
703 logical :: flag(ixi^s,1:
nw)
710 tmp2=wlp(ixo^s,
xi_)/wlp(ixo^s,
lfac_)**2.0d0
712 call srhd_get_csound2_prim_eos(ixo^
l,tmp1,tmp2,tmp3,csound2)
713 vidim(ixo^s) = wlp(ixo^s,
mom(idim))/wlp(ixo^s,
lfac_)
714 v2(ixo^s) = 1.0d0-1.0d0/wlp(ixo^s,
lfac_)**2
715 call srhd_get_cmax_loc(ixi^
l,ixo^
l,vidim,csound2,v2,cmaxl,cminl)
720 tmp2=wrp(ixo^s,
xi_)/wrp(ixo^s,
lfac_)**2.0d0
722 call srhd_get_csound2_prim_eos(ixo^
l,tmp1,tmp2,tmp3,csound2)
723 vidim(ixo^s) = wrp(ixo^s,
mom(idim))/wrp(ixo^s,
lfac_)
724 v2(ixo^s) = 1.0d0-1.0d0/wrp(ixo^s,
lfac_)**2
725 call srhd_get_cmax_loc(ixi^
l,ixo^
l,vidim,csound2,v2,cmaxr,cminr)
727 if(
present(cmin))
then
729 cmax(ixo^s,1)=max(cmaxl(ixo^s),cmaxr(ixo^s))
730 cmin(ixo^s,1)=min(cminl(ixo^s),cminr(ixo^s))
733 cmaxl(ixo^s)=max(cmaxl(ixo^s),dabs(cminl(ixo^s)))
734 cmaxr(ixo^s)=max(cmaxr(ixo^s),dabs(cminr(ixo^s)))
735 cmax(ixo^s,1)=max(cmaxl(ixo^s),cmaxr(ixo^s))
743 tmp1=wmean(ixo^s,
xi_)/wmean(ixo^s,
lfac_)**2.0d0
744 call srhd_get_csound2_rhoh(wmean,x,ixi^
l,ixo^
l,tmp1,csound2)
745 vidim(ixo^s) = wmean(ixo^s,
mom(idim))/wmean(ixo^s,
xi_)
746 v2(ixo^s)=1.0d0-1.0d0/wmean(ixo^s,
lfac_)**2
747 call srhd_get_cmax_loc(ixi^
l,ixo^
l,vidim,csound2,v2,cmaxl,cminl)
748 if(
present(cmin))
then
749 cmax(ixo^s,1)=cmaxl(ixo^s)
750 cmin(ixo^s,1)=cminl(ixo^s)
752 cmax(ixo^s,1)=max(cmaxl(ixo^s),dabs(cminl(ixo^s)))
760 tmp1=wmean(ixo^s,
rho_)
761 tmp2=wmean(ixo^s,
xi_)/wmean(ixo^s,
lfac_)**2.0d0
763 call srhd_get_csound2_prim_eos(ixo^
l,tmp1,tmp2,tmp3,csound2)
764 vidim(ixo^s) = wmean(ixo^s,
mom(idim))/wmean(ixo^s,
lfac_)
765 v2(ixo^s) = 1.0d0-1.0d0/wmean(ixo^s,
lfac_)**2
766 call srhd_get_cmax_loc(ixi^
l,ixo^
l,vidim,csound2,v2,cmaxl,cminl)
767 if(
present(cmin))
then
768 cmax(ixo^s,1)=cmaxl(ixo^s)
769 cmin(ixo^s,1)=cminl(ixo^s)
771 cmax(ixo^s,1)=max(cmaxl(ixo^s),dabs(cminl(ixo^s)))
775 end subroutine srhd_get_cbounds
778 subroutine srhd_get_flux(wC,wP,x,ixI^L,ixO^L,idim,f)
780 integer,
intent(in) :: ixi^
l, ixo^
l, idim
782 double precision,
intent(in) :: wc(ixi^s,nw)
784 double precision,
intent(in) :: wp(ixi^s,nw)
785 double precision,
intent(in) :: x(ixi^s,1:
ndim)
786 double precision,
intent(out) :: f(ixi^s,nwflux)
788 double precision :: pth(ixi^s)
789 double precision :: v(ixi^s,1:
ndir)
792 pth(ixo^s)=wp(ixo^s,
p_)
794 v(ixo^s,idir) = wp(ixo^s,
mom(idir))/wp(ixo^s,
lfac_)
798 f(ixo^s,
d_)=v(ixo^s,idim)*wc(ixo^s,
rho_)
808 f(ixo^s,
mom(idir))= v(ixo^s,idim)*wc(ixo^s,
mom(idir))
810 f(ixo^s,
mom(idim))=pth(ixo^s)+f(ixo^s,
mom(idim))
814 f(ixo^s,
e_)=v(ixo^s,idim)*(wc(ixo^s,
e_) + pth(ixo^s))
816 end subroutine srhd_get_flux
819 subroutine srhd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
823 integer,
intent(in) :: ixi^
l, ixo^
l
824 double precision,
intent(in) :: qdt, dtfactor, x(ixi^s, 1:
ndim)
825 double precision,
intent(inout) :: wct(ixi^s, 1:nw), wprim(ixi^s, 1:nw), w(ixi^s, 1:nw)
827 double precision :: pth(ixi^s),
source(ixi^s), v(ixi^s,1:
ndir)
828 integer :: idir, h1x^
l{^nooned, h2x^
l}
830 double precision :: exp_factor(ixi^s), del_exp_factor(ixi^s), exp_factor_primitive(ixi^s)
841 source(ixo^s) =
source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
852 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
853 source(ixo^s) = -wct(ixo^s, mphi_) * wprim(ixo^s,
mom(
r_))
854 w(ixo^s, mphi_) = w(ixo^s, mphi_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
856 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s,
r_)
861 h1x^
l=ixo^
l-
kr(1,^
d); {^nooned h2x^
l=ixo^
l-
kr(2,^
d);}
866 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
867 *(
block%surfaceC(ixo^s, 1) -
block%surfaceC(h1x^s, 1)) &
868 /
block%dvolume(ixo^s)
874 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt *
source(ixo^s) / x(ixo^s, 1)
878 source(ixo^s) = pth(ixo^s) * x(ixo^s, 1) &
879 * (
block%surfaceC(ixo^s, 2) -
block%surfaceC(h2x^s, 2)) &
880 /
block%dvolume(ixo^s)
885 w(ixo^s,
mom(2)) = w(ixo^s,
mom(2)) + qdt *
source(ixo^s) / x(ixo^s, 1)
889 source(ixo^s) = -(wct(ixo^s,
mom(3)) * wprim(ixo^s,
mom(1))) &
890 - (wct(ixo^s,
mom(3)) * wprim(ixo^s,
mom(2))) / dtan(x(ixo^s, 2))
891 w(ixo^s,
mom(3)) = w(ixo^s,
mom(3)) + qdt *
source(ixo^s) / x(ixo^s, 1)
896 end subroutine srhd_add_source_geom
899 subroutine srhd_handle_small_values(primitive, w, x, ixI^L, ixO^L, subname)
902 logical,
intent(in) :: primitive
903 integer,
intent(in) :: ixi^
l,ixo^
l
904 double precision,
intent(inout) :: w(ixi^s,1:nw)
905 double precision,
intent(in) :: x(ixi^s,1:
ndim)
906 character(len=*),
intent(in) :: subname
909 logical :: flag(ixi^s,1:nw),flagall(ixi^s)
917 flagall(ixo^s)=(flag(ixo^s,
rho_).or.flag(ixo^s,
e_))
919 where(flagall(ixo^s))
931 where(flagall(ixo^s)) w(ixo^s,
e_) =
small_e
951 if(.not.primitive)
then
953 write(*,*)
"handle_small_values default: note reporting conservatives!"
959 end subroutine srhd_handle_small_values
969 integer,
intent(in) :: ixi^
l,ixo^
l
970 double precision,
intent(in) :: w(ixi^s, 1:nw)
971 logical,
intent(in) :: varconserve
972 double precision,
intent(out) :: geff(ixi^s)
974 double precision,
dimension(ixO^S) :: pth,rho,e_th,e
977 if (varconserve)
then
978 pth(ixo^s)=w(ixo^s,
xi_)-w(ixo^s,
e_)-w(ixo^s,
d_)
979 rho(ixo^s)=w(ixo^s,
d_)/w(ixo^s,
lfac_)
981 e = e_th+dsqrt(e_th**2+rho**2)
983 (one-(rho(ixo^s)/e(ixo^s))**2)
987 e = e_th+dsqrt(e_th**2+w(ixo^s,
rho_)**2)
989 (one-(w(ixo^s,
rho_)/e(ixo^s))**2)
998 subroutine srhd_get_smallvalues_eos
1002 double precision :: lsmalle,lsmallp,lsmallrho
1009 call mpistop(
"must set finite values small-density/pressure for small value treatments")
1019 gamma_1*lsmallrho*(lsmallrho/lsmalle))
1028 end subroutine srhd_get_smallvalues_eos
1033 integer,
intent(in) :: ixo^
l
1034 double precision,
intent(in) :: rho(ixo^s),p(ixo^s)
1035 double precision,
intent(out) :: rhoh(ixo^s)
1037 double precision,
dimension(ixO^S) :: e_th,e
1042 e = e_th+dsqrt(e_th**2+rho**2)
1051 {
do ix^db= ixo^lim^db\}
1053 write(*,*)
"local pressure and density",p(ix^
d),rho(ix^
d)
1054 write(*,*)
"Error: small value of enthalpy rho*h=",rhoh(ix^
d),&
1055 " encountered when call srhd_get_enthalpy_eos"
1056 call mpistop(
'enthalpy below small_xi: stop (may need to turn on fixes)')
1062 {
do ix^db= ixo^lim^db\}
1073 subroutine srhd_get_pressure_eos(ixI^L,ixO^L,rho,rhoh,p,E)
1075 integer,
intent(in) :: ixi^
l, ixo^
l
1076 double precision,
intent(in) :: rho(ixo^s),rhoh(ixo^s)
1077 double precision,
intent(out) :: p(ixi^s)
1078 double precision,
intent(out) :: e(ixo^s)
1082 e = (rhoh+dsqrt(rhoh**2+(
srhd_gamma**2-one)*rho**2)) &
1084 p(ixo^s) = half*
gamma_1* (e-rho*(rho/e))
1090 {
do ix^db= ixo^lim^db\}
1092 write(*,*)
"local enthalpy rho*h and density rho",rhoh(ix^
d),rho(ix^
d)
1093 if(
srhd_eos)
write(*,*)
'E, rho^2/E, difference', &
1094 e(ix^
d),rho(ix^
d)**2/e(ix^
d),e(ix^
d)-rho(ix^
d)**2/e(ix^
d)
1095 write(*,*)
"Error: small value of gas pressure",p(ix^
d),&
1096 " encountered when call srhd_get_pressure_eos"
1097 call mpistop(
'pressure below small_pressure: stop (may need to turn on fixes)')
1103 {
do ix^db= ixo^lim^db\}
1111 end subroutine srhd_get_pressure_eos
1115 subroutine srhd_get_csound2_eos(ixI^L,ixO^L,rho,rhoh,csound2)
1117 integer,
intent(in) :: ixi^
l, ixo^
l
1118 double precision,
intent(in) :: rho(ixo^s),rhoh(ixo^s)
1119 double precision,
intent(out) :: csound2(ixo^s)
1121 double precision :: p(ixi^s)
1122 double precision :: e(ixo^s)
1125 call srhd_get_pressure_eos(ixi^
l,ixo^
l,rho,rhoh,p,e)
1128 +
gamma_1*(rho(ixo^s)/e(ixo^s))**2))&
1129 /(2.0d0*rhoh(ixo^s))
1131 csound2(ixo^s)=
srhd_gamma*p(ixo^s)/rhoh(ixo^s)
1135 {
do ix^db= ixo^lim^db\}
1136 if(csound2(ix^
d)>=1.0d0.or.csound2(ix^
d)<=0.0d0)
then
1137 write(*,*)
"sound speed error with p - rho - rhoh",p(ix^
d),rhoh(ix^
d),rho(ix^
d)
1138 if(
srhd_eos)
write(*,*)
'and E', e(ix^
d)
1139 write(*,*)
"Error: value of csound2",csound2(ix^
d),&
1140 " encountered when call srhd_get_csound2_eos"
1141 call mpistop(
'sound speed stop (may need to turn on fixes)')
1147 {
do ix^db= ixo^lim^db\}
1148 if(csound2(ix^
d)>=1.0d0)
then
1149 csound2(ix^
d)=1.0d0-1.0d0/
lfacmax**2
1151 if(csound2(ix^
d)<=0.0d0)
then
1157 end subroutine srhd_get_csound2_eos
1161 subroutine srhd_get_csound2_prim_eos(ixO^L,rho,rhoh,p,csound2)
1163 integer,
intent(in) :: ixo^
l
1164 double precision,
intent(in) :: rho(ixo^s),rhoh(ixo^s),p(ixo^s)
1165 double precision,
intent(out) :: csound2(ixo^s)
1167 double precision :: e(ixo^s)
1181 {
do ix^db= ixo^lim^db\}
1182 if(csound2(ix^
d)>=1.0d0.or.csound2(ix^
d)<=0.0d0)
then
1183 write(*,*)
"sound speed error with p - rho - rhoh",p(ix^
d),rhoh(ix^
d),rho(ix^
d)
1184 if(
srhd_eos)
write(*,*)
'and E', e(ix^
d)
1185 write(*,*)
"Error: value of csound2",csound2(ix^
d),&
1186 " encountered when call srhd_get_csound2_prim_eos"
1187 call mpistop(
'sound speed stop (may need to turn on fixes)')
1193 {
do ix^db= ixo^lim^db\}
1194 if(csound2(ix^
d)>=1.0d0)
then
1195 csound2(ix^
d)=1.0d0-1.0d0/
lfacmax**2
1197 if(csound2(ix^
d)<=0.0d0)
then
1203 end subroutine srhd_get_csound2_prim_eos
1206 subroutine con2prim_eos(lfac,xi,myd,myssqr,mytau,ierror)
1209 double precision,
intent(in) :: myd, myssqr, mytau
1210 double precision,
intent(inout) :: lfac, xi
1211 integer,
intent(inout) :: ierror
1214 double precision:: f,df,lfacl
1218 d = myd;
ssqr = myssqr;
tau = mytau;
1224 call funcd_eos(xi,f,df,lfacl,
d,
ssqr,
tau,ierror)
1225 if (ierror == 0 .and. dabs(f/df)<
absaccnr)
then
1233 write(*,*)
'entering con2prim_eos with xi=',xi
1242 call con2primhydro_eos(lfac,xi,
d,
ssqr,
tau,ierror)
1244 end subroutine con2prim_eos
1246 subroutine funcd_eos(xi,f,df,mylfac,d,ssqr,tau,ierror)
1247 double precision,
intent(in) :: xi,d,ssqr,tau
1248 double precision,
intent(out) :: f,df,mylfac
1249 integer,
intent(inout) :: ierror
1252 double precision :: dlfac
1253 double precision :: vsqr,p,dpdxi
1259 mylfac = one/dsqrt(one-vsqr)
1260 dlfac = -mylfac**3*ssqr/(xi**3)
1262 call funcpressure_eos(xi,mylfac,d,dlfac,p,dpdxi)
1273 end subroutine funcd_eos
1276 subroutine con2primhydro_eos(lfac,xi,d,sqrs,tau,ierror)
1277 double precision,
intent(out) :: xi,lfac
1278 double precision,
intent(in) :: d,sqrs,tau
1279 integer,
intent(inout) :: ierror
1282 integer :: ni,niiter,nit,n2it,ni3
1283 double precision :: pcurrent,pnew
1284 double precision :: er,er1,ff,df,dp,v2
1285 double precision :: pmin,lfac2inv,plabs,prabs,pprev
1286 double precision :: s2overcubeg2rh
1287 double precision :: xicurrent,h,dhdp
1288 double precision :: oldff1,oldff2,nff
1289 double precision :: pleft,pright
1304 pmin=dsqrt(sqrs)/(one-
dmaxvel)-tau-d
1305 plabs=max(
minp,pmin)
1324 pcurrent=half*(pcurrent+pprev)
1332 xicurrent=tau+d+pcurrent
1344 v2=sqrs/xicurrent**2
1346 if(lfac2inv>zero)
then
1347 lfac=one/dsqrt(lfac2inv)
1359 s2overcubeg2rh=sqrs/(xicurrent**3)
1361 call funcenthalpy_eos(pcurrent,lfac2inv,d,sqrs,xicurrent,&
1362 s2overcubeg2rh,h,dhdp)
1364 ff=-xicurrent*lfac2inv + h
1365 df=- two*sqrs/xicurrent**2 + dhdp - lfac2inv
1367 if (ff*df==zero)
then
1377 if (ff*df>zero)
then
1380 pnew=max(pnew,plabs)
1384 pnew=min(pnew,prabs)
1390 er=two*dabs(dp)/(pnew+pcurrent)
1396 if((dabs(oldff2-ff) < 1.0d-8 .or. niiter >=
maxitnr-
maxitnr/20).and.&
1397 ff * oldff1 < zero .and. dabs(ff)>
absaccnr)
then
1400 if(n2it<=3) pcurrent=half*(pnew+pcurrent)
1404 pcurrent=half*(pleft+pright)
1407 xicurrent=tau+d+pcurrent
1408 v2=sqrs/xicurrent**2
1411 if(lfac2inv>zero)
then
1412 lfac=one/dsqrt(lfac2inv)
1420 call bisection_enthalpy_eos(pnew,lfac2inv,d,xicurrent,h)
1421 nff=-xicurrent*lfac2inv + h
1424 if(ff * nff < zero)
then
1430 pcurrent=half*(pleft+pright)
1434 if(2.0d0*dabs(pleft-pright)/(pleft+pright)<
absaccnr &
1468 if(pcurrent<
minp)
then
1475 v2=sqrs/xicurrent**2
1477 if(lfac2inv>zero)
then
1478 lfac=one/dsqrt(lfac2inv)
1484 end subroutine con2primhydro_eos
1488 subroutine funcpressure_eos(xicurrent,lfac,d,dlfacdxi,p,dpdxi)
1490 double precision,
intent(in) :: xicurrent,lfac,d,dlfacdxi
1491 double precision,
intent(out) :: p,dpdxi
1493 double precision :: rho,h,e,dhdxi,rhotoe
1494 double precision :: dpdchi,dedxi
1497 h=xicurrent/(lfac**2)
1499 e = (h+dsqrt(h**2+(
srhd_gamma**2-one)*rho**2)) &
1503 p = half*
gamma_1*(e-rho*rhotoe)
1505 dhdxi = one/(lfac**2)-2.0d0*xicurrent/(lfac**2)*dlfacdxi/lfac
1507 dedxi=(dhdxi+(h*dhdxi-(
srhd_gamma**2-one)*rho**2*dlfacdxi/lfac)&
1512 dpdxi=half*
gamma_1*(2.0d0*rho*rhotoe*dlfacdxi/lfac+&
1513 (one+rhotoe**2)*dedxi)
1515 end subroutine funcpressure_eos
1519 subroutine funcenthalpy_eos(pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p,h,dhdp)
1521 double precision,
intent(in) :: pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p
1522 double precision,
intent(out):: h,dhdp
1525 double precision:: rho,e_th,e,de_thdp,dedp
1527 rho=d*dsqrt(lfac2inv)
1529 e = (e_th + dsqrt(e_th**2+rho**2))
1535 dedp = de_thdp * (one+e_th/dsqrt(e_th**2+rho**2))&
1536 + d**2*dv2d2p/dsqrt(e_th**2+rho**2)
1539 gamma_1*(rho*(rho/e))*(-2.0d0*dv2d2p/lfac2inv+dedp/e))
1540 end subroutine funcenthalpy_eos
1544 subroutine bisection_enthalpy_eos(pcurrent,lfac2inv,d,xicurrent,h)
1546 double precision,
intent(in) :: pcurrent,lfac2inv,d,xicurrent
1547 double precision,
intent(out):: h
1550 double precision:: rho,e_th,e
1552 rho=d*dsqrt(lfac2inv)
1554 e = (e_th + dsqrt(e_th**2+rho**2))
1559 end subroutine bisection_enthalpy_eos
1562 subroutine con2prim(lfac,xi,myd,myssqr,mytau,ierror)
1565 double precision,
intent(in) :: myd, myssqr, mytau
1566 double precision,
intent(inout) :: lfac, xi
1567 integer,
intent(inout) :: ierror
1570 double precision:: f,df,lfacl
1574 d = myd;
ssqr = myssqr;
tau = mytau;
1580 call funcd(xi,f,df,lfacl,
d,
ssqr,
tau,ierror)
1581 if (ierror == 0 .and. dabs(f/df)<
absaccnr)
then
1589 write(*,*)
'entering con2prim with xi=',xi
1598 call con2primhydro(lfac,xi,
d,
ssqr,
tau,ierror)
1600 end subroutine con2prim
1602 subroutine funcd(xi,f,df,mylfac,d,ssqr,tau,ierror)
1603 double precision,
intent(in) :: xi,d,ssqr,tau
1604 double precision,
intent(out) :: f,df,mylfac
1605 integer,
intent(inout) :: ierror
1608 double precision :: dlfac
1609 double precision :: vsqr,p,dpdxi
1615 mylfac = one/dsqrt(one-vsqr)
1616 dlfac = -mylfac**3*ssqr/(xi**3)
1618 call funcpressure(xi,mylfac,d,dlfac,p,dpdxi)
1629 end subroutine funcd
1632 subroutine con2primhydro(lfac,xi,d,sqrs,tau,ierror)
1633 double precision,
intent(out) :: xi,lfac
1634 double precision,
intent(in) :: d,sqrs,tau
1635 integer,
intent(inout) :: ierror
1638 integer :: ni,niiter,nit,n2it,ni3
1639 double precision :: pcurrent,pnew
1640 double precision :: er,er1,ff,df,dp,v2
1641 double precision :: pmin,lfac2inv,plabs,prabs,pprev
1642 double precision :: s2overcubeg2rh
1643 double precision :: xicurrent,h,dhdp
1644 double precision :: oldff1,oldff2,nff
1645 double precision :: pleft,pright
1660 pmin=dsqrt(sqrs)/(one-
dmaxvel)-tau-d
1661 plabs=max(
minp,pmin)
1680 pcurrent=half*(pcurrent+pprev)
1688 xicurrent=tau+d+pcurrent
1700 v2=sqrs/xicurrent**2
1702 if(lfac2inv>zero)
then
1703 lfac=one/dsqrt(lfac2inv)
1715 s2overcubeg2rh=sqrs/(xicurrent**3)
1717 call funcenthalpy(pcurrent,lfac2inv,d,sqrs,xicurrent,&
1718 s2overcubeg2rh,h,dhdp)
1720 ff=-xicurrent*lfac2inv + h
1721 df=- two*sqrs/xicurrent**2 + dhdp - lfac2inv
1723 if (ff*df==zero)
then
1733 if (ff*df>zero)
then
1736 pnew=max(pnew,plabs)
1740 pnew=min(pnew,prabs)
1746 er=two*dabs(dp)/(pnew+pcurrent)
1752 if((dabs(oldff2-ff) < 1.0d-8 .or. niiter >=
maxitnr-
maxitnr/20).and.&
1753 ff * oldff1 < zero .and. dabs(ff)>
absaccnr)
then
1756 if(n2it<=3) pcurrent=half*(pnew+pcurrent)
1760 pcurrent=half*(pleft+pright)
1763 xicurrent=tau+d+pcurrent
1764 v2=sqrs/xicurrent**2
1767 if(lfac2inv>zero)
then
1768 lfac=one/dsqrt(lfac2inv)
1776 call bisection_enthalpy(pnew,lfac2inv,d,xicurrent,h)
1777 nff=-xicurrent*lfac2inv + h
1780 if(ff * nff < zero)
then
1786 pcurrent=half*(pleft+pright)
1790 if(2.0d0*dabs(pleft-pright)/(pleft+pright)<
absaccnr &
1824 if(pcurrent<
minp)
then
1831 v2=sqrs/xicurrent**2
1833 if(lfac2inv>zero)
then
1834 lfac=one/dsqrt(lfac2inv)
1840 end subroutine con2primhydro
1844 subroutine funcpressure(xicurrent,lfac,d,dlfacdxi,p,dpdxi)
1846 double precision,
intent(in) :: xicurrent,lfac,d,dlfacdxi
1847 double precision,
intent(out) :: p,dpdxi
1849 double precision :: rho,h,e,dhdxi,rhotoe
1850 double precision :: dpdchi,dedxi
1853 h=xicurrent/(lfac**2)
1858 dpdxi = dpdchi * one/lfac**2
1860 if (dlfacdxi /= 0.0d0) &
1861 dpdxi = dpdxi + dpdchi * ((d*lfac-2.0d0*xicurrent)/lfac**3) * dlfacdxi
1863 end subroutine funcpressure
1867 subroutine funcenthalpy(pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p,h,dhdp)
1869 double precision,
intent(in) :: pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p
1870 double precision,
intent(out):: h,dhdp
1873 double precision:: rho,e_th,e,de_thdp,dedp
1875 rho=d*dsqrt(lfac2inv)
1878 end subroutine funcenthalpy
1882 subroutine bisection_enthalpy(pcurrent,lfac2inv,d,xicurrent,h)
1884 double precision,
intent(in) :: pcurrent,lfac2inv,d,xicurrent
1885 double precision,
intent(out):: h
1888 double precision:: rho,e_th,e
1890 rho=d*dsqrt(lfac2inv)
1894 end subroutine bisection_enthalpy
Calculate w(iw)=w(iw)+qdt*SOURCE[wCT,qtC,x] within ixO for all indices iw=iwmin......
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter spherical
integer, parameter cylindrical
integer, parameter cartesian_expansion
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 small_pressure
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.
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.
logical use_particles
Use particles module or not.
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
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision unit_velocity
Physical scaling factor for velocity.
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 small_density
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
integer boundspeed
bound (left/min and right.max) speed of Riemann fan
logical check_small_values
check and optionally fix unphysical small values (density, gas pressure)
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
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)
character(len=20), public small_values_method
How to handle small values.
Special Relativistic Hydrodynamics (with EOS) physics module.
subroutine, public srhd_get_auxiliary(ixil, ixol, w, x)
Compute auxiliary variables lfac and xi from a conservative state using srhd_con2prim to calculate en...
double precision, public tolernr
integer, public, protected srhd_n_tracer
Number of tracer species.
double precision, public inv_gamma_1
subroutine, public srhd_get_geff_eos(w, ixil, ixol, varconserve, geff)
calculate effective gamma
integer, dimension(:), allocatable, public, protected tracer
Indices of the tracers.
double precision, public absaccnr
double precision, public small_e
The smallest allowed energy.
integer, public, protected lfac_
Index of the Lorentz factor.
subroutine, public srhd_get_enthalpy_eos(ixol, rho, p, rhoh)
Compute the enthalpy rho*h from rho and pressure p.
subroutine, public srhd_get_auxiliary_prim(ixil, ixol, w)
Set auxiliary variables lfac and xi from a primitive state only used when handle_small_values average...
double precision, public gamma_to_gamma_1
integer, public maxitnr
parameters for NR in con2prim
subroutine, public srhd_phys_init()
Initialize the module.
double precision, public gamma_1
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
double precision, public small_xi
The smallest allowed inertia.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine, public srhd_get_pthermal(w, x, ixil, ixol, pth)
Calculate thermal pressure p within ixO^L must follow after update conservative with auxiliaries.
double precision, public minp
logical, public srhd_eos
Whether synge eos is used.
double precision, public lfacmax
integer, public, protected d_
subroutine, public srhd_check_w(primitive, ixil, ixol, w, flag)
Returns logical argument flag T where values are not ok.
subroutine, public srhd_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
subroutine, public srhd_check_params
integer, public, protected p_
Index of the gas pressure should equal e_.
subroutine, public srhd_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
logical, public, protected srhd_particles
Whether particles module is added.
double precision, public smalltau
double precision, public dmaxvel
integer, public, protected e_
Index of the energy density.
double precision, public minrho
double precision, public smallxi
integer, public, protected xi_
Index of the inertia.
double precision, public srhd_gamma
The adiabatic index and derived values.
integer, public, protected rho_
Index of the density (in the w array)
subroutine, public srhd_get_csound2(w, x, ixil, ixol, csound2)
Calculate the square of the thermal sound speed csound2 within ixO^L. here computed from conservative...
Module with all the methods that users can customize in AMRVAC.
procedure(set_surface), pointer usr_set_surface
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
integer nwflux
Number of flux variables.