Initial Models

Here we briefly describe the various routines for generating initial models for the main MAESTROeX problems and how the initial model is used to initialize both the base state and the full Cartesian state.

Note

When run, MAESTROeX outputs info that’s useful to check, including dr of base state and dr of the input file, which should be the same at the finest level; and Maximum HSE Error, which should be some small number.

Creating the Model Data from Raw Data

We have found that for the best performance, the MAESTROeX initialization procedure should be given model data with the same resolution as the base state resolution, Δr.

  • For planar problems, Δr=Δx.

  • For multilevel planar problems, we use Δr equal to Δx at the finest resolution.

  • For spherical problems we set Δr=Δx/drdxfac.

We generate our initial model either analytically or from raw data, ρraw,Traw,praw, and Xraw.

Here is the raw data file for each test problem:

(95)reacting_bubblenoneit is generated analyticallytest_convectnoneit is generated analyticallywdconvectmodel_6.e8.rawspherical_heatnoneit is generated analytically

We use a C++ program to interpolate the raw data, yielding the model data, ρmodel,Tmodel,pmodel, and Xmodel. These initial model routines generally uses an iterative procedure to modify the model data so that it is thermodynamically consistent with the MAESTROeX equation of state (EOS), and also satisfies our chosen hydrostatic equilibrium (HSE) discretization,

(96)prpr1Δr=ρr+ρr12gr12,

to a user defined tolerance. The shared initial_models git repo (https://github.com/AMReX-Astro/initial_models) has the routines used for creating initial models. Here are the names of the initial_models directories used for different problems:

problem name

initial_models routine

reacting_bubble

test2

test_convect

test2

wdconvect

none

spherical_heat

spherical

The model data is not generated at run-time—it must be generated in advance of running any MAESTROeX examples. The inputs file should point to the file containing the model data.

Creating the Initial Data from the Model Data

In base_state.f90, using the subroutine init_base_state (which is actually a terrible name since we are computing 1D initial data, which is not quite the same thing as the base satte), we set the initial data ρinit,Tinit,pinit and Xinit equal to the model data. Then, we set pinit,hinit=p,h(ρinit,Tinit,Xinit) through the EOS. Note that we overwrite the existing value of pinit but this should not change the value since we called pmodel=pmodel(ρmodel,Tmodel,Xmodel) at the end of the initial model generation routines in § 1.

For multilevel planar problems, we generate initial data at the non-finest levels by using linear interpolation from the two nearest model points (see the figure below). Note that the initial data is not in HSE nor is it thermodynamically consistent on non-finest cells. We will enforce these on the base state and full state later.

_images/multilevel_init.png

Fig. 11 Multilevel Initialization. The data from the initial model is represented by the dots on the right. The initial data at the finest level is represented by the circles. The initial data at non-finest levels is represented by the squares. We copy data from the dots directly into the circles. We use linear interpolation using the two nearest points to copy data from the dots to the squares.

We deal with the edge of the star by tracking the first coarse cell in which the density falls below base_cutoff_density. We note the radius of this cell center, and the values of ρ, ρh, Xk, p, and T in this cells. Then, at every level, if current cell-center is above this radius, we set the state equal to this stored state. This ensures a consistent treatment of the edge of the star at all levels.

Creating the Base State and Full State

Given pinit,ρinit,Tinit, and Xinit, this section describes how the base state (ρ0 and p0) and full state (ρ,h,X, and T) are computed. The base state is, general, not simply a direct copy of the initial model data, since we require that ρ0=ρ. Additionally, we require thebase state to be HSE according to Eq.96, and that the full state is thermodynamically consistent with p0. Overall we do:

  1. Fill ρinit,hinit,Xinit, and Tinit onto a multifab to obtain the full state ρ,h,X, and T.

  2. If perturb_model = T, a user-defined perturbation is added. This routine should make sure that the EOS is called so that there is some sense of thermodynamic consistency.

  3. Set ρ0=ρ.

  4. Compute p0 using enforce_HSE.

  5. Compute T,h=T,h(ρ,p0,Xk).

  6. Set (ρh)0=(ρh).

  7. Compute T. Note that we only use T as a diagnostic and as a seed for EOS calls.

Now ρ0=ρ, the base state is in HSE, and the full state is thermodynamically consistent with p0.

Coarse-Fine enforce_HSE Discretization

When integrating the HSE discretization upward, we must use a different differencing procedure at coarse-fine interfaces. The figure below shows the transition from coarse (level l1) to fine (level l), with the zone center indices noted.

_images/ctof.png

Fig. 12 A coarse-fine interface in the 1-d base state

To find the zone-centered pressure in the first fine zone, prl, from the zone-centered pressure in the coarse zone just below the coarse-fine interface, pr/21l1, we integrate in 2 steps. We allow for a spatially changing gravitational acceleration, for complete generality.

First we integrate up to the coarse-fine interface from the coarse-cell center as:

(97)pr12lpr/21l1Δrl1/2=ρr12l+ρr/21l12gr12l+gr/21l12

We can rewrite this as an expression for the pressure at the coarse-fine interface:

(98)pr12l=pr/21l1+Δrl18(ρr12l+ρr/21l1)(gr12l+gr/21l1).

Next we integrate up from the coarse-fine interface to the fine-cell center:

(99)prlpr12lΔrl/2=ρrl+ρr12l2grl+gr12l2

We can rewrite this as an expression for the pressure at the fine-cell center:

(100)prl=pr12l+Δrl8(ρrl+ρr12l)(grl+gr12l).

Combining Eq.98 and Eq.100 gives

(101)prl=pr/21l1+Δrl18(ρr12l+ρr/21l1)(gr12l+gr/21l1)+Δrl8(ρrl+ρr12l)(grl+gr12l).

We can simplify using

(102)Δrl1=2Δrl,

and by interpolating the cell-centered densities to the coarse-fine interface as:

(103)ρr12l=23ρrl+13ρr/21l1.

Because we carry both the cell- and edge-centered gravitational accelerations, we do not need to interpolate g to the interface. Simplifying, we have

(104)prl=pr/21l1+Δrl4(23ρrl+43ρr/21l1)(gr12l+gr/21l1)+Δrl8(53ρrl+13ρr/21l1)(grl+gr12l).

Finally, we note for constant g, this simplifies to:

(105)prl=pr/21l1+3Δrlg4(ρr/21l1+ρrl).

When integrating across a fine-coarse interface (see the figure below), the proceduce is similar.

_images/ftoc.png

Fig. 13 A fine-coarse interface in the 1-d base state

The expression for general gravity becomes:

(106)p(r+1)/2l1=prl+Δrl4(23ρrl+43ρ(r+1)/2l1)(g(r+1)/212l1+g(r+1)/2l1)+Δrl8(53ρrl+13ρ(r+1)/2l1)(grl+g(r+1)/212l1).

and for spatially-constant gravity, it simplifies to:

(107)p(r+1)/2l1=prl+3Δrlg4(ρrl+ρ(r+1)/2l1).