mod_afld Module Reference

Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynamics simulations using mod_rhd Based on Turner and stone 2001. See [1]Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R., “Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion”, Astronomy and Astrophysics, vol. 657, 2022. doi:10.1051/0004-6361/202141023. For more information. More...

Functions/Subroutines
subroutine	afld_params_read (files)
	public methods these are called in mod_rhd_phys More...
subroutine, public	afld_init (He_abundance, rhd_radiation_diffusion, afld_gamma)
	Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variables to w-array, flux, kappa, eddington Tensor Lambda and R ... More...
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 the radiation force More...
subroutine, public	afld_radforce_get_dt (w, ixIL, ixOL, dtnew, dxD, x)
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 the energy exchange between gas and radiation More...
subroutine, public	afld_get_opacity (w, x, ixIL, ixOL, fld_kappa)
	Sets the opacity in the w-array by calling mod_opal_opacity. More...
subroutine, public	afld_get_fluxlimiter (w, x, ixIL, ixOL, fld_lambda, fld_R)
	Calculate fld flux limiter This subroutine calculates flux limiter lambda using the prescription stored in fld_fluxlimiter. It also calculates the ratio of radiation scaleheight and mean free path. More...
subroutine, public	afld_get_radflux (w, x, ixIL, ixOL, rad_flux)
	Calculate Radiation Flux stores radiation flux in w-array. More...
subroutine	afld_get_eddington (w, x, ixIL, ixOL, eddington_tensor)
	Calculate Eddington-tensor Stores Eddington-tensor in w-array. More...
subroutine, public	afld_get_radpress (w, x, ixIL, ixOL, rad_pressure)
	Calculate Radiation Pressure Returns Radiation Pressure as tensor. More...
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 multigrid library. More...
subroutine	evaluate_e_rad_mg (qtC, psa)
	inplace update of psa==>F_im(psa) More...
subroutine	put_diffterm_onegrid (ixIL, ixOL, w)
	inplace update of psa==>F_im(psa) More...
subroutine	afld_get_diffcoef_central (w, wCT, x, ixIL, ixOL)
	Calculates cell-centered diffusion coefficient to be used in multigrid. More...
subroutine	update_diffcoeff (psa)
subroutine	afld_smooth_diffcoef (w, ixIL, ixOL, idim)
	Use running average on Diffusion coefficient. More...
subroutine	energy_interaction (w, wCT, x, ixIL, ixOL)
	This subroutine calculates the radiation heating, radiation cooling and photon tiring using an implicit scheme. These sourceterms are applied using the root-finding of a 4th order polynomial This routine loops over every cell in the domain and computes the coefficients of the polynomials in every cell. More...
subroutine	bisection_method (e_gas, E_rad, c0, c1)
	Find the root of the 4th degree polynomial using the bisection method. More...
subroutine	newton_method (e_gas, E_rad, c0, c1)
	Find the root of the 4th degree polynomial using the Newton-Ralphson method. More...
subroutine	halley_method (e_gas, E_rad, c0, c1)
	Find the root of the 4th degree polynomial using the Halley method. More...
double precision function	polynomial_bisection (e_gas, c0, c1)
	Evaluate polynomial at argument e_gas. More...
double precision function	dpolynomial_bisection (e_gas, c0, c1)
	Evaluate first derivative of polynomial at argument e_gas. More...
double precision function	ddpolynomial_bisection (e_gas, c0, c1)
	Evaluate second derivative of polynomial at argument e_gas. More...
subroutine	gradiento (q, x, ixIL, ixOL, idir, gradq, n)
	Calculate gradient of a scalar q within ixL in direction idir difference with gradient is gradq(ixO^S), NOT gradq(ixI^S) More...

Variables
logical	fld_eint_split = .false.
	source split for energy interact and radforce: More...
logical	fld_radforce_split = .false.
double precision, public	fld_kappa0 = 0.34d0
	Opacity per unit of unit_density. More...
double precision, public	afld_mu = 0.6d0
	mean particle mass More...
double precision, public	fld_bisect_tol = 1.d-4
	Tolerance for bisection method for Energy sourceterms This is a percentage of the minimum of gas- and radiation energy. More...
double precision, public	fld_diff_tol = 1.d-4
	Tolerance for adi method for radiative Energy diffusion. More...
double precision, public	diff_crit
	Number for splitting the diffusion module. More...
integer, public	i_test
	Index for testvariable. More...
character(len=8), dimension(:), allocatable	fld_opacity_law
	Use constant Opacity? More...
character(len=50)	fld_opal_table = 'Y09800'
character(len=16)	fld_fluxlimiter = 'Pomraning'
	Diffusion limit lambda = 0.33. More...
integer, dimension(:), allocatable, public	i_diff_mg
	Index for Diffusion coeficients. More...
character(len=8)	afld_diff_scheme = 'mg'
	Which method to solve diffusion part. More...
character(len=8)	fld_interaction_method = 'Halley'
	Which method to find the root for the energy interaction polynomial. More...
logical	fld_diff_testcase = .false.
	Set Diffusion coefficient to unity. More...
logical	diff_coef_filter = .false.
	Take running average for Diffusion coefficient. More...
integer	size_d_filter = 1
logical	flux_lim_filter = .false.
	Take a running average over the fluxlimiter. More...
integer	size_l_filter = 1
logical	lineforce_opacities = .false.
	Use or dont use lineforce opacities. More...
logical	diffcrash_resume = .true.
	Resume run when multigrid returns error. More...
double precision	dt_diff = 0.d0
	running timestep for diffusion solver, initialised as zero More...

Detailed Description

Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynamics simulations using mod_rhd Based on Turner and stone 2001. See [1]Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R., “Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion”, Astronomy and Astrophysics, vol. 657, 2022. doi:10.1051/0004-6361/202141023. For more information.

Function/Subroutine Documentation

afld_get_diffcoef_central()

subroutine mod_afld::afld_get_diffcoef_central	(	double precision, dimension(ixi^s, 1:nw), intent(inout)	w,
		double precision, dimension(ixi^s, 1:nw), intent(in)	wCT,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L
	)

Calculates cell-centered diffusion coefficient to be used in multigrid.

calculate diffusion coefficient

Definition at line 1262 of file mod_afld.t.

Here is the call graph for this function:

afld_get_eddington()

subroutine mod_afld::afld_get_eddington	(	double precision, dimension(ixi^s, 1:nw), intent(in)	w,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixo^s,1:ndim,1:ndim), intent(out)	eddington_tensor
	)

Calculate Eddington-tensor Stores Eddington-tensor in w-array.

Calculate R everywhere |grad E|/(rho kappa E)

Calculate radiation pressure P = (lambda + lambda^2 R^2)*E

Definition at line 639 of file mod_afld.t.

Here is the call graph for this function:

afld_get_fluxlimiter()

subroutine, public mod_afld::afld_get_fluxlimiter	(	double precision, dimension(ixi^s, 1:nw), intent(in)	w,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixo^s,1:ndim), intent(out)	fld_lambda,
		double precision, dimension(ixo^s,1:ndim), intent(out)	fld_R
	)

Calculate fld flux limiter This subroutine calculates flux limiter lambda using the prescription stored in fld_fluxlimiter. It also calculates the ratio of radiation scaleheight and mean free path.

Calculate R everywhere |grad E|/(rho kappa E)

Calculate the flux limiter, lambda

Calculate R everywhere |grad E|/(rho kappa E)

Calculate the flux limiter, lambda Levermore and Pomraning: lambda = (2 + R)/(6 + 3R + R^2)

Calculate R everywhere |grad E|/(rho kappa E)

Calculate the flux limiter, lambda Minerbo:

Definition at line 460 of file mod_afld.t.

Here is the call graph for this function:

afld_get_opacity()

subroutine, public mod_afld::afld_get_opacity	(	double precision, dimension(ixi^s, 1:nw), intent(in)	w,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixo^s,1:ndim), intent(out)	fld_kappa
	)

Sets the opacity in the w-array by calling mod_opal_opacity.

Take lower value of rho in domain

Opacity bump

Take lower value of rho in domain

Hard limit on kappa

Limit kappa through T

Definition at line 370 of file mod_afld.t.

Here is the call graph for this function:

afld_get_radflux()

subroutine, public mod_afld::afld_get_radflux	(	double precision, dimension(ixi^s, 1:nw), intent(in)	w,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixo^s, 1:ndim), intent(out)	rad_flux
	)

Calculate Radiation Flux stores radiation flux in w-array.

Calculate the Flux using the fld closure relation F = -c*lambda/(kappa*rho) *grad E

Definition at line 606 of file mod_afld.t.

Here is the call graph for this function:

afld_get_radpress()

subroutine, public mod_afld::afld_get_radpress	(	double precision, dimension(ixi^s, 1:nw), intent(in)	w,
		double precision, dimension(ixi^s, 1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixo^s,1:ndim,1:ndim), intent(out)	rad_pressure
	)

Calculate Radiation Pressure Returns Radiation Pressure as tensor.

Definition at line 715 of file mod_afld.t.

Here is the call graph for this function:

afld_init()

subroutine, public mod_afld::afld_init	(	double precision, intent(in)	He_abundance,
		logical, intent(in)	rhd_radiation_diffusion,
		double precision, intent(in)	afld_gamma
	)

Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variables to w-array, flux, kappa, eddington Tensor Lambda and R ...

allocate boundary conditions

read par files

Introduce test variable globally

DELETE WHEN DONE

Need mean molecular weight

set rhd_gamma

Definition at line 125 of file mod_afld.t.

Here is the call graph for this function:

afld_params_read()

subroutine mod_afld::afld_params_read

(

character(len=*), dimension(:), intent(in)

files

)

public methods these are called in mod_rhd_phys

Reading in fld-list parameters from .par file

Definition at line 99 of file mod_afld.t.

afld_radforce_get_dt()

subroutine, public mod_afld::afld_radforce_get_dt	(	double precision, dimension(ixi^s,1:nw), intent(in)	w,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, intent(inout)	dtnew,
		double precision, intent(in)	dx,
		double precision, intent(in)	D,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x
	)

Definition at line 315 of file mod_afld.t.

Here is the call graph for this function:

afld_smooth_diffcoef()

subroutine mod_afld::afld_smooth_diffcoef	(	double precision, dimension(ixi^s, 1:nw), intent(inout)	w,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		integer, intent(in)	idim
	)

Use running average on Diffusion coefficient.

Definition at line 1328 of file mod_afld.t.

bisection_method()

subroutine mod_afld::bisection_method	(	double precision, intent(inout)	e_gas,
		double precision, intent(in)	E_rad,
		double precision, intent(in)	c0,
		double precision, intent(in)	c1
	)

Find the root of the 4th degree polynomial using the bisection method.

Definition at line 1494 of file mod_afld.t.

Here is the call graph for this function:

ddpolynomial_bisection()

double precision function mod_afld::ddpolynomial_bisection	(	double precision, intent(in)	e_gas,
		double precision, intent(in)	c0,
		double precision, intent(in)	c1
	)

Evaluate second derivative of polynomial at argument e_gas.

Definition at line 1663 of file mod_afld.t.

diffuse_e_rad_mg()

subroutine mod_afld::diffuse_e_rad_mg	(	double precision, intent(in)	dtfactor,
		double precision, intent(in)	qdt,
		double precision, intent(in)	qtC,
		type(state), dimension(max_blocks), target	psa,
		type(state), dimension(max_blocks), target	psb
	)

Calling all subroutines to perform the multigrid method Communicates rad_e and diff_coeff to multigrid library.

Set diffusion timestep, add previous timestep if mg did not converge:

replace call set_rhs(mg, -1/dt, 0.0_dp)

Reset dt_diff when diffusion worked out

Definition at line 1071 of file mod_afld.t.

Here is the call graph for this function:

dpolynomial_bisection()

double precision function mod_afld::dpolynomial_bisection	(	double precision, intent(in)	e_gas,
		double precision, intent(in)	c0,
		double precision, intent(in)	c1
	)

Evaluate first derivative of polynomial at argument e_gas.

Definition at line 1652 of file mod_afld.t.

energy_interaction()

subroutine mod_afld::energy_interaction	(	double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:nw), intent(in)	wCT,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L
	)

This subroutine calculates the radiation heating, radiation cooling and photon tiring using an implicit scheme. These sourceterms are applied using the root-finding of a 4th order polynomial This routine loops over every cell in the domain and computes the coefficients of the polynomials in every cell.

e_gas is the INTERNAL ENERGY without KINETIC ENERGY

Calculate coefficients for polynomial

Loop over every cell for rootfinding method

advance E_rad

new w = w + dt f(wCT) e_gas,E_rad = wCT + dt f(wCT) dt f(wCT) = e_gas,E_rad - wCT new w = w + e_gas,E_rad - wCT

WAIT A SECOND?! DOES THIS EVEN WORK WITH THESE KIND OF IMPLICIT METHODS? NOT QUITE SURE TO BE HONEST IS IT POSSIBLE TO SHUT DOWN SOURCE SPLITTING FOR JUST THIS TERM? FIX BY PERFORMING Energy_interaction on (w,w,...)

Update gas-energy in w, internal + kinetic

Beginning of module substracted wCT Ekin So now add wCT Ekin

Update rad-energy in w

Definition at line 1371 of file mod_afld.t.

Here is the call graph for this function:

evaluate_e_rad_mg()

subroutine mod_afld::evaluate_e_rad_mg	(	double precision, intent(in)	qtC,
		type(state), dimension(max_blocks), target	psa
	)

inplace update of psa==>F_im(psa)

Definition at line 1205 of file mod_afld.t.

Here is the call graph for this function:

get_afld_energy_interact()

subroutine, public mod_afld::get_afld_energy_interact	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixi^s,1:nw), intent(in)	wCT,
		double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		logical, intent(in)	energy,
		logical, intent(in)	qsourcesplit,
		logical, intent(inout)	active
	)

w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO This subroutine handles the energy exchange between gas and radiation

Calculate and add sourceterms

Add energy sourceterms

Definition at line 345 of file mod_afld.t.

Here is the call graph for this function:

get_afld_rad_force()

subroutine, public mod_afld::get_afld_rad_force	(	double precision, intent(in)	qdt,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixi^s,1:nw), intent(in)	wCT,
		double precision, dimension(ixi^s,1:nw), intent(inout)	w,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		logical, intent(in)	energy,
		logical, intent(in)	qsourcesplit,
		logical, intent(inout)	active
	)

w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO This subroutine handles the radiation force

Calculate and add sourceterms

Radiation force = kappa*rho/c *Flux = lambda gradE

Momentum equation source term

Energy equation source term (kinetic energy)

Photon tiring calculate tensor div_v !$OMP PARALLEL DO

!$OMP END PARALLEL DO

VARIABLE NAMES DIV ARE ACTUALLY GRADIENTS

eq 34 Turner and stone (Only 2D)

Definition at line 206 of file mod_afld.t.

Here is the call graph for this function:

gradiento()

subroutine mod_afld::gradiento	(	double precision, dimension(ixi^s), intent(in)	q,
		double precision, dimension(ixi^s,1:ndim), intent(in)	x,
		integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		integer, intent(in)	idir,
		double precision, dimension(ixo^s), intent(out)	gradq,
		integer, intent(in)	n
	)

Calculate gradient of a scalar q within ixL in direction idir difference with gradient is gradq(ixO^S), NOT gradq(ixI^S)

Using higher order derivatives with wider stencil according to: https://en.wikipedia.org/wiki/Finite_difference_coefficient

coef 2/3

coef -1/12

divide by dx

coef 3/4

coef -3/20

coef 1/60

divide by dx

Definition at line 1675 of file mod_afld.t.

halley_method()

subroutine mod_afld::halley_method	(	double precision, intent(inout)	e_gas,
		double precision, intent(in)	E_rad,
		double precision, intent(in)	c0,
		double precision, intent(in)	c1
	)

Find the root of the 4th degree polynomial using the Halley method.

Compare error with dx = dx/dy dy

Definition at line 1604 of file mod_afld.t.

Here is the call graph for this function:

newton_method()

subroutine mod_afld::newton_method	(	double precision, intent(inout)	e_gas,
		double precision, intent(in)	E_rad,
		double precision, intent(in)	c0,
		double precision, intent(in)	c1
	)

Find the root of the 4th degree polynomial using the Newton-Ralphson method.

Compare error with dx = dx/dy dy

Definition at line 1570 of file mod_afld.t.

Here is the call graph for this function:

polynomial_bisection()

double precision function mod_afld::polynomial_bisection	(	double precision, intent(in)	e_gas,
		double precision, intent(in)	c0,
		double precision, intent(in)	c1
	)

Evaluate polynomial at argument e_gas.

Definition at line 1641 of file mod_afld.t.

put_diffterm_onegrid()

subroutine mod_afld::put_diffterm_onegrid	(	integer, intent(in)	ixI,
		integer, intent(in)	L,
		integer, intent(in)	ixO,
			L,
		double precision, dimension(ixi^s, 1:nw), intent(inout)	w
	)

inplace update of psa==>F_im(psa)

Definition at line 1233 of file mod_afld.t.

update_diffcoeff()

subroutine mod_afld::update_diffcoeff

(

type(state), dimension(max_blocks), target

psa

)

Definition at line 1306 of file mod_afld.t.

Here is the call graph for this function:

Variable Documentation

afld_diff_scheme

character(len=8) mod_afld::afld_diff_scheme = 'mg'

Which method to solve diffusion part.

Definition at line 53 of file mod_afld.t.

afld_mu

double precision, public mod_afld::afld_mu = 0.6d0

mean particle mass

Definition at line 22 of file mod_afld.t.

diff_coef_filter

logical mod_afld::diff_coef_filter = .false.

Take running average for Diffusion coefficient.

Definition at line 62 of file mod_afld.t.

diff_crit

double precision, public mod_afld::diff_crit

Number for splitting the diffusion module.

Definition at line 32 of file mod_afld.t.

diffcrash_resume

logical mod_afld::diffcrash_resume = .true.

Resume run when multigrid returns error.

Definition at line 73 of file mod_afld.t.

dt_diff

double precision mod_afld::dt_diff = 0.d0

running timestep for diffusion solver, initialised as zero

Definition at line 79 of file mod_afld.t.

fld_bisect_tol

double precision, public mod_afld::fld_bisect_tol = 1.d-4

Tolerance for bisection method for Energy sourceterms This is a percentage of the minimum of gas- and radiation energy.

Definition at line 26 of file mod_afld.t.

fld_diff_testcase

logical mod_afld::fld_diff_testcase = .false.

Set Diffusion coefficient to unity.

Definition at line 59 of file mod_afld.t.

fld_diff_tol

double precision, public mod_afld::fld_diff_tol = 1.d-4

Tolerance for adi method for radiative Energy diffusion.

Definition at line 29 of file mod_afld.t.

fld_eint_split

logical mod_afld::fld_eint_split = .false.

source split for energy interact and radforce:

Definition at line 15 of file mod_afld.t.

fld_fluxlimiter

character(len=16) mod_afld::fld_fluxlimiter = 'Pomraning'

Diffusion limit lambda = 0.33.

Definition at line 45 of file mod_afld.t.

fld_interaction_method

character(len=8) mod_afld::fld_interaction_method = 'Halley'

Which method to find the root for the energy interaction polynomial.

Definition at line 56 of file mod_afld.t.

fld_kappa0

double precision, public mod_afld::fld_kappa0 = 0.34d0

Opacity per unit of unit_density.

Definition at line 19 of file mod_afld.t.

fld_opacity_law

character(len=8), dimension(:), allocatable mod_afld::fld_opacity_law

Use constant Opacity?

Definition at line 41 of file mod_afld.t.

fld_opal_table

character(len=50) mod_afld::fld_opal_table = 'Y09800'

Definition at line 42 of file mod_afld.t.

fld_radforce_split

logical mod_afld::fld_radforce_split = .false.

Definition at line 16 of file mod_afld.t.

flux_lim_filter

logical mod_afld::flux_lim_filter = .false.

Take a running average over the fluxlimiter.

Definition at line 66 of file mod_afld.t.

i_diff_mg

integer, dimension(:), allocatable, public mod_afld::i_diff_mg

Index for Diffusion coeficients.

Definition at line 50 of file mod_afld.t.

i_test

integer, public mod_afld::i_test

Index for testvariable.

DELETE WHEN DONE

Definition at line 38 of file mod_afld.t.

lineforce_opacities

logical mod_afld::lineforce_opacities = .false.

Use or dont use lineforce opacities.

Definition at line 70 of file mod_afld.t.

size_d_filter

integer mod_afld::size_d_filter = 1

Definition at line 63 of file mod_afld.t.

size_l_filter

integer mod_afld::size_l_filter = 1

Definition at line 67 of file mod_afld.t.