Nicolas Moens Module for including flux limited diffusion (FLD)-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	fld_params_read (files)
	public methods these are called in mod_rhd_phys

subroutine, public	fld_init (he_abundance, radiation_diffusion, energy_interact, r_gamma)
	Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variables to w-array, flux, kappa, eddington Tensor Lambda and R ...

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 the radiation force

subroutine, public	fld_radforce_get_dt (w, ixil, ixol, dtnew, dxd, x)

subroutine, public	fld_get_opacity (w, x, ixil, ixol, fld_kappa)
	Sets the opacity in the w-array by calling mod_opal_opacity.

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 stored in fld_fluxlimiter. It also calculates the ratio of radiation scaleheight and mean free path.

subroutine, public	fld_get_radflux (w, x, ixil, ixol, rad_flux, nth)
	Calculate Radiation Flux stores radiation flux in w-array.

subroutine	fld_get_eddington (w, x, ixil, ixol, eddington_tensor, nth)
	Calculate Eddington-tensor Stores Eddington-tensor in w-array.

subroutine, public	fld_get_radpress (w, x, ixil, ixol, rad_pressure, nth)
	Calculate Radiation Pressure Returns Radiation Pressure as tensor.

subroutine	fld_implicit_update (dtfactor, qdt, qtc, psa, psb)

subroutine	energy_interaction (w, x, ixil, ixol, dtfactor, qdt)
	Energy interaction and photon tiring.

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.

subroutine	evaluate_e_rad_mg (qtc, psa)
	inplace update of psa==>F_im(psa)

subroutine	put_diffterm_onegrid (ixil, ixol, w)
	inplace update of psa==>F_im(psa)

subroutine	fld_get_diffcoef_central (w, wct, x, ixil, ixol)
	Calculates cell-centered diffusion coefficient to be used in multigrid.

subroutine	update_diffcoeff (psa)

subroutine	fld_smooth_diffcoef (w, ixil, ixol)
	Use running average on Diffusion coefficient.

subroutine	bisection_method (e_gas, e_rad, c0, c1)
	Find the root of the 4th degree polynomial using the bisection method.

subroutine	newton_method (e_gas, e_rad, c0, c1)
	Find the root of the 4th degree polynomial using the Newton-Ralphson method.

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.

double precision function	dpolynomial_bisection (e_gas, c0, c1)
	Evaluate first derivative of polynomial at argument e_gas.

double precision function	ddpolynomial_bisection (e_gas, c0, c1)
	Evaluate second derivative of polynomial at argument e_gas.

Variables
logical	fld_eint_split = .false.
	source split for energy interact and radforce:

logical	fld_radforce_split = .false.

double precision, public	fld_kappa0 = 0.34d0
	Opacity per unit of unit_density.

double precision, public	fld_mu = 0.6d0
	mean particle mass

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.

double precision, public	fld_diff_tol = 1.d-4
	Tolerance for adi method for radiative Energy diffusion.

double precision, public	diff_crit
	Number for splitting the diffusion module.

character(len=8)	fld_opacity_law = 'const'
	Use constant Opacity?

character(len=50)	fld_opal_table = 'Y09800'

character(len=16)	fld_fluxlimiter = 'Pomraning'
	Diffusion limit lambda = 0.33.

integer	i_diff_mg
	diffusion coefficient for multigrid method

character(len=8)	fld_diff_scheme = 'mg'
	Which method to solve diffusion part.

character(len=8)	fld_interaction_method = 'Halley'
	Which method to find the root for the energy interaction polynomial.

logical	diff_coef_filter = .false.
	Take running average for Diffusion coefficient.

integer	size_d_filter = 1

logical	flux_lim_filter = .false.
	Take a running average over the fluxlimiter.

integer	size_l_filter = 1

logical	lineforce_opacities = .false.
	Use or dont use lineforce opacities.

logical	diffcrash_resume = .true.
	Resume run when multigrid returns error.

integer, dimension(:), allocatable, public	i_opf
	Index for Flux weighted opacities.

double precision	dt_diff = 0.d0
	running timestep for diffusion solver, initialised as zero

Detailed Description

Nicolas Moens Module for including flux limited diffusion (FLD)-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

◆ bisection_method()

subroutine mod_fld::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 789 of file mod_fld.t.

Here is the call graph for this function:

◆ ddpolynomial_bisection()

double precision function mod_fld::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 932 of file mod_fld.t.

◆ diffuse_e_rad_mg()

subroutine mod_fld::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.

Parameters

psa	Advance psa=psb+dtfactorqdtF_im(psa)

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

Reset dt_diff when diffusion worked out

Definition at line 593 of file mod_fld.t.

Here is the call graph for this function:

◆ dpolynomial_bisection()

double precision function mod_fld::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 922 of file mod_fld.t.

◆ energy_interaction()

subroutine mod_fld::energy_interaction	(	double precision, dimension(ixi^s,1:nw), intent(inout)	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, intent(in)	dtfactor,
		double precision, intent(in)	qdt
	)

Energy interaction and photon tiring.

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)

e_gas is the INTERNAL ENERGY without KINETIC ENERGY

Coefficients for the polynomial in Moens et al. 2022, eq 37.

Loop over every cell for rootfinding method

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 492 of file mod_fld.t.

Here is the call graph for this function:

◆ evaluate_e_rad_mg()

subroutine mod_fld::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 676 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_get_diffcoef_central()

subroutine mod_fld::fld_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

so that w(i_diff_mg) in ghostcells won't accumulate to extreme values

Definition at line 718 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_get_eddington()

subroutine mod_fld::fld_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(ixi^s,1:ndim,1:ndim), intent(out)	eddington_tensor,
		integer, intent(in)	nth
	)

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 402 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_get_fluxlimiter()

subroutine, public mod_fld::fld_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), intent(out)	fld_lambda,
		double precision, dimension(ixo^s), intent(out)	fld_r,
		integer, intent(in)	nth
	)

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.

optically thick limit

optically thin limit 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 274 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_get_opacity()

subroutine, public mod_fld::fld_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), 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

Definition at line 205 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_get_radflux()

subroutine, public mod_fld::fld_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,
		integer, intent(in)	nth
	)

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 375 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_get_radpress()

subroutine, public mod_fld::fld_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,
		integer, intent(in)	nth
	)

Calculate Radiation Pressure Returns Radiation Pressure as tensor.

Definition at line 455 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_implicit_update()

subroutine mod_fld::fld_implicit_update	(	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
	)

Definition at line 472 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_init()

subroutine, public mod_fld::fld_init	(	double precision, intent(in)	he_abundance,
		logical, intent(in)	radiation_diffusion,
		logical, intent(in)	energy_interact,
		double precision, intent(in)	r_gamma
	)

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

read par files

Set lineforce opacities as variable

need to check

Need mean molecular weight

set gamma

Read in opacity table if necesary

Definition at line 93 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_params_read()

subroutine mod_fld::fld_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 67 of file mod_fld.t.

◆ fld_radforce_get_dt()

subroutine, public mod_fld::fld_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 180 of file mod_fld.t.

Here is the call graph for this function:

◆ fld_smooth_diffcoef()

subroutine mod_fld::fld_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
	)

Use running average on Diffusion coefficient.

Definition at line 762 of file mod_fld.t.

◆ get_fld_rad_force()

subroutine, public mod_fld::get_fld_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)

Definition at line 144 of file mod_fld.t.

Here is the call graph for this function:

◆ halley_method()

subroutine mod_fld::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 881 of file mod_fld.t.

Here is the call graph for this function:

◆ newton_method()

subroutine mod_fld::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 851 of file mod_fld.t.

Here is the call graph for this function:

◆ polynomial_bisection()

double precision function mod_fld::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 912 of file mod_fld.t.

◆ put_diffterm_onegrid()

subroutine mod_fld::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 697 of file mod_fld.t.

◆ update_diffcoeff()

subroutine mod_fld::update_diffcoeff ( type(state), dimension(max_blocks), target psa )

Definition at line 746 of file mod_fld.t.

Here is the call graph for this function:

Variable Documentation

◆ diff_coef_filter

logical mod_fld::diff_coef_filter = .false.

Take running average for Diffusion coefficient.

Definition at line 39 of file mod_fld.t.

◆ diff_crit

double precision, public mod_fld::diff_crit

Number for splitting the diffusion module.

Definition at line 26 of file mod_fld.t.

◆ diffcrash_resume

logical mod_fld::diffcrash_resume = .true.

Resume run when multigrid returns error.

Definition at line 47 of file mod_fld.t.

◆ dt_diff

double precision mod_fld::dt_diff = 0.d0

running timestep for diffusion solver, initialised as zero

Definition at line 53 of file mod_fld.t.

◆ fld_bisect_tol

double precision, public mod_fld::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 22 of file mod_fld.t.

◆ fld_diff_scheme

character(len=8) mod_fld::fld_diff_scheme = 'mg'

Which method to solve diffusion part.

Definition at line 35 of file mod_fld.t.

◆ fld_diff_tol

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

Tolerance for adi method for radiative Energy diffusion.

Definition at line 24 of file mod_fld.t.

◆ fld_eint_split

logical mod_fld::fld_eint_split = .false.

source split for energy interact and radforce:

Definition at line 14 of file mod_fld.t.

◆ fld_fluxlimiter

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

Diffusion limit lambda = 0.33.

Definition at line 31 of file mod_fld.t.

◆ fld_interaction_method

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

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

Definition at line 37 of file mod_fld.t.

◆ fld_kappa0

double precision, public mod_fld::fld_kappa0 = 0.34d0

Opacity per unit of unit_density.

Definition at line 17 of file mod_fld.t.

◆ fld_mu

double precision, public mod_fld::fld_mu = 0.6d0

mean particle mass

Definition at line 19 of file mod_fld.t.

◆ fld_opacity_law

character(len=8) mod_fld::fld_opacity_law = 'const'

Use constant Opacity?

Definition at line 28 of file mod_fld.t.

◆ fld_opal_table

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

Definition at line 29 of file mod_fld.t.

◆ fld_radforce_split

logical mod_fld::fld_radforce_split = .false.

Definition at line 15 of file mod_fld.t.

◆ flux_lim_filter

logical mod_fld::flux_lim_filter = .false.

Take a running average over the fluxlimiter.

Definition at line 42 of file mod_fld.t.

◆ i_diff_mg

integer mod_fld::i_diff_mg

diffusion coefficient for multigrid method

Definition at line 33 of file mod_fld.t.

◆ i_opf

integer, dimension(:), allocatable, public mod_fld::i_opf

Index for Flux weighted opacities.

Definition at line 49 of file mod_fld.t.

◆ lineforce_opacities

logical mod_fld::lineforce_opacities = .false.

Use or dont use lineforce opacities.

Definition at line 45 of file mod_fld.t.

◆ size_d_filter

integer mod_fld::size_d_filter = 1

Definition at line 40 of file mod_fld.t.

◆ size_l_filter

integer mod_fld::size_l_filter = 1

Definition at line 43 of file mod_fld.t.

Functions/Subroutines

Variables

Detailed Description

Function/Subroutine Documentation

◆ bisection_method()

◆ ddpolynomial_bisection()

◆ diffuse_e_rad_mg()

◆ dpolynomial_bisection()

◆ energy_interaction()

◆ evaluate_e_rad_mg()

◆ fld_get_diffcoef_central()

◆ fld_get_eddington()

◆ fld_get_fluxlimiter()

◆ fld_get_opacity()

◆ fld_get_radflux()

◆ fld_get_radpress()

◆ fld_implicit_update()

◆ fld_init()

◆ fld_params_read()

◆ fld_radforce_get_dt()

◆ fld_smooth_diffcoef()

◆ get_fld_rad_force()

◆ halley_method()

◆ newton_method()

◆ polynomial_bisection()

◆ put_diffterm_onegrid()

◆ update_diffcoeff()

Variable Documentation

◆ diff_coef_filter

◆ diff_crit

◆ diffcrash_resume

◆ dt_diff

◆ fld_bisect_tol

◆ fld_diff_scheme

◆ fld_diff_tol

◆ fld_eint_split

◆ fld_fluxlimiter

◆ fld_interaction_method

◆ fld_kappa0

◆ fld_mu

◆ fld_opacity_law

◆ fld_opal_table

◆ fld_radforce_split

◆ flux_lim_filter

◆ i_diff_mg

◆ i_opf

◆ lineforce_opacities

◆ size_d_filter

◆ size_l_filter