overview.dox

/** \page overview Overview

ESTER uses Newton's method to solve the set of stellar PDEs.

Solving the PDE \f$ L(u) = 0 \tag{1} \f$ with Newton's method is done
by refining the solution (last iterate) \f$ u_n\f$, by first solving:

\f$ 
J_L(\delta u) = -L(u_n) \tag{2}
\f$
for \f$\delta u\f$ and then updating the solution with \f$ u = u_n +
\delta u \f$

In order to write equations in a simple way, the ESTER code uses a
specific formalism, which can be easily understood with an example.

Let us take Poisson equation :

\f$ \Delta \phi = \pi_c \rho,
\qquad {\rm with}\quad  \pi_c = \frac{4 \pi G \rho_c^2 R^2}{p_c} \tag{3}\f$

The evaluation of the Newton increment demands the solution of :
\f$
\Delta \delta\phi - \rho\delta\pi_c  - \pi_c\delta\rho = -(\Delta \phi -
\pi_c \rho)_n
\tag{4}\f$




This is equivalent to the code:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.py}
lap_phi=lap(phi);    
lap_phi.add(op, "eq_Phi", "Phi");
lap_phi.add(op, "eq_Phi", "r");

op->add_d("eq_Phi", "rho", -pi_c*ones(nr, nth));
op->add_d("eq_Phi", "pi_c", -rho);

rhs1 = -lap_phi.eval()+pi_c*rho;
rhs2=-lap_phi.eval();

rhs=zeros(nr+nex,nth);
rhs.setblock(0,nr-1,0,-1,rhs1);
rhs.setblock(nr,nr+nex-1,0,-1,rhs2);

op->set_rhs("eq_Phi",rhs);

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Line 1 defines the symbolic expression for \f$ \Delta \phi \f$ (see
\ref symbolic).

Line 2 says that the equation associated with lap_phi is tagged "eq_Phi"
and that it depends on the variable  \f$\phi\f$.

Line 3 says that equation "eq_Phi" also depends on \f$ r \f$.

Lines 5 and 6 add the remaining terms to the Jacobian matrix. We see
the coefficient \f$ -\pi_c \f$ of \f$\delta\rho \f$ and \f$-\rho\f$
of \f$\delta\pi_c\f$. This completes the definition of the left-hand
side of the equation.

Line 8 prepares the RHS of the equation for the domains within the star.

Line 9 prepares the RHS of the equation for the domain outside the star.

Line 11 initializes the RHS-matrix

line 12 inserts rhs1 into the block describing radial points from 0 to
nr-1, and for all latitudinal points "0,-1" (-1 is C++ coding for the
last point).

Line 13 inserts rhs2 into the block for the outer (vacuum) domain.

Line 15 inserts the RHS into equation "eq_Phi"

\section var_sec Variables in ESTER

Variable            | Name in the code | Description 
------------------- | ---------------- | -----------
\f$ \phi \f$        | Phi              | gravitational potential
\f$ p \f$           | p                | pressure
\f$ \log p \f$      | log_p            | log of the pressure
\f$ \pi_c \f$       | \pi_c            | \f$\frac{4\pi G \rho^2_c R^2}{p_c}\f$
\f$ T \f$           | T                | temperature
\f$ \log T \f$      | log_T            | log of the temperature
\f$ \Lambda \f$     | Lambda           | \f$ \frac{\rho_c R^2}{T_c} \f$
\f$ \eta \f$        | eta              | array of the polar radii of the domains
\f$ \Delta\eta \f$  | deta             | array deta(i)=eta(i+1)-eta(i)
\f$ R_i \f$         | Ri               | array \f$R_i(\theta_k)\f$
\f$ \Delta R_i \f$  | dRi              | array \f$R_{i+1}(\theta_k)-R_i(\theta_k)\f$
\f$ r \f$           | r                | radial distance (spherical coordinate)
\f$ r_\zeta \f$     | rz               | \f$ \frac{\partial r}{\partial \zeta} \f$
\f$ \gamma \f$      | gamma            | metric element \f$\sqrt{1+r^2_\theta/r^2}\f$
\f$ \Omega \f$      | Omega            | angular velocity (equator)
\f$ \log\rho_c \f$ | log_rhoc          | log of central density
\f$ \log p_c \f$    | log_pc           | log of central pressure
\f$ \log T_c \f$    | log_Tc           | log of central temperature
\f$ \log R  \f$     | log_R            | log of polar radius
\f$ m \f$           | m                | mass
\f$ p_s \f$         | ps               | polar pressure
\f$ T_s \f$         | Ts               | polar temperature
\f$ lum \f$         | lum              | luminosity
\f$ Frad \f$        | Frad             | normal radiative flux   \f$\mathbf{E}^\zeta\cdot\mathbf{F}_{\rm rad}\f$ (not normalized)
\f$ T_{eff} \f$     | Teff             | effective temperature at star's surface
\f$ g_{sup} \f$     | gsup             | gravity at star's surface
\f$ \omega \f$      | w                | angular velocity
\f$ \Psi \f$        | G                | stream function (meridional circulation)
\f$ \rho \f$        | rho              | density
\f$ \xi \f$         | opa.xi           | radial conductivity
\f$ \kappa \f$      | opa.k            | opacity
\f$ \epsilon \f$    | nuc.eps          | nuclear power
\f$ s \f$           | s                | entropy
\f$ \mu \f$         | mu               | dynamic shear viscosity


\section eqlist_sec Equations solved by ESTER

Equation                                                                | Name
----------------------------------------------------------------------- | --------
\f$ \Delta \phi - \pi_c \rho = 0 \f$       | Poisson - for the gravitational potential
\f$ \nabla p + \rho \nabla \phi - \rho s \Omega^2 \hat{s} = 0 \f$       | momentum (steady and inviscid) - for the pressure
\f$ \hat{\varphi} \cdot \frac{\nabla P\times\nabla\rho}{\rho^2} - s\frac{\partial \Omega^2}{\partial z} = 0 \f$         | vorticity - for the differential rotation \f$\Omega(r,\theta)\f$
\f$ \nabla\cdot (\rho s^2 \Omega \mathbf{V})-\nabla\cdot (\mu s^2\nabla \Omega) = 0 \f$  | transport of angular momentum - for the meridional circulation
\f$ -\frac{\nabla\cdot (\xi \nabla T))}{\xi} + \frac{\Lambda \rho \epsilon}{\xi}= 0 \f$ | Heat - for the temperature

Where:
- \f$ s = r \sin\theta \f$
- \f$ \hat{s} = \nabla s\f$
- \f$ \hat{\varphi}\f$ is the azimuthal unit vector
- \f$ \mathbf{V} = \frac{\nabla\times (\Psi \hat{\varphi})}{\rho} \f$ is the meridional velocity

\section algo_sec Code Overview

The Newton iteration is performed in star2d::solve.
The set of equation and boundary conditions are written in the solver in
functions:
- star2d::solve_definitions
- star2d::solve_poisson
- star2d::solve_mov
- star2d::solve_temp
- star2d::solve_dim
- star2d::solve_map
- star2d::solve_Omega
- star2d::solve_atm
- star2d::solve_gsup
- star2d::solve_Teff

Function solver::solve stores all the terms registered in the operator in the
previous steps into the operator matrix (see solver::wrap, solver::unwrap
and solver::create).
And then solves the matrix and update the solution the solution fields.

*/