19 double precision,
public ::
fld_mu = 0.6d0
49 integer,
allocatable,
public ::
i_opf(:)
51 double precision,
private,
protected :: fld_gamma
70 character(len=*),
intent(in) :: files(:)
80 open(
unitpar, file=trim(files(n)), status=
"old")
81 read(
unitpar, fld_list,
end=111)
93 subroutine fld_init(He_abundance, radiation_diffusion, energy_interact, r_gamma)
100 double precision,
intent(in) :: he_abundance, r_gamma
101 logical,
intent(in) :: radiation_diffusion, energy_interact
102 double precision :: sigma_thomson
104 character(len=1) :: ind_1
105 character(len=1) :: ind_2
106 character(len=2) :: cmp_f
107 character(len=5) :: cmp_e
115 write(ind_1,
'(I1)') idir
120 if(radiation_diffusion .or. energy_interact)
then
126 mg%operator_type = mg_vhelmholtz
131 fld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
137 sigma_thomson = 6.6524585
d-25
138 fld_kappa0 = sigma_thomson/const_mp * (1.+2.*he_abundance)/(1.+4.*he_abundance)
145 energy,qsourcesplit,active)
150 integer,
intent(in) :: ixi^
l, ixo^
l
151 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
152 double precision,
intent(in) :: wct(ixi^s,1:nw)
153 double precision,
intent(inout) :: w(ixi^s,1:nw)
154 logical,
intent(in) :: energy,qsourcesplit
155 logical,
intent(inout) :: active
156 integer :: idir, jdir
157 double precision :: radiation_forcect(ixo^s,1:
ndim)
158 double precision :: rad_fluxct(ixo^s,1:
ndim)
159 double precision :: grad_e(ixi^s)
160 double precision :: lambda(ixo^s),fld_r(ixo^s)
169 radiation_forcect(ixo^s,idir) = -lambda(ixo^s)*grad_e(ixo^s)
171 w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir)) &
172 + qdt*radiation_forcect(ixo^s,idir)
174 w(ixo^s,iw_e) = w(ixo^s,iw_e) &
175 + qdt*wct(ixo^s,iw_mom(idir))/wct(ixo^s,iw_rho)*radiation_forcect(ixo^s,idir)
184 integer,
intent(in) :: ixi^
l, ixo^
l
185 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
186 double precision,
intent(inout) :: dtnew
187 double precision :: radiation_force(ixo^s,1:
ndim)
188 double precision :: dxinv(1:
ndim), max_grad
189 double precision :: lambda(ixo^s),fld_r(ixo^s),grad_e(ixi^s)
192 ^
d&dxinv(^
d)=one/
dx^
d;
196 radiation_force(ixo^s,idir) = -lambda(ixo^s)*grad_e(ixo^s)
197 max_grad = maxval(abs(radiation_force(ixo^s,idir)))
198 max_grad = max(max_grad, epsilon(1.0d0))
199 dtnew = min(dtnew,
courantpar / sqrt(max_grad * dxinv(idir)))
211 integer,
intent(in) :: ixi^
l, ixo^
l
212 double precision,
intent(in) :: w(ixi^s, 1:nw)
213 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
214 double precision,
intent(out) :: fld_kappa(ixo^s)
215 double precision :: temp(ixi^s), pth(ixi^s), a2(ixo^s)
216 double precision :: rho0,temp0,n,sigma_b
217 double precision :: akram, bkram
218 double precision :: vth(ixo^s), gradv(ixi^s), eta(ixo^s), t(ixo^s)
219 integer :: i,j,ix^
d, idir
236 temp(ixo^s)=temp(ixo^s)/w(ixo^s,iw_rho)
244 akram = 13.1351597305
245 bkram = -4.5182188206
247 * (1.d0+10.d0**akram*w(ixo^s,iw_rho)*
unit_density*(a2(ixo^s)/1.d12)**bkram)
248 {
do ix^
d = ixomin^
d,ixomax^
d\ }
254 {
do ix^
d=ixomin^
d,ixomax^
d\ }
262 call mpistop(
"special opacity not defined")
266 call mpistop(
"Doesn't know opacity law")
278 integer,
intent(in) :: ixi^
l,ixo^
l,nth
279 double precision,
intent(in) :: w(ixi^s,1:nw)
280 double precision,
intent(in) :: x(ixi^s,1:
ndim)
281 double precision,
intent(out) :: fld_r(ixo^s),fld_lambda(ixo^s)
282 integer :: idir, ix^
d
284 double precision :: kappa(ixo^s)
285 double precision :: normgrad2(ixo^s)
286 double precision :: grad_r_e(ixi^s)
287 double precision :: rad_e(ixi^s)
288 double precision :: tmp_l(ixi^s), filtered_l(ixi^s)
293 fld_lambda(ixo^s) = 1.d0/3.d0
299 normgrad2(ixo^s) = zero
300 rad_e(ixi^s) = w(ixi^s, iw_r_e)
302 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
303 normgrad2(ixo^s) = normgrad2(ixo^s)+grad_r_e(ixo^s)**2
306 fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
308 fld_lambda(ixo^s) = one/fld_r(ixo^s)
312 normgrad2(ixo^s) = zero
313 rad_e(ixi^s) = w(ixi^s, iw_r_e)
315 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
316 normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_e(ixo^s)**2
319 fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
322 fld_lambda(ixo^s) = (2.d0+fld_r(ixo^s))/(6.d0+3*fld_r(ixo^s)+fld_r(ixo^s)**2.d0)
326 normgrad2(ixo^s) = zero
327 rad_e(ixi^s) = w(ixi^s, iw_r_e)
329 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
330 normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_e(ixo^s)**2
333 fld_r(ixo^s) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
336 {
do ix^
d = ixomin^
d,ixomax^
d\ }
337 if(fld_r(ix^
d) .lt. 3.d0/2.d0)
then
338 fld_lambda(ix^
d) = 2.d0/(3.d0+dsqrt(9.d0+12.d0*fld_r(ix^
d)**2.d0))
340 fld_lambda(ix^
d) = 1.d0/(1.d0+fld_r(ix^
d)+dsqrt(1.d0+2.d0*fld_r(ix^
d)))
345 call mpistop(
"special fluxlimiter not defined")
349 call mpistop(
'Fluxlimiter unknown')
353 if(
size_l_filter .lt. 1)
call mpistop(
"D filter of size < 1 makes no sense")
355 tmp_l(ixo^s) = fld_lambda(ixo^s)
356 filtered_l(ixo^s) = zero
360 filtered_l(ix^
d) = filtered_l(ix^
d) &
361 + tmp_l(ix^
d+filter*
kr(idir,^
d)) &
362 + tmp_l(ix^
d-filter*
kr(idir,^
d))
369 fld_lambda(ixo^s) = tmp_l(ixo^s)
378 integer,
intent(in) :: ixi^
l, ixo^
l, nth
379 double precision,
intent(in) :: w(ixi^s, 1:nw)
380 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
381 double precision,
intent(out) :: rad_flux(ixo^s, 1:
ndim)
382 double precision :: cc
383 double precision :: grad_r_e(ixi^s)
384 double precision :: rad_e(ixi^s)
385 double precision :: kappa(ixo^s), lambda(ixo^s), fld_r(ixo^s)
386 integer :: ix^
d, idir
389 rad_e(ixi^s) = w(ixi^s, iw_r_e)
395 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
396 rad_flux(ixo^s, idir) = -cc*lambda(ixo^s)/(kappa(ixo^s)*w(ixo^s,iw_rho))*grad_r_e(ixo^s)
405 integer,
intent(in) :: ixI^L, ixO^L, nth
406 double precision,
intent(in) :: w(ixI^S, 1:nw)
407 double precision,
intent(in) :: x(ixI^S, 1:ndim)
408 double precision,
intent(out) :: eddington_tensor(ixI^S,1:ndim,1:ndim)
409 double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
410 double precision :: normgrad2(ixO^S), f(ixO^S)
411 double precision :: grad_r_e(ixI^S,1:ndim), rad_e(ixI^S)
412 double precision :: lambda(ixO^S), fld_R(ixO^S)
413 integer :: i,j, idir,jdir
417 normgrad2(ixo^s) = zero
418 rad_e(ixi^s) = w(ixi^s, iw_r_e)
420 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e(ixi^s,idir),nth)
421 normgrad2(ixo^s)=normgrad2(ixo^s)+grad_r_e(ixo^s,idir)**2
426 f(ixo^s) = lambda(ixo^s)+lambda(ixo^s)**two*fld_r(ixo^s)**two
427 f(ixo^s) = half*(one-f(ixo^s))+half*(3.d0*f(ixo^s)-one)
429 eddington_tensor(ixo^s,1,1) = f(ixo^s)
433 eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s))
437 if(idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
438 tnsr2(ixo^s,idir,jdir) = half*(3.d0*f(ixo^s)-one) &
439 * grad_r_e(ixo^s,idir)*grad_r_e(ixo^s,jdir)/max(normgrad2(ixo^s),1.e-8)
444 where((tnsr2(ixo^s,idir,jdir) .eq. tnsr2(ixo^s,idir,jdir)) &
445 .and. (normgrad2(ixo^s) .gt. smalldouble))
446 eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir)+tnsr2(ixo^s,idir,jdir)
457 integer,
intent(in) :: ixi^
l, ixo^
l, nth
458 double precision,
intent(in) :: w(ixi^s, 1:nw)
459 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
460 double precision,
intent(out):: rad_pressure(ixo^s,1:
ndim,1:
ndim)
462 double precision :: eddington_tensor(ixi^s,1:
ndim,1:
ndim)
467 rad_pressure(ixo^s,i,j) = eddington_tensor(ixo^s,i,j)*w(ixo^s,iw_r_e)
474 type(state),
target :: psa(max_blocks)
475 type(state),
target :: psb(max_blocks)
476 double precision,
intent(in) :: qdt
477 double precision,
intent(in) :: qtC
478 double precision,
intent(in) :: dtfactor
479 integer :: ixO^L,iigrid,igrid
483 ixo^l=ixg^
ll^lsubnghostcells;
485 do iigrid=1,igridstail; igrid=igrids(iigrid);
497 integer,
intent(in) :: ixI^L,ixO^L
498 double precision,
intent(in) :: x(ixI^S,1:ndim)
499 double precision,
intent(in) :: dtfactor, qdt
500 double precision,
intent(inout) :: w(ixI^S,1:nw)
501 double precision,
dimension(ixO^S) :: a1,a2,a3,c0,c1,kappa
502 double precision,
dimension(ixI^S) :: e_gas,E_rad,vel,grad_v,nabla_vP
503 double precision,
dimension(ixI^S,1:ndim,1:ndim) :: div_v,edd
504 double precision :: sigma_b,cc
505 integer :: i,j,idir,jdir,ix^D
512 vel(ixi^s) = w(ixi^s,iw_mom(jdir))/w(ixi^s,iw_rho)
513 call gradient(vel,ixi^l,ixo^l,idir,grad_v)
514 div_v(ixo^s,idir,jdir) = grad_v(ixo^s)
521 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)
525 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
526 + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
527 + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
528 + div_v(ixo^s,2,2)*edd(ixo^s,2,2)
531 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
532 + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
533 + div_v(ixo^s,1,3)*edd(ixo^s,1,3) &
534 + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
535 + div_v(ixo^s,2,2)*edd(ixo^s,2,2) &
536 + div_v(ixo^s,2,3)*edd(ixo^s,2,3) &
537 + div_v(ixo^s,3,1)*edd(ixo^s,3,1) &
538 + div_v(ixo^s,3,2)*edd(ixo^s,3,2) &
539 + div_v(ixo^s,3,3)*edd(ixo^s,3,3)
545 e_gas(ixo^s) = w(ixo^s,iw_e)-half*sum(w(ixo^s,iw_mom(:))**2,dim=ndim+1)/w(ixo^s,iw_rho)
546 if(
allocated(iw_mag))
then
547 e_gas(ixo^s) = e_gas(ixo^s)-half*sum(w(ixo^s,iw_mag(:))**2,dim=ndim+1)
549 {
do ix^d = ixomin^d,ixomax^d\ }
552 e_rad(ixo^s) = w(ixo^s,iw_r_e)
559 a1(ixo^s) = 4.d0*kappa(ixo^s)*sigma_b*(fld_gamma-one)**4/w(ixo^s,iw_rho)**3*dtfactor*qdt
560 a2(ixo^s) = cc*kappa(ixo^s)*w(ixo^s,iw_rho)*dtfactor*qdt
561 a3(ixo^s) = nabla_vp(ixo^s)*dtfactor*qdt
562 c0(ixo^s) = ((one+a2(ixo^s)+a3(ixo^s))*e_gas(ixo^s)+a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)/(one+a3(ixo^s))
563 c1(ixo^s) = (one+a2(ixo^s)+a3(ixo^s))/a1(ixo^s)/(one+a3(ixo^s))
566 {
do ix^d = ixomin^d,ixomax^d\}
571 call newton_method(e_gas(ix^d),e_rad(ix^d),c0(ix^d),c1(ix^d))
573 call halley_method(e_gas(ix^d),e_rad(ix^d),c0(ix^d),c1(ix^d))
575 call mpistop(
'root-method not known')
578 e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0+e_rad(ixo^s))/(one+a2(ixo^s)+a3(ixo^s))
580 w(ixo^s,iw_e) = e_gas(ixo^s)
583 w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(w(ixo^s,iw_mom(:))**2,dim=ndim+1)/w(ixo^s,iw_rho)
584 if(
allocated(iw_mag))
then
585 w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(w(ixo^s,iw_mag(:))**2,dim=ndim+1)
588 w(ixo^s,iw_r_e) = e_rad(ixo^s)
597 type(state),
target :: psa(max_blocks)
598 type(state),
target :: psb(max_blocks)
599 double precision,
intent(in) :: qdt
600 double precision,
intent(in) :: qtC
601 double precision,
intent(in) :: dtfactor
603 double precision :: res, max_residual, lambda
604 integer,
parameter :: max_its = 100
605 integer :: iw_to,iw_from
606 integer :: iigrid, igrid, id
607 integer :: nc, lvl, i
609 double precision :: fac, facD
618 print *,
'skipping implicit solve: dt too small!'
624 mg%operator_type = mg_vhelmholtz
625 mg%smoother_type = mg_smoother_gs
626 call mg_set_methods(
mg)
627 if(.not.
mg%is_allocated)
call mpistop(
"multigrid tree not allocated yet")
628 lambda = 1.d0/(dtfactor *
dt_diff)
629 call vhelmholtz_set_lambda(lambda)
637 call mg_restrict(
mg, mg_iveps)
638 call mg_fill_ghost_cells(
mg, mg_iveps)
640 call mg_fas_fmg(
mg, .true., max_res=res)
642 if(res < max_residual)
exit
643 call mg_fas_vcycle(
mg, max_res=res)
645 if(res .le. 0.d0)
then
647 if (
mg%my_rank == 0) &
648 write(*,*)
it,
' resiudal zero ', res
651 if(
mg%my_rank == 0)
then
653 error stop
"Diffusion residual to zero"
679 type(state),
target :: psa(max_blocks)
680 double precision,
intent(in) :: qtC
681 integer :: iigrid, igrid, level
687 do iigrid=1,igridstail; igrid=igrids(iigrid);
693 call getbc(qtc,0.d0,psa,1,nwflux+nwaux)
699 integer,
intent(in) :: ixI^L, ixO^L
700 double precision,
intent(inout) :: w(ixI^S, 1:nw)
701 double precision :: gradE(ixO^S),divF(ixO^S)
702 double precision :: divF_h(ixO^S),divF_j(ixO^S)
703 double precision :: diff_term(ixO^S)
704 integer :: idir, jxO^L, hxO^L
708 hxo^l=ixo^l-
kr(idir,^
d);
709 jxo^l=ixo^l+
kr(idir,^
d);
712 divf(ixo^s) = divf(ixo^s) + 2.d0*(divf_h(ixo^s) + divf_j(ixo^s))/
dxlevel(idir)**2
714 w(ixo^s,iw_r_e) = divf(ixo^s)
722 integer,
intent(in) :: ixI^L, ixO^L
723 double precision,
intent(inout) :: w(ixI^S,1:nw)
724 double precision,
intent(in) :: wCT(ixI^S,1:nw)
725 double precision,
intent(in) :: x(ixI^S,1:ndim)
726 integer :: idir,i,j,ix^D
727 double precision :: cc
728 double precision :: kappa(ixO^S),lambda(ixO^S),fld_R(ixO^S)
729 double precision :: max_D(ixI^S),grad_r_e(ixI^S),rad_e(ixI^S)
736 w(ixo^s,
i_diff_mg) = cc*lambda(ixo^s)/(kappa(ixo^s)*wct(ixo^s,iw_rho))
748 type(state),
target :: psa(max_blocks)
749 integer :: iigrid, igrid, level
754 do iigrid=1,igridstail; igrid=igrids(iigrid);
764 integer,
intent(in) :: ixI^L, ixO^L
765 double precision,
intent(inout) :: w(ixI^S, 1:nw)
766 double precision :: tmp_D(ixI^S), filtered_D(ixI^S)
767 integer :: ix^D, filter, idir
769 if(
size_d_filter .lt. 1)
call mpistop(
"D filter of size < 1 makes no sense")
772 filtered_d(ixo^s) = zero
776 filtered_d(ix^d) = filtered_d(ix^d) &
777 + tmp_d(ix^d+filter*
kr(idir,^d)) &
778 + tmp_d(ix^d-filter*
kr(idir,^d))
791 double precision,
intent(in) :: c0, c1
792 double precision,
intent(in) :: E_rad
793 double precision,
intent(inout) :: e_gas
794 double precision :: bisect_a, bisect_b, bisect_c
795 integer :: n, max_its
800 bisect_b = min(abs(c0/c1),abs(c0)**(1.d0/4.d0))+smalldouble
801 do while(abs(bisect_b-bisect_a) .ge.
fld_bisect_tol*min(e_gas,e_rad))
802 bisect_c = (bisect_a + bisect_b)/two
804 if(n .gt. max_its)
then
806 call mpistop(
'No convergece in bisection scheme')
829 call mpistop(
"Problem with fld bisection method")
837 print*,
"IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
847 2435 e_gas = (bisect_a + bisect_b)/two
853 double precision,
intent(in) :: c0, c1
854 double precision,
intent(in) :: E_rad
855 double precision,
intent(inout) :: e_gas
856 double precision :: xval, yval, der, deltax
872 print*,
'skip to bisection algorithm'
883 double precision,
intent(in) :: c0, c1
884 double precision,
intent(in) :: E_rad
885 double precision,
intent(inout) :: e_gas
886 double precision :: xval, yval, der, dder, deltax
900 deltax = -two*yval*der/(two*der**2 - yval*dder)
914 double precision,
intent(in) :: e_gas
915 double precision,
intent(in) :: c0, c1
916 double precision :: val
918 val = e_gas**4.d0 + c1*e_gas - c0
924 double precision,
intent(in) :: e_gas
925 double precision,
intent(in) :: c0, c1
926 double precision :: der
928 der = 4.d0*e_gas**3.d0 + c1
934 double precision,
intent(in) :: e_gas
935 double precision,
intent(in) :: c0, c1
936 double precision :: dder
938 dder = 4.d0*3.d0*e_gas**2.d0
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
Nicolas Moens Module for including flux limited diffusion (FLD)-approximation in Radiation-hydrodynam...
subroutine put_diffterm_onegrid(ixil, ixol, w)
inplace update of psa==>F_im(psa)
subroutine, public fld_get_opacity(w, x, ixil, ixol, fld_kappa)
Sets the opacity in the w-array by calling mod_opal_opacity.
double precision, public fld_mu
mean particle mass
double precision, public fld_bisect_tol
Tolerance for bisection method for Energy sourceterms This is a percentage of the minimum of gas- and...
subroutine energy_interaction(w, x, ixil, ixol, dtfactor, qdt)
Energy interaction and photon tiring.
double precision, public fld_diff_tol
Tolerance for adi method for radiative Energy diffusion.
subroutine fld_params_read(files)
public methods these are called in mod_rhd_phys
logical diffcrash_resume
Resume run when multigrid returns error.
subroutine update_diffcoeff(psa)
integer, dimension(:), allocatable, public i_opf
Index for Flux weighted opacities.
character(len=8) fld_opacity_law
Use constant Opacity?
subroutine, public get_fld_rad_force(qdt, ixil, ixol, wct, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO This subroutine handles th...
double precision function ddpolynomial_bisection(e_gas, c0, c1)
Evaluate second derivative of polynomial at argument e_gas.
subroutine evaluate_e_rad_mg(qtc, psa)
inplace update of psa==>F_im(psa)
subroutine, public fld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
subroutine, public fld_init(he_abundance, radiation_diffusion, energy_interact, r_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
character(len=16) fld_fluxlimiter
Diffusion limit lambda = 0.33.
subroutine fld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
subroutine bisection_method(e_gas, e_rad, c0, c1)
Find the root of the 4th degree polynomial using the bisection method.
character(len=50) fld_opal_table
double precision, public fld_kappa0
Opacity per unit of unit_density.
subroutine, public fld_get_fluxlimiter(w, x, ixil, ixol, fld_lambda, fld_r, nth)
Calculate fld flux limiter This subroutine calculates flux limiter lambda using the prescription stor...
character(len=8) fld_diff_scheme
Which method to solve diffusion part.
logical lineforce_opacities
Use or dont use lineforce opacities.
double precision function polynomial_bisection(e_gas, c0, c1)
Evaluate polynomial at argument e_gas.
logical fld_eint_split
source split for energy interact and radforce:
subroutine halley_method(e_gas, e_rad, c0, c1)
Find the root of the 4th degree polynomial using the Halley method.
character(len=8) fld_interaction_method
Which method to find the root for the energy interaction polynomial.
subroutine newton_method(e_gas, e_rad, c0, c1)
Find the root of the 4th degree polynomial using the Newton-Ralphson method.
logical fld_radforce_split
subroutine, public fld_get_radflux(w, x, ixil, ixol, rad_flux, nth)
Calculate Radiation Flux stores radiation flux in w-array.
double precision function dpolynomial_bisection(e_gas, c0, c1)
Evaluate first derivative of polynomial at argument e_gas.
subroutine fld_smooth_diffcoef(w, ixil, ixol)
Use running average on Diffusion coefficient.
subroutine diffuse_e_rad_mg(dtfactor, qdt, qtc, psa, psb)
Calling all subroutines to perform the multigrid method Communicates rad_e and diff_coeff to multigri...
subroutine, public fld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
integer i_diff_mg
diffusion coefficient for multigrid method
logical flux_lim_filter
Take a running average over the fluxlimiter.
subroutine fld_get_eddington(w, x, ixil, ixol, eddington_tensor, nth)
Calculate Eddington-tensor Stores Eddington-tensor in w-array.
logical diff_coef_filter
Take running average for Diffusion coefficient.
double precision dt_diff
running timestep for diffusion solver, initialised as zero
double precision, public diff_crit
Number for splitting the diffusion module.
subroutine fld_implicit_update(dtfactor, qdt, qtc, psa, psb)
Module with basic grid data structures.
Module with geometry-related routines (e.g., divergence, curl)
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision small_pressure
double precision unit_density
Physical scaling factor for density.
double precision unit_opacity
Physical scaling factor for Opacity.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer it
Number of time steps taken.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
double precision courantpar
The Courant (CFL) number used for the simulation.
double precision unit_velocity
Physical scaling factor for velocity.
logical time_advance
do time evolving
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc, factor, state_from)
Copy a variable to the multigrid tree, including a layer of ghost cells.
subroutine mg_copy_from_tree_gc(iw_from, iw_to, state_to)
Copy from multigrid tree with one layer of ghost cells. Corner ghost cells are not used/set.
This module reads in Rosseland-mean opacities from OPAL tables. Table opacity values are given in bas...
subroutine, public init_opal_table(tablename)
This subroutine is called when the FLD radiation module is initialised. The OPAL tables for different...
subroutine, public set_opal_opacity(rho, temp, kappa)
This subroutine calculates the opacity for a given temperature-density structure. Opacities are read ...
This module defines the procedures of a physics module. It contains function pointers for the various...
procedure(sub_get_tgas), pointer phys_get_tgas
procedure(sub_evaluate_implicit), pointer phys_evaluate_implicit
procedure(sub_get_pthermal), pointer phys_get_pthermal
procedure(sub_implicit_update), pointer phys_implicit_update
Module with all the methods that users can customize in AMRVAC.
procedure(special_opacity), pointer usr_special_opacity
procedure(special_diffcoef), pointer usr_special_diffcoef
procedure(special_fluxlimiter), pointer usr_special_fluxlimiter
procedure(special_opacity_qdot), pointer usr_special_opacity_qdot
integer function var_set_extravar(name_cons, name_prim, ix)
Set extra variable in w, which is not advected and has no boundary conditions. This has to be done af...
The data structure that contains information about a tree node/grid block.