18 double precision,
public ::
afld_mu = 0.6d0
50 double precision,
private,
protected :: afld_gamma
69 character(len=*),
intent(in) :: files(:)
79 open(
unitpar, file=trim(files(n)), status=
"old")
80 read(
unitpar, fld_list,
end=111)
92 subroutine afld_init(He_abundance, rhd_radiation_diffusion, afld_gamma)
98 double precision,
intent(in) :: he_abundance, afld_gamma
99 logical,
intent(in) :: rhd_radiation_diffusion
100 double precision :: sigma_thomson
102 character(len=1) :: ind_1
103 character(len=1) :: ind_2
104 character(len=2) :: cmp_f
105 character(len=5) :: cmp_e
108 call mpistop(
"Anisotropic in 1D... really?")
118 if (rhd_radiation_diffusion)
then
124 mg%operator_type = mg_ahelmholtz
129 write(ind_1,
'(I1)') idir
134 afld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
140 energy,qsourcesplit,active)
145 integer,
intent(in) :: ixi^
l, ixo^
l
146 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
147 double precision,
intent(in) :: wct(ixi^s,1:nw)
148 double precision,
intent(inout) :: w(ixi^s,1:nw)
149 logical,
intent(in) :: energy,qsourcesplit
150 logical,
intent(inout) :: active
151 double precision :: radiation_forcect(ixo^s,1:
ndim)
152 double precision :: kappact(ixo^s,1:
ndim)
153 double precision :: rad_fluxct(ixo^s,1:
ndim)
154 double precision :: div_v(ixi^s,1:
ndim,1:
ndim)
155 double precision :: edd(ixo^s,1:
ndim,1:
ndim)
156 double precision :: nabla_vp(ixo^s)
157 double precision :: vel(ixi^s), grad_v(ixi^s), grad0_v(ixo^s)
158 double precision :: grad_e(ixi^s)
159 integer :: idir, jdir
160 double precision :: afld_r(ixo^s,1:
ndim), lambda(ixo^s,1:
ndim)
169 radiation_forcect(ixo^s,idir) = -lambda(ixo^s,idir)*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)
182 vel(ixi^s) = wct(ixi^s,iw_mom(jdir))/wct(ixi^s,iw_rho)
184 div_v(ixo^s,idir,jdir) = grad_v(ixo^s)
191 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)
195 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
196 + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
197 + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
198 + div_v(ixo^s,2,2)*edd(ixo^s,2,2)
201 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
202 + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
203 + div_v(ixo^s,1,3)*edd(ixo^s,1,3) &
204 + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
205 + div_v(ixo^s,2,2)*edd(ixo^s,2,2) &
206 + div_v(ixo^s,2,3)*edd(ixo^s,2,3) &
207 + div_v(ixo^s,3,1)*edd(ixo^s,3,1) &
208 + div_v(ixo^s,3,2)*edd(ixo^s,3,2) &
209 + div_v(ixo^s,3,3)*edd(ixo^s,3,3)
211 w(ixo^s,iw_r_e) = w(ixo^s,iw_r_e) &
212 - qdt * nabla_vp(ixo^s)*wct(ixo^s,iw_r_e)
220 integer,
intent(in) :: ixi^
l, ixo^
l
221 double precision,
intent(in) ::
dx^
d, x(ixi^s,1:
ndim), w(ixi^s,1:nw)
222 double precision,
intent(inout) :: dtnew
224 double precision :: radiation_force(ixo^s,1:
ndim)
225 double precision :: lambda(ixo^s,1:
ndim), grad_e(ixi^s)
226 double precision :: dxinv(1:
ndim), afld_r(ixo^s,1:
ndim)
227 double precision :: max_grad
229 ^
d&dxinv(^
d)=one/
dx^
d;
233 radiation_force(ixo^s,idir) = -lambda(ixo^s,idir)*grad_e(ixo^s)
234 max_grad = maxval(abs(radiation_force(ixo^s,idir)))
235 max_grad = max(max_grad, epsilon(1.0d0))
236 dtnew = min(dtnew,
courantpar/sqrt(max_grad*dxinv(idir)))
243 energy,qsourcesplit,active)
248 integer,
intent(in) :: ixi^
l, ixo^
l
249 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
250 double precision,
intent(in) :: wct(ixi^s,1:nw)
251 double precision,
intent(inout) :: w(ixi^s,1:nw)
252 logical,
intent(in) :: energy,qsourcesplit
253 logical,
intent(inout) :: active
271 integer,
intent(in) :: ixi^
l, ixo^
l
272 double precision,
intent(in) :: w(ixi^s, 1:nw)
273 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
274 double precision,
intent(out) :: fld_kappa(ixo^s,1:
ndim)
275 integer :: i,j,ix^
d, idir
276 double precision :: temp(ixi^s), pth(ixi^s), a2(ixo^s)
277 double precision :: rho0,temp0,n,sigma_b
278 double precision :: akram, bkram
279 double precision :: vth(ixo^s), gradv(ixi^s), eta(ixo^s), t(ixo^s)
295 fld_kappa(ixo^s,idir) =
fld_kappa0/
unit_opacity*(one + n*dexp(-one/sigma_b*(dlog(w(ixo^s,iw_rho)/rho0))**two))
298 temp(ixo^s) = temp(ixo^s)/w(ixo^s,iw_rho)
306 akram = 13.1351597305d0
307 bkram = -4.5182188206d0
309 * (1.d0+10.d0**akram*w(ixo^s,iw_rho)*
unit_density*(a2(ixo^s)/1.d12)**bkram)
310 {
do ix^
d=ixomin^
d,ixomax^
d\ }
316 {
do ix^
d=ixomin^
d,ixomax^
d\ }
324 call mpistop(
"special opacity not defined")
328 call mpistop(
"Doesn't know opacity law")
341 integer,
intent(in) :: ixi^
l, ixo^
l, nth
342 double precision,
intent(in) :: w(ixi^s, 1:nw)
343 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
344 double precision,
intent(out) :: fld_r(ixo^s,1:
ndim), fld_lambda(ixo^s,1:
ndim)
345 integer :: idir, ix^
d
346 integer :: filter, idim
347 double precision :: kappa(ixo^s,1:
ndim)
348 double precision :: normgrad2(ixo^s)
349 double precision :: grad_r_e(ixi^s)
350 double precision :: rad_e(ixi^s)
351 double precision :: tmp_l(ixi^s), filtered_l(ixi^s)
355 fld_lambda(ixo^s,1:
ndim) = 1.d0/3.d0
356 fld_r(ixo^s,1:
ndim) = zero
360 rad_e(ixi^s) = w(ixi^s, iw_r_e)
363 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
364 normgrad2(ixo^s) = grad_r_e(ixo^s)**2
365 fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
367 fld_lambda(ixo^s,idir) = one/fld_r(ixo^s,idir)
372 normgrad2(ixo^s) = zero
373 rad_e(ixi^s) = w(ixi^s,iw_r_e)
376 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
377 normgrad2(ixo^s) = grad_r_e(ixo^s)**2
378 fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
381 fld_lambda(ixo^s,idir) = (2.d0+fld_r(ixo^s,idir))/(6.d0+3*fld_r(ixo^s,idir)+fld_r(ixo^s,idir)**2.d0)
386 normgrad2(ixo^s) = zero
387 rad_e(ixi^s) = w(ixi^s,iw_r_e)
390 call gradient(rad_e,ixi^
l,ixo^
l,idir,grad_r_e,nth)
391 normgrad2(ixo^s) = grad_r_e(ixo^s)**2
392 fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
395 {
do ix^
d=ixomin^
d,ixomax^
d\ }
396 if(fld_r(ix^
d,idir) .lt. 3.d0/2.d0)
then
397 fld_lambda(ix^
d,idir) = 2.d0/(3.d0 + dsqrt(9.d0 + 12.d0*fld_r(ix^
d,idir)**2.d0))
399 fld_lambda(ix^
d,idir) = 1.d0/(1.d0 + fld_r(ix^
d,idir) + dsqrt(1.d0 + 2.d0*fld_r(ix^
d,idir)))
405 call mpistop(
"special fluxlimiter not defined")
409 call mpistop(
'Fluxlimiter unknown')
413 if(
size_l_filter .lt. 1)
call mpistop(
"D filter of size < 1 makes no sense")
416 tmp_l(ixo^s) = fld_lambda(ixo^s,idir)
417 filtered_l(ixo^s) = zero
421 filtered_l(ix^
d) = filtered_l(ix^
d) &
422 + tmp_l(ix^
d+filter*
kr(idim,^
d)) &
423 + tmp_l(ix^
d-filter*
kr(idim,^
d))
431 fld_lambda(ixo^s,idir) = tmp_l(ixo^s)
441 integer,
intent(in) :: ixi^
l, ixo^
l, nth
442 double precision,
intent(in) :: w(ixi^s, 1:nw)
443 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
444 double precision,
intent(out) :: rad_flux(ixo^s, 1:
ndim)
445 integer :: ix^
d, idir
446 double precision :: cc
447 double precision :: grad_r_e(ixi^s)
448 double precision :: rad_e(ixi^s)
449 double precision :: kappa(ixo^s,1:
ndim), lambda(ixo^s,1:
ndim), fld_r(ixo^s,1:
ndim)
452 rad_e(ixi^s) = w(ixi^s, iw_r_e)
457 rad_flux(ixo^s, idir) = -cc*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*w(ixo^s,iw_rho))*grad_r_e(ixo^s)
466 integer,
intent(in) :: ixI^L, ixO^L, nth
467 double precision,
intent(in) :: w(ixI^S, 1:nw)
468 double precision,
intent(in) :: x(ixI^S, 1:ndim)
469 double precision,
intent(out) :: eddington_tensor(ixO^S,1:ndim,1:ndim)
470 double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
471 double precision :: normgrad2(ixO^S), f(ixO^S,1:ndim)
472 double precision :: grad_r_e(ixI^S,1:ndim), rad_e(ixI^S)
473 double precision :: lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
474 integer :: i,j, idir,jdir
479 normgrad2(ixo^s) = zero
480 rad_e(ixi^s) = w(ixi^s, iw_r_e)
482 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e(ixi^s,idir),nth)
483 normgrad2(ixo^s) = normgrad2(ixo^s)+grad_r_e(ixo^s,idir)**2
486 f(ixo^s,idir) = lambda(ixo^s,idir)+lambda(ixo^s,idir)**two*fld_r(ixo^s,idir)**two
487 f(ixo^s,idir) = half*(one-f(ixo^s,idir))+half*(3.d0*f(ixo^s,idir)-one)
491 eddington_tensor(ixo^s,1,1) = f(ixo^s,1)
495 eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s,idir))
499 if(idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
500 tnsr2(ixo^s,idir,jdir) = half*(3.d0*f(ixo^s,idir) - 1)&
501 * grad_r_e(ixo^s,idir)*grad_r_e(ixo^s,jdir)/normgrad2(ixo^s)
506 where((tnsr2(ixo^s,idir,jdir) .eq. tnsr2(ixo^s,idir,jdir)) &
507 .and. (normgrad2(ixo^s) .gt. smalldouble))
508 eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir)+tnsr2(ixo^s,idir,jdir)
519 integer,
intent(in) :: ixi^
l, ixo^
l, nth
520 double precision,
intent(in) :: w(ixi^s, 1:nw)
521 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
522 double precision,
intent(out):: rad_pressure(ixo^s,1:
ndim,1:
ndim)
524 double precision :: eddington_tensor(ixo^s,1:
ndim,1:
ndim)
529 rad_pressure(ixo^s,i,j) = eddington_tensor(ixo^s,i,j)*w(ixo^s,iw_r_e)
542 type(state),
target :: psa(max_blocks)
543 type(state),
target :: psb(max_blocks)
544 double precision,
intent(in) :: qdt
545 double precision,
intent(in) :: qtC
546 double precision,
intent(in) :: dtfactor
547 integer,
parameter :: max_its = 50
549 integer :: iw_to,iw_from
550 integer :: iigrid, igrid, id
551 integer :: nc, lvl, i
553 double precision :: fac, facD
554 double precision :: res, max_residual, lambda
560 if (qdt <
dtmin)
then
562 print *,
'skipping implicit solve: dt too small!'
569 mg%operator_type = mg_ahelmholtz
570 mg%smoother_type = mg_smoother_gsrb
571 call mg_set_methods(
mg)
572 if (.not.
mg%is_allocated)
call mpistop(
"multigrid tree not allocated yet")
573 lambda = 1.d0/(dtfactor *
dt_diff)
574 call ahelmholtz_set_lambda(lambda)
582 call mg_restrict(
mg, mg_iveps)
583 call mg_fill_ghost_cells(
mg, mg_iveps)
590 call mg_restrict(
mg, mg_iveps1)
591 call mg_fill_ghost_cells(
mg, mg_iveps1)
593 call mg_restrict(
mg, mg_iveps2)
594 call mg_fill_ghost_cells(
mg, mg_iveps2)
602 call mg_fas_fmg(
mg, .true., max_res=res)
605 if (res < max_residual)
exit
606 call mg_fas_vcycle(
mg, max_res=res)
608 if (res .le. 0.d0)
then
610 if (
mg%my_rank == 0) &
611 write(*,*)
it,
' resiudal zero ', res
614 if (
mg%my_rank == 0)
then
616 error stop
"Diffusion residual to zero"
619 if(n == max_its + 1)
then
621 if(
mg%my_rank == 0) &
622 write(*,*)
it,
' resiudal high ', res
625 if(
mg%my_rank == 0)
then
626 print *,
"Did you specify boundary conditions correctly?"
627 print *,
"Or is the variation in diffusion too large?"
629 print *,
mg%bc(1, mg_iphi)%bc_value,
mg%bc(2, mg_iphi)%bc_value
631 error stop
"diffusion_solve: no convergence"
643 type(state),
target :: psa(max_blocks)
644 double precision,
intent(in) :: qtC
645 integer :: iigrid, igrid, level
651 do iigrid=1,igridstail; igrid=igrids(iigrid);
658 call getbc(qtc,0.d0,psa,1,nwflux+nwaux)
664 integer,
intent(in) :: ixI^L, ixO^L
665 double precision,
intent(inout) :: w(ixI^S, 1:nw)
666 double precision :: gradE(ixO^S),divF(ixO^S)
667 double precision :: divF_h(ixO^S),divF_j(ixO^S)
668 double precision :: diff_term(ixO^S)
669 integer :: idir, jxO^L, hxO^L
673 hxo^l=ixo^l-
kr(idir,^
d);
674 jxo^l=ixo^l+
kr(idir,^
d);
677 divf(ixo^s) = divf(ixo^s) + 2.d0*(divf_h(ixo^s) + divf_j(ixo^s))/
dxlevel(idir)**2
679 w(ixo^s,iw_r_e) = divf(ixo^s)
687 integer,
intent(in) :: ixI^L, ixO^L
688 double precision,
intent(inout) :: w(ixI^S, 1:nw)
689 double precision,
intent(in) :: wCT(ixI^S, 1:nw)
690 double precision,
intent(in) :: x(ixI^S, 1:ndim)
691 double precision :: kappa(ixO^S,1:ndim), lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
692 double precision :: max_D(ixI^S), grad_r_e(ixI^S), rad_e(ixI^S)
693 double precision :: cc
694 integer :: idir,i,j,ix^D
704 w(ixo^s,
i_diff_mg(idir)) = cc*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*wct(ixo^s,iw_rho))
705 where(w(ixo^s,
i_diff_mg(idir)) .lt. 0.d0)
719 type(state),
target :: psa(max_blocks)
720 integer :: iigrid, igrid, level
725 do iigrid=1,igridstail; igrid=igrids(iigrid);
735 integer,
intent(in) :: ixI^L, ixO^L, idir
736 double precision,
intent(inout) :: w(ixI^S, 1:nw)
737 double precision :: tmp_D(ixI^S), filtered_D(ixI^S)
738 integer :: ix^D, filter, idim
740 if(
size_d_filter .lt. 1)
call mpistop(
"D filter of size < 1 makes no sense")
743 filtered_d(ixo^s) = zero
747 filtered_d(ix^d) = filtered_d(ix^d) &
748 + tmp_d(ix^d+filter*
kr(idim,^d)) &
749 + tmp_d(ix^d-filter*
kr(idim,^d))
769 integer,
intent(in) :: ixI^L, ixO^L
770 double precision,
intent(in) :: x(ixI^S,1:ndim)
771 double precision,
intent(in) :: wCT(ixI^S,1:nw)
772 double precision,
intent(inout) :: w(ixI^S,1:nw)
773 double precision :: a1(ixO^S), a2(ixO^S)
774 double precision :: c0(ixO^S), c1(ixO^S)
775 double precision :: e_gas(ixO^S), E_rad(ixO^S)
776 double precision :: kappa(ixO^S)
777 double precision :: sigma_b,cc
783 e_gas(ixo^s) = wct(ixo^s,iw_e)-half*sum(wct(ixo^s, iw_mom(:))**2, dim=ndim+1)/wct(ixo^s, iw_rho)
784 if(
allocated(iw_mag))
then
785 e_gas(ixo^s) = e_gas(ixo^s)-half*sum(wct(ixo^s,iw_mag(:))**2,dim=ndim+1)
787 {
do ix^d=ixomin^d,ixomax^d\ }
790 e_rad(ixo^s) = wct(ixo^s,iw_r_e)
797 a1(ixo^s) = 4.d0*sigma_b/cc*(afld_gamma-one)**4.d0/wct(ixo^s,iw_rho)**4.d0
798 a2(ixo^s) = e_gas(ixo^s) + e_rad(ixo^s)
799 c0(ixo^s) = a2(ixo^s)/a1(ixo^s)
800 c1(ixo^s) = 1.d0/a1(ixo^s)
803 a1(ixo^s) = 4.d0*kappa(ixo^s)*sigma_b*(afld_gamma-one)**4.d0/wct(ixo^s,iw_rho)**3.d0*
dt
805 c0(ixo^s) = ((one+a2(ixo^s))*e_gas(ixo^s)+a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)
806 c1(ixo^s) = (one+a2(ixo^s))/a1(ixo^s)
809 {
do ix^d=ixomin^d,ixomax^d\ }
814 call newton_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
816 call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
818 call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
820 call mpistop(
'root-method not known')
825 e_rad(ixo^s) = a2(ixo^s)-e_gas(ixo^s)
828 e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0+e_rad(ixo^s))/(one + a2(ixo^s))
831 w(ixo^s,iw_e) = e_gas(ixo^s)
834 w(ixo^s,iw_e) = w(ixo^s,iw_e) + half*sum(wct(ixo^s, iw_mom(:))**2,dim=ndim+1)/wct(ixo^s, iw_rho)
835 if(
allocated(iw_mag))
then
836 w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(wct(ixo^s, iw_mag(:))**2,dim=ndim+1)
839 w(ixo^s,iw_r_e) = e_rad(ixo^s)
845 double precision,
intent(in) :: c0, c1
846 double precision,
intent(in) :: E_rad
847 double precision,
intent(inout) :: e_gas
848 double precision :: bisect_a, bisect_b, bisect_c
849 integer :: n, max_its
854 bisect_b = 1.2d0*max(abs(c0/c1),abs(c0)**(1.d0/4.d0))
855 do while(abs(bisect_b-bisect_a) .ge.
fld_bisect_tol*min(e_gas,e_rad))
856 bisect_c = (bisect_a + bisect_b)/two
858 if(n .gt. max_its)
then
860 call mpistop(
'No convergece in bisection scheme')
883 call mpistop(
"Problem with fld bisection method")
891 print*,
"IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
901 2435 e_gas = (bisect_a + bisect_b)/two
907 double precision,
intent(in) :: c0, c1
908 double precision,
intent(in) :: E_rad
909 double precision,
intent(inout) :: e_gas
910 double precision :: xval, yval, der, deltax
925 if (ii .gt. 1d3)
then
936 double precision,
intent(in) :: c0, c1
937 double precision,
intent(in) :: E_rad
938 double precision,
intent(inout) :: e_gas
939 double precision :: xval, yval, der, dder, deltax
953 deltax = -two*yval*der/(two*der**2 - yval*dder)
956 if (ii .gt. 1d3)
then
967 double precision,
intent(in) :: e_gas
968 double precision,
intent(in) :: c0, c1
969 double precision :: val
971 val = e_gas**4.d0 + c1*e_gas - c0
977 double precision,
intent(in) :: e_gas
978 double precision,
intent(in) :: c0, c1
979 double precision :: der
981 der = 4.d0*e_gas**3.d0 + c1
987 double precision,
intent(in) :: e_gas
988 double precision,
intent(in) :: c0, c1
989 double precision :: dder
991 dder = 4.d0*3.d0*e_gas**2.d0
Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynami...
logical lineforce_opacities
Use or dont use lineforce opacities.
subroutine afld_smooth_diffcoef(w, ixil, ixol, idir)
Use running average on Diffusion coefficient.
character(len=8) fld_interaction_method
Which method to find the root for the energy interaction polynomial.
logical fld_eint_split
source split for energy interact and radforce:
subroutine afld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
logical flux_lim_filter
Take a running average over the fluxlimiter.
subroutine, public afld_get_radflux(w, x, ixil, ixol, rad_flux, nth)
Calculate Radiation Flux Calculate the Flux using the fld closure relation F = -c*lambda/(kappa*rho) ...
double precision, public fld_bisect_tol
Tolerance for bisection method for Energy sourceterms This is a percentage of the minimum of gas- and...
subroutine, public get_afld_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, public fld_kappa0
Opacity per unit of unit_density.
subroutine, public afld_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...
double precision, public fld_diff_tol
Tolerance for adi method for radiative Energy diffusion.
character(len=16) fld_fluxlimiter
Diffusion limit lambda = 0.33.
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 afld_init(he_abundance, rhd_radiation_diffusion, afld_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
double precision function ddpolynomial_bisection(e_gas, c0, c1)
Evaluate second derivative of polynomial at argument e_gas.
subroutine bisection_method(e_gas, e_rad, c0, c1)
Find the root of the 4th degree polynomial using the bisection method.
logical fld_diff_testcase
Set Diffusion coefficient to unity.
logical diffcrash_resume
Resume run when multigrid returns error.
double precision function dpolynomial_bisection(e_gas, c0, c1)
Evaluate first derivative of polynomial at argument e_gas.
character(len=8), dimension(:), allocatable fld_opacity_law
Use constant Opacity?
subroutine, public afld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
subroutine evaluate_e_rad_mg(qtc, psa)
inplace update of psa==>F_im(psa)
subroutine energy_interaction(w, wct, x, ixil, ixol)
This subroutine calculates the radiation heating, radiation cooling and photon tiring using an implic...
subroutine, public afld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
subroutine update_diffcoeff(psa)
subroutine halley_method(e_gas, e_rad, c0, c1)
Find the root of the 4th degree polynomial using the Halley method.
double precision function polynomial_bisection(e_gas, c0, c1)
Evaluate polynomial at argument e_gas.
subroutine newton_method(e_gas, e_rad, c0, c1)
Find the root of the 4th degree polynomial using the Newton-Ralphson method.
character(len=50) fld_opal_table
character(len=8) afld_diff_scheme
Which method to solve diffusion part.
subroutine, public get_afld_energy_interact(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, public diff_crit
Number for splitting the diffusion module.
logical fld_radforce_split
subroutine, public afld_get_opacity(w, x, ixil, ixol, fld_kappa)
Sets the opacity in the w-array by calling mod_opal_opacity.
logical diff_coef_filter
Take running average for Diffusion coefficient.
double precision, public afld_mu
mean particle mass
integer, dimension(:), allocatable i_diff_mg
subroutine put_diffterm_onegrid(ixil, ixol, w)
inplace update of psa==>F_im(psa)
subroutine afld_params_read(files)
public methods these are called in mod_rhd_phys
double precision dt_diff
running timestep for diffusion solver, initialised as zero
subroutine afld_get_eddington(w, x, ixil, ixol, eddington_tensor, nth)
Calculate Eddington-tensor Stores Eddington-tensor in w-array.
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
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 dt
global time step
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_set_mg_bounds), pointer phys_set_mg_bounds
procedure(sub_implicit_update), pointer phys_implicit_update
Module with all the methods that users can customize in AMRVAC.
procedure(special_diffcoef), pointer usr_special_diffcoef
procedure(special_fluxlimiter), pointer usr_special_fluxlimiter
procedure(special_opacity_qdot), pointer usr_special_opacity_qdot
procedure(special_aniso_opacity), pointer usr_special_aniso_opacity
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.