Merge pull request #4 from orlox/LABS

PR 14th november last

docs/make.jl

@@ -28,9 +28,8 @@ makedocs(sitename="Stellar Structure and Evolution",
                                               "Exercises"=>"6_convection_problems.md"],
                 "Nucleosynthesis (07/11/23)" => ["Notes" => "7_nucleo1.md",
                                               "Exercises"=>"7_nucleo1_problems.md"],
-                "LAB1 (14/11/2023)" => [
+                "LAB1 (14/11/2023)" => ["Stellar evolution models" => "8_SE_codes.md",
                                         "Instructions" => "8_instructions.md",
-                                        "Stellar evolution models" => "8_SE_codes.md",
                                         "Questions"=>"8_questions.md"],
                 ])
 

docs/src/8_SE_codes.md

@@ -14,13 +14,13 @@ Let's collect and summarize the system of differential equations derived in the
 
 3. **ENERGY EQUATION**
 
-   $T\dfrac{\partial s}{\partial t}=-\dfrac{\partial L}{\partial m}+\epsilon_{\mathrm{nuc}}\hspace{8.1cm}$  
+   $\dfrac{\partial L}{\partial m}=\epsilon_{\mathrm{nuc}}-c_P\dfrac{\partial T}{\partial t}+\dfrac{\delta}{\rho}\dfrac{\partial P}{\partial t}\hspace{8.1cm}$  
 
 4. **ENERGY TRANSPORT EQUATION**
 
    $\dfrac{\partial T}{\partial m}=-\dfrac{Gm}{4\pi r^4}\dfrac{T}{P}\nabla\hspace{6.56cm}$    
 
-   with $\hspace{0.25cm}\nabla=\nabla_{\mathrm{rad}}\hspace{0.25cm}\mathrm{if}\hspace{0.25cm}\nabla_{\mathrm{rad}}\leq\nabla_{\mathrm{ad}}\:, \nabla=\nabla_{\mathrm{ad}}+\Delta\nabla\hspace{0.25cm}\mathrm{if}\hspace{0.25cm}\nabla_{\mathrm{rad}}>\nabla_{\mathrm{ad}}$
+   with $\hspace{0.25cm}\nabla=\nabla_{\mathrm{rad}}\hspace{0.25cm}\mathrm{if}\hspace{0.25cm}\nabla_{\mathrm{rad}}\leq\nabla_{\mathrm{ad}}\:, \nabla=\nabla_{\mathrm{ad}}\hspace{0.25cm}\mathrm{if}\hspace{0.25cm}\nabla_{\mathrm{rad}}>\nabla_{\mathrm{ad}}$
 
 5. **ABUNDANCE EQUATION**
 
@@ -53,17 +53,19 @@ with $Q_{ij}$ being the $Q$-value of the reaction between species $i$ with $j$,
 $$\mathrm{EOS}\hspace{1cm}\rho=\rho(P,T,X_i)$$
 
 #### Assumed microphysics ingredients
-These are again characteristics of the stellar material, which are describing the physics we put in our model to reach the desired degree of complexity: the Rosseland mean opacity $\kappa$, the nuclear reaction rates $r_{i,j}$ for the chosen network of isotopes, the thermodynamic gradients $\nabla_{\mathrm{rad}}$ and $\nabla_{\mathrm{ad}}$, the specific entropy $s$. As studied in the lectures, we assume that these are known as a function of $(P,T,X_i)$. Summarizing:
+These are again characteristics of the stellar material, which are describing the physics we put in our model to reach the desired degree of complexity: the Rosseland mean opacity $\kappa$, the nuclear reaction rates $r_{i,j}$ for the chosen network of isotopes, the specific entropy $s$. As studied in the lectures, we assume that these are known as a function of $(P,T,X_i)$:
 
 $$\kappa=\kappa(P,T,X_i)$$
 $$r_{jk}=r_{jk}(P,T,X_i)\hspace{0.5cm}\Rightarrow\hspace{0.5cm}\epsilon_{\mathrm{nuc}}=\epsilon_{\mathrm{nuc}}(P,T,X_i)$$
-$$\nabla_{\mathrm{rad}}=\nabla_{\mathrm{rad}}(P,T,X_i)\:,\hspace{1cm}\nabla_{\mathrm{ad}}=\nabla_{\mathrm{ad}}(P,T,X_i)$$
-$$s=s(P,T,X_i)$$
+
+Lastly, the thermodynamic gradients $\nabla_{\mathrm{rad}}$ and $\nabla_{\mathrm{ad}}$ are also assumed to be known functions of $(P,T,X_i)$, once the EOS is known:
+
+$$\nabla_{\mathrm{rad}}=\nabla_{\mathrm{rad}}(P,T,X_i,)\:,\hspace{1cm}\nabla_{\mathrm{ad}}=\nabla_{\mathrm{ad}}(P,T,X_i)$$
 
 If you now look at the symbols that are left, you can count the following: $r(m,t)$, $P(t,m)$, $T(m,t)$, $L(m,t)$, to be added to the $I$ mass fractions $X_i(m,t)$. We have therefore $I+4$ unknowns for $I+4$ equations.
 
 #### Initial conditions $t=t_0$
-By looking at the quantities that show time derivatives in the system of equations, we see that the radius $r(m,t)$ appears as first and second derivative with respect to time; specific entropy $s(m,t)$ and mass fractions $X_i(m,t)$ appear as first derivatives with respect to time. Therefore, we need to know the initial values of $r(m,t_0)$, $\dot{r}(m,t_0)$, $s(m,t_0)$, $X_i(m,t_0)$, $\forall\: m$.
+By looking at the quantities that show time derivatives in the system of equations, we see that the radius $r(m,t)$ appears as first and second derivative with respect to time; specific entropy $s(m,t)$ (if you recast the time derivatives of $P$ and $T$ into -$T\partial s/\partial t$) and mass fractions $X_i(m,t)$ appear as first derivatives with respect to time. Therefore, we need to know the initial values of $r(m,t_0)$, $\dot{r}(m,t_0)$, $s(m,t_0)$, $X_i(m,t_0)$, $\forall\: m$.
 
 ### What is a stellar model?
 The $I+4$ non-linear partial differential equations must be solved simultaneously as a function of the two independent variables $(m,t)$, namely a solution to the above set of equations will span the interval $0\leq m\leq M$, with $M$ depending on mass loss, and $t\geq t_0$, with $t_0$ being an initial time. Notice that we will only find physically relevant solutions once we specify
@@ -73,88 +75,16 @@ The $I+4$ non-linear partial differential equations must be solved simultaneousl
 
 **DEF**: A solution $r(m,t)$, $P(t,m)$, $T(m,t)$, $L(m,t)$, $X_i(m,t)$ for which all the above is satisfied is called **a stellar model**.
 
-## Boundary conditions
-TODO
-### Center $m=0$
-
-### Surface $m=M$
-
-## Henyey scheme
-For realistic functions no analytic solutions are possible, so that one should search for numerical solutions to solve the system of equations. Here we shall describe the ***Henyey numerical scheme***, a simplified method to solve differential equations with boundary conditions at both ends of the solution interval. The method is based on 1) discretizing the equations, subsequently 2) linearizing them and lastly 3) find a solution with iterations. The trial solution for the whole interval is gradually improved upon in consecutive iterations until the required degree of accuracy is reached, in a Newton-Raphson method fashion.
-
-#### Standard models  
-In absence of mixing terms (i.e. $\partial X_i/\partial m=0$), we can first solve the system of equations 1-4 with a given $X_i(m)$, then update $X_i$ by applying equation 5 with a small timestep $\Delta t$. Then again solve 1-4 with the updated $X_i$, and so on. For illustrative reasons, we will limit ourselves to use this assumption and focus on how to solve the equations 1-4 simultaneously. Additionally, we will assume full hydrostatic and thermal equilibrium:
-* *Hydrostatic Equilibrium* $\hspace{1cm}\ddot{r}=0\hspace{0.5cm}\Rightarrow\hspace{0.5cm}0=-\dfrac{Gm}{r^2}-4\pi r^2\dfrac{\partial P}{\partial m}$
-
-* *Thermal Equilibirum* $\hspace{1cm}\dot{s}=\dot{T}=\dot{P}=0\hspace{0.5cm}\Rightarrow\hspace{0.5cm}0=-\dfrac{\partial L}{\partial m}+\epsilon_{\mathrm{nuc}}$
-
-At each given $X_i(m)$, we need to solve for 4 unknown functions $r$, $P$, $T$, $L$ in the interval $0\leq m\leq M$. And as for initial conditions, we only need to know $s(m,t_0)$ (actually just $\nabla(m,t_0)$), and then $\dot{s}=0$.
-
-The method under these assumptions is particularly suitable to compute stellar models during the core-hydrogen burning stage, and the outcome gives us what we call **standard models**. In reality more complex, state-of-the-art stellar evolution codes, like MESA, solve the full set with mixing terms, and are able to treat all the stages of the evolution.
-
-#### From PDEs to ODEs
-With the above simplifications, we got rid of all the time derivatives in the equations 1-4. Therefore, we are not dealing anymore with a system of partial differential equations, but with a system of **ordinary** differential equations. AS a result, we can write the 4 equations in this way
-
-$$\dfrac{dy_i}{dm}=f_i(y_1,y_2,y_3,y_4)\hspace{1cm}\mathrm{with}\hspace{0.5cm}i=1,2,3,4$$
-where we have called $y_1=r$, $y_2=P$, $y_3=T$ and $y_4=L$.
-
 ### Discretization
-Let us now discretize the equations 1-4: this is just replacing them by differential equations for a finite mass-interval $\Delta m=[m^j,\:m^{j+1}]$. At each end of the finite mass interval $\Delta m$, the variables $r$, $P$, $T$, $L$ will assume the values $y_1^j,\: y_1^{j+1}$, ... $y_4^j,\: y_4^{j+1}$. The functions $f_i$ have to be evaluated in an average argument, which is indicated with $y_i^{j+1/2}$, often being just a geometric mean between $y_i^{j}$ and $y_i^{j+1}$. 
-We can now write our discretized system 1-4 of ordinary differential equations in this way:
-
-$$A_i^j=0$$
-
-$$\mathrm{with}\hspace{0.5cm}A_i^j\equiv \dfrac{y_i^j-y_i^{j+1}}{m^j-m^{j+1}}-f_i(y_1^{j+1/2},\:...\:,y_4^{j+1/2})\:,\hspace{1cm}i=1,2,3,4\:.$$
-
-#### Surface $j=1$
-Let us write in the same fashion also the boundary conditions at $m=0$ and $m=M$. If the boundary conditions for $P$ and $T$ at the surface are fixed at the mass coordinate $m^1$, namely at index $j=1$, we can write
-
-$$B_i=0$$
-
-$$\mathrm{with}\hspace{0.5cm}B_1\equiv y_2^1-\pi(y_1^1,y_4^1)\:,\hspace{0.5cm}B_2\equiv y_3^1-\theta(y_1^1,y_4^1)\:,\hspace{1cm}i=1,2\:.$$
-
-#### Core $j=K$
-Lastly, we need to impose the boundary conditions in the center $m=0$. Notice that if we want to know the solution $r(m)$, $P(m)$, $T(m)$, $L(m)$ in $K$ grid points, we are diving the star in $K-1$ intervals: every discretization point will have an $m^j$ with $j=1,\:...\: K$. The upper boundary is located at $j=1$, and the core at $j=K$. We can write
-
-$$C_i\left(y_1^{K-1},\: ... \:,y_4^{K-1},y_2^K,y_3^K\right)=0\:,\hspace{1cm}i=1,2,3,4\:,$$
-
-where we have made use of the fact that $P_c\simeq P+C_2(\rho_c=\rho_c(r))\Leftrightarrow y_2^K=y_2^{K-1}+C_2(y_1^{K-1})$ and $T_c\simeq T+C_3(\kappa_c,\epsilon_{\mathrm{nuc},c},\rho_c)\Leftrightarrow y_3^{K}=y_3^{K-1}+C_3(y_1^{K-1},\: ... \:,y_4^{K-1})$. The conditions $y_1^K=0$ and $y_4^K=0$ are incorporated. 
-
-#### The complete discretized system
-1. $A_i^j=0\:,\hspace{0.5cm}i=1,\:...\:,4\:,\hspace{0.25cm}j=1,\:...\:,K-2$
-2. $B_i=0\:,\hspace{0.5cm}i=1,2$
-3. $C_i=0\:,\hspace{0.5cm}i=1,\:...\:,4$
-
-and you have $4(K-2)+2+4=4K-2$ equations, which correspond to $4K-2$ unknowns (as we want to know the values of $y_i^j$ at every $j=1,\:...\:,K$, with $y_1^K=y_4^K=0$).
+TODO
 
 #### Convergence
 These difference equations represent an **approximation** of the full differential equations, the accuracy of which can be improved by reducing $\Delta t$ and $\Delta m$. A good approach is to choose $m_j$ for each $j$ such, that all variables change by less than a predefined upper limit between points $j$ and $j-1$. *Spatial convergence* of the stellar model is achieved when the change in variables is reduced such that the numerical solution no longer depends on the $\Delta m$. Similarly, *temporal convergence* is achieved when the change in the variables is reduced till independence on $\Delta t$.
 
 ### Linearization
-Starting from the above system 1-3, the solving procedure is now very similar to a Newton-Raphson algorithm in linear algebra. Let us start from a guessed initial solution $\left(y_i^j\right)_{1}$ for $i=1,2,3,4$ and $j=1,\:...\:,K$. This first guess won't be a perfect solution for the system 1-3, i.e. when we use them as arguments in the functions $A_i^j$, $B_i$ and $C_i$ we will find
-
-$$A_i^j(1)\neq 0\:,\hspace{0.25cm} B_i(1)\neq0\:,\hspace{0.25cm} C_i(1)\neq 0\:.$$
-
-Let us now look for corrections $\delta y_i^j$ for all variables at all mesh points such that the second approximation $\left(y_i^j\right)_{2}$ will make the functions $A_i^j(2)$, $B_i(2)$ and $C_i(2)$ disappear, i.e.
-
-$$A_i^j(1)+\delta A_i^j=0\:,\hspace{0.25cm}B_i(1)+\delta B_i=0\:,\hspace{0.25cm}C_i(1)+\delta C_i=0\:,$$
-
-$$\mathrm{with}\hspace{0.5cm}\left(y_i^j\right)_{2}=\left(y_i^j\right)_{1}+\delta y_i^j\:.$$
-
-For small enough corrections, we may expand the $\delta A_i^j$, $\delta B_i$ and $\delta C_i$ in terms of powers of the corrections $\delta y_i^j$, and keep only the linear terms in this expansion; for example,
-
-$$\delta B_1\approx \left.\dfrac{\partial B_1}{\partial y_1^1}\right|_1\delta y_1^1+\left.\dfrac{\partial B_1}{\partial y_2^1}\right|_1\delta y_2^1+\left.\dfrac{\partial B_1}{\partial y_3^1}\right|_1\delta y_3^1+\left.\dfrac{\partial B_1}{\partial y_4^1}\right|_1\delta y_4^1\:,$$
-where the partial derivatives are to be evaluated at the first approximation arguments $\left(y_i^j\right)_1$. We will have in the end a system that looks like this:
+TODO
 
-1. $\left.\dfrac{\partial A_i^j}{\partial y_1^j}\right|_1\delta y_1^j+...+\left.\dfrac{\partial A_i^j}{\partial y_4^j}\right|_1\delta y_4^j+\left.\dfrac{\partial A_i^j}{\partial y_1^{j+1}}\right|_1\delta y_1^{j+1}+...+\left.\dfrac{\partial A_i^j}{\partial y_4^{j+1}}\right|_1\delta y_4^{j+1}=\left.-A_i^j\right|_1$
-2. $\left.\dfrac{\partial B_i}{\partial y_1^1}\right|_1\delta y_1^1+...+\left.\dfrac{\partial B_i}{\partial y_4^1}\right|_1\delta y_4^1=\left.-B_i\right|_1$
-3. $\left.\dfrac{\partial C_i^{K-1}}{\partial y_1^{K-1}}\right|_1\delta y_1^{K-1}+...+\left.\dfrac{\partial C_i^{K-1}}{\partial y_4^{K-1}}\right|_1\delta y_4^{K-1}+\left.\dfrac{\partial C_i^{K}}{\partial y_2^{K}}\right|_1\delta y_2^{K}+\left.\dfrac{\partial C_i^{K}}{\partial y_3^{K}}\right|_1\delta y_3^{K}=\left.-C_i\right|_1$
-
-Notice that we know all the partial derivatives and the RHS of the equations, since they are evaluated at $\left(y_i^j\right)_1$. This is a **linear** (inhomogeneous) system that needs to be solved for the small corrections $\delta y_i^j$. The matrix containing all the partial derivatives is called *Henyey matrix*, and when its determinant is not zero, the system can be solved.
+**NB** There's no guarantee for the iteration procedure for improving the approximations  to always converge! In fact if the chosen approximation is too far from the solution, then the required corrections are so large that one cannot neglect the second-order terms, and the linearization is not appropriate anymore.
 
-### Iterations
-When the determinant of the Henyey matrix is not zero, we can obtain the corrections $\delta y_i^j$ and use them to compute a better approximation of the solution, i.e. $\left(y_i^j\right)_2=\left(y_i^j\right)_1+\delta y_i^j$. When using these second approximations as arguments, we will generally still find $A_i^j(2)\neq 0$, $B_i(2)\neq0$ and $C_i(2)\neq 0$. This is because the corrections were calculated from the linearized equations 1-3, while the full discretized equations are non-linear, and in any case the accuracy is limited, as discussed when talking about convergence.
 
-In general, we will need more consecutive iterations. The approximate solution can be improved until either the absolute values of all corrections $\delta y_i^j$, or the absolute values of all right-hand sides in 1-3, drop below a chosen limit, i.e. the required accuracy.
 
-**NB** There's no guarantee for the iteration procedure for improving the approximations  to always converge! In fact if the chosen approximation is too far from the solution, then the required corrections are so large that one cannot neglect the second-order terms, and the linearization is not appropriate anymore.

docs/src/8_instructions.md

@@ -1,14 +1,14 @@
 # Instructions
 
-Modules for Experiments in Stellar Astrophysics (MESA) [^Jermyn2022] is a state of the art, open-source 1D stellar evolution code. MESA is built to allow users to run experiments in stellar evolution, with a huge variety of possibilities for input physics and customization. You are strongly encouraged to give a look at the [documentation](https://docs.mesastar.org/en/release-r21.12.1/) to discover more about the software. We will specifically use the MESAstar module to evolve a single star.
+Modules for Experiments in Stellar Astrophysics (MESA) [^Paxton2011] is a state of the art, open-source 1D stellar evolution code. MESA is built to allow users to run experiments in stellar evolution, with a huge variety of possibilities for input physics and customization. You are strongly encouraged to give a look at the [documentation](https://docs.mesastar.org/en/release-r21.12.1/) to discover more about the software. We will specifically use the MESAstar module to evolve a single star.
 
 ## Lab intro
 In this lab you will learn how to evolve a star with MESA and how to interpret its outputs in terms of the theory of stellar structure and evolution you have seen so far. In particular
 
 - **FIRST HOUR**: see [Stellar evolution models](https://orlox.github.io/stars_2023_2024/dev/). Theoretical overview of the software. The basic numerical scheme used to solve the set of stellar structure and evolution equations; the boundary conditions for the differential equations; the software infrastructure and usage of inlists, output files.
 
 - **SECOND HOUR**: set up your simulation and save the outputs in your local machines / USB memory / personal folder in the system, as you will need to use it at home to post-process the information and write a report. This time at the lab is mainly meant for you to run the simulation while you have someone that can help you (me); if, at the end, you still have some time, you can start answering the questions of the report.  
-  NB: you will have to run **two** simulations of two stars with different masses, and save **both** those outputs. Be sure to do that before you leave the session.
+  NB: you will have to run **two** simulations of two stars with different masses, and save **both** those outputs. Be sure to do that before you leave the session. 
 
 
 ## The report
@@ -140,9 +140,9 @@ teff = matrix["log_Teff"]
 ## The assigned masses
 Due to the number of available computers in the lab, you all will be paired during the session. This means that every student (A) will have a mate (B) to look at the simulations with. Each team (student A + student B) of students will have to produce, **during the session**, a total of two simulations (simulation 1 + simulation 2), and save the outputs for later use (at home). These two simulations will serve the following purposes:
 
-- **Simulation 1**  will evolve a star of initial mass $M$. Student A will have to answer all the questions aside from **Massive vs Low Mass evolution** by taking into account exclusively the outcome of simulation 1. Student B will make use of the outcome of simulation 1 just to answer the aforesaid question.
+- **Simulation 1**  will evolve a star of initial mass $M$. Student A will have to answer all the questions, aside from **Massive vs Low Mass evolution** and a small part of the question on **Variations of the EoS: Radiation Pressure and Degeneracy**, by taking into account exclusively the outcome of simulation 1. Student B will make use of the outcome of simulation 1 just to answer the aforesaid question.
 
-- **Simulation 2** will evolve a star of initial mass $M'\neq M$. Student B will have to answer all the questions aside from  **Massive vs Low Mass evolution** by taking into account exclusively the outcome of simulation 2. Student A will make use of the outcome of simulation 2 just to answer the aforesaid question.
+- **Simulation 2** will evolve a star of initial mass $M'\neq M$. Student B will have to answer all the questions, aside from  **Massive vs Low Mass evolution** and a small part of the question on **Variations of the EoS: Radiation Pressure and Degeneracy**, by taking into account exclusively the outcome of simulation 2. Student A will make use of the outcome of simulation 2 just to answer the aforesaid question.
 
 Please find below the list of pairs of masses, together with the team (student A + student B) number. I'll give you the team number once you're paired in the lab and make sure that no pairs of masses is picked by two different teams.
 
@@ -159,6 +159,6 @@ Please find below the list of pairs of masses, together with the team (student A
 |    **9**    |         8        |         0.8         |
 
 
-[^Jermyn2022]:
-    The Astrophysical Journal Supplement Series, Volume 265, Issue 1, id.15, 38 pp.
-    https://ui.adsabs.harvard.edu/abs/2023ApJS..265...15J/abstract
+[^Paxton2011]:
+    The Astrophysical Journal Supplement, Volume 192, Issue 1, article id. 3, 35 pp. (2011).
+    https://ui.adsabs.harvard.edu/abs/2011ApJS..192....3P/abstract