Nektar::SolverUtils::DriverParareal Class Reference

Base class for the development of solvers. More...

#include <DriverParareal.h>

Inheritance diagram for Nektar::SolverUtils::DriverParareal:

Static Public Member Functions
static DriverSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className
	Name of the class. More...

Protected Member Functions
SOLVER_UTILS_EXPORT	DriverParareal (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Constructor. More...
virtual SOLVER_UTILS_EXPORT	~DriverParareal ()
	Destructor. More...
virtual SOLVER_UTILS_EXPORT void	v_InitObject (std::ostream &out=std::cout) override
	Second-stage initialisation. More...
virtual SOLVER_UTILS_EXPORT void	v_Execute (std::ostream &out=std::cout) override
	Virtual function for solve implementation. More...
void	SetPararealSessionFile (void)
	Set the Parareal (coarse solver) session file. More...
void	RunCoarseSolve (const NekDouble time, const int nstep, const int iter, const Array< OneD, const Array< OneD, NekDouble >> &input, Array< OneD, Array< OneD, NekDouble >> &output)
void	RunFineSolve (const NekDouble time, const int nstep, const int iter, const Array< OneD, const Array< OneD, NekDouble >> &input, Array< OneD, Array< OneD, NekDouble >> &output)
Protected Member Functions inherited from Nektar::SolverUtils::Driver
	Driver (const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
	Initialises EquationSystem class members. More...
virtual SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	v_GetRealEvl (void)
virtual SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	v_GetImagEvl (void)

Protected Attributes
FieldUtils::Interpolator< Array< OneD, MultiRegions::ExpListSharedPtr > >	m_interp
LibUtilities::SessionReaderSharedPtr	m_sessionCoarse
	Parareal (coarse solver) session reader object. More...
SpatialDomains::MeshGraphSharedPtr	m_graphCoarse
	Parareal (coarse solver) MeshGraph object. More...
NekDouble	m_fineTimeStep
	Timestep for fine solver. More...
NekDouble	m_coarseTimeStep
	Timestep for coarse solver. More...
NekDouble	m_totalTime
	Total time integration interval. More...
NekDouble	m_chunkTime
	Time for chunks. More...
NekDouble	m_coarseSolveFactor = 100.0
	Coarse solver time factor. More...
int	m_fineSteps = 1
	Number of steps for the fine solver. More...
int	m_coarseSteps = 1
	Number of steps for the coarse solver. More...
int	m_numChunks = 1
	Number of time chunks. More...
int	m_chunkRank = 0
	Rank in time. More...
int	m_pararealIterMax = 0
	Maximum number of parareal iteration. More...
bool	m_exactSolution = 0
	Using exact solution to compute error norms. More...
NekDouble	m_pararealToler = 1e-15
	Parareal tolerance. More...
Protected Attributes inherited from Nektar::SolverUtils::Driver
LibUtilities::CommSharedPtr	m_comm
	Communication object. More...
LibUtilities::SessionReaderSharedPtr	m_session
	Session reader object. More...
LibUtilities::SessionReaderSharedPtr	session_LinNS
	I the Coupling between SFD and arnoldi. More...
SpatialDomains::MeshGraphSharedPtr	m_graph
	MeshGraph object. More...
Array< OneD, EquationSystemSharedPtr >	m_equ
	Equation system to solve. More...
int	m_nequ
	number of equations More...
enum EvolutionOperatorType	m_EvolutionOperator
	Evolution Operator. More...

Static Protected Attributes
static std::string	driverLookupId
Static Protected Attributes inherited from Nektar::SolverUtils::Driver
static std::string	evolutionOperatorLookupIds []
static std::string	evolutionOperatorDef
static std::string	driverDefault

Friends
class	MemoryManager< DriverParareal >

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual	~Driver ()
	Destructor. More...
SOLVER_UTILS_EXPORT void	InitObject (std::ostream &out=std::cout)
	Initialise Object. More...
SOLVER_UTILS_EXPORT void	Execute (std::ostream &out=std::cout)
	Execute driver. More...
SOLVER_UTILS_EXPORT Array< OneD, EquationSystemSharedPtr >	GetEqu ()
SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	GetRealEvl (void)
SOLVER_UTILS_EXPORT Array< OneD, NekDouble >	GetImagEvl (void)

Detailed Description

Base class for the development of solvers.

Definition at line 47 of file DriverParareal.h.

Constructor & Destructor Documentation

DriverParareal()

Nektar::SolverUtils::DriverParareal::DriverParareal ( const LibUtilities::SessionReaderSharedPtr  pSession, const SpatialDomains::MeshGraphSharedPtr  pGraph  )

protected

Constructor.

Definition at line 54 of file DriverParareal.cpp.

    : Driver(pSession, pGraph)
{
}

~DriverParareal()

Nektar::SolverUtils::DriverParareal::~DriverParareal ( )

protectedvirtual

Destructor.

Definition at line 64 of file DriverParareal.cpp.

{
}

Member Function Documentation

create()

static DriverSharedPtr Nektar::SolverUtils::DriverParareal::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const SpatialDomains::MeshGraphSharedPtr &  pGraph  )

inlinestatic

Creates an instance of this class.

Definition at line 53 of file DriverParareal.h.

    {
        DriverSharedPtr p =
            MemoryManager<DriverParareal>::AllocateSharedPtr(pSession, pGraph);
        p->InitObject();
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

RunCoarseSolve()

void Nektar::SolverUtils::DriverParareal::RunCoarseSolve ( const NekDouble  time, const int  nstep, const int  iter, const Array< OneD, const Array< OneD, NekDouble >> &  input, Array< OneD, Array< OneD, NekDouble >> &  output  )

protected

Definition at line 279 of file DriverParareal.cpp.

{
    // Output coarse timestep.
    if (m_chunkRank == iter && m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "RUNNING COARSE SOLVE: dt = " << m_coarseTimeStep
                  << " nsteps = " << m_coarseSteps / m_numChunks << std::endl
                  << std::flush;
    }
 
    // Set to coarse timestep.
    m_equ[1]->SetTime(time);
    m_equ[1]->SetSteps(nstep);
 
    // Copy initial condition from input.
    if (m_equ[0]->GetNpoints() == m_equ[1]->GetNpoints())
    {
        // Interpolation not necessary, directly copy data.
        for (int i = 0; i < m_equ[1]->GetNvariables(); ++i)
        {
            m_equ[1]->CopyToPhysField(i, input[i]);
        }
    }
    else
    {
        // Copy data to fine solver and interpolate from fine to coarse (WIP).
        for (int i = 0; i < m_equ[0]->GetNvariables(); ++i)
        {
            m_equ[0]->CopyToPhysField(i, input[i]);
        }
        m_interp.InterpExp1ToExp2(m_equ[0]->UpdateFields(),
                                  m_equ[1]->UpdateFields());
    }
 
    // Solve equations.
    m_equ[1]->DoSolve();
 
    // Copy solution to output.
    if (m_equ[0]->GetNpoints() == m_equ[1]->GetNpoints())
    {
        // Interpolation not necessary, directly copy data.
        for (int i = 0; i < m_equ[1]->GetNvariables(); ++i)
        {
            m_equ[1]->CopyFromPhysField(i, output[i]);
        }
    }
    else
    {
        // Copy data from coarse solver and interpolate coarse to fine (WIP).
        m_interp.InterpExp1ToExp2(m_equ[1]->UpdateFields(),
                                  m_equ[0]->UpdateFields());
        for (int i = 0; i < m_equ[1]->GetNvariables(); ++i)
        {
            m_equ[0]->CopyFromPhysField(i, output[i]);
        }
    }
}

References m_chunkRank, m_coarseSteps, m_coarseTimeStep, Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, m_interp, and m_numChunks.

Referenced by v_Execute().

RunFineSolve()

void Nektar::SolverUtils::DriverParareal::RunFineSolve ( const NekDouble  time, const int  nstep, const int  iter, const Array< OneD, const Array< OneD, NekDouble >> &  input, Array< OneD, Array< OneD, NekDouble >> &  output  )

protected

Definition at line 340 of file DriverParareal.cpp.

{
    // Output fine timestep.
    if (m_chunkRank == iter && m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "RUNNING FINE SOLVE: dt = " << m_fineTimeStep
                  << " nsteps = " << m_fineSteps / m_numChunks << std::endl
                  << std::flush;
    }
 
    // Number of checkpoint by chunk.
    int nChkPts =
        m_session->GetParameter("IO_CheckSteps")
            ? m_fineSteps /
                  int(m_session->GetParameter("IO_CheckSteps") * m_numChunks)
            : 1;
 
    // Parareal iteration number.
    int nIter = m_equ[0]->GetPararealIterationNumber();
 
    // Set to fine timestep.
    m_equ[0]->SetTime(time);
    m_equ[0]->SetSteps(nstep);
 
    // Reinitialize check point number for each parareal iteration.
    m_equ[0]->SetCheckpointNumber(m_chunkRank * nChkPts + 1);
 
    // Update parareal iteration number.
    m_equ[0]->SetPararealIterationNumber(++nIter);
 
    // Copy initial condition from input.
    for (int i = 0; i < m_equ[0]->GetNvariables(); ++i)
    {
        m_equ[0]->CopyToPhysField(i, input[i]);
    }
 
    // Solve equations.
    m_equ[0]->DoSolve();
 
    // Copy solution to output.
    for (int i = 0; i < m_equ[0]->GetNvariables(); ++i)
    {
        m_equ[0]->CopyFromPhysField(i, output[i]);
    }
}

References m_chunkRank, Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, m_fineSteps, m_fineTimeStep, m_numChunks, and Nektar::SolverUtils::Driver::m_session.

Referenced by v_Execute().

SetPararealSessionFile()

void Nektar::SolverUtils::DriverParareal::SetPararealSessionFile ( void  )

protected

Set the Parareal (coarse solver) session file.

Definition at line 180 of file DriverParareal.cpp.

{
    // Get the coarse solver session file.
    string meshFile;
    string coarseSolverFile;
    vector<string> coarseSolverFilenames;
    bool opt         = (m_session->GetFilenames()[0].substr(
                    m_session->GetFilenames()[0].size() - 3) == "opt");
    meshFile         = m_session->GetFilenames()[0 + opt];
    coarseSolverFile = m_session->GetFilenames().size() > 1 + opt
                           ? m_session->GetFilenames()[1 + opt]
                           : m_session->GetFilenames()[0 + opt];
    coarseSolverFile = coarseSolverFile.substr(0, coarseSolverFile.size() - 4);
    coarseSolverFile += "_coarseSolver.xml";
    std::ifstream f(coarseSolverFile);
 
    if (f.good())
    {
        // if _coarseSolver.xml exit, read session file
        if (m_session->GetFilenames().size() > 1 + opt)
        {
            coarseSolverFilenames.push_back(meshFile);
        }
        coarseSolverFilenames.push_back(coarseSolverFile);
    }
    else
    {
        // if _coarseSolver.xml does not exit, use original session file
        coarseSolverFilenames.push_back(m_session->GetFilenames()[0 + opt]);
        if (m_session->GetFilenames().size() > 1 + opt)
        {
            coarseSolverFilenames.push_back(m_session->GetFilenames()[1 + opt]);
        }
    }
 
    // Define argument for the coarse parareal solver.
    int npx = m_session->DefinesCmdLineArgument("npx")
                  ? m_session->GetCmdLineArgument<int>("npx")
                  : 1;
    int npy = m_session->DefinesCmdLineArgument("npy")
                  ? m_session->GetCmdLineArgument<int>("npy")
                  : 1;
    int npz = m_session->DefinesCmdLineArgument("npz")
                  ? m_session->GetCmdLineArgument<int>("npz")
                  : 1;
    int nsz = m_session->DefinesCmdLineArgument("nsz")
                  ? m_session->GetCmdLineArgument<int>("nsz")
                  : 1;
    int npt = m_session->DefinesCmdLineArgument("npt")
                  ? m_session->GetCmdLineArgument<int>("npt")
                  : 1;
 
    // Convert into string.
    std::string npx_string = std::to_string(npx);
    std::string npy_string = std::to_string(npy);
    std::string npz_string = std::to_string(npz);
    std::string nsz_string = std::to_string(nsz);
    std::string npt_string = std::to_string(npt);
 
    char *argv[] = {const_cast<char *>("Solver"), // this is just a place holder
                    const_cast<char *>("--npx"),
                    const_cast<char *>(npx_string.c_str()),
                    const_cast<char *>("--npy"),
                    const_cast<char *>(npy_string.c_str()),
                    const_cast<char *>("--npz"),
                    const_cast<char *>(npz_string.c_str()),
                    const_cast<char *>("--nsz"),
                    const_cast<char *>(nsz_string.c_str()),
                    const_cast<char *>("--npt"),
                    const_cast<char *>(npt_string.c_str()),
                    nullptr};
 
    // Set session for coarse solver.
    m_sessionCoarse = LibUtilities::SessionReader::CreateInstance(
        11, argv, coarseSolverFilenames, m_session->GetComm());
 
    // Set graph for coarse solver.
    m_graphCoarse = SpatialDomains::MeshGraph::Read(m_sessionCoarse);
 
    // Set BndRegionOrdering (necessary for DG with periodic BC) FIXME
    m_graphCoarse->SetBndRegionOrdering(m_graph->GetBndRegionOrdering());
 
    // Set CompositeOrdering (necessary for DG with periodic BC) FIXME
    m_graphCoarse->SetCompositeOrdering(m_graph->GetCompositeOrdering());
 
    // If a coarse solver session file is not specified, use
    // m_coarseSolveFactor to determine the timestep of the coarse solver
    // (default value is 100.0).
    if (!f.good())
    {
        double timeStep =
            m_session->GetParameter("TimeStep") * m_coarseSolveFactor;
        int numSteps =
            m_session->GetParameter("NumSteps") / m_coarseSolveFactor;
        m_sessionCoarse->SetParameter("TimeStep", timeStep);
        m_sessionCoarse->SetParameter("NumSteps", numSteps);
    }
}

References Nektar::LibUtilities::SessionReader::CreateInstance(), m_coarseSolveFactor, Nektar::SolverUtils::Driver::m_graph, m_graphCoarse, Nektar::SolverUtils::Driver::m_session, m_sessionCoarse, and Nektar::SpatialDomains::MeshGraph::Read().

Referenced by v_InitObject().

v_Execute()

void Nektar::SolverUtils::DriverParareal::v_Execute ( std::ostream &  out = std::cout )

overrideprotectedvirtual

Virtual function for solve implementation.

Implements Nektar::SolverUtils::Driver.

Definition at line 389 of file DriverParareal.cpp.

{
    Nektar::LibUtilities::Timer timer;
    NekDouble CPUtime = 0.0;
 
    m_numChunks = m_session->GetComm()->GetTimeComm()->GetSize();
    m_chunkRank = m_session->GetComm()->GetTimeComm()->GetRank();
 
    // Set parameters from session file.
    m_pararealToler   = m_session->DefinesParameter("PararealToler")
                            ? m_session->GetParameter("PararealToler")
                            : m_pararealToler;
    m_pararealIterMax = m_session->DefinesParameter("PararealIterMax")
                            ? m_session->GetParameter("PararealIterMax")
                            : m_numChunks;
    m_fineTimeStep    = m_session->GetParameter("TimeStep");
    m_coarseTimeStep  = m_sessionCoarse->GetParameter("TimeStep");
    m_fineSteps       = m_session->GetParameter("NumSteps");
    m_coarseSteps     = m_sessionCoarse->GetParameter("NumSteps");
    m_totalTime       = m_fineTimeStep * m_fineSteps;
    m_chunkTime       = m_totalTime / m_numChunks;
    m_exactSolution   = m_session->DefinesParameter("ExactSolution")
                            ? m_session->GetParameter("ExactSolution")
                            : m_exactSolution;
 
    // Turnoff I/O for coarse solver.
    m_equ[1]->SetInfoSteps(0);
    m_equ[1]->SetCheckpointSteps(0);
 
    // Check time step inputs.
    ASSERTL0(
        m_fineSteps % m_numChunks == 0,
        "Total number of fine step should be divisible by number of chunks.");
 
    ASSERTL0(
        m_coarseSteps % m_numChunks == 0,
        "Total number of coarse step should be divisible by number of chunks.");
 
    ASSERTL0(fabs(m_coarseTimeStep * m_coarseSteps -
                  m_fineTimeStep * m_fineSteps) < 1e-12,
             "Fine and coarse total computational times do not match");
 
    // Check I/O inputs.
    if (m_session->GetParameter("IO_InfoSteps"))
    {
        ASSERTL0(m_fineSteps % int(m_session->GetParameter("IO_InfoSteps") *
                                   m_numChunks) ==
                     0,
                 "number of IO_InfoSteps should divide number of fine steps "
                 "per time chunk");
    }
    if (m_session->GetParameter("IO_CheckSteps"))
    {
        ASSERTL0(m_fineSteps % int(m_session->GetParameter("IO_CheckSteps") *
                                   m_numChunks) ==
                     0,
                 "number of IO_CheckSteps should divide number of fine steps "
                 "per time chunk");
    }
 
    // Fine solver summary.
    if (m_chunkRank == 0 && m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "========================================================="
                     "=============="
                  << std::endl
                  << std::flush;
        std::cout << "========================= FINE PROPAGATOR INFO "
                     "========================"
                  << std::endl
                  << std::flush;
        m_equ[0]->PrintSummary(out);
 
        std::cout << std::endl << std::flush;
    }
 
    // Coarse solver summary.
    if (m_chunkRank == 0 && m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "========================================================="
                     "=============="
                  << std::endl
                  << std::flush;
        std::cout << "======================== COARSE PROPAGATOR INFO "
                     "======================="
                  << std::endl
                  << std::flush;
        m_equ[1]->PrintSummary(out);
    }
 
    timer.Start();
 
    // Allocate storage for parareal solver.
    int nPts = m_equ[0]->GetNpoints();
    int nVar = m_equ[0]->GetNvariables();
    Array<OneD, Array<OneD, NekDouble>> solution(nVar);
    Array<OneD, Array<OneD, NekDouble>> solutionCoarsePrev(nVar);
    Array<OneD, Array<OneD, NekDouble>> solutionCoarseCurr(nVar);
    Array<OneD, Array<OneD, NekDouble>> ic(nVar);
    Array<OneD, Array<OneD, NekDouble>> exactsoln(nVar);
    for (int i = 0; i < nVar; ++i)
    {
        solution[i]           = Array<OneD, NekDouble>(nPts, 0.0);
        solutionCoarsePrev[i] = Array<OneD, NekDouble>(nPts, 0.0);
        solutionCoarseCurr[i] = Array<OneD, NekDouble>(nPts, 0.0);
        ic[i]                 = Array<OneD, NekDouble>(nPts, 0.0);
        exactsoln[i]          = Array<OneD, NekDouble>(nPts, 0.0);
    }
 
    // Initialize fine solver.
    m_equ[0]->DoInitialise();
 
    // Initialize coarse solver.
    m_equ[1]->SetUseInitialCondition(false);
    m_equ[1]->DoInitialise();
 
    // Get initial conditions.
    for (int i = 0; i < nVar; ++i)
    {
        m_equ[0]->CopyFromPhysField(i, ic[i]);
    }
 
    // Run coarse solution, G(y_j^k) to get initial conditions.
    if (m_chunkRank == 0 && m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "** INITIAL CONDITION **" << std::endl << std::flush;
    }
    LibUtilities::CommSharedPtr tComm = m_session->GetComm()->GetTimeComm();
    int recvProc                      = m_chunkRank - 1;
    int sendProc                      = m_chunkRank + 1;
 
    // Iterate to improve on the approximation
    // For the moment, we use the maximum number of iterations
    // We can add a tolerance threshold later.
    // We start the iteration with the current approximation stored in
    // the 'solution' array.
    int k               = 0;
    int kmax            = 0;
    int convergenceCurr = false;
    int convergencePrev = (m_chunkRank == 0);
    NekDouble vL2Error  = 0.0;
    Array<OneD, NekDouble> vL2Errors(nVar, 0.0);
    Array<OneD, NekDouble> vLinfErrors(nVar, 0.0);
 
    // To compute the initial conditions, the coarse solver is run serially
    // for each time chunk to avoid communication time between processors.
    if (m_chunkRank > 0)
    {
        RunCoarseSolve(0.0, m_chunkRank * m_coarseSteps / m_numChunks, -1, ic,
                       ic);
    }
 
    // Compute initial coarse solution.
    if (m_chunkRank == 0 && m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "** ITERATION " << 0 << " **" << std::endl << std::flush;
    }
 
    RunCoarseSolve(m_chunkRank * m_chunkTime, m_coarseSteps / m_numChunks, 0,
                   ic, solutionCoarsePrev);
 
    for (int i = 0; i < nVar; ++i)
    {
        for (int q = 0; q < nPts; ++q)
        {
            solution[i][q] = solutionCoarsePrev[i][q];
        }
    }
 
    if (m_exactSolution)
    {
        // Evaluate exact solution.
        for (int i = 0; i < nVar; ++i)
        {
            m_equ[0]->EvaluateExactSolution(i, exactsoln[i],
                                            (m_chunkRank + 1) * m_chunkTime);
        }
    }
 
    timer.Stop();
    CPUtime += timer.Elapsed().count();
    if (m_chunkRank == m_numChunks - 1 &&
        m_comm->GetSpaceComm()->GetRank() == 0)
    {
        std::cout << "-------------------------------------------" << std::endl
                  << std::flush;
        std::cout << "Total Computation Time = " << CPUtime << "s" << std::endl
                  << std::flush;
        std::cout << "-------------------------------------------" << std::endl
                  << std::flush;
    }
    for (int i = 0; i < nVar; ++i)
    {
        vL2Errors[i]   = m_equ[1]->L2Error(i, exactsoln[i], 1);
        vLinfErrors[i] = m_equ[1]->LinfError(i, exactsoln[i]);
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "L2 error (variable " << m_equ[1]->GetVariable(i)
                      << ") : " << vL2Errors[i] << std::endl
                      << std::flush;
            std::cout << "Linf error (variable " << m_equ[1]->GetVariable(i)
                      << ") : " << vLinfErrors[i] << std::endl
                      << std::flush;
        }
    }
    timer.Start();
 
    // Start Parareal iteration loop.
    while (k < m_pararealIterMax && !convergenceCurr)
    {
        if (m_chunkRank == min(k, m_numChunks - 1) &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "** ITERATION " << k + 1 << " **" << std::endl
                      << std::flush;
        }
 
        // Use previous parareal solution as "exact solution".
        if (!m_exactSolution)
        {
            for (int i = 0; i < nVar; ++i)
            {
                for (int q = 0; q < nPts; ++q)
                {
                    exactsoln[i][q] = solution[i][q];
                }
            }
        }
 
        // Calculate fine solution, F(y_j^k)
        // Again no communication necessary.
        RunFineSolve(m_chunkRank * m_chunkTime, m_fineSteps / m_numChunks, k,
                     ic, solution);
 
        // Calculate coarse solve correction G(y_j^{k+1})
        // These are dependent on the previous time slice, so need to
        // compute serially.
        if (m_chunkRank > 0 && !convergencePrev)
        {
            // All time slices, apart from the first, receive their initial
            // state from the previous time slice.
            tComm->Recv(recvProc, convergencePrev);
            for (int i = 0; i < nVar; ++i)
            {
                tComm->Recv(recvProc, ic[i]);
            }
        }
 
        // Run the coarse solver
        RunCoarseSolve(m_chunkRank * m_chunkTime, m_coarseSteps / m_numChunks,
                       k, ic, solutionCoarseCurr);
 
        // Calculate the new approximation y_{j+1}^{k+1}
        // This is calculated point-wise.
        for (int i = 0; i < nVar; ++i)
        {
            for (int q = 0; q < nPts; ++q)
            {
                solution[i][q] +=
                    solutionCoarseCurr[i][q] - solutionCoarsePrev[i][q];
 
                solutionCoarsePrev[i][q] = solutionCoarseCurr[i][q];
            }
        }
 
        // Compute L2 error for each time chunk.
        vL2Error = 0.0;
        for (int i = 0; i < nVar; ++i)
        {
            // Copy the new approximation back.
            m_equ[0]->CopyToPhysField(i, solution[i]);
 
            vL2Error = max(vL2Error, m_equ[0]->L2Error(i, exactsoln[i], 1));
        }
 
        // Check convergence of L2 error for each time chunk.
        if ((vL2Error < m_pararealToler && convergencePrev) || m_chunkRank == k)
        {
            convergenceCurr = true;
        }
 
        // All but the last time slice should communicate the solution to
        // the next time slice. This will become the initial condition for
        // the next slice.
        if (m_chunkRank < m_numChunks - 1)
        {
            tComm->Send(sendProc, convergenceCurr);
            for (int i = 0; i < nVar; ++i)
            {
                tComm->Send(sendProc, solution[i]);
            }
        }
 
        if (m_exactSolution)
        {
            // Evaluate exact solution.
            for (int i = 0; i < nVar; ++i)
            {
                m_equ[0]->EvaluateExactSolution(
                    i, exactsoln[i], (m_chunkRank + 1) * m_chunkTime);
            }
        }
 
        timer.Stop();
        CPUtime += timer.Elapsed().count();
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Total Computation Time = " << CPUtime << "s"
                      << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
        for (int i = 0; i < nVar; ++i)
        {
            vL2Errors[i]   = m_equ[0]->L2Error(i, exactsoln[i], 1);
            vLinfErrors[i] = m_equ[0]->LinfError(i, exactsoln[i]);
            if (m_chunkRank == m_numChunks - 1 &&
                m_comm->GetSpaceComm()->GetRank() == 0)
            {
                std::cout << "L2 error (variable " << m_equ[0]->GetVariable(i)
                          << ") : " << vL2Errors[i] << std::endl
                          << std::flush;
                std::cout << "Linf error (variable " << m_equ[0]->GetVariable(i)
                          << ") : " << vLinfErrors[i] << std::endl
                          << std::flush;
            }
        }
        timer.Start();
 
        // Increment counter.
        k++;
        kmax = k;
    }
    timer.Stop();
 
    // If already converged, simply copy previous computed solution.
    m_comm->GetTimeComm()->AllReduce(kmax, Nektar::LibUtilities::ReduceMax);
    for (; k < kmax; k++)
    {
        if (m_comm->GetSpaceComm()->GetRank() == 0 &&
            m_session->GetParameter("IO_CheckSteps"))
        {
            std::string olddir = m_equ[0]->GetSessionName() + "_" +
                                 boost::lexical_cast<std::string>(k) + ".pit";
            std::string outdir = m_equ[0]->GetSessionName() + "_" +
                                 boost::lexical_cast<std::string>(k + 1) +
                                 ".pit";
 
            int nChkPts =
                m_fineSteps /
                int(m_session->GetParameter("IO_CheckSteps") * m_numChunks);
            for (int i = 0; i < nChkPts; i++)
            {
                // Old file name
                std::string oldname =
                    olddir + "/" + m_equ[0]->GetSessionName() + "_" +
                    boost::lexical_cast<std::string>(m_chunkRank * nChkPts + i +
                                                     1) +
                    ".chk";
 
                // New file name
                std::string outname =
                    outdir + "/" + m_equ[0]->GetSessionName() + "_" +
                    boost::lexical_cast<std::string>(m_chunkRank * nChkPts + i +
                                                     1) +
                    ".chk";
 
                // Remove file if already existing
                fs::remove_all(outname);
 
                // Copy from previous converged solution
                fs::copy(oldname, outname);
            }
        }
    }
 
    // Update solution before printing restart solution.
    for (int i = 0; i < nVar; ++i)
    {
        m_equ[0]->CopyToPhysField(i, solution[i]);
        m_equ[0]->UpdateFields()[i]->FwdTransLocalElmt(
            m_equ[0]->UpdateFields()[i]->GetPhys(),
            m_equ[0]->UpdateFields()[i]->UpdateCoeffs());
    }
 
    // Print restart solution.
    m_equ[0]->Output();
 
    // Wait for all processors to finish their writing activities
    m_comm->Block();
 
    // Print total computational time.
    if (m_chunkRank == m_numChunks - 1)
    {
        if (m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Total Computation Time = " << CPUtime << "s"
                      << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
 
        // Print solution errors.
        for (int i = 0; i < nVar; ++i)
        {
            // Copy the new approximation back.
            m_equ[0]->CopyToPhysField(i, solution[i]);
 
            // Evaluate exact solution.
            m_equ[0]->EvaluateExactSolution(i, exactsoln[i], m_totalTime);
 
            // Evaluate error norms.
            NekDouble vL2Error   = m_equ[0]->L2Error(i, exactsoln[i]);
            NekDouble vLinfError = m_equ[0]->LinfError(i, exactsoln[i]);
 
            if (m_comm->GetSpaceComm()->GetRank() == 0)
            {
                std::cout << "L 2 error (variable " << m_equ[0]->GetVariable(i)
                          << ") : " << vL2Error << std::endl
                          << std::flush;
                std::cout << "L inf error (variable "
                          << m_equ[0]->GetVariable(i) << ") : " << vLinfError
                          << std::endl
                          << std::flush;
            }
        }
    }
 
    // Speed-up analysis
    if (true)
    {
        int count = 20;
 
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "PARAREAL SPEED-UP ANALYSIS" << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
 
        // Mean communication time.
        NekDouble commTime = 0.0;
        timer.Start();
        for (int n = 0; n < count; n++)
        {
            if (m_chunkRank > 0)
            {
                for (int i = 0; i < nVar; ++i)
                {
                    tComm->Recv(recvProc, ic[i]);
                }
            }
 
            if (m_chunkRank < m_numChunks - 1)
            {
                for (int i = 0; i < nVar; ++i)
                {
                    tComm->Send(sendProc, solution[i]);
                }
            }
        }
        timer.Stop();
        commTime = timer.Elapsed().count() / count;
 
        // Print communication time.
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Mean Communication Time = " << commTime << "s"
                      << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
 
        // Mean restriction time.
        NekDouble restTime = 0.0;
        if (m_equ[0]->GetNpoints() != m_equ[1]->GetNpoints())
        {
            timer.Start();
            for (int n = 0; n < count; n++)
            {
                m_interp.InterpExp1ToExp2(m_equ[0]->UpdateFields(),
                                          m_equ[1]->UpdateFields());
            }
            timer.Stop();
            restTime = timer.Elapsed().count() / count;
 
            // Print restriction time.
            if (m_chunkRank == m_numChunks - 1 &&
                m_comm->GetSpaceComm()->GetRank() == 0)
            {
                std::cout << "-------------------------------------------"
                          << std::endl
                          << std::flush;
                std::cout << "Mean Restriction Time = " << restTime << "s"
                          << std::endl
                          << std::flush;
                std::cout << "-------------------------------------------"
                          << std::endl
                          << std::flush;
            }
        }
 
        // Mean interpolation time.
        NekDouble interTime = 0.0;
        if (m_equ[0]->GetNpoints() != m_equ[1]->GetNpoints())
        {
            timer.Start();
            for (int n = 0; n < count; n++)
            {
                m_interp.InterpExp1ToExp2(m_equ[1]->UpdateFields(),
                                          m_equ[0]->UpdateFields());
            }
            timer.Stop();
            interTime = timer.Elapsed().count() / count;
 
            // Print restriction time.
            if (m_chunkRank == m_numChunks - 1 &&
                m_comm->GetSpaceComm()->GetRank() == 0)
            {
                std::cout << "-------------------------------------------"
                          << std::endl
                          << std::flush;
                std::cout << "Mean Interpolation Time = " << interTime << "s"
                          << std::endl
                          << std::flush;
                std::cout << "-------------------------------------------"
                          << std::endl
                          << std::flush;
            }
        }
 
        // Mean coarse solver time.
        NekDouble coarseSolveTime = 0.0;
        timer.Start();
        for (int n = 0; n < count; n++)
        {
            RunCoarseSolve(0.0, 100, -1, ic, solution);
        }
        timer.Stop();
        coarseSolveTime = 0.01 * timer.Elapsed().count() / count *
                              (m_coarseSteps / m_numChunks) -
                          restTime - interTime;
 
        // Print restriction time.
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Mean Coarse Solve Time = " << coarseSolveTime << "s"
                      << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
 
        // Fine solver time.
        NekDouble fineSolveTime = 0.0;
        timer.Start();
        // Turnoff I/O for fine solver.
        m_equ[0]->SetInfoSteps(0);
        m_equ[0]->SetCheckpointSteps(0);
        for (int n = 0; n < count; n++)
        {
            RunFineSolve(0.0, 100, -1, ic, solution);
        }
        timer.Stop();
        fineSolveTime = 0.01 * timer.Elapsed().count() / count *
                        (m_fineSteps / m_numChunks);
 
        // Print fine solve time.
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Mean Fine Solve Time = " << fineSolveTime << "s"
                      << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
 
        // Print speedup time.
        if (m_chunkRank == m_numChunks - 1 &&
            m_comm->GetSpaceComm()->GetRank() == 0)
        {
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Maximum Speed-up" << std::endl << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            for (int k = 1; k <= m_numChunks; k++)
            {
                NekDouble ratio      = double(k) / m_numChunks;
                NekDouble ratioSolve = coarseSolveTime / fineSolveTime;
                NekDouble speedup = 1.0 / ((1.0 + ratio) * ratioSolve + ratio);
                std::cout << "Speed-up (" << k << ") = " << speedup << std::endl
                          << std::flush;
            }
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Speed-up with comm." << std::endl << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            for (int k = 1; k <= m_numChunks; k++)
            {
                NekDouble ratio      = double(k) / m_numChunks;
                NekDouble ratioComm  = commTime / fineSolveTime;
                NekDouble ratioSolve = coarseSolveTime / fineSolveTime;
                NekDouble speedup = 1.0 / ((1.0 + ratio) * ratioSolve + ratio +
                                           ratioComm / m_numChunks);
                std::cout << "Speed-up (" << k << ") = " << speedup << std::endl
                          << std::flush;
            }
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            std::cout << "Speed-up with comm. and interp." << std::endl
                      << std::flush;
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
            for (int k = 1; k <= m_numChunks; k++)
            {
                NekDouble ratio     = double(k) / m_numChunks;
                NekDouble ratioComm = commTime / fineSolveTime;
                NekDouble ratioSolve =
                    (coarseSolveTime + restTime + interTime) / fineSolveTime;
                NekDouble speedup =
                    1.0 / ((1.0 + ratio) * ratioSolve + ratio +
                           ratioComm * k * (2 * m_numChunks - k - 1) / 2.0 /
                               m_numChunks);
                std::cout << "Speed-up (" << k << ") = " << speedup << std::endl
                          << std::flush;
            }
            std::cout << "-------------------------------------------"
                      << std::endl
                      << std::flush;
        }
    }
}

References ASSERTL0, CellMLToNektar.pycml::copy(), Nektar::LibUtilities::Timer::Elapsed(), m_chunkRank, m_chunkTime, m_coarseSteps, m_coarseTimeStep, Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, m_exactSolution, m_fineSteps, m_fineTimeStep, m_interp, m_numChunks, m_pararealIterMax, m_pararealToler, Nektar::SolverUtils::Driver::m_session, m_sessionCoarse, m_totalTime, Nektar::LibUtilities::ReduceMax, RunCoarseSolve(), RunFineSolve(), Nektar::LibUtilities::Timer::Start(), and Nektar::LibUtilities::Timer::Stop().

v_InitObject()

void Nektar::SolverUtils::DriverParareal::v_InitObject ( std::ostream &  out = std::cout )

overrideprotectedvirtual

Second-stage initialisation.

Todo:

At the moment this is Navier-Stokes specific - generalise?

Reimplemented from Nektar::SolverUtils::Driver.

Definition at line 71 of file DriverParareal.cpp.

{
    try
    {
        // Retrieve the equation system to solve.
        ASSERTL0(m_session->DefinesSolverInfo("EqType"),
                 "EqType SolverInfo tag must be defined.");
        std::string vEquation = m_session->GetSolverInfo("EqType");
        if (m_session->DefinesSolverInfo("SolverType"))
        {
            vEquation = m_session->GetSolverInfo("SolverType");
        }
 
        // Check such a module exists for this equation.
        ASSERTL0(
            GetEquationSystemFactory().ModuleExists(vEquation),
            "EquationSystem '" + vEquation +
                "' is not defined.\n"
                "Ensure equation name is correct and module is compiled.\n");
 
        // Retrieve the type of evolution operator to use
        /// @todo At the moment this is Navier-Stokes specific - generalise?
        m_EvolutionOperator =
            m_session->GetSolverInfoAsEnum<EvolutionOperatorType>(
                "EvolutionOperator");
 
        m_nequ = 2;
 
        m_equ = Array<OneD, EquationSystemSharedPtr>(m_nequ);
 
        // Set the AdvectiveType tag and create EquationSystem objects.
        switch (m_EvolutionOperator)
        {
            case eNonlinear:
                // Set coarse parareal session file.
                SetPararealSessionFile();
 
                // Set fine parareal solver.
                m_session->SetTag("AdvectiveType", "Convective");
                m_session->SetTag("PararealSolver", "FineSolver");
                m_equ[0] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_session, m_graph);
 
                // Set coarse parareal solver.
                m_sessionCoarse->SetTag("AdvectiveType", "Convective");
                m_sessionCoarse->SetTag("PararealSolver", "CoarseSolver");
                m_equ[1] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_sessionCoarse, m_graphCoarse);
                break;
            case eDirect:
                // Set coarse parareal session file.
                SetPararealSessionFile();
 
                // Set fine parareal solver.
                m_session->SetTag("AdvectiveType", "Linearised");
                m_session->SetTag("PararealSolver", "FineSolver");
                m_equ[0] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_session, m_graph);
 
                // Set coarse parareal solver.
                m_sessionCoarse->SetTag("AdvectiveType", "Linearised");
                m_sessionCoarse->SetTag("PararealSolver", "CoarseSolver");
                m_equ[1] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_sessionCoarse, m_graphCoarse);
                break;
            case eAdjoint:
                // Set coarse parareal session file.
                SetPararealSessionFile();
 
                // Set fine parareal solver.
                m_session->SetTag("AdvectiveType", "Adjoint");
                m_session->SetTag("PararealSolver", "FineSolver");
                m_equ[0] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_session, m_graph);
 
                // Set coarse parareal solver.
                m_sessionCoarse->SetTag("AdvectiveType", "Adjoint");
                m_sessionCoarse->SetTag("PararealSolver", "CoarseSolver");
                m_equ[1] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_sessionCoarse, m_graphCoarse);
                break;
            case eSkewSymmetric:
                // Set coarse parareal session file.
                SetPararealSessionFile();
 
                // Set fine parareal solver.
                m_session->SetTag("AdvectiveType", "SkewSymmetric");
                m_session->SetTag("PararealSolver", "FineSolver");
                m_equ[0] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_session, m_graph);
 
                // Set coarse parareal solver.
                m_sessionCoarse->SetTag("AdvectiveType", "SkewSymmetric");
                m_sessionCoarse->SetTag("PararealSolver", "CoarseSolver");
                m_equ[1] = GetEquationSystemFactory().CreateInstance(
                    vEquation, m_sessionCoarse, m_graphCoarse);
                break;
            default:
                ASSERTL0(false, "Unrecognised evolution operator.");
        }
    }
    catch (int e)
    {
        ASSERTL0(e == -1, "No such class class defined.");
        out << "An error occurred during driver initialisation." << endl;
    }
}

References ASSERTL0, Nektar::LibUtilities::NekFactory< tKey, tBase, tParam >::CreateInstance(), Nektar::SolverUtils::eAdjoint, Nektar::SolverUtils::eDirect, Nektar::SolverUtils::eNonlinear, Nektar::SolverUtils::eSkewSymmetric, Nektar::SolverUtils::GetEquationSystemFactory(), Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_EvolutionOperator, Nektar::SolverUtils::Driver::m_graph, m_graphCoarse, Nektar::SolverUtils::Driver::m_nequ, Nektar::SolverUtils::Driver::m_session, m_sessionCoarse, and SetPararealSessionFile().

Friends And Related Function Documentation

MemoryManager< DriverParareal >

friend class MemoryManager< DriverParareal >

friend

Definition at line 1 of file DriverParareal.h.

Member Data Documentation

className

string Nektar::SolverUtils::DriverParareal::className

static

Initial value:

= GetDriverFactory().RegisterCreatorFunction(
    "Parareal", DriverParareal::create)

Name of the class.

Definition at line 64 of file DriverParareal.h.

driverLookupId

string Nektar::SolverUtils::DriverParareal::driverLookupId

staticprotected

Initial value:

=
    LibUtilities::SessionReader::RegisterEnumValue("Driver", "Parareal", 0)

Definition at line 140 of file DriverParareal.h.

m_chunkRank

int Nektar::SolverUtils::DriverParareal::m_chunkRank = 0

protected

Rank in time.

Definition at line 102 of file DriverParareal.h.

Referenced by RunCoarseSolve(), RunFineSolve(), and v_Execute().

m_chunkTime

NekDouble Nektar::SolverUtils::DriverParareal::m_chunkTime

protected

Time for chunks.

Definition at line 87 of file DriverParareal.h.

Referenced by v_Execute().

m_coarseSolveFactor

NekDouble Nektar::SolverUtils::DriverParareal::m_coarseSolveFactor = 100.0

protected

Coarse solver time factor.

Definition at line 90 of file DriverParareal.h.

Referenced by SetPararealSessionFile().

m_coarseSteps

int Nektar::SolverUtils::DriverParareal::m_coarseSteps = 1

protected

Number of steps for the coarse solver.

Definition at line 96 of file DriverParareal.h.

Referenced by RunCoarseSolve(), and v_Execute().

m_coarseTimeStep

NekDouble Nektar::SolverUtils::DriverParareal::m_coarseTimeStep

protected

Timestep for coarse solver.

Definition at line 81 of file DriverParareal.h.

Referenced by RunCoarseSolve(), and v_Execute().

m_exactSolution

bool Nektar::SolverUtils::DriverParareal::m_exactSolution = 0

protected

Using exact solution to compute error norms.

Definition at line 108 of file DriverParareal.h.

Referenced by v_Execute().

m_fineSteps

int Nektar::SolverUtils::DriverParareal::m_fineSteps = 1

protected

Number of steps for the fine solver.

Definition at line 93 of file DriverParareal.h.

Referenced by RunFineSolve(), and v_Execute().

m_fineTimeStep

NekDouble Nektar::SolverUtils::DriverParareal::m_fineTimeStep

protected

Timestep for fine solver.

Definition at line 78 of file DriverParareal.h.

Referenced by RunFineSolve(), and v_Execute().

m_graphCoarse

SpatialDomains::MeshGraphSharedPtr Nektar::SolverUtils::DriverParareal::m_graphCoarse

protected

Parareal (coarse solver) MeshGraph object.

Definition at line 75 of file DriverParareal.h.

Referenced by SetPararealSessionFile(), and v_InitObject().

m_interp

FieldUtils::Interpolator<Array<OneD, MultiRegions::ExpListSharedPtr> > Nektar::SolverUtils::DriverParareal::m_interp

protected

Definition at line 69 of file DriverParareal.h.

Referenced by RunCoarseSolve(), and v_Execute().

m_numChunks

int Nektar::SolverUtils::DriverParareal::m_numChunks = 1

protected

Number of time chunks.

Definition at line 99 of file DriverParareal.h.

Referenced by RunCoarseSolve(), RunFineSolve(), and v_Execute().

m_pararealIterMax

int Nektar::SolverUtils::DriverParareal::m_pararealIterMax = 0

protected

Maximum number of parareal iteration.

Definition at line 105 of file DriverParareal.h.

Referenced by v_Execute().

m_pararealToler

NekDouble Nektar::SolverUtils::DriverParareal::m_pararealToler = 1e-15

protected

Parareal tolerance.

Definition at line 111 of file DriverParareal.h.

Referenced by v_Execute().

m_sessionCoarse

LibUtilities::SessionReaderSharedPtr Nektar::SolverUtils::DriverParareal::m_sessionCoarse

protected

Parareal (coarse solver) session reader object.

Definition at line 72 of file DriverParareal.h.

Referenced by SetPararealSessionFile(), v_Execute(), and v_InitObject().

m_totalTime

NekDouble Nektar::SolverUtils::DriverParareal::m_totalTime

protected

Total time integration interval.

Definition at line 84 of file DriverParareal.h.

Referenced by v_Execute().