MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_srhd_phys.t
Go to the documentation of this file.
1!> Special Relativistic Hydrodynamics (with EOS) physics module
5 use mod_comm_lib, only: mpistop
6 implicit none
7 private
8
9 !> Whether particles module is added
10 logical, public, protected :: srhd_particles = .false.
11
12 !> Number of tracer species
13 integer, public,protected :: srhd_n_tracer = 0
14
15 !> Index of the density (in the w array)
16 integer, public,protected :: rho_
17 integer, public,protected :: d_
18
19 !> Indices of the momentum density
20 integer, allocatable, public, protected :: mom(:)
21
22 !> Indices of the tracers
23 integer, allocatable, public, protected :: tracer(:)
24
25 !> Index of the energy density
26 integer, public,protected :: e_
27 !> Index of the gas pressure should equal e_
28 integer, public,protected :: p_
29
30 !> Index of the Lorentz factor
31 integer, public,protected :: lfac_
32
33 !> Index of the inertia
34 integer, public,protected :: xi_
35
36 !> Whether synge eos is used
37 logical, public :: srhd_eos = .false.
38
39 !> The adiabatic index and derived values
40 double precision, public :: srhd_gamma = 5.d0/3.0d0
41 double precision, public :: gamma_1,inv_gamma_1,gamma_to_gamma_1
42
43 !> The smallest allowed energy
44 double precision, public :: small_e
45 !> The smallest allowed inertia
46 double precision, public :: small_xi
47
48 !> Helium abundance over Hydrogen
49 double precision, public, protected :: he_abundance=0.0d0
50
51 !> parameters for NR in con2prim
52 integer, public :: maxitnr = 100
53 double precision, public :: absaccnr = 1.0d-8
54 double precision, public :: tolernr = 1.0d-9
55 double precision, public :: dmaxvel = 1.0d-7
56 double precision, public :: lfacmax
57 double precision, public :: minp, minrho, smalltau, smallxi
58
59 ! Public methods
60 public :: srhd_phys_init
61 public :: srhd_get_pthermal
62 public :: srhd_get_auxiliary
64 public :: srhd_get_csound2
65 public :: srhd_to_conserved
66 public :: srhd_to_primitive
67 public :: srhd_check_params
68 public :: srhd_check_w
69 public :: srhd_get_geff_eos
70 public :: srhd_get_enthalpy_eos
71contains
72
73 !> Read this module's parameters from a file
74 subroutine srhd_read_params(files)
76 character(len=*), intent(in) :: files(:)
77 integer :: n
78
79 namelist /srhd_list/ srhd_n_tracer, srhd_eos, srhd_gamma, &
82
83 do n = 1, size(files)
84 open(unitpar, file=trim(files(n)), status="old")
85 read(unitpar, srhd_list, end=111)
86111 close(unitpar)
87 end do
88
89 end subroutine srhd_read_params
90
91 !> Write this module's parameters to a snapshot
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
99 integer :: er
100
101 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
102
103 names(1) = "gamma"
104 values(1) = srhd_gamma
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)
107
108 end subroutine srhd_write_info
109
110 !> Initialize the module
111 subroutine srhd_phys_init()
114 integer :: itr,idir
115
116 call srhd_read_params(par_files)
117
118 physics_type = "srhd"
119 phys_energy = .true.
120 phys_total_energy = .true.
121 phys_gamma = srhd_gamma
122
123 ! unused physics options
124 phys_internal_e = .false.
125 phys_partial_ionization=.false.
126 phys_trac=.false.
127
129
130 ! note: number_species is 1 for srhd
131 allocate(start_indices(number_species),stop_indices(number_species))
132 ! set the index of the first flux variable for species 1
133 start_indices(1)=1
134
135 ! Determine flux variables
136 rho_ = var_set_rho()
137 d_=rho_
138
139 allocate(mom(ndir))
140 mom(:) = var_set_momentum(ndir)
141
142 ! Set index of energy variable
143 e_ = var_set_energy()
144 p_ = e_
145
146 ! derive units from basic units
147 call srhd_physical_units()
148
149 allocate(tracer(srhd_n_tracer))
150
151 ! Set starting index of tracers
152 do itr = 1, srhd_n_tracer
153 tracer(itr) = var_set_fluxvar("trc", "trp", itr, need_bc=.false.)
154 end do
155
156 ! Set index for auxiliary variables
157 ! MUST be after the possible tracers (which have fluxes)
158 xi_ = var_set_auxvar('xi','xi')
159 lfac_= var_set_auxvar('lfac','lfac')
160
161 ! set number of variables which need update ghostcells
162 nwgc=nwflux+nwaux
163
164 ! set the index of the last flux variable for species 1
165 stop_indices(1)=nwflux
166
167 ! Number of variables need reconstruction in w
168 nw_recon=nwflux
169
170 ! Check whether custom flux types have been defined
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")
176 end if
177
178 nvector = 1 ! No. vector vars
179 allocate(iw_vector(nvector))
180 iw_vector(1) = mom(1) - 1
181
182 ! dummy for now, no extra source terms precoded
183 phys_add_source => srhd_add_source
184 phys_get_dt => srhd_get_dt
185 ! copied in from HD/MHD, for certain limiters
186 phys_get_a2max => srhd_get_a2max
187
188 ! actual srhd routines
189 phys_check_params => srhd_check_params
190 phys_check_w => srhd_check_w
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
195 phys_to_conserved => srhd_to_conserved
196 phys_to_primitive => srhd_to_primitive
197 phys_get_pthermal => srhd_get_pthermal
198 phys_get_auxiliary => srhd_get_auxiliary
199 phys_get_auxiliary_prim => srhd_get_auxiliary_prim
200 phys_get_pthermal => srhd_get_pthermal
201 phys_get_v => srhd_get_v
202 phys_write_info => srhd_write_info
203 phys_handle_small_values => srhd_handle_small_values
204
205 ! Initialize particles module
206 if (srhd_particles) then
207 call particles_init()
208 end if
209
210 end subroutine srhd_phys_init
211
214
215 if (srhd_gamma <= 0.0d0 .or. srhd_gamma == 1.0d0) &
216 call mpistop ("Error: srhd_gamma <= 0 or srhd_gamma == 1")
217 ! additional useful values
218 gamma_1=srhd_gamma-1.0d0
219 inv_gamma_1=1.0d0/gamma_1
221
222 ! the following sets small_e and small_xi from small_density/small_pressure
223 ! according to the srhd_eos used
224 call srhd_get_smallvalues_eos
225
226 if(mype==0)then
227 write(*,*)'------------------------------------------------------------'
228 write(*,*)'Using EOS set via srhd_eos=',srhd_eos
229 write(*,*)'Maximal lorentz factor (via dmaxvel) is=',lfacmax
230 write(*,*)'Use fixes set through check/fix small values:', check_small_values,fix_small_values
231 write(*,*)'Controlled with small pressure/density:', small_pressure,small_density
232 write(*,*)'Derived small values: xi and e ',small_xi,small_e
233 write(*,*)'------------------------------------------------------------'
234 endif
235
236 end subroutine srhd_check_params
237
238 subroutine srhd_physical_units
240 double precision :: mp,kb
241
242 ! Derive scaling units
243 if(si_unit) then
244 mp=mp_si
245 kb=kb_si
246 unit_velocity=c_si
247 else
248 mp=mp_cgs
249 kb=kb_cgs
250 unit_velocity=const_c
251 end if
252 if(unit_numberdensity*unit_length<=0.0d0)then
253 call mpistop("Abort: must set positive values for unit length and numberdensity")
254 endif
255 ! we assume user sets: unit_numberdensity, unit_length, He_abundance
256 ! then together with light speed c, all units fixed
262
263 end subroutine srhd_physical_units
264
265 !> Returns logical argument flag T where values are not ok
266 subroutine srhd_check_w(primitive, ixI^L, ixO^L, w, flag)
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)
272
273 flag=.false.
274
275 ! NOTE: we should not check or use nwaux variables here
276 if(primitive) then
277 where(w(ixo^s,rho_) < small_density) flag(ixo^s,rho_) = .true.
278 where(w(ixo^s,p_) < small_pressure) flag(ixo^s,p_) = .true.
279 else
280 where(w(ixo^s,d_) < small_density) flag(ixo^s,d_) = .true.
281 where(w(ixo^s,e_) < small_e) flag(ixo^s,e_) = .true.
282 endif
283
284 end subroutine srhd_check_w
285
286 !> Returns logical argument flag T where auxiliary values are not ok
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)
292
293 flag=.false.
294
295 where(w(ixo^s,xi_) < small_xi) flag(ixo^s,xi_) = .true.
296 where(w(ixo^s,lfac_) < one) flag(ixo^s,lfac_) = .true.
297
298 if(any(flag(ixo^s,xi_)))then
299 write(*,*)'auxiliary xi too low: abort'
300 call mpistop('auxiliary check failed')
301 endif
302 if(any(flag(ixo^s,lfac_)))then
303 write(*,*)'auxiliary lfac too low: abort'
304 call mpistop('auxiliary check failed')
305 endif
306
307 end subroutine srhd_check_w_aux
308
309 !> Set auxiliary variables lfac and xi from a primitive state
310 !> only used when handle_small_values average on primitives
311 subroutine srhd_get_auxiliary_prim(ixI^L,ixO^L,w)
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
316
317 ! assume four-velocity in momentum vector (i.e. lfac*v)
318 rho(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
319 w(ixo^s,lfac_) = dsqrt(1.0d0+rho(ixo^s))
320
321 rho(ixo^s)=w(ixo^s,rho_)
322 pth(ixo^s)=w(ixo^s,p_)
323 ! compute rho*h (enthalpy h) from density-pressure
324 call srhd_get_enthalpy_eos(ixo^l,rho,pth,rhoh)
325
326 ! fill auxiliary variable xi= lfac^2 rhoh
327 w(ixo^s,xi_) = w(ixo^s,lfac_)**2.0d0*rhoh(ixo^s)
328
329 end subroutine srhd_get_auxiliary_prim
330
331 !> Compute auxiliary variables lfac and xi from a conservative state
332 !> using srhd_con2prim to calculate enthalpy and lorentz factor
333 subroutine srhd_get_auxiliary(ixI^L,ixO^L,w,x)
335 implicit none
336
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)
340
341 integer :: ix^d,ierror,idir
342 integer :: flag_error(ixo^s)
343 double precision :: ssqr
344
345 if(srhd_eos)then
346 {do ix^db=ixomin^db,ixomax^db\}
347 ierror=0
348 ssqr=0.0d0
349 do idir=1,ndir
350 ssqr= ssqr+w(ix^d,mom(idir))**2
351 enddo
352 if(w(ix^d,d_)<small_density)then
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)
357 if(check_small_values) call mpistop('small density on entry con2prim')
359 endif
360 if(w(ix^d,e_)<small_e)then
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)
365 if(check_small_values) call mpistop('small energy on entry con2prim')
366 if(fix_small_values) w(ix^d,e_)=small_e
367 endif
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
371 {enddo^d&\}
372 else
373 {do ix^db=ixomin^db,ixomax^db\}
374 ierror=0
375 ssqr=0.0d0
376 do idir=1,ndir
377 ssqr= ssqr+w(ix^d,mom(idir))**2
378 enddo
379 if(w(ix^d,d_)<small_density)then
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)
384 if(check_small_values) call mpistop('small density on entry con2prim')
386 endif
387 if(w(ix^d,e_)<small_e)then
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)
392 if(check_small_values) call mpistop('small energy on entry con2prim')
393 if(fix_small_values) w(ix^d,e_)=small_e
394 endif
395 call con2prim(w(ix^d,lfac_),w(ix^d,xi_), &
396 w(ix^d,d_),ssqr,w(ix^d,e_),ierror)
397 flag_error(ix^d) = ierror
398 {enddo^d&\}
399 endif
400
401 if(check_small_values)then
402 if(any(flag_error(ixo^s)/=0))then
403 print *,flag_error(ixo^s)
404 call mpistop('Problem when getting auxiliaries')
405 ! call srhd_handle_small_values(.false.,w,x,ixI^L,ixO^L,'srhd_get_auxiliary')
406 end if
407 end if
408
409 end subroutine srhd_get_auxiliary
410
411 !> Transform primitive variables into conservative ones
412 subroutine srhd_to_conserved(ixI^L, ixO^L, w, x)
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)
417 integer :: idir,itr
418 double precision, dimension(ixO^S) :: rhoh,rho,pth
419
420 ! assume four-velocity in momentum vector (i.e. lfac*v)
421 ! use rhoh slot for temporary array
422 rhoh(ixo^s) = sum(w(ixo^s, mom(:))**2, dim=ndim+1)
423 w(ixo^s,lfac_) = dsqrt(1.0d0+rhoh(ixo^s))
424
425 rho(ixo^s)=w(ixo^s,rho_)
426 pth(ixo^s)=w(ixo^s,p_)
427 ! compute rho*h (enthalpy h) from density-pressure
428 call srhd_get_enthalpy_eos(ixo^l,rho,pth,rhoh)
429
430 ! compute rhoh*lfac (recycle rhoh)
431 rhoh(ixo^s)= rhoh(ixo^s)*w(ixo^s,lfac_)
432 ! fill auxiliary variable xi= lfac^2 rhoh
433 w(ixo^s,xi_) = w(ixo^s,lfac_)*rhoh(ixo^s)
434
435 ! set conservative density: d = lfac * rho
436 w(ixo^s,d_)=w(ixo^s,lfac_)*rho(ixo^s)
437
438 ! Convert four-velocity (lfac*v) to momentum (xi*v=[rho*h*lfac^2]*v)
439 do idir = 1, ndir
440 w(ixo^s, mom(idir)) = rhoh(ixo^s)*w(ixo^s, mom(idir))
441 end do
442
443 ! set tau = xi-p-d energy variable
444 w(ixo^s,e_) = w(ixo^s,xi_)-pth(ixo^s)-w(ixo^s,d_)
445
446 do itr=1,srhd_n_tracer
447 w(ixo^s,tracer(itr)) = w(ixo^s,d_)*w(ixo^s,tracer(itr))
448 end do
449
450 end subroutine srhd_to_conserved
451
452 !> Transform conservative variables into primitive ones
453 subroutine srhd_to_primitive(ixI^L, ixO^L, w, x)
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)
458 integer :: idir,itr
459 double precision, dimension(ixO^S) :: rho,rhoh,e
460 double precision, dimension(ixI^S) :: pth
461 character(len=30) :: subname_loc
462
463 ! get auxiliary variables lfac and xi from conserved set
464 call srhd_get_auxiliary(ixi^l,ixo^l,w,x)
465
466 ! from d to rho (d=rho*lfac)
467 rho(ixo^s) = w(ixo^s,d_)/w(ixo^s,lfac_)
468
469 ! compute pressure
470 ! deduce rho*h from xi/lfac^2
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)
473
474 w(ixo^s,rho_)=rho(ixo^s)
475 ! from xi*v to U=lfac*v (four-velocity)
476 do idir=1,ndir
477 w(ixo^s,mom(idir)) = w(ixo^s,lfac_)*w(ixo^s,mom(idir))&
478 /w(ixo^s,xi_)
479 end do
480 w(ixo^s,p_)=pth(ixo^s)
481
482 do itr=1,srhd_n_tracer
483 w(ixo^s,tracer(itr)) = w(ixo^s,tracer(itr)) &
484 /(rho(ixo^s)*w(ixo^s,lfac_))
485 end do
486
487 end subroutine srhd_to_primitive
488
489 !> Calculate v vector from conservatives
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)
495 integer :: idir
496
497 ! get v from xi*v
498 do idir=1,ndir
499 v(ixo^s,idir) = w(ixo^s, mom(idir))/w(ixo^s,xi_)
500 end do
501
502 end subroutine srhd_get_v
503
504 !> Calculate the square of the thermal sound speed csound2 within ixO^L.
505 !> here computed from conservative set WITH ADDED rho*h
506 !> local version: does not do con2prim
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)
513
514 double precision :: rho(ixo^s)
515
516 rho=w(ixo^s,d_)/w(ixo^s,lfac_)
517 call srhd_get_csound2_eos(ixi^l,ixo^l,rho,rhoh,csound2)
518
519 end subroutine srhd_get_csound2_rhoh
520
521 !> Calculate the square of the thermal sound speed csound2 within ixO^L.
522 !> here computed from conservative set and uses con2prim!!!
523 !> public version!
524 subroutine srhd_get_csound2(w,x,ixI^L,ixO^L,csound2)
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)
530
531 double precision :: rho(ixo^s),rhoh(ixo^s)
532
533 ! get auxiliary variables lfac and xi from conserved set
534 call srhd_get_auxiliary(ixi^l,ixo^l,w,x)
535 ! quantify rho, rho*h
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)
539
540 end subroutine srhd_get_csound2
541
542 !> Calculate thermal pressure p within ixO^L
543 !> must follow after update conservative with auxiliaries
544 subroutine srhd_get_pthermal(w, x, ixI^L, ixO^L, pth)
547
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)
552
553 integer :: iw, ix^d
554 double precision :: rho(ixo^s),rhoh(ixo^s),e(ixo^s)
555
556 ! quantify rho, rho*h, tau and get pthermal
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)
560
561 end subroutine srhd_get_pthermal
562
563 !> dummy addsource subroutine
564 ! w[iw]= w[iw]+qdt*S[wCT, qtC, x] where S is the source based on wCT within ixO
565 subroutine srhd_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
567
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
574
575 end subroutine srhd_add_source
576
577 !> dummy get_dt subroutine
578 subroutine srhd_get_dt(w, ixI^L, ixO^L, dtnew, dx^D, x)
580
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
585
586 dtnew = bigdouble
587
588 end subroutine srhd_get_dt
589
590 !> Calculate cmax_idim within ixO^L
591 !> used especially for setdt CFL limit
592 subroutine srhd_get_cmax(w, x, ixI^L, ixO^L, idim, cmax)
594
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)
598
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
602
603 logical :: flag(ixi^s,1:nw)
604
605 !!call srhd_check_w_aux(ixI^L, ixO^L, w, flag)
606
607 ! input w is in primitive form TODO use it
608 wc=w
609 call srhd_to_conserved(ixi^l, ixo^l, wc, x)
610 ! auxiliaries are filled here
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))
622 ! Limit by speed of light
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)
627 ! now take extremal value only for dt limit
628 cmax(ixo^s) = max(dabs(cmax(ixo^s)),dabs(cmin(ixo^s)))
629
630 end subroutine srhd_get_cmax
631
632 subroutine srhd_get_a2max(w,x,ixI^L,ixO^L,a2max)
634
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
640
641 a2=zero
642 do i = 1,ndim
643 !> 4th order
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
651 end do
652
653 end subroutine srhd_get_a2max
654
655 !> local version for recycling code when computing cmax-cmin
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)
664
665 double precision, dimension(ixI^S):: tmp1,tmp2
666
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))
673 ! Limit by speed of light
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))
677 ! Limit by speed of light
678 cmin(ixo^s) = max(cmin(ixo^s), - 1.0d0)
679 cmin(ixo^s) = min(cmin(ixo^s), 1.0d0)
680
681 end subroutine srhd_get_cmax_loc
682
683 !> Estimating bounds for the minimum and maximum signal velocities
684 !> here we will not use Hspeed at all (one species only)
685 subroutine srhd_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
687 use mod_variables
688
689 integer, intent(in) :: ixi^l, ixo^l, idim
690 ! conservative left and right status
691 double precision, intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
692 ! primitive left and right status
693 double precision, intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
694 double precision, intent(in) :: x(ixi^s, 1:ndim)
695 double precision, intent(inout) :: cmax(ixi^s,1:number_species)
696 double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
697 double precision, intent(in) :: hspeed(ixi^s,1:number_species)
698
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
702
703 logical :: flag(ixi^s,1:nw)
704
705 select case(boundspeed)
706 case(1) ! we do left-right first and take maximals
707 !!call srhd_check_w(.true.,ixI^L, ixO^L, wLp, flag)
708 !!call srhd_check_w_aux(ixI^L, ixO^L, wLp, flag)
709 tmp1=wlp(ixo^s,rho_)
710 tmp2=wlp(ixo^s,xi_)/wlp(ixo^s,lfac_)**2.0d0
711 tmp3=wlp(ixo^s,p_)
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)
716
717 !!call srhd_check_w(.true.,ixI^L, ixO^L, wRp, flag)
718 !!call srhd_check_w_aux(ixI^L, ixO^L, wRp, flag)
719 tmp1=wrp(ixo^s,rho_)
720 tmp2=wrp(ixo^s,xi_)/wrp(ixo^s,lfac_)**2.0d0
721 tmp3=wrp(ixo^s,p_)
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)
726
727 if(present(cmin))then
728 ! for HLL
729 cmax(ixo^s,1)=max(cmaxl(ixo^s),cmaxr(ixo^s))
730 cmin(ixo^s,1)=min(cminl(ixo^s),cminr(ixo^s))
731 else
732 ! for TVDLF
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))
736 endif
737 case(2) ! this is cmaxmean from conservatives
738 ! here we do arithmetic mean of conservative vars
739 wmean(ixo^s,1:nwflux)=0.5d0*(wlc(ixo^s,1:nwflux)+wrc(ixo^s,1:nwflux))
740 ! get auxiliary variables
741 call srhd_get_auxiliary(ixi^l,ixo^l,wmean,x)
742 ! here tmp1 is rhoh
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)
751 else
752 cmax(ixo^s,1)=max(cmaxl(ixo^s),dabs(cminl(ixo^s)))
753 endif
754 case(3) ! this is cmaxmean from primitives
755 ! here we do arithmetic mean of primitive vars
756 wmean(ixo^s,1:nwflux)=0.5d0*(wlp(ixo^s,1:nwflux)+wrp(ixo^s,1:nwflux))
757 ! get auxiliary variables for wmean (primitive array)
758 call srhd_get_auxiliary_prim(ixi^l,ixo^l,wmean)
759 ! here tmp1 is rhoh
760 tmp1=wmean(ixo^s,rho_)
761 tmp2=wmean(ixo^s,xi_)/wmean(ixo^s,lfac_)**2.0d0
762 tmp3=wmean(ixo^s,p_)
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)
770 else
771 cmax(ixo^s,1)=max(cmaxl(ixo^s),dabs(cminl(ixo^s)))
772 endif
773 end select
774
775 end subroutine srhd_get_cbounds
776
777 !> Calculate fluxes within ixO^L.
778 subroutine srhd_get_flux(wC,wP,x,ixI^L,ixO^L,idim,f)
780 integer, intent(in) :: ixi^l, ixo^l, idim
781 ! conservative w
782 double precision, intent(in) :: wc(ixi^s,nw)
783 ! primitive w
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)
787
788 double precision :: pth(ixi^s)
789 double precision :: v(ixi^s,1:ndir)
790 integer :: iw,idir
791
792 pth(ixo^s)=wp(ixo^s,p_)
793 do idir=1,ndir
794 v(ixo^s,idir) = wp(ixo^s, mom(idir))/wp(ixo^s,lfac_)
795 end do
796
797 ! Get flux of density d, namely D*v
798 f(ixo^s,d_)=v(ixo^s,idim)*wc(ixo^s,rho_)
799
800 ! Get flux of tracer
801 do iw=1,srhd_n_tracer
802 f(ixo^s,tracer(iw))=v(ixo^s,idim)*wc(ixo^s,tracer(iw))
803 end do
804
805 ! Get flux of momentum
806 ! f_i[m_k]=v_i*m_k [+pth if i==k]
807 do idir=1,ndir
808 f(ixo^s,mom(idir))= v(ixo^s,idim)*wc(ixo^s,mom(idir))
809 end do
810 f(ixo^s,mom(idim))=pth(ixo^s)+f(ixo^s,mom(idim))
811
812 ! Get flux of energy
813 ! f_i[e]=v_i*e+v_i*pth
814 f(ixo^s,e_)=v(ixo^s,idim)*(wc(ixo^s,e_) + pth(ixo^s))
815
816 end subroutine srhd_get_flux
817
818 !> Add geometrical source terms to w
819 subroutine srhd_add_source_geom(qdt, dtfactor, ixI^L, ixO^L, wCT, wprim, w, x)
822 use mod_geometry
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)
826
827 double precision :: pth(ixi^s), source(ixi^s), v(ixi^s,1:ndir)
828 integer :: idir, h1x^l{^nooned, h2x^l}
829 integer :: mr_,mphi_ ! Polar var. names
830 double precision :: exp_factor(ixi^s), del_exp_factor(ixi^s), exp_factor_primitive(ixi^s)
831
832 select case (coordinate)
833
835 !the user provides the functions of exp_factor and del_exp_factor
836 if(associated(usr_set_surface)) &
837 call usr_set_surface(ixi^l,x,block%dx,exp_factor,del_exp_factor,exp_factor_primitive)
838 ! get auxiliary variables lfac and xi from conserved set
839 call srhd_get_auxiliary(ixi^l,ixo^l,wct,x)
840 call srhd_get_pthermal(wct, x, ixi^l, ixo^l, source)
841 source(ixo^s) = source(ixo^s)*del_exp_factor(ixo^s)/exp_factor(ixo^s)
842 w(ixo^s,mom(1)) = w(ixo^s,mom(1)) + qdt*source(ixo^s)
843
844 case (cylindrical)
845 mr_ = mom(r_)
846 ! get auxiliary variables lfac and xi from conserved set
847 call srhd_get_auxiliary(ixi^l,ixo^l,wct,x)
848 call srhd_get_pthermal(wct, x, ixi^l, ixo^l, source)
849 if (phi_ > 0) then
850 mphi_ = mom(phi_)
851 source(ixo^s) = source(ixo^s) + wct(ixo^s, mphi_)*wprim(ixo^s,mom(phi_))
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_)
855 else
856 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, r_)
857 end if
858
859 case (spherical)
860 mr_ = mom(r_)
861 h1x^l=ixo^l-kr(1,^d); {^nooned h2x^l=ixo^l-kr(2,^d);}
862 ! s[mr]=((stheta*vtheta+sphi*vphi)+2*p)/r
863 ! get auxiliary variables lfac and xi from conserved set
864 call srhd_get_auxiliary(ixi^l,ixo^l,wct,x)
865 call srhd_get_pthermal(wct, x, ixi^l, ixo^l, pth)
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)
869 if (ndir > 1) then
870 do idir = 2, ndir
871 source(ixo^s) = source(ixo^s) + wct(ixo^s, mom(idir))*wprim(ixo^s,mom(idir))
872 end do
873 end if
874 w(ixo^s, mr_) = w(ixo^s, mr_) + qdt * source(ixo^s) / x(ixo^s, 1)
875
876 {^nooned
877 ! s[mtheta]=-(stheta*vr)/r+cot(theta)*(sphi*vphi+p)/r
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)
881 if (ndir == 3) then
882 source(ixo^s) = source(ixo^s) + (wct(ixo^s, mom(3))*v(ixo^s,ndir)) / dtan(x(ixo^s, 2))
883 end if
884 source(ixo^s) = source(ixo^s) - wct(ixo^s, mom(2)) * wprim(ixo^s, mom(1))
885 w(ixo^s, mom(2)) = w(ixo^s, mom(2)) + qdt * source(ixo^s) / x(ixo^s, 1)
886
887 if (ndir == 3) then
888 ! s[mphi]=-(sphi*vr)/r-cot(theta)*(sphi*vtheta)/r
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)
892 end if
893 }
894 end select
895
896 end subroutine srhd_add_source_geom
897
898 !> handles bootstrapping
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
907
908 integer :: n,idir
909 logical :: flag(ixi^s,1:nw),flagall(ixi^s)
910
911 call srhd_check_w(primitive, ixi^l, ixo^l, w, flag)
912
913 if (any(flag)) then
914 select case (small_values_method)
915 case ("replace")
916 ! any faulty cell is replaced by physical lower limit
917 flagall(ixo^s)=(flag(ixo^s,rho_).or.flag(ixo^s,e_))
918
919 where(flagall(ixo^s))
920 ! D or rho: no difference primitive-conservative
921 w(ixo^s,rho_) = small_density
922 !w(ixO^S,lfac_)= 1.0d0
923 !w(ixO^S,xi_) = small_xi
924 endwhere
925 !do idir = 1, ndir
926 ! where(flagall(ixO^S)) w(ixO^S, mom(idir)) = 0.0d0
927 !end do
928 if(primitive) then
929 where(flagall(ixo^s)) w(ixo^s, p_) = small_pressure
930 else
931 where(flagall(ixo^s)) w(ixo^s, e_) = small_e
932 endif
933
934 case ("average")
935 ! note: in small_values_average we use
936 ! small_values_fix_iw(1:nw) and small_values_daverage
937 ! when fails, may use small_pressure/small_density
938 if(primitive)then
939 ! averaging for all primitive fields (p, lfac*v, tau))
940 call small_values_average(ixi^l, ixo^l, w, x, flag)
941 ! update the auxiliaries from primitives
942 call srhd_get_auxiliary_prim(ixi^l,ixo^l,w)
943 else
944 ! do averaging of density d
945 call small_values_average(ixi^l, ixo^l, w, x, flag, d_)
946 ! do averaging of energy tau
947 call small_values_average(ixi^l, ixo^l, w, x, flag, e_)
948 ! and now hope for the best....
949 endif
950 case default
951 if(.not.primitive) then
952 ! note that we throw error here, which assumes w is primitive
953 write(*,*) "handle_small_values default: note reporting conservatives!"
954 end if
955 call small_values_error(w, x, ixi^l, ixo^l, flag, subname)
956 end select
957 end if
958
959 end subroutine srhd_handle_small_values
960
961 !> calculate effective gamma
962 subroutine srhd_get_geff_eos(w,ixI^L,ixO^L,varconserve,Geff)
963 use mod_global_parameters, only: nw
964 !================== IMPORTANT ==================!
965 !This subroutine is used with conserved variables in w when varconserve=T
966 !This subroutine is used with primitive variables in w when varconserve=F
967 ! both cases assume updated auxiliary variables xi_ en lfac_
968 !===============================================!
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)
973
974 double precision, dimension(ixO^S) :: pth,rho,e_th,e
975
976 if (srhd_eos) then
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_)
980 e_th = pth*inv_gamma_1
981 e = e_th+dsqrt(e_th**2+rho**2)
982 geff(ixo^s) = srhd_gamma-half*gamma_1 * &
983 (one-(rho(ixo^s)/e(ixo^s))**2)
984 else
985 ! primitives available
986 e_th = w(ixo^s,p_)*inv_gamma_1
987 e = e_th+dsqrt(e_th**2+w(ixo^s,rho_)**2)
988 geff(ixo^s) = srhd_gamma-half*gamma_1 * &
989 (one-(w(ixo^s,rho_)/e(ixo^s))**2)
990 end if
991 else
992 geff(ixo^s) = srhd_gamma
993 endif
994
995 end subroutine srhd_get_geff_eos
996
997 !> Compute the small value limits
998 subroutine srhd_get_smallvalues_eos
1000 implicit none
1001 ! local small values
1002 double precision :: lsmalle,lsmallp,lsmallrho
1003
1004 ! the maximal allowed Lorentz factor
1005 lfacmax=one/dsqrt(one-(one-dmaxvel)**2)
1008 if(small_density*small_pressure<=0.0d0)then
1009 call mpistop("must set finite values small-density/pressure for small value treatments")
1010 endif
1011 if(srhd_eos)then
1012 lsmallp=(one+10.d0*small_pressure)*small_pressure
1013 lsmallrho=(one+10.d0*small_density)*small_density
1014 !!Lsmallp=small_pressure
1015 !!Lsmallrho=small_density
1016 lsmalle=lsmallp*inv_gamma_1+&
1017 dsqrt((lsmallp*inv_gamma_1)**2+lsmallrho**2)
1018 small_xi=half*((srhd_gamma+one)*lsmalle-&
1019 gamma_1*lsmallrho*(lsmallrho/lsmalle))
1020 small_e=small_xi-lsmallp-lsmallrho
1021 else
1024 endif
1027
1028 end subroutine srhd_get_smallvalues_eos
1029
1030 !> Compute the enthalpy rho*h from rho and pressure p
1031 subroutine srhd_get_enthalpy_eos(ixO^L,rho,p,rhoh)
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)
1036
1037 double precision, dimension(ixO^S) :: e_th,e
1038 integer :: ix^d
1039
1040 if(srhd_eos) then
1041 e_th = p*inv_gamma_1
1042 e = e_th+dsqrt(e_th**2+rho**2)
1043 ! writing rho/E on purpose, for numerics
1044 rhoh = half*((srhd_gamma+one)*e &
1045 - gamma_1*rho*(rho/e))
1046 else
1047 rhoh = rho+gamma_to_gamma_1*p
1048 end if
1049
1050 if (check_small_values) then
1051 {do ix^db= ixo^lim^db\}
1052 if(rhoh(ix^d)<small_xi) then
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)')
1057 end if
1058 {enddo^d&\}
1059 end if
1060
1061 if (fix_small_values) then
1062 {do ix^db= ixo^lim^db\}
1063 if(rhoh(ix^d)<small_xi) then
1064 rhoh(ix^d)=small_xi
1065 endif
1066 {enddo^d&\}
1067 endif
1068
1069 end subroutine srhd_get_enthalpy_eos
1070
1071 !> Calculate thermal pressure p from density rho and enthalpy rho*h
1072 !> will provide p (and E if srhd_eos)
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)
1079 integer :: ix^d
1080
1081 if(srhd_eos) then
1082 e = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-one)*rho**2)) &
1083 /(srhd_gamma+one)
1084 p(ixo^s) = half*gamma_1* (e-rho*(rho/e))
1085 else
1086 p(ixo^s) = (rhoh-rho)/gamma_to_gamma_1
1087 end if
1088
1089 if (check_small_values) then
1090 {do ix^db= ixo^lim^db\}
1091 if(p(ix^d)<small_pressure) then
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)')
1098 end if
1099 {enddo^d&\}
1100 end if
1101
1102 if (fix_small_values) then
1103 {do ix^db= ixo^lim^db\}
1104 if(p(ix^d)<small_pressure) then
1105 p(ix^d)=small_pressure
1106 if(srhd_eos)e(ix^d)=max(small_e,e(ix^d))
1107 endif
1108 {enddo^d&\}
1109 endif
1110
1111 end subroutine srhd_get_pressure_eos
1112
1113 !> Calculate the square of the thermal sound speed csound2 within ixO^L.
1114 !> available rho - rho*h
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)
1120
1121 double precision :: p(ixi^s)
1122 double precision :: e(ixo^s)
1123 integer :: ix^d
1124
1125 call srhd_get_pressure_eos(ixi^l,ixo^l,rho,rhoh,p,e)
1126 if(srhd_eos) then
1127 csound2(ixo^s)=(p(ixo^s)*((srhd_gamma+one)&
1128 +gamma_1*(rho(ixo^s)/e(ixo^s))**2))&
1129 /(2.0d0*rhoh(ixo^s))
1130 else
1131 csound2(ixo^s)=srhd_gamma*p(ixo^s)/rhoh(ixo^s)
1132 end if
1133
1134 if (check_small_values) then
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)')
1142 end if
1143 {enddo^d&\}
1144 end if
1145
1146 if (fix_small_values) then
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
1150 endif
1151 if(csound2(ix^d)<=0.0d0) then
1153 endif
1154 {enddo^d&\}
1155 endif
1156
1157 end subroutine srhd_get_csound2_eos
1158
1159 !> Calculate the square of the thermal sound speed csound2 within ixO^L.
1160 !> available rho - rho*h - p
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)
1166
1167 double precision :: e(ixo^s)
1168 integer :: ix^d
1169
1170 if(srhd_eos) then
1171 e = (rhoh+dsqrt(rhoh**2+(srhd_gamma**2-one)&
1172 *rho**2))/(srhd_gamma+one)
1173 csound2(ixo^s)=(p*((srhd_gamma+one)&
1174 +gamma_1*(rho/e)**2))&
1175 /(2.0d0*rhoh)
1176 else
1177 csound2(ixo^s)=srhd_gamma*p/rhoh
1178 end if
1179
1180 if (check_small_values) then
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)')
1188 end if
1189 {enddo^d&\}
1190 end if
1191
1192 if (fix_small_values) then
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
1196 endif
1197 if(csound2(ix^d)<=0.0d0) then
1199 endif
1200 {enddo^d&\}
1201 endif
1202
1203 end subroutine srhd_get_csound2_prim_eos
1204
1205 !> con2prim: (D,S**2,tau) --> compute auxiliaries lfac and xi
1206 subroutine con2prim_eos(lfac,xi,myd,myssqr,mytau,ierror)
1208
1209 double precision, intent(in) :: myd, myssqr, mytau
1210 double precision, intent(inout) :: lfac, xi
1211 integer, intent(inout) :: ierror
1212
1213 ! .. local ..
1214 double precision:: f,df,lfacl
1215 !------------------------------------------------------------------
1216
1217 ! Save the input-state in mod_con2prim_vars
1218 d = myd; ssqr = myssqr; tau = mytau;
1219
1220 ierror=0
1221
1222 ! Check if guess is close enough: gives f,df,lfacl
1223 if(xi>smallxi)then
1224 call funcd_eos(xi,f,df,lfacl,d,ssqr,tau,ierror)
1225 if (ierror == 0 .and. dabs(f/df)<absaccnr) then
1226 xi = xi - f/df
1227 lfac = lfacl
1228 return
1229 else
1230 ierror = 0
1231 end if
1232 else
1233 write(*,*)'entering con2prim_eos with xi=',xi
1234 end if
1235
1236 ! ierror=1 : must have D>=minrho, tau>=smalltau
1237 !if(d<minrho .or. tau<smalltau)then
1238 ! ierror=1
1239 ! return
1240 !endif
1241
1242 call con2primhydro_eos(lfac,xi,d,ssqr,tau,ierror)
1243
1244 end subroutine con2prim_eos
1245
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
1250
1251 ! .. local ..
1252 double precision :: dlfac
1253 double precision :: vsqr,p,dpdxi
1254 !-----------------------------------------------------------------
1255
1256 vsqr = ssqr/xi**2
1257
1258 if (vsqr<one) then
1259 mylfac = one/dsqrt(one-vsqr)
1260 dlfac = -mylfac**3*ssqr/(xi**3)
1261 !===== Pressure, calculate using EOS =====!
1262 call funcpressure_eos(xi,mylfac,d,dlfac,p,dpdxi)
1263 !=========================================!
1264 f = xi-tau-d-p
1265 df = one-dpdxi
1266 else
1267 ! print *,'Erroneous input to funcd since vsqr=',vsqr,' >=1'
1268 ! print *,'input values d, ssqr, tau:',d,ssqr,tau
1269 ierror =6
1270 return
1271 end if
1272
1273 end subroutine funcd_eos
1274
1275 !> SRHD iteration solves for p via NR, and then gives xi as output
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
1280
1281 ! .. local ..
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
1290 !---------------------------------------------------------------------
1291
1292 ierror=0
1293 ! ierror=0 : ok
1294 ! we already checked D>=minrho, tau>=smalltau (ierror=1)
1295 !
1296 ! ierror<>0
1297 !
1298 ! ierror=2 : maxitnr reached without convergence
1299 ! ierror=3 : final pressure value < smallp or xi<smallxi during iteration
1300 ! ierror=4 : final v^2=1 hence problem as lfac=1/0
1301 ! ierror=5 : nonmonotonic function f (as df=0)
1302
1303 ! left and right brackets for p-range
1304 pmin=dsqrt(sqrs)/(one-dmaxvel)-tau-d
1305 plabs=max(minp,pmin)
1306 prabs=1.0d99
1307 ! start value from input
1308 pcurrent=plabs
1309
1310 er1=one
1311 pprev=pcurrent
1312
1313 ! Fudge Parameters
1314 oldff1=1.0d7 ! High number
1315 oldff2=1.0d9 ! High number bigger then oldff1
1316 n2it = 0
1317 nit = 0
1318
1319 loopnr: do ni=1,maxitnr
1320 nit = nit + 1
1321 !=== Relax NR iteration accuracy=======!
1322 if(nit>maxitnr/4)then
1323 ! mix pressure value for convergence
1324 pcurrent=half*(pcurrent+pprev)
1325 ! relax accuracy requirement
1326 er1=10.0d0*er1
1327 nit = nit - maxitnr/10
1328 endif
1329 !=======================================!
1330
1331 niiter=ni
1332 xicurrent=tau+d+pcurrent
1333
1334 if(xicurrent<smallxi) then
1335 ! print *,'stop: too small xi iterate:',xicurrent
1336 ! print *,'for pressure iterate p',pcurrent
1337 ! print *,'pressure bracket pLabs pRabs',pLabs,pRabs
1338 ! print *,'iteration number:',ni
1339 ! print *,'values for d,s,tau,s2:',d,sqrs,tau,sqrs
1340 ierror=3
1341 return
1342 endif
1343
1344 v2=sqrs/xicurrent**2
1345 lfac2inv=one - v2
1346 if(lfac2inv>zero) then
1347 lfac=one/dsqrt(lfac2inv)
1348 else
1349 ! print *,'stop: negative or zero factor 1-v2:',lfac2inv
1350 ! print *,'for pressure iterate p',pcurrent
1351 ! print *,'absolute pressure bracket pLabs pRabs',pLabs,pRabs
1352 ! print *,'iteration number:',ni
1353 ! print *,'values for d,s,tau,s2:',d,sqrs,tau,sqrs
1354 ! print *,'values for v2,xi:',v2,xicurrent
1355 ierror=4
1356 return
1357 endif
1358
1359 s2overcubeg2rh=sqrs/(xicurrent**3)
1360 !== calculation done using the EOS ==!
1361 call funcenthalpy_eos(pcurrent,lfac2inv,d,sqrs,xicurrent,&
1362 s2overcubeg2rh,h,dhdp)
1363 !=======================================!
1364 ff=-xicurrent*lfac2inv + h
1365 df=- two*sqrs/xicurrent**2 + dhdp - lfac2inv
1366
1367 if (ff*df==zero) then
1368 if (ff==zero) then
1369 exit ! zero found
1370 else
1371 ! print *,'stop: df becomes zero, non-monotonic f(p)!'
1372 ierror=5
1373 return
1374 endif
1375 else
1376 pnew=pcurrent-ff/df
1377 if (ff*df>zero) then
1378 ! pressure iterate has decreased
1379 ! restrict to left
1380 pnew=max(pnew,plabs)
1381 else ! ff*df<0
1382 ! pressure iterate has increased
1383 ! restrict to right
1384 pnew=min(pnew,prabs)
1385 endif
1386 endif
1387
1388 !===============================================!
1389 dp=pcurrent-pnew
1390 er=two*dabs(dp)/(pnew+pcurrent)
1391 if(((er<tolernr*er1).or.(dabs(dp)<absaccnr))) exit loopnr
1392 !===============================================!
1393
1394 ! For very small values of pressure, NR algorithm is not efficient to
1395 ! find root, use Euler algorithm to find precise value of pressure
1396 if((dabs(oldff2-ff) < 1.0d-8 .or. niiter >= maxitnr-maxitnr/20).and.&
1397 ff * oldff1 < zero .and. dabs(ff)>absaccnr)then
1398
1399 n2it=n2it+1
1400 if(n2it<=3) pcurrent=half*(pnew+pcurrent)
1401 if(n2it>3)then
1402 pright =pcurrent
1403 pleft=pprev
1404 pcurrent=half*(pleft+pright)
1405 dicho: do ni3=1,maxitnr
1406 !===================!
1407 xicurrent=tau+d+pcurrent
1408 v2=sqrs/xicurrent**2
1409 lfac2inv=one - v2
1410
1411 if(lfac2inv>zero)then
1412 lfac=one/dsqrt(lfac2inv)
1413 else
1414 ierror=4
1415 return
1416 endif
1417 !===================!
1418
1419 !== calculation done using the EOS ==!
1420 call bisection_enthalpy_eos(pnew,lfac2inv,d,xicurrent,h)
1421 nff=-xicurrent*lfac2inv + h
1422 !=======================================!
1423 !==== Iterate ====!
1424 if(ff * nff < zero)then
1425 pleft=pcurrent
1426 else
1427 pright=pcurrent
1428 endif
1429
1430 pcurrent=half*(pleft+pright)
1431 !==================!
1432
1433 !=== The iteration converged ===!
1434 if(2.0d0*dabs(pleft-pright)/(pleft+pright)< absaccnr &
1435 .or. dabs(ff)<absaccnr)then
1436 pnew=pcurrent
1437 exit loopnr
1438 endif
1439 !==============================!
1440
1441 !==============================!
1442
1443 !=== conserve the last value of Nff ===!
1444 ff=nff
1445 !======================================!
1446 enddo dicho
1447 endif
1448
1449 else
1450 !====== There is no problems, continue the NR iteration ======!
1451 pprev=pcurrent
1452 pcurrent=pnew
1453 !=============================================================!
1454 endif
1455
1456
1457 !=== keep the values of the 2 last ff ===!
1458 oldff2=oldff1
1459 oldff1=ff
1460 !========================================!
1461 enddo loopnr
1462
1463 if(niiter==maxitnr)then
1464 ierror=2
1465 return
1466 endif
1467
1468 if(pcurrent<minp) then
1469 ierror=3
1470 return
1471 endif
1472
1473 !------------------------------!
1474 xi=tau+d+pcurrent
1475 v2=sqrs/xicurrent**2
1476 lfac2inv=one - v2
1477 if(lfac2inv>zero) then
1478 lfac=one/dsqrt(lfac2inv)
1479 else
1480 ierror=4
1481 return
1482 endif
1483
1484 end subroutine con2primhydro_eos
1485
1486 !> pointwise evaluations used in con2prim
1487 !> compute pointwise value for pressure p and dpdxi
1488 subroutine funcpressure_eos(xicurrent,lfac,d,dlfacdxi,p,dpdxi)
1489
1490 double precision, intent(in) :: xicurrent,lfac,d,dlfacdxi
1491 double precision, intent(out) :: p,dpdxi
1492 ! .. local ..
1493 double precision :: rho,h,e,dhdxi,rhotoe
1494 double precision :: dpdchi,dedxi
1495
1496 ! rhoh here called h
1497 h=xicurrent/(lfac**2)
1498 rho=d/lfac
1499 e = (h+dsqrt(h**2+(srhd_gamma**2-one)*rho**2)) &
1500 /(srhd_gamma+one)
1501 ! output pressure
1502 rhotoe = rho/e
1503 p = half*gamma_1*(e-rho*rhotoe)
1504
1505 dhdxi = one/(lfac**2)-2.0d0*xicurrent/(lfac**2)*dlfacdxi/lfac
1506
1507 dedxi=(dhdxi+(h*dhdxi-(srhd_gamma**2-one)*rho**2*dlfacdxi/lfac)&
1508 /dsqrt(h**2+(srhd_gamma**2-one)*rho**2))&
1509 /(srhd_gamma+one)
1510
1511 ! output pressure derivative to xi
1512 dpdxi=half*gamma_1*(2.0d0*rho*rhotoe*dlfacdxi/lfac+&
1513 (one+rhotoe**2)*dedxi)
1514
1515 end subroutine funcpressure_eos
1516
1517 !> pointwise evaluations used in con2prim
1518 !> returns enthalpy rho*h (h) and derivative d(rho*h)/dp (dhdp)
1519 subroutine funcenthalpy_eos(pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p,h,dhdp)
1520
1521 double precision, intent(in) :: pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p
1522 double precision, intent(out):: h,dhdp
1523
1524 ! local
1525 double precision:: rho,e_th,e,de_thdp,dedp
1526
1527 rho=d*dsqrt(lfac2inv)
1528 e_th = pcurrent*inv_gamma_1
1529 e = (e_th + dsqrt(e_th**2+rho**2))
1530 !== Enthalpy ==!
1531 h = half*((srhd_gamma+one)*e-gamma_1*rho*(rho/e))
1532 !=== Derivative of thermal energy ===!
1533 de_thdp = one*inv_gamma_1
1534 !=== Derivative of internal energy ===!
1535 dedp = de_thdp * (one+e_th/dsqrt(e_th**2+rho**2))&
1536 + d**2*dv2d2p/dsqrt(e_th**2+rho**2)
1537 !====== Derivative of Enthalpy ======!
1538 dhdp = half*((srhd_gamma+one)*dedp + &
1539 gamma_1*(rho*(rho/e))*(-2.0d0*dv2d2p/lfac2inv+dedp/e))
1540 end subroutine funcenthalpy_eos
1541
1542 !> pointwise evaluations used in con2prim
1543 !> returns enthalpy rho*h (h)
1544 subroutine bisection_enthalpy_eos(pcurrent,lfac2inv,d,xicurrent,h)
1545
1546 double precision, intent(in) :: pcurrent,lfac2inv,d,xicurrent
1547 double precision, intent(out):: h
1548
1549 ! local
1550 double precision:: rho,e_th,e
1551
1552 rho=d*dsqrt(lfac2inv)
1553 e_th = pcurrent*inv_gamma_1
1554 e = (e_th + dsqrt(e_th**2+rho**2))
1555 !== Enthalpy ==!
1556 h = half*((srhd_gamma+one)*e-gamma_1*rho*(rho/e))
1557
1558 return
1559 end subroutine bisection_enthalpy_eos
1560
1561 !> con2prim: (D,S**2,tau) --> compute auxiliaries lfac and xi
1562 subroutine con2prim(lfac,xi,myd,myssqr,mytau,ierror)
1564
1565 double precision, intent(in) :: myd, myssqr, mytau
1566 double precision, intent(inout) :: lfac, xi
1567 integer, intent(inout) :: ierror
1568
1569 ! .. local ..
1570 double precision:: f,df,lfacl
1571 !------------------------------------------------------------------
1572
1573 ! Save the input-state in mod_con2prim_vars
1574 d = myd; ssqr = myssqr; tau = mytau;
1575
1576 ierror=0
1577
1578 ! Check if guess is close enough: gives f,df,lfacl
1579 if(xi>smallxi)then
1580 call funcd(xi,f,df,lfacl,d,ssqr,tau,ierror)
1581 if (ierror == 0 .and. dabs(f/df)<absaccnr) then
1582 xi = xi - f/df
1583 lfac = lfacl
1584 return
1585 else
1586 ierror = 0
1587 end if
1588 else
1589 write(*,*) 'entering con2prim with xi=',xi
1590 end if
1591
1592 ! ierror=1 : must have D>=minrho, tau>=smalltau
1593 !if(d<minrho .or. tau<smalltau)then
1594 ! ierror=1
1595 ! return
1596 !endif
1597
1598 call con2primhydro(lfac,xi,d,ssqr,tau,ierror)
1599
1600 end subroutine con2prim
1601
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
1606
1607 ! .. local ..
1608 double precision :: dlfac
1609 double precision :: vsqr,p,dpdxi
1610 !-----------------------------------------------------------------
1611
1612 vsqr = ssqr/xi**2
1613
1614 if (vsqr<one) then
1615 mylfac = one/dsqrt(one-vsqr)
1616 dlfac = -mylfac**3*ssqr/(xi**3)
1617 !===== Pressure, calculate using EOS =====!
1618 call funcpressure(xi,mylfac,d,dlfac,p,dpdxi)
1619 !=========================================!
1620 f = xi-tau-d-p
1621 df = one-dpdxi
1622 else
1623 ! print *,'Erroneous input to funcd since vsqr=',vsqr,' >=1'
1624 ! print *,'input values d, ssqr, tau:',d,ssqr,tau
1625 ierror =6
1626 return
1627 end if
1628
1629 end subroutine funcd
1630
1631 !> SRHD iteration solves for p via NR, and then gives xi as output
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
1636
1637 ! .. local ..
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
1646 !---------------------------------------------------------------------
1647
1648 ierror=0
1649 ! ierror=0 : ok
1650 ! we already checked D>=minrho, tau>=smalltau (ierror=1)
1651 !
1652 ! ierror<>0
1653 !
1654 ! ierror=2 : maxitnr reached without convergence
1655 ! ierror=3 : final pressure value < smallp or xi<smallxi during iteration
1656 ! ierror=4 : final v^2=1 hence problem as lfac=1/0
1657 ! ierror=5 : nonmonotonic function f (as df=0)
1658
1659 ! left and right brackets for p-range
1660 pmin=dsqrt(sqrs)/(one-dmaxvel)-tau-d
1661 plabs=max(minp,pmin)
1662 prabs=1.0d99
1663 ! start value from input
1664 pcurrent=plabs
1665
1666 er1=one
1667 pprev=pcurrent
1668
1669 ! Fudge Parameters
1670 oldff1=1.0d7 ! High number
1671 oldff2=1.0d9 ! High number bigger then oldff1
1672 n2it = 0
1673 nit = 0
1674
1675 loopnr: do ni=1,maxitnr
1676 nit = nit + 1
1677 !=== Relax NR iteration accuracy=======!
1678 if(nit>maxitnr/4)then
1679 ! mix pressure value for convergence
1680 pcurrent=half*(pcurrent+pprev)
1681 ! relax accuracy requirement
1682 er1=10.0d0*er1
1683 nit = nit - maxitnr/10
1684 endif
1685 !=======================================!
1686
1687 niiter=ni
1688 xicurrent=tau+d+pcurrent
1689
1690 if(xicurrent<smallxi) then
1691 ! print *,'stop: too small xi iterate:',xicurrent
1692 ! print *,'for pressure iterate p',pcurrent
1693 ! print *,'pressure bracket pLabs pRabs',pLabs,pRabs
1694 ! print *,'iteration number:',ni
1695 ! print *,'values for d,s,tau,s2:',d,sqrs,tau,sqrs
1696 ierror=3
1697 return
1698 endif
1699
1700 v2=sqrs/xicurrent**2
1701 lfac2inv=one - v2
1702 if(lfac2inv>zero) then
1703 lfac=one/dsqrt(lfac2inv)
1704 else
1705 ! print *,'stop: negative or zero factor 1-v2:',lfac2inv
1706 ! print *,'for pressure iterate p',pcurrent
1707 ! print *,'absolute pressure bracket pLabs pRabs',pLabs,pRabs
1708 ! print *,'iteration number:',ni
1709 ! print *,'values for d,s,tau,s2:',d,sqrs,tau,sqrs
1710 ! print *,'values for v2,xi:',v2,xicurrent
1711 ierror=4
1712 return
1713 endif
1714
1715 s2overcubeg2rh=sqrs/(xicurrent**3)
1716 !== calculation done using the EOS ==!
1717 call funcenthalpy(pcurrent,lfac2inv,d,sqrs,xicurrent,&
1718 s2overcubeg2rh,h,dhdp)
1719 !=======================================!
1720 ff=-xicurrent*lfac2inv + h
1721 df=- two*sqrs/xicurrent**2 + dhdp - lfac2inv
1722
1723 if (ff*df==zero) then
1724 if (ff==zero) then
1725 exit ! zero found
1726 else
1727 ! print *,'stop: df becomes zero, non-monotonic f(p)!'
1728 ierror=5
1729 return
1730 endif
1731 else
1732 pnew=pcurrent-ff/df
1733 if (ff*df>zero) then
1734 ! pressure iterate has decreased
1735 ! restrict to left
1736 pnew=max(pnew,plabs)
1737 else ! ff*df<0
1738 ! pressure iterate has increased
1739 ! restrict to right
1740 pnew=min(pnew,prabs)
1741 endif
1742 endif
1743
1744 !===============================================!
1745 dp=pcurrent-pnew
1746 er=two*dabs(dp)/(pnew+pcurrent)
1747 if(((er<tolernr*er1).or.(dabs(dp)<absaccnr))) exit loopnr
1748 !===============================================!
1749
1750 ! For very small values of pressure, NR algorithm is not efficient to
1751 ! find root, use Euler algorithm to find precise value of pressure
1752 if((dabs(oldff2-ff) < 1.0d-8 .or. niiter >= maxitnr-maxitnr/20).and.&
1753 ff * oldff1 < zero .and. dabs(ff)>absaccnr)then
1754
1755 n2it=n2it+1
1756 if(n2it<=3) pcurrent=half*(pnew+pcurrent)
1757 if(n2it>3)then
1758 pright =pcurrent
1759 pleft=pprev
1760 pcurrent=half*(pleft+pright)
1761 dicho: do ni3=1,maxitnr
1762 !===================!
1763 xicurrent=tau+d+pcurrent
1764 v2=sqrs/xicurrent**2
1765 lfac2inv=one - v2
1766
1767 if(lfac2inv>zero)then
1768 lfac=one/dsqrt(lfac2inv)
1769 else
1770 ierror=4
1771 return
1772 endif
1773 !===================!
1774
1775 !== calculation done using the EOS ==!
1776 call bisection_enthalpy(pnew,lfac2inv,d,xicurrent,h)
1777 nff=-xicurrent*lfac2inv + h
1778 !=======================================!
1779 !==== Iterate ====!
1780 if(ff * nff < zero)then
1781 pleft=pcurrent
1782 else
1783 pright=pcurrent
1784 endif
1785
1786 pcurrent=half*(pleft+pright)
1787 !==================!
1788
1789 !=== The iteration converged ===!
1790 if(2.0d0*dabs(pleft-pright)/(pleft+pright)< absaccnr &
1791 .or. dabs(ff)<absaccnr)then
1792 pnew=pcurrent
1793 exit loopnr
1794 endif
1795 !==============================!
1796
1797 !==============================!
1798
1799 !=== conserve the last value of Nff ===!
1800 ff=nff
1801 !======================================!
1802 enddo dicho
1803 endif
1804
1805 else
1806 !====== There is no problems, continue the NR iteration ======!
1807 pprev=pcurrent
1808 pcurrent=pnew
1809 !=============================================================!
1810 endif
1811
1812
1813 !=== keep the values of the 2 last ff ===!
1814 oldff2=oldff1
1815 oldff1=ff
1816 !========================================!
1817 enddo loopnr
1818
1819 if(niiter==maxitnr)then
1820 ierror=2
1821 return
1822 endif
1823
1824 if(pcurrent<minp) then
1825 ierror=3
1826 return
1827 endif
1828
1829 !------------------------------!
1830 xi=tau+d+pcurrent
1831 v2=sqrs/xicurrent**2
1832 lfac2inv=one - v2
1833 if(lfac2inv>zero) then
1834 lfac=one/dsqrt(lfac2inv)
1835 else
1836 ierror=4
1837 return
1838 endif
1839
1840 end subroutine con2primhydro
1841
1842 !> pointwise evaluations used in con2prim
1843 !> compute pointwise value for pressure p and dpdxi
1844 subroutine funcpressure(xicurrent,lfac,d,dlfacdxi,p,dpdxi)
1845
1846 double precision, intent(in) :: xicurrent,lfac,d,dlfacdxi
1847 double precision, intent(out) :: p,dpdxi
1848 ! .. local ..
1849 double precision :: rho,h,e,dhdxi,rhotoe
1850 double precision :: dpdchi,dedxi
1851
1852 ! rhoh here called h
1853 h=xicurrent/(lfac**2)
1854 rho=d/lfac
1855 ! output pressure
1856 p = (h - rho)/gamma_to_gamma_1
1857 dpdchi = one/gamma_to_gamma_1
1858 dpdxi = dpdchi * one/lfac**2
1859 ! zero case dlfacdxi implies zero velocity (ssqr=0)
1860 if (dlfacdxi /= 0.0d0) &
1861 dpdxi = dpdxi + dpdchi * ((d*lfac-2.0d0*xicurrent)/lfac**3) * dlfacdxi
1862
1863 end subroutine funcpressure
1864
1865 !> pointwise evaluations used in con2prim
1866 !> returns enthalpy rho*h (h) and derivative d(rho*h)/dp (dhdp)
1867 subroutine funcenthalpy(pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p,h,dhdp)
1868
1869 double precision, intent(in) :: pcurrent,lfac2inv,d,sqrs,xicurrent,dv2d2p
1870 double precision, intent(out):: h,dhdp
1871
1872 ! local
1873 double precision:: rho,e_th,e,de_thdp,dedp
1874
1875 rho=d*dsqrt(lfac2inv)
1876 h = rho + gamma_to_gamma_1 * pcurrent
1877 dhdp = gamma_to_gamma_1 + d/dsqrt(lfac2inv)*sqrs/xicurrent**3
1878 end subroutine funcenthalpy
1879
1880 !> pointwise evaluations used in con2prim
1881 !> returns enthalpy rho*h (h)
1882 subroutine bisection_enthalpy(pcurrent,lfac2inv,d,xicurrent,h)
1883
1884 double precision, intent(in) :: pcurrent,lfac2inv,d,xicurrent
1885 double precision, intent(out):: h
1886
1887 ! local
1888 double precision:: rho,e_th,e
1889
1890 rho=d*dsqrt(lfac2inv)
1891 h = rho + gamma_to_gamma_1 * pcurrent
1892
1893 return
1894 end subroutine bisection_enthalpy
1895
1896end module mod_srhd_phys
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.
double precision ssqr
double precision d
double precision tau
Module for physical and numeric constants.
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
integer coordinate
Definition mod_geometry.t:7
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 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
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...
Definition mod_physics.t:4
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.