Bimetric 3+1 toolkit for spherical symmetry
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Functions
Maximal slicing

Functions

void BimetricEvolve::maximalSlice_2 (Int m, Int N, Real gAlp_at_N)
 Finds the maximal slice using a second-order accurate finite difference scheme.
void BimetricEvolve::maximalSlice_4 (Int m, Int N, Real gAlp_at_N)
 Finds the maximal slice using a fourth-order accurate finite difference scheme.
void BimetricEvolve::maximalSlice_2opt (Int m, Int N, Real gAlp_at_N)
 Finds the maximal slice using a second-order accurate finite difference scheme.
void BimetricEvolve::maximalSlice_PostSteps (Int m, Int N)
 Additional steps for gAlp and gDAlp (like fixing boundary conditions.)

Function Documentation

void BimetricEvolve::maximalSlice_2 ( Int  m,
Int  N,
Real  gAlp_at_N 
) [private]

Finds the maximal slice using a second-order accurate finite difference scheme.

Todo:
fixme: Should be k_g/2 in the above.

Removing 1/2, however, speeds up the collapse.

Todo:
fixme: Current implementation uses rho + J of the perfect fluid (should use bimetric).

The maximal slicing involves solving a differential equation: 

        y''(x) = P(x) y'(x) + Q(x) y(x)

with the Neumann BC on the left y´(0) = 0, and the Dirichlet BC on the right y(r_1) = 1.

The solution is found using a linear finite-difference method. The algorithm comprises two major steps:

  • Prepare the band-diagonal linear system corresponding to a discretized DE.
  • Solve the linear system (using Crout's method).

The coefficients for the discretized DE are generated in the Mathematica notebook: Maximal slicing, FD method.nb. The linear system is solved used Bandec.

Parameters:
mThe time slice
NThe radial coordinate of the right boundary
gAlp_at_NThe right boundary condition

Definition at line 45 of file maximalSlice.h.

{
    const Int offset = nGhost - 1;  N = N + 1;
    const Real h = delta_r;
    MatReal A( N, 3 );
    VecReal b( N );

    /////////////////////////////////////////////////////////////////////////////////////

    OMP_parallel_for( Int i = 0; i < N; ++i )
    {
        A[i][0] = -2         - h * PP( offset + i );
        A[i][1] =  4 + 2 * h * h * QQ( offset + i );
        A[i][2] = -2         + h * PP( offset + i );
        b[i]    =  0;
    }

    // Override the default row values

    A[0][0] = 0;
    A[0][2] = -4;

    A[N-1][2] = 0;
    b[N-1] = gAlp_at_N * ( 2 - h * PP( offset + (N-1) ) ) ;

    /////////////////////////////////////////////////////////////////////////////////////

    VecReal_O x( N );
    BandLUDecomposition( A, 1, 1 ).solve( b, x );

    // Retrieve the results from x

    gAlp( m, offset + N ) = gAlp_at_N;

    for( Int i = 0; i < N; ++i ) {
        gAlp( m, offset + i ) = x[i];
    }

    maximalSlice_PostSteps( m, N );
}
void BimetricEvolve::maximalSlice_2opt ( Int  m,
Int  N,
Real  gAlp_at_N 
) [private]

Finds the maximal slice using a second-order accurate finite difference scheme.

This is a slightly optimized version where the differential equation is solved using Algorithm 11.3 from Burden & Faires, Numerical Analysis, 2016 (found in the section Finite-Difference Methods for Linear Problems on p. 700).

See also:
maximalSlice_2
Parameters:
mThe time slice
NThe radial coordinate of the right boundary
gAlp_at_NThe right boundary condition

Definition at line 171 of file maximalSlice.h.

{
    const Int offset = nGhost - 1; N = N + 1;
    const Real h = delta_r;

    // First, allocate memory used to store the intermediate variables.
    //
    static Real* workArea = NULL;
    static Real *a, *b, *c, *d; // Used in steps 1-3
    static Real *l, *u, *z, *w; // Used in steps 4-8
    if ( workArea == NULL ) {
        size_t chunkSize = nLen + 2 * nGhost;
        workArea = new Real[ 8 * chunkSize ];
        a = workArea;      b = a + chunkSize; c = b + chunkSize; d = c + chunkSize;
        l = d + chunkSize; u = l + chunkSize; z = u + chunkSize; w = z + chunkSize;
    }

    // Note: a[] = diagonal, b[] = upper band, c[] = lower band,
    //       and d[] is on the right hand side of the equation A x = d

    /////////////////////////////////////////////////////////////////////////////////////
    // Prepare the tridiagonal linear system corresponding to DE.

    // Steps 1-3

    OMP_parallel_for( Int i = 0; i < N; ++i )
    {
        c[i] = -2         - h * PP( offset + i );
        a[i] =  4 + 2 * h * h * QQ( offset + i );
        b[i] = -2         + h * PP( offset + i );
        d[i] =  0;
    }

    c[0] =  0;
    b[0] = -4;  // Set by Neumann BC: w{-1} = w_{1}

    b[N-1] =   0;
    d[N-1] = ( 2 - h * PP( offset + N-1 ) ) * gAlp_at_N;

    /////////////////////////////////////////////////////////////////////////////////////
    // Solve the linear system using Crout's factorization.

    // Step 4
    //
    l[0] = a[0];
    u[0] = b[0] / a[0];
    z[0] = d[0] / l[0];

    // Step 5
    //
    for( Int i = 0; i < N - 1; ++i )
    {
        l[i] = a[i] - c[i] * u[i-1];
        u[i] = b[i] / l[i];
        z[i] = ( d[i] - c[i] * z[i-1] ) / l[i];
    }

    // Step 6-7
    //
    Int i = N - 1;
    l[i] = a[i] - c[i] * u[i-1];
    gAlp( m, offset + i ) = w[i] = z[i] = ( d[i] - c[i] * z[i-1] ) / l[i];

    // Step 8
    //
    for( Int i = N - 2; i >= 1; --i ) {
        gAlp( m, offset + i ) = w[i] = z[i] - u[i] * w[i+1];
    }

    gAlp( m, offset + N ) = w[N] = gAlp_at_N;

    maximalSlice_PostSteps( m, N );
}
void BimetricEvolve::maximalSlice_4 ( Int  m,
Int  N,
Real  gAlp_at_N 
) [private]

Finds the maximal slice using a fourth-order accurate finite difference scheme.

See also:
maximalSlice_2
Parameters:
mThe time slice
NThe radial coordinate of the right boundary
gAlp_at_NThe right boundary condition

Definition at line 94 of file maximalSlice.h.

{
    const Int offset = nGhost;
    const Real h = delta_r;
    MatReal A( N, 7 );
    VecReal b( N );

    /////////////////////////////////////////////////////////////////////////////////////

    OMP_parallel_for( Int i = 0; i < N; ++i )
    {
        A[i][0] =   0;
        A[i][1] =  -1          - h * PP( offset + i );
        A[i][2] =  16      + 8 * h * PP( offset + i );
        A[i][3] = -30 - 12 * h * h * QQ( offset + i );
        A[i][4] =  16      - 8 * h * PP( offset + i );
        A[i][5] =  -1        +   h * PP( offset + i );
        A[i][6] =   0;
        b[i]    =   0;
    }

    // Override the default row values

    Int i = 0;
    A[i][1] =   0;
    A[i][2] =   0;
    A[i][3] = -20 + 10 * h * PP( offset + i ) - 12 * h * h * QQ( offset + i );
    A[i][4] =  17 - 15 * h * PP( offset + i );
    A[i][5] =   4 +  6 * h * PP( offset + i );
    A[i][6] =  -1 -      h * PP( offset + i );

    i = 1;
    A[i][1] = 0;
    A[i][3] = -31 - h * PP( offset + i ) - 12 * h * h * QQ( offset + i );

    i = N - 2;
    A[i][5] = 0;
    b[i] = gAlp_at_N * ( 1 - 8 * h * PP( offset + i ) ) ;

    i = N - 1;
    A[i][0] =  -1 +      h * PP( offset + i );
    A[i][1] =   4 -  6 * h * PP( offset + i );
    A[i][2] =   6 + 18 * h * PP( offset + i );
    A[i][3] = -20 - 10 * h * PP( offset + i ) - 12 * h * h * QQ( offset + i );
    A[i][4] =   0;
    A[i][5] =   0;
    b[i] = gAlp_at_N * ( -11 + 3 * h * PP( offset + i ) ) ;

    /////////////////////////////////////////////////////////////////////////////////////

    VecReal_O x( N );
    BandLUDecomposition( A, 3, 3 ).solve( b, x );

    // Retrieve the results from x

    gAlp( m, offset + N ) = gAlp_at_N;

    for( Int i = 0; i < N; ++i ) {
        gAlp( m, offset + i ) = x[i];
    }

    maximalSlice_PostSteps( m, N );
}
void BimetricEvolve::maximalSlice_PostSteps ( Int  m,
Int  N 
) [private]

Additional steps for gAlp and gDAlp (like fixing boundary conditions.)

Todo:
fixme: add +TINY in the case if zero

Definition at line 251 of file maximalSlice.h.

{
    // Fix the ghost cells on the left (using parity BC)
    //
    for( Int i = 0; i < nGhost; ++i ) {
        gAlp( m, nGhost - 1 - i  ) = gAlp( m, nGhost + i );
    }

    // Determine cells on the right by solving the differential equation.
    //
    Real h = delta_r;
    Int overlap = 10;
    for( Int n = nGhost + N - overlap; n < nGhost + nLen + nGhost - 1; ++n )
    {
        gAlp( m, n+1 ) =
            ( ( 4  + 2 * h * h * QQ(n) ) * gAlp(m,n) - ( 2 + h * PP(n) ) * gAlp(m,n-1) )
            / ( 2 - h * PP(n) );
    }

    #if 0
    // Alternative smoothing:
    //
    OMP_parallel_for( Int n = nGhost; n < nTotal - CFDS_ORDER/2; ++n ) {
        gAlp_r ( m, n ) = GF_r( gAlp,  m, n );
    }
    smoothenGF( m, fld::gAlp_r, fld::gdbg, fld::gAlp_r, -1 );
    OMP_parallel_for( Int n = nGhost; n < nTotal - CFDS_ORDER/2; ++n ) {
        gDAlp(m,n) = gAlp_r(m,n) / gAlp(m,n);
    }
    smoothenGF( m, fld::gDAlp, fld::fdbg, fld::gDAlp, -1 );
    return;
    #endif

    // Calculate gDAlp
    //
    OMP_parallel_for( Int n = nGhost; n < nTotal - CFDS_ORDER/2; ++n )
    {
        gAlp_r ( m, n ) = GF_r( gAlp,  m, n );
        /// @todo fixme:  add +TINY in the case if zero
        gDAlp(m,n) = gAlp_r(m,n) / gAlp(m,n);
    }

    // The accuracy is very low at low r. Assume that a few first cells are linear,
    // then apply a six-point cubic spline to smooth the region around r = 0.
    //
    cubicSplineShmooth( m, fld::gDAlp, lin2n, cub2n );

    // Fix the left boundary using parity BC
    //
    for( Int i = 0; i < nGhost; ++i ) {
        gDAlp( m, nGhost - 1 - i  ) = - gDAlp( m, nGhost + i );
    }
}