Integration

Functions
void	BimetricEvolve::calculateDerivedVariables (Int m, Int n)
	Calculate the derived variables R, Lt, and pfv.
void	BimetricEvolve::determineGaugeFunctions (Int m)
	Determine gAlp (e.g.
void	BimetricEvolve::applySommerfeldBC (Int m, Int n, Int m_old, Int gf, Real background=0.0, Int fall_off=1)
	Sommerfeld outgoing (radiative) wave condition on a variable.
virtual void	BimetricEvolve::integStep_Prepare (Int m)
	Calculate the dependent variables from the prime state variables which are needed for the integration.
virtual void	BimetricEvolve::integStep_CalcEvolutionRHS (Int m)
	Prepare the right-hand side of the evolution equations.
virtual void	BimetricEvolve::integStep_Finalize (Int mNext, Int mPrev)
	Perform the additional steps after each integration step.
virtual bool	BimetricEvolve::integStep_Diagnostics (Int m, Int chkNaNs_nFrom, Int chkNaNs_nTo)
	At each checkpoint calculate diagnostics variables and check for eventual NaNs.
virtual void	BimetricEvolve::applyLeftBoundaryCondition (Int m)
	Fix the state variables at the left boundary.
virtual void	BimetricEvolve::applyRightBoundaryCondition (Int m)
	Fix the state variables at the right boundary.

---

Function Documentation

void BimetricEvolve::applyLeftBoundaryCondition

(

Int

m

)

[private, virtual]

Fix the state variables at the left boundary.

Here we impose the parity conditions at the left boundary.

Implements IntegFace.

Definition at line 666 of file bim-solver.cpp.

{
    for( Int i = 0; i < nGhost; ++i )
    {
        Int n  = nGhost - i - 1;
        Int nR = nGhost + i;

        /// Here we impose the parity conditions at the left boundary.

        t     (m,n) = t      (m,nR);
        r     (m,n) = -r     (m,nR);

        p     (m,n) = -p     (m,nR);

        gA    (m,n) = gA     (m,nR);      fA   (m,n) = fA    (m,nR);
        gB    (m,n) = gB     (m,nR);      fB   (m,n) = fB    (m,nR);
        gK    (m,n) = gK     (m,nR);      fK   (m,n) = fK    (m,nR);
        gKD   (m,n) = gKD    (m,nR);      fKD  (m,n) = fKD   (m,nR);
        gDA   (m,n) = -gDA   (m,nR);      fDA  (m,n) = -fDA  (m,nR);
        gDB   (m,n) = -gDB   (m,nR);      fDB  (m,n) = -fDB  (m,nR);
        gSig  (m,n) = -gSig  (m,nR);      fSig (m,n) = -fSig (m,nR);

        pfD   (m,n) = pfD    (m,nR);
        pfS   (m,n) = -pfS   (m,nR);
        pftau (m,n) = pftau  (m,nR);

        q     (m,n) = -q     (m,nR);
        gAlp  (m,n) = gAlp   (m,nR);
        fAlp  (m,n) = fAlp   (m,nR);
        gDAlp (m,n) = -gDAlp (m,nR);

        Lt    (m,n) = Lt     (m,nR);
        R     (m,n) = R      (m,nR);
        pfv   (m,n) = -pfv   (m,nR);
    }
}

void BimetricEvolve::applyRightBoundaryCondition

(

Int

m

)

[private, virtual]

Fix the state variables at the right boundary.

Implements IntegFace.

Definition at line 703 of file bim-solver.cpp.

{
    for( Int n = nGhost + nLen; n < nTotal; ++n )
    {
        t(m,n) = t(m,n-1);
        r(m,n) = r(m,n-1) + delta_r;

        extrapolate_R( fld::p,     m, n );   // Extrapolate [n] from [n-1], [n-2], ...

        extrapolate_R( fld::gA,    m, n );      extrapolate_R( fld::fA,    m, n );
        extrapolate_R( fld::gB,    m, n );      extrapolate_R( fld::fB,    m, n );
        extrapolate_R( fld::gK,    m, n );      extrapolate_R( fld::fK,    m, n );
        extrapolate_R( fld::gKD,   m, n );      extrapolate_R( fld::fKD,   m, n );
        extrapolate_R( fld::gDA,   m, n );      extrapolate_R( fld::fDA,   m, n );
        extrapolate_R( fld::gDB,   m, n );      extrapolate_R( fld::fDB,   m, n );
        extrapolate_R( fld::gSig,  m, n );      extrapolate_R( fld::fSig,  m, n );

        extrapolate_R( fld::pfD,   m, n );
        extrapolate_R( fld::pfS,   m, n );
        extrapolate_R( fld::pftau, m, n );

        extrapolate_R( fld::q,     m, n );
        extrapolate_R( fld::gAlp,  m, n );
        extrapolate_R( fld::gDAlp, m, n );
        extrapolate_R( fld::fAlp,  m, n );

        calculateDerivedVariables( m, n ); // Calculate R, Lt, and pfv.

        extrapolate_R( fld::pfv, m, n ); // Nevertheless, we extrapolate pfv!
    }
}

void BimetricEvolve::applySommerfeldBC	(	Int	m,
		Int	n,
		Int	m_old,
		Int	gf,
		Real	background = 0.0,
		Int	fall_off = 1
	)		[private]

Sommerfeld outgoing (radiative) wave condition on a variable.

Parameters:

m	The time slice
n	The spatial position
m_old	The earlier time
gf	A grid function index
background	Background value; 1.0 for ALPHA
fall_off	Fall-off power exponent; 2 for VET, BET

Definition at line 736 of file bim-solver.cpp.

{
    const Real wave_speed = 1.0;
    const Real Courant = wave_speed * delta_t / delta_r;

    for( Int n = nGhost + nLen + 1; n < nTotal - 1; ++n )
    {
        const Real r0     = r( m, n-2 );
        const Real r1     = r( m, n-1 );
        const Real r2     = r( m, n   );
        const Real r_last = Courant * r1 + ( 1.0 - Courant ) * r2;

        Real factor = r_last / r2;
        if( fall_off == 2 ) {
            factor *= factor;
        }

        const Real P0   = GF( gf, m, n-2 );
        const Real P1   = GF( gf, m, n-1 );
        const Real P2   = GF( gf, m, n );
        const Real P01  = ( ( r_last - r1 ) * P0  + ( r0 - r_last ) * P1  ) / ( r0 - r1 );
        const Real P12  = ( ( r_last - r2 ) * P1  + ( r1 - r_last ) * P2  ) / ( r1 - r2 );
        const Real P012 = ( ( r_last - r2 ) * P01 + ( r0 - r_last ) * P12 ) / ( r0 - r2 );

        GF( gf, m1, n ) = factor * ( P012 - background ) + background;
    }
}

void BimetricEvolve::calculateDerivedVariables	(	Int	m,
		Int	n
	)		[inline, private]

Calculate the derived variables R, Lt, and pfv.

Note:

TINY_Real is added to the denominator of pfv to avoid dividing by zero.

Definition at line 521 of file bim-solver.cpp.

    {
        R(m,n)   = fB(m,n) / gB(m,n);
        Lt(m,n)  = sqrt( 1 + p(m,n) * p(m,n) );
        pfv(m,n) = pfS(m,n) == 0 ? 0 : eq_pf_v(m,n);
    }

void BimetricEvolve::determineGaugeFunctions

(

Int

m

)

[private]

Determine gAlp (e.g.

by maximal slicing). Then find fAlp and also calculate the derivatives of the lapse related grid-functions.

• Maximal slicing (optional) -- this will calculate gAlp and gDAlp

Todo:

Go through the order of evaluation for the gauge grid functions

Todo:

GF_r and GF_rr requires fixed ghost zones

Definition at line 781 of file bim-solver.cpp.

{
    /////////////////////////////////////////////////////////////////////////////////////
    /// - Maximal slicing (optional) -- this will calculate gAlp and gDAlp

    Int sliceBC = smooth >= 1 ? round(40/delta_r) : nLen/2;
    switch( slicing )
    {
        case SLICE_MS2OPT: maximalSlice_2opt ( m, sliceBC, 1 );  break;
        case SLICE_MS2:    maximalSlice_2    ( m, sliceBC, 1 );  break;
        case SLICE_MS4:    maximalSlice_4    ( m, sliceBC, 1 );  break;
    }

#if 1
    #define KDT
#else
    #define KDT 0*
    if( true )
    {
        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
        {
            gW   ( m, n ) = eq_gW   ( m, n );
            fW   ( m, n ) = eq_fW   ( m, n );
            cW   ( m, n ) = eq_cW   ( m, n );
            fAlp ( m, n ) = eq_fAlp ( m, n );

            gDAlp_r( m, n ) = 0;
            fAlp_r ( m, n ) = 0;
            fAlp_rr( m, n ) = 0;
            fDAlp  ( m, n ) = 0;
            fDAlp_r( m, n ) = 0;
        }
    } else
#endif

    /// @todo Go through the order of evaluation for the gauge grid functions
    /// @todo GF_r and GF_rr requires fixed ghost zones
    ///
    if ( isGR () || slicing == SLICE_CONSTGF )
    {
        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
        {
            gAlp_r ( m, n ) = GF_r( gAlp,  m, n );
            gDAlp_r( m, n ) = GF_r( gDAlp, m, n );
            fAlp_r ( m, n ) = GF_r( fAlp,  m, n );
            fAlp_rr( m, n ) = GF_rr( fAlp, m, n );
            fDAlp  ( m, n ) = fAlp_r(m,n) / fAlp(m,n);
            fDAlp_r( m, n ) = ( fAlp_rr(m,n) - fAlp_r(m,n) * fAlp_r(m,n) / fAlp(m,n) )
                              / fAlp(m,n);
        }
    }
    else // if not GR
    {
        // gAlp and gDAlp must the BC fixed at this point
        //
        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) {
            gAlp_r ( m, n ) = GF_r( gAlp,  m, n );
            gDAlp_r( m, n ) = GF_r( gDAlp, m, n );
        }

        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) {
            cW   ( m, n ) = eq_cW   ( m, n );
            gW   ( m, n ) = eq_gW   ( m, n );
            fW   ( m, n ) = eq_fW   ( m, n );
            fAlp ( m, n ) = eq_fAlp ( m, n );
        }

        if( smooth ) {
            smoothenGF( m, fld::fAlp, fld::tmp, fld::fAlp, 1 );
        }
        else {
            applyBoundaryConditions( m, fld::fAlp, +1 );
        }

        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) {
            fAlp_r( m, n ) = GF_r( fAlp, m, n );
        }

        if( smooth ) {
            smoothenGF( m, fld::fAlp_r, fld::tmp, fld::fAlp_r, -1 );
        }
        else {
            applyBoundaryConditions( m, fld::fAlp_r, -1 );
        }

        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n ) {
            fAlp_rr( m, n ) = GF_r( fAlp_r, m, n );
        }

        if( smooth ) {
            smoothenGF( m, fld::fAlp_rr, fld::tmp, fld::fAlp_rr, 1 );
        }
        else {
            applyBoundaryConditions( m, fld::fAlp_r, -1 );
        }

        OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
        {
            fDAlp  ( m, n ) = fAlp_r(m,n) / fAlp(m,n);
            fDAlp_r( m, n ) = ( fAlp_rr(m,n) - fAlp_r(m,n) * fAlp_r(m,n) / fAlp(m,n) )
                               / fAlp(m,n);
        }
    }
}

void BimetricEvolve::integStep_CalcEvolutionRHS

(

Int

m

)

[private, virtual]

Prepare the right-hand side of the evolution equations.

• First, calculate the values of the spatial derivatives

  Todo:

  fixme: Severely lost accuracy in all second derivatives

• Calculate the gauge (e.g., do maximal slicing)

• Calculate the intermediate variables that do not depend on derivatives.

• Calculate the variables that depend on the spatial derivatives

• Smoothen the time derivatives inside the grid zone near the outer boundary

  Todo:

  Smoothen the time derivatives inside the left grid zone

Implements IntegFace.

Definition at line 886 of file bim-solver.cpp.

{
    /////////////////////////////////////////////////////////////////////////////////////
    /// - First, calculate the values of the spatial derivatives
    /// @todo fixme: Severely lost accuracy in all second derivatives

    OMP_parallel_for( Int n = CFDS_ORDER/2; n < nTotal - CFDS_ORDER/2; ++n )
    {
        gA_r   (m,n) = GF_r( gA,    m,n );    fA_r   (m,n) = GF_r( fA,    m,n );
        gB_r   (m,n) = GF_r( gB,    m,n );    fB_r   (m,n) = GF_r( fB,    m,n );
        gK_r   (m,n) = GF_r( gK,    m,n );    fK_r   (m,n) = GF_r( fK,    m,n );
        gKD_r  (m,n) = GF_r( gKD,   m,n );    fKD_r  (m,n) = GF_r( fKD,   m,n );
        gDA_r  (m,n) = GF_r( gDA,   m,n );    fDA_r  (m,n) = GF_r( fDA,   m,n );
        gDB_r  (m,n) = GF_r( gDB,   m,n );    fDB_r  (m,n) = GF_r( fDB,   m,n );

        p_r    (m,n) = GF_r( p,     m,n );    q_r    (m,n) = GF_r( q,     m,n );
        p_rr   (m,n) = GF_rr( p,    m,n );    q_rr   (m,n) = GF_rr( q,    m,n );
        eq_pr_r(m,n) = GF_r( eq_pr, m,n );    eq_qr_r(m,n) = GF_r( eq_qr, m,n );

        pfD_r  (m,n) = GF_r( pfD,   m,n );
        pfS_r  (m,n) = GF_r( pfS,   m,n );
        pftau_r(m,n) = GF_r( pftau, m,n );

        // pfv_r(m,n) should be last as it depends on pfD_r, pfS_r, and pftau_r
        //
        pfv_r(m,n) = pfS(m,n) == 0 ? 0 : eq_pf_v_r( m, n );
    }

    if( smooth >= 2 )
    {
        smoothenGF( m, fld::p_r,      fld::tmp, fld::p_r,      +1 );
        smoothenGF( m, fld::p_rr,     fld::tmp, fld::p_rr,     -1 );
        smoothenGF( m, fld::eq_pr_r,  fld::tmp, fld::eq_pr_r,  +1 );

        smoothenGF( m, fld::gA_r,  fld::tmp, fld::gA_r,  -1 );
        smoothenGF( m, fld::gB_r,  fld::tmp, fld::gB_r,  -1 );
        smoothenGF( m, fld::gK_r,  fld::tmp, fld::gK_r,  -1 );
        smoothenGF( m, fld::gKD_r, fld::tmp, fld::gKD_r, -1 );
        smoothenGF( m, fld::gDA_r, fld::tmp, fld::gDA_r, +1 );
        smoothenGF( m, fld::gDB_r, fld::tmp, fld::gDB_r, +1 );

        smoothenGF( m, fld::fA_r,  fld::tmp, fld::fA_r,  -1 );
        smoothenGF( m, fld::fB_r,  fld::tmp, fld::fB_r,  -1 );
        smoothenGF( m, fld::fK_r,  fld::tmp, fld::fK_r,  -1 );
        smoothenGF( m, fld::fKD_r, fld::tmp, fld::fKD_r, -1 );
        smoothenGF( m, fld::fDA_r, fld::tmp, fld::fDA_r, +1 );
        smoothenGF( m, fld::fDB_r, fld::tmp, fld::fDB_r, +1 );
    }

    /// - Calculate the gauge (e.g., do maximal slicing)
    ///
    determineGaugeFunctions( m );

    /////////////////////////////////////////////////////////////////////////////////////
    /// - Calculate the intermediate variables that do not depend on derivatives.

    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        // The sources (dependent on both p and the primary dynamical fields)
        //
        g_rho (m,n) = eq_g_rho (m,n);      f_rho (m,n) = eq_f_rho (m,n);
        g_JK1 (m,n) = eq_g_JK1 (m,n);      f_JK1 (m,n) = eq_f_JK1 (m,n);
        g_JK2 (m,n) = eq_g_JK2 (m,n);      f_JK2 (m,n) = eq_f_JK2 (m,n);
    }

    /////////////////////////////////////////////////////////////////////////////////////
    /// - Calculate the variables that depend on the spatial derivatives

    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        // The time evolution right-hand side
        //
        p_t     (m,n) = isGR() ? 0 : KDT eq_p_t(m,n); // + eq_C_bimConsLaw (m,n);

        gA_t    (m,n) = eq_gA_t    (m,n);       fA_t    (m,n) = eq_fA_t    (m,n);
        gB_t    (m,n) = eq_gB_t    (m,n);       fB_t    (m,n) = eq_fB_t    (m,n);
        gK_t    (m,n) = eq_gK_t    (m,n);       fK_t    (m,n) = eq_fK_t    (m,n);
        gKD_t   (m,n) = KDT eq_gKD_t   (m,n);   fKD_t   (m,n) = KDT eq_fKD_t   (m,n);
        gDA_t   (m,n) = KDT eq_gDA_t   (m,n);   fDA_t   (m,n) = KDT eq_fDA_t   (m,n);
        gDB_t   (m,n) = KDT eq_gDB_t   (m,n);   fDB_t   (m,n) = KDT eq_fDB_t   (m,n);
        gSig_t  (m,n) = KDT eq_gSig_t  (m,n);   fSig_t  (m,n) = KDT eq_fSig_t  (m,n);

        pfD_t   (m,n) = eq_pf_D_t  (m,n);
        pfS_t   (m,n) = eq_pf_S_t  (m,n);
        pftau_t (m,n) = eq_pf_tau_t(m,n);

        if( slicing == SLICE_BM )
        {
            gAlp_t  (m,n) = eq_BM_gAlp_t  (m,n);
            gDAlp_t (m,n) = eq_BM_gDAlp_t (m,n);
        }
    }

    if( smooth >= 2 )
    {
        // smoothenGF2( m, fld::p_t,   fld::tmp, fld::p_t,   -1 );
        // smoothenGF( m, fld::gA_t,  fld::tmp, fld::gA_t,  +1 );
        // smoothenGF( m, fld::gB_t,  fld::tmp, fld::gB_t,  +1 );
        smoothenGF2( m, fld::gK_t,  fld::tmp, fld::gK_t,  +1 );
        smoothenGF2( m, fld::gKD_t, fld::tmp, fld::gKD_t, +1 );
        // smoothenGF2( m, fld::gDA_t, fld::tmp, fld::gDA_t, -1 );
        // smoothenGF2( m, fld::gDB_t, fld::tmp, fld::gDB_t, -1 );
        // smoothenGF( m, fld::fA_t,  fld::tmp, fld::fA_t,  +1 );
        // smoothenGF( m, fld::fB_t,  fld::tmp, fld::fB_t,  +1 );
        smoothenGF2( m, fld::fK_t,  fld::tmp, fld::fK_t,  +1 );
        smoothenGF2( m, fld::fKD_t, fld::tmp, fld::fKD_t, +1 );
        // smoothenGF2( m, fld::fDA_t, fld::tmp, fld::fDA_t, -1 );
        // smoothenGF2( m, fld::fDB_t, fld::tmp, fld::fDB_t, -1 );
    }

    /////////////////////////////////////////////////////////////////////////////////////
    /// - Smoothen the time derivatives inside the grid zone near the outer boundary
    /// @todo Smoothen the time derivatives inside the left grid zone

    if ( gridDriver->isLastInRank() )
    {
        for( Int n = nGhost + nLen - CFDS_ORDER; n < nGhost + nLen; ++n )
        {
            for( auto e: fld::bimEvolvedGF ) {
                extrapolate_lin( e.f_t, m, n );
            }
        }
    }
}

bool BimetricEvolve::integStep_Diagnostics	(	Int	m,
		Int	chkNaNs_nFrom,
		Int	chkNaNs_nTo
	)		[private, virtual]

At each checkpoint calculate diagnostics variables and check for eventual NaNs.

Find horizons, calculate constraints and similar.

Returns 'false' as the indicator to abort the integration.

Also check for NaNs.

• Apparent horizon finder (dependent only on the primary dynamical fields)

• The constraint violations

• Calculate p in an alternative way using the algebraic equations

• Check for NaNs in gA in the given zone.

Implements IntegFace.

Definition at line 1021 of file bim-solver.cpp.

{
    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        /// - Apparent horizon finder (dependent only on the primary dynamical fields)
        ///
        gHorz (m,n) = eq_gHorz (m,n);
        fHorz (m,n) = eq_fHorz (m,n);

        /// - The constraint violations
        ///
        C_1 (m,n) = eq_C_1 (m,n);
        C_2 (m,n) = eq_C_2 (m,n);
        C_3 (m,n) = eq_C_3 (m,n);
        C_4 (m,n) = eq_C_bimConsLaw (m,n);

        // C_1 (m,n) = abs( eq_C_1 (m,n) / (1e-10 + g_rho (m,n) ) );
        // C_2 (m,n) = abs( eq_C_2 (m,n) / (1e-10 + f_rho (m,n) ) );
        // C_3 (m,n) = eq_C_3 (m,n);
        // C_4 (m,n) = abs( eq_C_bimConsLaw (m,n) / (1e-10 + p_t (m,n) ) );

        // C_1 (m,n) = abs( eq_C_1 (m,n) );
        // C_2 (m,n) = abs( eq_C_2 (m,n) );
        // C_3 (m,n) = abs( eq_C_3 (m,n) );
        // C_4 (m,n) = abs( eq_C_bimConsLaw (m,n) );

        /// - Calculate p in an alternative way using the algebraic equations
        ///
        p_g (m,n) = eq_p_g (m,n);
        p_f (m,n) = eq_p_f (m,n);
    }

    if ( chkNaNs_nTo < 0 ) { // CheckNaNs is disabled
        return true;
    }

    /// - Check for NaNs in gA in the given zone.
    ///
    Int nFrom = nGhost + std::min( nLen, std::max( Int(0), chkNaNs_nFrom - nOffset ) );
    Int nTo   = nGhost + std::min( nLen, std::max( Int(0), chkNaNs_nTo   - nOffset ) );

    // if( false )
    for( Int n = nFrom; n < nTo; ++n )
    {
        if( std::isnan( gA( m, n ) ) ) {
            std::cerr << "*** Detected gA NaN at t = " << t(m,n)
                      << ", r = " << r(m,n) << std::endl;
            return false;
        }
    }
    return true;
}

void BimetricEvolve::integStep_Finalize	(	Int	mNext,
		Int	mPrev
	)		[private, virtual]

Perform the additional steps after each integration step.

Post-evolution smoothing and calculation of diagnostics.

Implements IntegFace.

Definition at line 1014 of file bim-solver.cpp.

{
}

void BimetricEvolve::integStep_Prepare

(

Int

m

)

[private, virtual]

Calculate the dependent variables from the prime state variables which are needed for the integration.

Implements IntegFace.

Definition at line 768 of file bim-solver.cpp.

{
    // Move gKD and fKD to 0 at r = 0
    // Real dgKD = gKD(m,nGhost)/2, dfKD = fKD(m,nGhost)/2;

    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        calculateDerivedVariables( m, n );
        // gKD(m,n) -= dgKD;
        // fKD(m,n) -= dfKD;
    }
}