grid/integrator.h

Go to the documentation of this file.

/**
 *  @file      integrator.h
 *  @brief     Time evolution using the Method of Lines (MoL).
 *  @authors   Mikica Kocic, Wilkinson & Reinsch, Press et al
 *  @copyright GNU General Public License (GPLv3).
 */

#include <csignal>

/////////////////////////////////////////////////////////////////////////////////////////
/** @defgroup g11 Numerical Methods                                                    */
/////////////////////////////////////////////////////////////////////////////////////////
                                                                                   /*@{*/
/** Time evolution using the Method of Lines (MoL).
 *
 *  The implemented methods are declared in MoL::knownMethods.
 *
 *  The integrator can evolve equations from several modules at the same time.
 *  All such modules should implement @ref IntegFace ans sign in to the integrator.
 *  The modules should also specify which grid functions are either evolved or
 *  kept constant by the integrator. (Those grid functions that are unknown to
 *  the integrator have to be maintained by the the modules themselves.)
 *
 *  The system provides the coordinates `t` and `r` as global grid functions
 *  (they are defined in @ref fld::sysIndex). The integrator keeps time `t` evolved and
 *  the spatial coordinate `r` constant.
 *
 *  The integrator handles process signals (e.g., Ctrl-C), which will cause
 *  a premature end of the integration.
 */
class MoL : GridUser
{
    std::vector<IntegFace*> eomList;       //!< Equations of motion that we evolve
    std::vector<Int> constantGF;           //!< Grid functions that are kept constant
    std::vector<fld::EvolvedBy> evolvedGF; //!< Grid functions that are integrated in time

    GridOutputWriter* output; //!< The integration output is directed here

    time_hr  rt_start;        //!< The realtime when the integration started
    Int      method;          //!< The integration method (Euler, ICN, MoL, ...)
    bool     running;         //!< Indicates if the integration should proceed
    Real     t_0;             //!< The initial `t`.
    Real     t_1;             //!< Integrate up to this `t`.
    Real     delta_t;         //!< The integration step
    Real     dissip;          //!< Kreiss-Oliger dissipation coefficient.
    Real     dissip_delta_r;  //!< `dissip / delta_r`
    Real     cur_t;           //!< Current time

    /** A list of the implemented integration methods.
     */
    static std::map<std::string,int> knownMethods;

    /** Keep the list of all created integrator objects (used by signalHandler).
     */
    static std::vector<MoL*> knownIntegrators;

    /** Forward the signal to all integrators (signaling them to quit).
     */
    static void signalHandler( int signum )
    {
        std::cout << std::endl;
        std::cerr << "*** Signal " << signum << " received" << std::endl;

        if( knownIntegrators.size () == 0 ) {
            exit( -1 );
        }
        for( auto i: knownIntegrators ) {
            i->quit ();
        }
    }

    /** Report the integration time and the estimated real-time to finish.
     */
    void reportIntegrationTime( Int m )
    {
        if ( mpiRank() > 1 ) {
            return; // only the first in rank can access the standard output
        }

        Real tnow = GF( fld::t, m, nGhost ); // current integration time

        auto elapsed = 1e-3 * std::chrono::duration_cast<std::chrono::milliseconds>
                ( std::chrono::high_resolution_clock::now () - rt_start ).count ();

        auto rt_factor = elapsed / ( tnow - t_0 );  // real time/integration time
        auto rt_total = round( ( t_1 - t_0 ) * rt_factor ); // projected total time
        auto rt_left = round( ( t_1 - tnow ) * rt_factor ); // estimated time to cover

        std::cout << "    t = " << std::setw(6) << std::setprecision(2) << std::fixed
             << tnow << std::setw(-1) << std::setprecision(-1) << std::defaultfloat
             << ",   real = " << round( elapsed )
             << " s   (" << round( 100 * ( tnow - t_0 ) / ( t_1 - t_0 ) )
             << "% of " << rt_total
             << " s),   ETA " << rt_left << " s"
             << "                        \r";
        std::flush( std::cout );
    }

    /** Executed at the beginning of each integration step.
     */
    void integStep_Begin( Int m )
    {
        /// - Prepare the intermediate variables.
        ///
        for( auto eom: eomList ) {
            eom->integStep_Prepare( m );
        }

        /// - Exchange boundaries with the neighboring MPI ranks and/or fix the
        ///   left/right boundary conditions (depending on our MPI rank).
        ///   Quit if MPI has failed or we received abort message.
        ///
        if( gridDriver->mpiSize() > 1
            && ! gridDriver->exchangeBoundaries( m ) )
        {
            running = false;
            return;
        }
        if ( gridDriver->isFirstInRank () ) {
            for( auto eom: eomList ) {
                eom->applyLeftBoundaryCondition( m );
            }
        }
        if ( gridDriver->isLastInRank () ) {
            for( auto eom: eomList ) {
                eom->applyRightBoundaryCondition( m );
            }
        }

        /// - Calculate the RHS of time evolution equations.
        ///
        for( auto eom: eomList ) {
            eom->integStep_CalcEvolutionRHS( m );
        }
    }

    /** Executed at the end of each integration step.
     *  Returns `false` if we need to break the integration loop.
     */
    void integStep_End( Int m1, Int m )
    {
        /// - Finalize the integration step (e.g., doing diagnostic and post
        ///   processing of EoM like calculating the constraint violations).
        ///
        for( auto eom: eomList ) {
            eom->integStep_Finalize( m1, m );
        }

        /// - Wait on finishing the exchange of the boundaries (when using MPI).
        ///
        if ( ! gridDriver->waitExchange () ) {
            quit( m ); // Quit if MPI has failed
        }
    }

    /** Executed at each checkpoint.
     */
    void integStep_Checkpoint( Int m )
    {
        /// - Report the elapsed (both integration and real) time.
        ///
        reportIntegrationTime( m );

        /// - Check for eventual NaNs in the central/lower region of the grid.
        ///
        for( auto eom: eomList )
        {
            if ( ! eom->integStep_Diagnostics( m,
                            output->get_nFrom (), output->get_nTo () ) ) {
                if ( running ) { // quit only once
                    quit( m );
                }
            }
        }

        /// - Send the grid slice to the output.
        ///
        output->write( m, cur_t, delta_t );
    }

    /** Evolves the state variables using the FT scheme.
     */
    void integStep_Euler( Int mNext, Int m );

    /** Evolves the sate variables using the ICN averaging.
     */
    void integStep_AvgICN( Int mNext, Int mMid, Int m );

    /** Time evolution using the iterated Crank-Nicolson (ICN).
     */
    void integrate_AvgICN( int ICN_steps );

    /** Prepares the state variables for the intermediate MoL steps.
     */
    void integStep_MoLInit( Int m1, Int m, Real beta_delta_t );

    /** The intermediate method of lines (MoL) steps.
     */
    void integStep_MoL( Int m1, Int m, Real alpha );

    /** Time evolution using the method of lines (MoL).
     */
    void integrate_MoL( const MoLDescriptor& MoL );

public:

    /** Constructs the MoL integrator as specified in the parameter file.
     */
    MoL( Parameters& params, UniformGrid& ug, GridOutputWriter& out )
        : GridUser( ug ), output( &out )
    {
        constantGF.reserve( 64 );
        evolvedGF.reserve( 128 );

        constantGF.push_back( fld::r ); // Radial coordinate is kept constant by default

        params.get( "integ.t_0",     t_0,      0.0 );
        params.get( "integ.t_1",     t_1,      1.0 );
        params.get( "integ.delta_t", delta_t,  0.1 );
        params.get( "integ.dissip",  dissip,   0.0 );

        std::string name = params.get( "integ.method", method, 0, knownMethods );

        // Fix the sign of delta_t
        //
        delta_t = t_1 < t_0 ? -abs(delta_t) : abs(delta_t);

        cur_t = t_0;  // The initial `t`
        running = false; // Stopped integration by default (will be enabled later on)

        // Calculate Kreiss-Oliger dissipation corrected delta r
        //
        dissip_delta_r = dissip / delta_r;

        // Expected output file data size (in bytes)
        //
        const Real sz = ( ( t_1 - t_0 ) / delta_t / output->get_mSkip() + 1 )
                           * output->recordSizeInBytes();

        slog << "Integrator:" << std::endl << std::endl
             << "    t_0 = " << t_0 << ",  t_1 = " << t_1 << ",  delta_t = " << delta_t
             << ",  method = " << name << " (#" << method << ")" << std::endl
             << "    Kreiss-Oliger dissipation = " << dissip << std::endl
             << "    Expected output data size = "
             << round( ( sz < 1e3 ? sz : sz < 1e6 ? sz/1e3 : sz/1e6 ) * 10 ) / 10
             << ( sz < 1e3 ? "" : sz < 1e6 ? " kB" : " MB" )
             << ",  GF count = " << output->GFCount ()
             << std::endl << std::endl;

        // Register to receive the signals from the system
        //
        knownIntegrators.push_back( this );
        if ( knownIntegrators.size () == 1 ) {
            signal( SIGINT, signalHandler );
            signal( SIGTERM, signalHandler );
        }
    }

    /** The destructor (here only responsible for cleaning up standard output).
     */
    ~MoL()
    {
        if( mpiRank() == 0 ) {
            std::cout << std::endl << std::endl;
        }
    }

    /** Moves the integration final time to the current time and marks
     *  the integration stopped.
     */
    void quit( Int m = -1 )
    {
        running = false; // Stop the integration

        t_1 = m < 0 ? cur_t : GF( fld::t, m, nGhost );

        if( mpiRank() == 0 ) std::cout << std::endl;
        slog << "*** Quitting at t = " << t_1 << std::endl;

        /// - In case of NaNs, signal the other MPI ranks that we are quitting.
        /// @todo Distribute the abort msg to all ranks instead of calling MPI_Abort.
        ///
        gridDriver->abortExchange ();
    }

    /** Adds the EoM to the list of all evolved subsystems.
     */
    void addToEvolution( IntegFace* eomInterface )
    {
        eomList.push_back( eomInterface );
    }

    /** The given fields will be kept constant by the integrator.
     */
    void keepConstant( const std::vector<Int>& gfs )
    {
        constantGF.insert( constantGF.end(), gfs.begin(), gfs.end() );
    }

    /** Add the given functions to the evolution list.
     */
    void keepEvolved( const std::vector<fld::EvolvedBy>& gfs )
    {
        evolvedGF.insert( evolvedGF.end(), gfs.begin(), gfs.end() );
    }

    /** Returns the time step.
     */
    Real dt () const {
        return delta_t;
    }

    /** Integrates the evolution equations in the given grid realm.
     */
    bool evolveEquations ()
    {
        #if 0
        if( delta_r != r(0,nGhost+1) - r(0,nGhost) ) {
            slog << "    Error: delta_r = " << delta_r << " != grid spacing = "
                 << r(0,nGhost+1) - r(0,nGhost)
                 << ", diff = " << delta_r - ( r(0,nGhost+1) - r(0,nGhost) )
                 << std::endl;
        }
        #endif

        output->writeHeader ();

        // Start timing the real-time
        rt_start = std::chrono::high_resolution_clock::now ();

        switch( method )
        {
            case 0:  integrate_AvgICN ( 0    );  break;
            case 1:  integrate_AvgICN ( 1    );  break;
            case 2:  integrate_AvgICN ( 2    );  break;
            case 3:  integrate_AvgICN ( 3    );  break;
            case 4:  integrate_MoL    ( ICN2 );  break;
            case 5:  integrate_MoL    ( ICN3 );  break;
            case 6:  integrate_MoL    ( RK1  );  break;
            case 7:  integrate_MoL    ( RK2  );  break;
            case 8:  integrate_MoL    ( RK3  );  break;
            case 9:  integrate_MoL    ( RK4  );  break;
        }

        return true;
    }
};

/////////////////////////////////////////////////////////////////////////////////////////
// Definitions of MoL static members
/////////////////////////////////////////////////////////////////////////////////////////

std::vector<MoL*> MoL::knownIntegrators;

std::map<std::string,int> MoL::knownMethods =
{
    { "Euler",   0 },  { "avgICN2", 1 },
    { "avgICN3", 2 },  { "avgICN4", 3 },
    { "ICN2",    4 },  { "ICN3",    5 },
    { "RK1",     6 },  { "RK2",     7 },
    { "RK3",     8 },  { "RK4",     9 }
};

/////////////////////////////////////////////////////////////////////////////////////////
// Definitions of MoL methods
/////////////////////////////////////////////////////////////////////////////////////////

void MoL::integStep_Euler( Int m1, Int m )
{
    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        GF( fld::t, m1, n ) = GF( fld::t, m, n ) + delta_t; // Evolve time

        for( auto f: constantGF )
        {
            GF( f, m1, n ) = GF( f, m, n );
        }

        for( auto e: evolvedGF )
        {
            GF( e.f, m1, n ) = GF( e.f, m, n ) + delta_t * GF( e.f_t, m, n )
                             - KreissOligerTerm( e.f, delta_t );
        }
    }
}

void MoL::integStep_AvgICN
(
    Int m2,  // final time
    Int m1,  // middle time
    Int m    // initial time
    )
{
    const Real half_dt = delta_t / 2;

    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        GF( fld::t, m2, n ) = GF( fld::t, m, n ) + delta_t; // Evolve time

        for( auto f: constantGF )
        {
            GF( f, m2, n ) = GF( f, m, n );
        }

        for( auto e: evolvedGF )
        {
            GF( e.f, m2, n ) = GF( e.f, m, n )
                    + half_dt * ( GF( e.f_t, m1, n) + GF( e.f_t, m, n ) )
                    - KreissOligerTerm( e.f, delta_t );
        }
    }
}

/** Integrates the evolution equations using the iterated Crank-Nicolson (ICN) method.
 *  If ICN_steps == 0, then the method effectively becomes FTCS.
 *  See @cite Press:2007nrec, @cite Baumgarte:2010nr
 */
void MoL::integrate_AvgICN( int ICN_steps )
{
    slog << "Employing Method of Lines, Iterative Crank-Nicolson (avg)"
         << ", N = " << ICN_steps << std::endl << std::endl;

    cur_t = t_0;  // The initial `t`
    running = true; // Enable integration

    for( Int n = 0; n < nTotal; ++n ) {
        GF( fld::t, 0, n ) = t_0;
    }

    Int  mStep = 0; // The step counter (used to filter checkpoints)
    while( running && abs(cur_t) < abs(t_1) )
    {
        for( Int m = 0; running && m < mLen && abs(cur_t) < abs(t_1); ++m )
        {
            Int mNext = m + 1 >= mLen ? 0 : m + 1;

            /////////////////////////////////////////////////////////////////////////////
            // Predictor. Use FTCS for the first step.
            // Note that this will be the only step, if ICN_steps == 0.
            //
            Int m_1 = ICN_steps >= 1 ? mLen + 1 : mNext;

            integStep_Begin ( m      );   if( ! running ) break;
            integStep_Euler ( m_1, m );
            integStep_End   ( m_1, m );

            /////////////////////////////////////////////////////////////////////////////
            // Corrector iterations. There will be no iterations, if ICN_steps == 0.
            //
            for( int i = 0; running && i < ICN_steps; ++i )
            {
                Int m_2 = ( i + 1 >= ICN_steps ? mNext : m_1 + 1 );

                integStep_Begin  ( m_1         );   if( ! running ) break;
                integStep_AvgICN ( m_2, m_1, m );
                integStep_End    ( m_2, m_1    );

                m_1 = m_2;
            }

            if( ( mStep++ % output->get_mSkip () ) == 0 ) {
                integStep_Checkpoint( m );
            }

            cur_t = GF( fld::t, m, nGhost ); // The current time step
        }
    }
}

void MoL::integStep_MoL( Int m1, Int m, Real ac )
{
    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        GF( fld::t, m1, n ) += ac * GF( fld::t, m, n ); // Evolve time

        for( auto e: evolvedGF )
        {
            GF( e.f, m1, n ) += ac * GF( e.f, m, n );
        }
    }
}

void MoL::integStep_MoLInit( Int m1, Int m, Real bc_delta_t )
{
    OMP_parallel_for( Int n = nGhost; n < nGhost + nLen; ++n )
    {
        // There is no explicit dependence on t, and t_t(m,n) = 1.
        // Also, r is constant along the `line`, i.e., r_t(m,n) = 0.
        //
        GF( fld::t, m1, n ) = bc_delta_t; // Evolve time

        for( auto f: constantGF )
        {
            GF( f, m1, n ) = GF( f, m, n );
        }

        for( auto e: evolvedGF )
        {
            GF( e.f, m1, n ) = bc_delta_t * GF( e.f_t, m, n )
                             - KreissOligerTerm( e.f, bc_delta_t );
        }
    }
}

/** The method of lines (MoL) separates the time integration from the rest of
 *  an evolution scheme. A `N`-step MoL method evolves the equation:
 *  <pre>
 *       partial_t Y = L[Y],  </pre>
 *
 *  using the following algorithm:
 *  <pre>
 *      Y^{(0)}   = Y^{m},
 *      Y^{(i+1)} = Sum_{j=0}^{i} ( alpha_{ij} Y^{(j)} ) +            <- integStep_MoL
 *                  Delta t beta_i L[Y^{(i)}],   for i = 0, ..., N-1, <- integStep_MoLInit
 *      Y^{m+1}   = Y^{(N)}.    </pre>
 *
 *  The variables `Y^{(i)}`, `i = 0, ..., N`, are intermediate.
 *
 *  The method is completely specified by `N`, and arrays alpha and beta
 *  in the structure MoLDescriptor.
 */
void MoL::integrate_MoL( const MoLDescriptor& MoL )
{
    slog << "Employing Method of Lines, " << MoL.name
         << ", N = " << MoL.N << std::endl << std::endl;

    cur_t = t_0;  // The initial `t`
    running = true; // Enable integration

    for( Int n = 0; n < nTotal; ++n ) {
        GF( fld::t, 0, n ) = t_0;
    }

    // The array mk[] contains redirections to grid rows (m's):
    //   mk[0]                   points to Y^{m}   == Y^{(0)},
    //   mk[1], ..., mk[Mol.N-1] point to the intermediate Y^{(1)}, ..., Y^{(N-1)},
    //   mk[MoL.N]               points to Y^{m+1} == Y^{(N)}.
    //
    int mk[ 10 ]; // int mk[ MoL.N + 1 ]; // Note that N is at most 8
    for( int i = 0; i <= MoL.N; ++i ) {
        mk[i] = mLen + i;
    }

    Int  mStep = 0; // The step counter
    while( running && abs(cur_t) < abs(t_1) )
    {
        for( Int m = 0; running && m < mLen && abs(cur_t) < abs(t_1); ++m )
        {
            mk[0] = m; // Y^{(m)}
            mk[MoL.N] = m + 1 >= mLen ? 0 : m + 1; // Y^{(m+1)}

            for( int i = 0; running && i < MoL.N; ++i )
            {
                integStep_Begin( mk[i] );   if( ! running ) break;

                integStep_MoLInit( mk[i+1], mk[i], MoL.beta[i] * delta_t );

                for( int j = 0; j <= i; j++ ) {
                    integStep_MoL( mk[i+1], mk[j], MoL.alpha[i][j] );
                }

                integStep_End( mk[i+1], mk[i] );
            }

            if( ( mStep++ % output->get_mSkip () ) == 0 ) {
                integStep_Checkpoint( m );
            }

            cur_t = GF( fld::t, m, nGhost ); // The current time step
        }
    }
}
                                                                                   /*@}*/
/////////////////////////////////////////////////////////////////////////////////////////