Our goal is to simulate a isothermal sphere undergoing gravitational collapse. We can write down the force balance equation for the sphere, which describes the competition between gravitational infall versus thermal pressure support.
$$ -\frac{1}{\rho}\nabla P -\nabla \Phi_g =0$$
Using Poisson's equation ($\nabla \Phi_g^2 =4\pi G \rho$) and the isothermal equation of state ($P = \rho a_T^2$), we can rearrange the equation to the form of the isothermal Lane-Emden equation: 


$$\frac{1}{\psi^2} \frac{d}{d\psi}\Bigg(\xi^2 \frac{d\psi}{d\xi}\Bigg)=e^{-\psi}$$
Using product rule, we can rewrite the Lane-Emden equation as: 
\begin{align}
\frac{d^2\psi}{d\xi^2}+\frac{2}{\xi}\frac{d\psi}{d\xi}=e^{-\psi}
\end{align}
with boundary conditions as: $\psi(0)=0$ and $\frac{d\psi}{d\xi}(0)=0$.

The Lane-Emden equation describes a family of solution based on the single parameter $\xi_{max}$, which specifies the cloud radius and determine the density contrast $\rho_c/\rho_0$ (ratio between central density and edge of the sphere). Analytically, the sphere is marginally stable (i.e. just at the verge of collpasing) when $\rho_c/\rho_0$=14.1 [1]. This marginally-stable state is called the Bonner-Ebert sphere. Numerically, we can see that a cloud radius of  $\xi$=6.451 corresponds to a density contrast of $\rho_c/\rho_0$=14.038. The singular isothermal sphere case corresponds to $\xi_{max}=\infty$. [1]

We want to solve the Lane-Emden equation to figure out what is the density distribution that yields a marginally-stable sphere, but the Lane-Emden solution can not be solved in closed form, so instead we need to do numerical intergration to get the solution. 


```python
%pylab inline
```

    Populating the interactive namespace from numpy and matplotlib


We start integrating at a small poisitive number (1e-6) to avoid div-by-0 error and integrate up to xi_max, with interval of $\Delta \xi$ = 0.01.


```python
xi_max = 6.451
from scipy import integrate
def solvr(Y, t):
    return [Y[1], exp(-Y[0])-2/t*Y[1]]
xi = np.arange(1e-6, xi_max, 0.01)
IC = [0, 0] # initial conditions for psi and psi prime
psi = integrate.odeint(solvr, IC, xi)
```

For a spherically symmetric cloud, $\rho(r)=\rho_c e^{-\psi}$, so we can use the numerically integrated results for $\psi$ to compute the density profile.


```python
rho = np.exp(-psi[:,0]) 
plt.plot(xi,psi[:,0],label ='$\psi$')
plt.plot(xi,rho,label =r'$\rho/\rho_c$') 
plt.xlabel(r"$\xi$",fontsize=15)
plt.xlim(0,xi_max)
from matplotlib.legend_handler import HandlerLine2D
plt.legend(loc='upper left',prop={'size':13},numpoints=1)
```

Note that since ``rmax`` is the only parameter that needs to be adjusted for different models of the Lane-Emden solution. To make sure that the boundary of the sphere is properly filled with the Bonnert-Ebert values, we extend the density profile in ``density.txt`` to $\xi$= 20. 


```python
xi = np.arange(1e-6, 20, 0.01) 
psi = integrate.odeint(solvr, IC, xi)
rho = np.exp(-psi[:,0]) 
np.savetxt("FLASH_4.3/density.txt",rho)
```

Inside ``Simulation_initBlock.F90``, we can read in the values of the array 
~~~fortran
  open(12,file="../density.txt")
  read(12,*,IOSTAT=io) dens_arr
  close(12)
~~~

Now we have a 1D array of non-dimensional density values from numerical integration, we need to assign each cell in the simulation box with the density values that corresponds to the right $\rho(r)$. Since we are working with a radial distribution on a cartesian grid, we don't have the exact radius value at each cell, so need to do linear interpolation to map the cells to somewhere between their closest density values. As explained by [Wikipedia](https://en.wikipedia.org/wiki/Linear_interpolation), linear interpolation basically weighs what the y values are  based on how close the x values lie to the previous and subseqeunt datapoints. In our case, 

$$ \rho(r) = \rho_0 + (\rho_1-\rho_0)\frac{r-r_0}{r_1-r_0}$$

Translating this into our Simulation_init.F90 and picking dr ($\Delta \xi$) as 0.01 to correspond with np.arange used in the numerical integration: 

    if (rc <= rmax) then
       rc0 = int(rc/dr)+1 !to prevent from hitting index 0 which is yields zero density
       rho0 =  rho_c*dens_arr(rc0,1)
       rho1 = rho_c*dens_arr(rc0+1,1)
       rc0 = rc0*dr
       rc1 = rc0+dr
       rhoZone = fattening_factor*(rho0+(rho1-rho0)*(rc-rc0)/(rc1-rc0) )!linear interpolation

Simulation_init.F90 reads in values up to ``rmax`` $\xi$= 16.90 for initialize the sphere, note that this has a much higher density contrast than the Bonner-Ebert sphere. This ensures that the sphere in our toy problem will actually collapse in a short time. We also multiply the sphere's density by a "fattening_factor" of 100 for the same reason.


```python
idx = np.where(np.isclose(xi,16.9))[0][0]
print "At xi = ", xi[idx]
print "Density is: ", rho[idx]*rho_c, "g/cm^3"
print "density contrast is ", 1./rho[idx]
```

    At xi =  16.900001
    Density is:  7.2918475616e-20 g/cm^3
    density contrast is  150.853400418


We can now run our simulation to check that the sphere is initialized with the right density distribution.


```python
cd ~/dlee/FLASH4.3/object/
```

    /global/u2/d/dorislee/dlee/FLASH4.3/object



```python
import yt
from matplotlib.legend_handler import HandlerLine2D
from mpl_toolkits.axes_grid1 import AxesGrid
from yt.mods import *
yt.mylog.setLevel(50)
```


```python
def compareplot1Dprofile(timestep):
    fattening_factor=100
    rho_c = 1.1e-19*fattening_factor
    pf= yt.load("sphere_hdf5_chk_{}".format(str(timestep).zfill(4)))
    sp = pf.sphere(pf.domain_center, (0.32,"pc"))
    rp = yt.create_profile(sp,'radius','density')
    plt.plot(rp.x.value*1.05e-17,rp["density"].in_units("g/cm**3").value/rho_c,label="Simulation")
    xi_max = 6.451
    def solvr(Y, t):
        return [Y[1], exp(-Y[0])-2/t*Y[1]]
    xi = np.arange(1e-6, xi_max, 0.01) #start at small poisitive number to avoid div-by-0
    psi = integrate.odeint(solvr, [0, 0], xi)
    rho = np.exp(-psi[:,0])  
    plt.plot(xi,rho,"--",label ='Integration',color="red",linewidth=5)
    plt.xlabel(r"$\xi$",fontsize=15)
    plt.ylabel(r"$\rho/\rho_c$",fontsize=15)
    plt.legend(loc='upper right',prop={'size':12},numpoints=1)
```


```python
compareplot1Dprofile(0)
```

## Reference: 

[1] Stahler, Steven William., and F. Palla. *The Formation of Stars*. pg. 243-248. Weinheim: Wiley-VCH, 2004.  
[2] Typanski,Nathan, *Solving a second-order ODE with NumPy and SciPy*. Online. [Link](http://www.nathantypanski.com/blog/2014-08-23-ode-solver-py.html)

e67a4071f51dbd0328a050745a36b869422838c5
ipynb
Jupyter Notebook
tutorial/3-Lane Emden Numerical Solution.ipynb
dorisjlee/astroSim-tutorial
0f0b07b4c4021c6e9f166b51e7dc8c3a1836920d
BSD-3-Clause
2016-07-24T17:28:34.000Z
2022-03-03T14:08:18.000Z
tutorial/3-Lane Emden Numerical Solution.ipynb
dorisjlee/astroSim-tutorial
0f0b07b4c4021c6e9f166b51e7dc8c3a1836920d
BSD-3-Clause
tutorial/3-Lane Emden Numerical Solution.ipynb
dorisjlee/astroSim-tutorial
0f0b07b4c4021c6e9f166b51e7dc8c3a1836920d
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES