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

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

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

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 :

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

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

This is equivalent to the code:

 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 \( \Delta \phi \) (see symbolic).

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

Line 3 says that equation "eq_Phi" also depends on \( r \).

Lines 5 and 6 add the remaining terms to the Jacobian matrix. We see the coefficient \( -\pi_c \) of \(\delta\rho \) and \(-\rho\) of \(\delta\pi_c\). 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"

Variables in ESTER

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

Equations solved by ESTER

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

Where:

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

Code Overview

The Newton iteration is performed in star2d::solve. The set of equation and boundary conditions are written in the solver in functions:

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.