Basic example: 1D diffusion equation
This example demonstrates how to solve a 1D radiative diffusion problem with pyFLD. Assuming constant density \(\rho\), opacity \(\kappa_\mathrm{R}=\kappa_\mathrm{P}=\kappa\), and flux limiter \(\lambda\), the coupled equation system reduces to a simpler form:
which, in steady state, further reduces to:
where \(E\) is the energy density and \(D\) is the diffusion coefficient. The solution to this equation is a straight line connecting the boundary values of \(E_\mathrm{r}\). This example will show how to set up and solve this problem with pyFLD and compare the results to the expected solution.
Note that to reach steady state, we will need to run the model for several diffusion timescales. For a box of size \(L\), the diffusion timescale is approximately \(t_\mathrm{diff} \sim L^2 / D\). In this example, with \(L = 1\,\mathrm{cm}\), \(\kappa = 1\,\mathrm{cm}^2/\mathrm{g}\), \(\rho = 1\,\mathrm{g}/\mathrm{cm}^3\), and \(\lambda = 1/3\), the diffusion coefficient is \(D = \lambda c / (\kappa \rho) \sim 10^{10}\,\mathrm{cm}^2/\mathrm{s}\), leading to a diffusion timescale of \(t_\mathrm{diff} \sim 10^{-10}\,\mathrm{s}\). Here, we will just run the model for an absurdly long time to ensure we reach steady state.
from astropy import units as u
import pyFLD as fld
import numpy as np
# grid setup
xi = fld.grid.get_domain_with_ghosts(1, 2, Nx=100) * u.cm
yi = fld.grid.get_domain_with_ghosts(0, 1, Nx=1) * u.cm
zi = fld.grid.get_domain_with_ghosts(0, 1, Nx=1) * u.cm
grid = fld.grid.Grid(xi, yi, zi, geometry='cartesian')
# fluid setup
Gas = fld.fluids.Gas(grid=grid)
E0 = 1 * u.erg / u.cm**3
rho = 1 * u.g/u.cm**3 * np.ones(grid.shape)
tmp = fld.tools.Er_to_T(E0) * np.ones(grid.shape)
Gas.set_density(rho)
Gas.set_temperature(tmp)
Gas.enroll_kappaR(lambda f: 1 * u.cm**2 / u.g)
Gas.enroll_kappaP(lambda f: 1 * u.cm**2 / u.g)
Gas.setup()
# radiative environment setup
Box = fld.FLD.RadiativeEnvironment([Gas], verbose=True, debug=False, constant_fluxlimiter=True)
# boundary conditions
val, ref = fld.grid.BOUNDARY_FIXEDVALUE, fld.grid.BOUNDARY_REFLECTIVE
bounds = [val, val, ref, ref, ref, ref]
vals = [E0, 2*E0, 0, 0, 0, 0]
Box.declare_boundaries(bounds, vals)
# time integration to steady state
for _ in range(5): Box.advance(dt=1e10 * u.s)
# visualization
import matplotlib.pyplot as plt
fig, ax = plt.subplots(ncols=2, figsize=(8, 3), dpi=120)
# compute expected solution
x = grid.x1[grid.act1].to_value('cm')
xstart, xend = grid.x1[grid.ibeg-1].to('cm'), grid.x1[grid.iend+1].to('cm')
Estart, Eend = Box.boundary_val[0], Box.boundary_val[1]
Eexp = Estart + (Eend-Estart)/(xend-xstart) * (x*u.cm - xstart)
Eexp = Eexp.to_value(E0)
Er_ = Box.Er[grid.active][0,0,:].to_value(E0) # [z=0, y=0, x=all] — only one active cell in y and z
ax[0].plot(x, Er_, label='pyFLD')
ax[0].plot(x, Eexp, label='expectation', ls=':', lw=2)
ax[1].plot(x, np.gradient(Er_, x, edge_order=2), label='pyFLD')
ax[1].plot(x, np.gradient(Eexp, x, edge_order=2), label='expectation', ls=':', lw=2)
ax[0].legend()
for axis in ax: axis.set_xlabel(r'$x$ [cm]')
ax[0].set_ylabel(r'$E_\mathrm{rad} / E_0$')
ax[1].set_ylabel(r'$\mathrm{d}E_\mathrm{rad}/\mathrm{d}x$ [$E_0$/cm]')
fig.tight_layout()
plt.show()