EquationSystem.cpp

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
// File: EquationSystem.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description:
//
///////////////////////////////////////////////////////////////////////////////
 
#include <FieldUtils/Interpolator.h>
#include <SolverUtils/EquationSystem.h>
 
#include <LibUtilities/BasicUtils/Equation.h>
#include <LocalRegions/MatrixKey.h>
#include <MultiRegions/ContField.h>
#include <MultiRegions/ContField3DHomogeneous1D.h>
#include <MultiRegions/ContField3DHomogeneous2D.h>
 
#include <MultiRegions/ExpList.h>
#include <MultiRegions/ExpList3DHomogeneous1D.h>
#include <MultiRegions/ExpList3DHomogeneous2D.h>
 
#include <SolverUtils/AdvectionSystem.h>
#include <SolverUtils/Diffusion/Diffusion.h>
 
#include <GlobalMapping/Mapping.h>
 
#include <boost/format.hpp>
 
#include <iostream>
#include <string>
 
using namespace std;
 
namespace Nektar::SolverUtils
{
 
std::string EquationSystem::equationSystemTypeLookupIds[2] = {
    LibUtilities::SessionReader::RegisterEnumValue("DEALIASING", "True", 0),
    LibUtilities::SessionReader::RegisterEnumValue("DEALIASING", "False", 1)};
 
std::string EquationSystem::projectionTypeLookupIds[7] = {
    LibUtilities::SessionReader::RegisterEnumValue("Projection", "Continuous",
                                                   MultiRegions::eGalerkin),
    LibUtilities::SessionReader::RegisterEnumValue("Projection", "CONTINUOUS",
                                                   MultiRegions::eGalerkin),
    LibUtilities::SessionReader::RegisterEnumValue("Projection", "Galerkin",
                                                   MultiRegions::eGalerkin),
    LibUtilities::SessionReader::RegisterEnumValue("Projection", "GALERKIN",
                                                   MultiRegions::eGalerkin),
    LibUtilities::SessionReader::RegisterEnumValue(
        "Projection", "DisContinuous", MultiRegions::eDiscontinuous),
    LibUtilities::SessionReader::RegisterEnumValue(
        "Projection", "Mixed_CG_Discontinuous",
        MultiRegions::eMixed_CG_Discontinuous),
    LibUtilities::SessionReader::RegisterEnumValue(
        "Projection", "MixedCGDG", MultiRegions::eMixed_CG_Discontinuous),
};
 
/**
 * @class EquationSystem
 *
 * This class is a base class for all solver implementations. It
 * provides the underlying generic functionality and interface for
 * solving equations.
 *
 * To solve a steady-state equation, create a derived class from this
 * class and reimplement the virtual functions to provide custom
 * implementation for the problem.
 *
 * To solve unsteady problems, derive from the UnsteadySystem class
 * instead which provides general time integration.
 */
EquationSystemFactory &GetEquationSystemFactory()
{
    static EquationSystemFactory instance;
    return instance;
}
 
/**
 * This constructor is protected as the objects of this class are never
 * instantiated directly.
 * @param   pSession The session reader holding problem parameters.
 */
EquationSystem::EquationSystem(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const SpatialDomains::MeshGraphSharedPtr &pGraph)
    : m_comm(pSession->GetComm()), m_session(pSession), m_graph(pGraph),
      m_lambda(0), m_infosteps(10),
      m_fieldMetaDataMap(LibUtilities::NullFieldMetaDataMap)
{
    // set up session names in fieldMetaDataMap
    const vector<std::string> filenames = m_session->GetFilenames();
 
    for (int i = 0; i < filenames.size(); ++i)
    {
        string sessionname = "SessionName";
        sessionname += std::to_string(i);
        m_fieldMetaDataMap[sessionname]  = filenames[i];
        m_fieldMetaDataMap["ChkFileNum"] = std::to_string(0);
    }
}
 
/**
 * @brief Destructor for class EquationSystem.
 */
EquationSystem::~EquationSystem()
{
    LibUtilities::NekManager<LocalRegions::MatrixKey, DNekScalMat,
                             LocalRegions::MatrixKey::opLess>::ClearManager();
    LibUtilities::NekManager<LocalRegions::MatrixKey, DNekScalBlkMat,
                             LocalRegions::MatrixKey::opLess>::ClearManager();
}
 
/**
 * @brief Initialisation object for EquationSystem.
 */
void EquationSystem::v_InitObject(bool DeclareFields)
{
    // Save the basename of input file name for output details
    m_sessionName = m_session->GetSessionName();
 
    // Instantiate a field reader/writer
    m_fld = LibUtilities::FieldIO::CreateDefault(m_session);
 
    // Also read and store the boundary conditions
    m_boundaryConditions =
        MemoryManager<SpatialDomains::BoundaryConditions>::AllocateSharedPtr(
            m_session, m_graph);
 
    // Set space dimension for use in class
    m_spacedim = m_graph->GetSpaceDimension();
 
    // Setting parameteres for homogenous problems
    m_HomoDirec          = 0;
    m_useFFT             = false;
    m_homogen_dealiasing = false;
    m_singleMode         = false;
    m_halfMode           = false;
    m_multipleModes      = false;
    m_HomogeneousType    = eNotHomogeneous;
 
    m_verbose = m_session->DefinesCmdLineArgument("verbose");
 
    if (m_session->DefinesSolverInfo("HOMOGENEOUS"))
    {
        std::string HomoStr = m_session->GetSolverInfo("HOMOGENEOUS");
        m_spacedim          = 3;
 
        if ((HomoStr == "HOMOGENEOUS1D") || (HomoStr == "Homogeneous1D") ||
            (HomoStr == "1D") || (HomoStr == "Homo1D"))
        {
            m_HomogeneousType = eHomogeneous1D;
            m_session->LoadParameter("LZ", m_LhomZ);
            m_HomoDirec = 1;
 
            if (m_session->DefinesSolverInfo("ModeType"))
            {
                m_session->MatchSolverInfo("ModeType", "SingleMode",
                                           m_singleMode, false);
                m_session->MatchSolverInfo("ModeType", "HalfMode", m_halfMode,
                                           false);
                m_session->MatchSolverInfo("ModeType", "MultipleModes",
                                           m_multipleModes, false);
            }
 
            // Stability Analysis flags
            if (m_session->DefinesSolverInfo("ModeType"))
            {
                if (m_singleMode)
                {
                    m_npointsZ = 2;
                }
                else if (m_halfMode)
                {
                    m_npointsZ = 1;
                }
                else if (m_multipleModes)
                {
                    m_npointsZ = m_session->GetParameter("HomModesZ");
                }
                else
                {
                    ASSERTL0(false, "SolverInfo ModeType not valid");
                }
            }
            else
            {
                m_npointsZ = m_session->GetParameter("HomModesZ");
            }
        }
 
        if ((HomoStr == "HOMOGENEOUS2D") || (HomoStr == "Homogeneous2D") ||
            (HomoStr == "2D") || (HomoStr == "Homo2D"))
        {
            m_HomogeneousType = eHomogeneous2D;
            m_session->LoadParameter("HomModesY", m_npointsY);
            m_session->LoadParameter("LY", m_LhomY);
            m_session->LoadParameter("HomModesZ", m_npointsZ);
            m_session->LoadParameter("LZ", m_LhomZ);
            m_HomoDirec = 2;
        }
 
        if ((HomoStr == "HOMOGENEOUS3D") || (HomoStr == "Homogeneous3D") ||
            (HomoStr == "3D") || (HomoStr == "Homo3D"))
        {
            m_HomogeneousType = eHomogeneous3D;
            m_session->LoadParameter("HomModesY", m_npointsY);
            m_session->LoadParameter("LY", m_LhomY);
            m_session->LoadParameter("HomModesZ", m_npointsZ);
            m_session->LoadParameter("LZ", m_LhomZ);
            m_HomoDirec = 2;
        }
 
        m_session->MatchSolverInfo("USEFFT", "FFTW", m_useFFT, false);
 
        m_session->MatchSolverInfo("DEALIASING", "True", m_homogen_dealiasing,
                                   false);
    }
    else
    {
        // set to default value so can use to identify 2d or 3D
        // (homogeneous) expansions
        m_npointsZ = 1;
    }
 
    m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "True",
                               m_specHP_dealiasing, false);
    if (m_specHP_dealiasing == false)
    {
        m_session->MatchSolverInfo("SPECTRALHPDEALIASING", "On",
                                   m_specHP_dealiasing, false);
    }
 
    // Options to determine type of projection from file or directly
    // from constructor
    if (m_session->DefinesSolverInfo("PROJECTION"))
    {
        std::string ProjectStr = m_session->GetSolverInfo("PROJECTION");
 
        if ((ProjectStr == "Continuous") || (ProjectStr == "Galerkin") ||
            (ProjectStr == "CONTINUOUS") || (ProjectStr == "GALERKIN"))
        {
            m_projectionType = MultiRegions::eGalerkin;
        }
        else if ((ProjectStr == "MixedCGDG") ||
                 (ProjectStr == "Mixed_CG_Discontinuous"))
        {
            m_projectionType = MultiRegions::eMixed_CG_Discontinuous;
        }
        else if (ProjectStr == "DisContinuous")
        {
            m_projectionType = MultiRegions::eDiscontinuous;
        }
        else
        {
            ASSERTL0(false, "PROJECTION value not recognised");
        }
    }
    else
    {
        cerr << "Projection type not specified in SOLVERINFO,"
                "defaulting to continuous Galerkin"
             << endl;
        m_projectionType = MultiRegions::eGalerkin;
    }
 
    // Enforce singularity check for some problems
    m_checkIfSystemSingular = v_GetSystemSingularChecks();
 
    int i;
    int nvariables              = m_session->GetVariables().size();
    bool DeclareCoeffPhysArrays = true;
 
    m_fields   = Array<OneD, MultiRegions::ExpListSharedPtr>(nvariables);
    m_spacedim = m_graph->GetSpaceDimension() + m_HomoDirec;
    m_expdim   = m_graph->GetMeshDimension();
 
    if (DeclareFields) // declare field if required
    {
        /// Continuous field
        if (m_projectionType == MultiRegions::eGalerkin ||
            m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
        {
            switch (m_expdim)
            {
                case 1:
                {
                    if (m_HomogeneousType == eHomogeneous2D ||
                        m_HomogeneousType == eHomogeneous3D)
                    {
                        const LibUtilities::PointsKey PkeyY(
                            m_npointsY, LibUtilities::eFourierEvenlySpaced);
                        const LibUtilities::BasisKey BkeyY(
                            LibUtilities::eFourier, m_npointsY, PkeyY);
                        const LibUtilities::PointsKey PkeyZ(
                            m_npointsZ, LibUtilities::eFourierEvenlySpaced);
                        const LibUtilities::BasisKey BkeyZ(
                            LibUtilities::eFourier, m_npointsZ, PkeyZ);
 
                        for (i = 0; i < m_fields.size(); i++)
                        {
                            m_fields[i] = MemoryManager<
                                MultiRegions::ContField3DHomogeneous2D>::
                                AllocateSharedPtr(m_session, BkeyY, BkeyZ,
                                                  m_LhomY, m_LhomZ, m_useFFT,
                                                  m_homogen_dealiasing, m_graph,
                                                  m_session->GetVariable(i));
                        }
                    }
                    else
                    {
                        for (i = 0; i < m_fields.size(); i++)
                        {
                            m_fields[i] =
                                MemoryManager<MultiRegions::ContField>::
                                    AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i));
                        }
                    }
                    break;
                }
                case 2:
                {
                    if (m_HomogeneousType == eHomogeneous1D)
                    {
                        // Fourier single mode stability analysis
                        if (m_singleMode)
                        {
                            const LibUtilities::PointsKey PkeyZ(
                                m_npointsZ,
                                LibUtilities::eFourierSingleModeSpaced);
 
                            const LibUtilities::BasisKey BkeyZ(
                                LibUtilities::eFourierSingleMode, m_npointsZ,
                                PkeyZ);
 
                            for (i = 0; i < m_fields.size(); i++)
                            {
                                m_fields[i] = MemoryManager<
                                    MultiRegions::ContField3DHomogeneous1D>::
                                    AllocateSharedPtr(
                                        m_session, BkeyZ, m_LhomZ, m_useFFT,
                                        m_homogen_dealiasing, m_graph,
                                        m_session->GetVariable(i),
                                        m_checkIfSystemSingular[i]);
                            }
                        }
                        // Half mode stability analysis
                        else if (m_halfMode)
                        {
                            const LibUtilities::PointsKey PkeyZ(
                                m_npointsZ,
                                LibUtilities::eFourierSingleModeSpaced);
 
                            const LibUtilities::BasisKey BkeyZR(
                                LibUtilities::eFourierHalfModeRe, m_npointsZ,
                                PkeyZ);
 
                            const LibUtilities::BasisKey BkeyZI(
                                LibUtilities::eFourierHalfModeIm, m_npointsZ,
                                PkeyZ);
 
                            for (i = 0; i < m_fields.size(); i++)
                            {
                                if (m_session->GetVariable(i).compare("w") == 0)
                                {
                                    m_fields[i] = MemoryManager<
                                        MultiRegions::
                                            ContField3DHomogeneous1D>::
                                        AllocateSharedPtr(
                                            m_session, BkeyZI, m_LhomZ,
                                            m_useFFT, m_homogen_dealiasing,
                                            m_graph, m_session->GetVariable(i),
                                            m_checkIfSystemSingular[i]);
                                }
                                else
                                {
                                    m_fields[i] = MemoryManager<
                                        MultiRegions::
                                            ContField3DHomogeneous1D>::
                                        AllocateSharedPtr(
                                            m_session, BkeyZR, m_LhomZ,
                                            m_useFFT, m_homogen_dealiasing,
                                            m_graph, m_session->GetVariable(i),
                                            m_checkIfSystemSingular[i]);
                                }
                            }
                        }
                        // Normal homogeneous 1D
                        else
                        {
                            const LibUtilities::PointsKey PkeyZ(
                                m_npointsZ, LibUtilities::eFourierEvenlySpaced);
                            const LibUtilities::BasisKey BkeyZ(
                                LibUtilities::eFourier, m_npointsZ, PkeyZ);
 
                            for (i = 0; i < m_fields.size(); i++)
                            {
                                m_fields[i] = MemoryManager<
                                    MultiRegions::ContField3DHomogeneous1D>::
                                    AllocateSharedPtr(
                                        m_session, BkeyZ, m_LhomZ, m_useFFT,
                                        m_homogen_dealiasing, m_graph,
                                        m_session->GetVariable(i),
                                        m_checkIfSystemSingular[i]);
                            }
                        }
                    }
                    else
                    {
                        i = 0;
                        MultiRegions::ContFieldSharedPtr firstfield;
                        firstfield = MemoryManager<MultiRegions::ContField>::
                            AllocateSharedPtr(m_session, m_graph,
                                              m_session->GetVariable(i),
                                              DeclareCoeffPhysArrays,
                                              m_checkIfSystemSingular[0]);
                        m_fields[0] = firstfield;
                        for (i = 1; i < m_fields.size(); i++)
                        {
                            if (m_graph->SameExpansionInfo(
                                    m_session->GetVariable(0),
                                    m_session->GetVariable(i)))
                            {
                                m_fields[i] =
                                    MemoryManager<MultiRegions::ContField>::
                                        AllocateSharedPtr(
                                            *firstfield, m_graph,
                                            m_session->GetVariable(i),
                                            DeclareCoeffPhysArrays,
                                            m_checkIfSystemSingular[i]);
                            }
                            else
                            {
                                m_fields[i] =
                                    MemoryManager<MultiRegions::ContField>::
                                        AllocateSharedPtr(
                                            m_session, m_graph,
                                            m_session->GetVariable(i),
                                            DeclareCoeffPhysArrays,
                                            m_checkIfSystemSingular[i]);
                            }
                        }
 
                        if (m_projectionType ==
                            MultiRegions::eMixed_CG_Discontinuous)
                        {
                            /// Setting up the normals
                            m_traceNormals =
                                Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
 
                            for (i = 0; i < m_spacedim; ++i)
                            {
                                m_traceNormals[i] =
                                    Array<OneD, NekDouble>(GetTraceNpoints());
                            }
 
                            m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
                        }
                    }
 
                    break;
                }
                case 3:
                {
                    i = 0;
                    MultiRegions::ContFieldSharedPtr firstfield =
                        MemoryManager<MultiRegions::ContField>::
                            AllocateSharedPtr(m_session, m_graph,
                                              m_session->GetVariable(i),
                                              DeclareCoeffPhysArrays,
                                              m_checkIfSystemSingular[i]);
 
                    m_fields[0] = firstfield;
                    for (i = 1; i < m_fields.size(); i++)
                    {
                        if (m_graph->SameExpansionInfo(
                                m_session->GetVariable(0),
                                m_session->GetVariable(i)))
                        {
                            m_fields[i] =
                                MemoryManager<MultiRegions::ContField>::
                                    AllocateSharedPtr(
                                        *firstfield, m_graph,
                                        m_session->GetVariable(i),
                                        DeclareCoeffPhysArrays,
                                        m_checkIfSystemSingular[i]);
                        }
                        else
                        {
                            m_fields[i] =
                                MemoryManager<MultiRegions::ContField>::
                                    AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i),
                                        DeclareCoeffPhysArrays,
                                        m_checkIfSystemSingular[i]);
                        }
                    }
 
                    if (m_projectionType ==
                        MultiRegions::eMixed_CG_Discontinuous)
                    {
                        /// Setting up the normals
                        m_traceNormals =
                            Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
                        for (i = 0; i < m_spacedim; ++i)
                        {
                            m_traceNormals[i] =
                                Array<OneD, NekDouble>(GetTraceNpoints());
                        }
 
                        m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
                        // Call the trace on all fields to ensure DG setup.
                        for (i = 1; i < m_fields.size(); ++i)
                        {
                            m_fields[i]->GetTrace();
                        }
                    }
                    break;
                }
                default:
                    ASSERTL0(false, "Expansion dimension not recognised");
                    break;
            }
        }
        // Discontinuous field
        else
        {
            switch (m_expdim)
            {
                case 1:
                {
                    if (m_HomogeneousType == eHomogeneous2D ||
                        m_HomogeneousType == eHomogeneous3D)
                    {
                        const LibUtilities::PointsKey PkeyY(
                            m_npointsY, LibUtilities::eFourierEvenlySpaced);
                        const LibUtilities::BasisKey BkeyY(
                            LibUtilities::eFourier, m_npointsY, PkeyY);
                        const LibUtilities::PointsKey PkeyZ(
                            m_npointsZ, LibUtilities::eFourierEvenlySpaced);
                        const LibUtilities::BasisKey BkeyZ(
                            LibUtilities::eFourier, m_npointsZ, PkeyZ);
 
                        for (i = 0; i < m_fields.size(); i++)
                        {
                            m_fields[i] = MemoryManager<
                                MultiRegions::DisContField3DHomogeneous2D>::
                                AllocateSharedPtr(m_session, BkeyY, BkeyZ,
                                                  m_LhomY, m_LhomZ, m_useFFT,
                                                  m_homogen_dealiasing, m_graph,
                                                  m_session->GetVariable(i));
                        }
                    }
                    else
                    {
                        for (i = 0; i < m_fields.size(); i++)
                        {
                            m_fields[i] =
                                MemoryManager<MultiRegions::DisContField>::
                                    AllocateSharedPtr(
                                        m_session, m_graph,
                                        m_session->GetVariable(i));
                        }
                    }
 
                    break;
                }
                case 2:
                {
                    if (m_HomogeneousType == eHomogeneous1D)
                    {
                        const LibUtilities::PointsKey PkeyZ(
                            m_npointsZ, LibUtilities::eFourierEvenlySpaced);
                        const LibUtilities::BasisKey BkeyZ(
                            LibUtilities::eFourier, m_npointsZ, PkeyZ);
 
                        for (i = 0; i < m_fields.size(); i++)
                        {
                            m_fields[i] = MemoryManager<
                                MultiRegions::DisContField3DHomogeneous1D>::
                                AllocateSharedPtr(m_session, BkeyZ, m_LhomZ,
                                                  m_useFFT,
                                                  m_homogen_dealiasing, m_graph,
                                                  m_session->GetVariable(i));
                        }
                    }
                    else
                    {
                        i = 0;
                        MultiRegions::DisContFieldSharedPtr firstfield;
                        firstfield = MemoryManager<MultiRegions::DisContField>::
                            AllocateSharedPtr(m_session, m_graph,
                                              m_session->GetVariable(i));
                        m_fields[0] = firstfield;
                        for (i = 1; i < m_fields.size(); i++)
                        {
                            if (m_graph->SameExpansionInfo(
                                    m_session->GetVariable(0),
                                    m_session->GetVariable(i)))
                            {
                                m_fields[i] =
                                    MemoryManager<MultiRegions::DisContField>::
                                        AllocateSharedPtr(
                                            *firstfield, m_graph,
                                            m_session->GetVariable(i));
                            }
                            else
                            {
                                m_fields[i] =
                                    MemoryManager<MultiRegions::DisContField>::
                                        AllocateSharedPtr(
                                            m_session, m_graph,
                                            m_session->GetVariable(i));
                            }
                        }
                    }
 
                    break;
                }
                case 3:
                {
                    if (m_HomogeneousType == eHomogeneous3D)
                    {
                        ASSERTL0(
                            false,
                            "3D fully periodic problems not implemented yet");
                    }
                    else
                    {
                        i = 0;
                        MultiRegions::DisContFieldSharedPtr firstfield =
                            MemoryManager<MultiRegions::DisContField>::
                                AllocateSharedPtr(m_session, m_graph,
                                                  m_session->GetVariable(i));
                        m_fields[0] = firstfield;
                        for (i = 1; i < m_fields.size(); i++)
                        {
                            if (m_graph->SameExpansionInfo(
                                    m_session->GetVariable(0),
                                    m_session->GetVariable(i)))
                            {
                                m_fields[i] =
                                    MemoryManager<MultiRegions::DisContField>::
                                        AllocateSharedPtr(
                                            *firstfield, m_graph,
                                            m_session->GetVariable(i));
                            }
                            else
                            {
                                m_fields[i] =
                                    MemoryManager<MultiRegions::DisContField>::
                                        AllocateSharedPtr(
                                            m_session, m_graph,
                                            m_session->GetVariable(i));
                            }
                        }
                    }
                    break;
                }
                default:
                    ASSERTL0(false, "Expansion dimension not recognised");
                    break;
            }
 
            // Setting up the normals
            m_traceNormals = Array<OneD, Array<OneD, NekDouble>>(m_spacedim);
 
            for (i = 0; i < m_spacedim; ++i)
            {
                m_traceNormals[i] =
                    Array<OneD, NekDouble>(GetTraceNpoints(), 0.0);
            }
 
            m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
        }
        // Zero all physical fields initially
        ZeroPhysFields();
    }
 
    // Set Default Parameter
    m_session->LoadParameter("Time", m_time, 0.0);
    m_session->LoadParameter("TimeStep", m_timestep, 0.0);
    m_session->LoadParameter("NumSteps", m_steps, 0);
    m_session->LoadParameter("IO_CheckSteps", m_checksteps, 0);
    m_session->LoadParameter("IO_CheckTime", m_checktime, 0.0);
    m_session->LoadParameter("FinTime", m_fintime, 0);
    m_session->LoadParameter("NumQuadPointsError", m_NumQuadPointsError, 0);
 
    // Check uniqueness of checkpoint output
    ASSERTL0((m_checktime == 0.0 && m_checksteps == 0) ||
                 (m_checktime > 0.0 && m_checksteps == 0) ||
                 (m_checktime == 0.0 && m_checksteps > 0),
             "Only one of IO_CheckTime and IO_CheckSteps "
             "should be set!");
    WARNINGL0(!m_session->DefinesParameter("IO_CheckSteps"),
              "The IO_CheckSteps parameter is now deprecated.\nPlease use the "
              "'Checkpoint' filter instead");
    WARNINGL0(!m_session->DefinesParameter("IO_CheckTime"),
              "The IO_CheckTime parameter is now deprecated.\nPlease use the "
              "'Checkpoint' filter instead");
 
    m_session->LoadParameter("TimeIncrementFactor", m_TimeIncrementFactor, 1.0);
 
    // Check for parallel-in-time
    if (m_comm->IsParallelInTime())
    {
        ASSERTL0(m_fintime == 0.0,
                 "Only specify NumSteps and TimeSteps for Parallel-in-Time. "
                 "FinTime should not be used! ");
    }
 
    m_nchk    = 0;
    m_iterPIT = 0;
}
 
/**
 *
 */
SessionFunctionSharedPtr EquationSystem::GetFunction(
    std::string name, const MultiRegions::ExpListSharedPtr &field, bool cache)
{
    MultiRegions::ExpListSharedPtr vField = field;
    if (!field)
    {
        vField = m_fields[0];
    }
 
    if (cache)
    {
        if ((m_sessionFunctions.find(name) == m_sessionFunctions.end()) ||
            (m_sessionFunctions[name]->GetSession() != m_session) ||
            (m_sessionFunctions[name]->GetExpansion() != vField))
        {
            m_sessionFunctions[name] =
                MemoryManager<SessionFunction>::AllocateSharedPtr(
                    m_session, vField, name, cache);
        }
 
        return m_sessionFunctions[name];
    }
    else
    {
        return SessionFunctionSharedPtr(
            new SessionFunction(m_session, vField, name, cache));
    }
}
 
/**
 * If boundary conditions are time-dependent, they will be evaluated at
 * the time specified.
 * @param   time            The time at which to evaluate the BCs
 */
void EquationSystem::SetBoundaryConditions(NekDouble time)
{
    std::string varName;
    int nvariables = m_fields.size();
    for (int i = 0; i < nvariables; ++i)
    {
        varName = m_session->GetVariable(i);
        m_fields[i]->EvaluateBoundaryConditions(time, varName);
    }
}
 
/**
 * Compute the error in the L2-norm.
 * @param   field           The field to compare.
 * @param   exactsoln       The exact solution to compare with.
 * @param   Normalised      Normalise L2-error.
 * @returns                 Error in the L2-norm.
 */
NekDouble EquationSystem::v_L2Error(unsigned int field,
                                    const Array<OneD, NekDouble> &exactsoln,
                                    bool Normalised)
{
    NekDouble L2error = -1.0;
 
    if (m_NumQuadPointsError == 0)
    {
        if (m_fields[field]->GetPhysState() == false)
        {
            m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
                                      m_fields[field]->UpdatePhys());
        }
 
        // Transform from (phys, wave) -> (phys, phys), for 2.5D
        bool physwavetrans = false;
        if (m_fields[field]->GetWaveSpace() == true)
        {
            m_fields[field]->HomogeneousBwdTrans(
                m_fields[field]->GetTotPoints(), m_fields[field]->GetPhys(),
                m_fields[field]->UpdatePhys());
            physwavetrans = true;
        }
 
        if (exactsoln.size())
        {
            L2error =
                m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
        }
        else if (m_session->DefinesFunction("ExactSolution"))
        {
            Array<OneD, NekDouble> exactsoln(m_fields[field]->GetNpoints());
 
            GetFunction("ExactSolution")
                ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
 
            L2error =
                m_fields[field]->L2(m_fields[field]->GetPhys(), exactsoln);
        }
        else
        {
            L2error = m_fields[field]->L2(m_fields[field]->GetPhys());
        }
 
        if (Normalised == true)
        {
            Array<OneD, NekDouble> one(m_fields[field]->GetNpoints(), 1.0);
 
            NekDouble Vol = m_fields[field]->Integral(one);
            L2error       = sqrt(L2error * L2error / Vol);
        }
 
        // Transform from (phys, phys) -> (phys, wave), for 2.5D
        if (physwavetrans == true)
        {
            m_fields[field]->HomogeneousFwdTrans(
                m_fields[field]->GetTotPoints(), m_fields[field]->GetPhys(),
                m_fields[field]->UpdatePhys());
            physwavetrans = false;
        }
    }
    else
    {
        Array<OneD, NekDouble> errNorms(3);
        errNorms = ErrorExtraPoints(field);
        L2error  = errNorms[0];
    }
    return L2error;
}
 
/**
 * Compute the error in the L_inf-norm
 * @param   field           The field to compare.
 * @param   exactsoln       The exact solution to compare with.
 * @returns                 Error in the L_inft-norm.
 */
NekDouble EquationSystem::v_LinfError(unsigned int field,
                                      const Array<OneD, NekDouble> &exactsoln)
{
    NekDouble Linferror = -1.0;
 
    if (m_NumQuadPointsError == 0)
    {
        if (m_fields[field]->GetPhysState() == false)
        {
            m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
                                      m_fields[field]->UpdatePhys());
        }
 
        if (exactsoln.size())
        {
            Linferror =
                m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
        }
        else if (m_session->DefinesFunction("ExactSolution"))
        {
            Array<OneD, NekDouble> exactsoln(m_fields[field]->GetNpoints());
 
            GetFunction("ExactSolution")
                ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
 
            Linferror =
                m_fields[field]->Linf(m_fields[field]->GetPhys(), exactsoln);
        }
        else
        {
            Linferror = m_fields[field]->Linf(m_fields[field]->GetPhys());
        }
    }
    else
    {
        Array<OneD, NekDouble> errNorms(3);
        errNorms  = ErrorExtraPoints(field);
        Linferror = errNorms[1];
    }
 
    return Linferror;
}
 
/**
 * Compute the error in the H1-norm.
 * @param   field           The field to compare.
 * @param   exactsoln       The exact solution to compare with.
 * @param   Normalised      Normalise H1-error.
 * @returns                 Error in the H1-norm.
 */
NekDouble EquationSystem::v_H1Error(unsigned int field,
                                    const Array<OneD, NekDouble> &exactsoln,
                                    bool Normalised)
{
    NekDouble error = -1.0;
 
    if (m_NumQuadPointsError == 0)
    {
        if (m_fields[field]->GetPhysState() == false)
        {
            m_fields[field]->BwdTrans(m_fields[field]->GetCoeffs(),
                                      m_fields[field]->UpdatePhys());
        }
 
        if (exactsoln.size())
        {
            error = m_fields[field]->H1(m_fields[field]->GetPhys(), exactsoln);
        }
        else if (m_session->DefinesFunction("ExactSolution"))
        {
            Array<OneD, NekDouble> exactsoln(m_fields[field]->GetNpoints());
 
            GetFunction("ExactSolution")
                ->Evaluate(m_session->GetVariable(field), exactsoln, m_time);
 
            error = m_fields[field]->H1(m_fields[field]->GetPhys(), exactsoln);
        }
        else
        {
            error = m_fields[field]->H1(m_fields[field]->GetPhys());
        }
 
        if (Normalised == true)
        {
            Array<OneD, NekDouble> one(m_fields[field]->GetNpoints(), 1.0);
 
            NekDouble Vol = m_fields[field]->Integral(one);
            error         = sqrt(error * error / Vol);
        }
    }
    else
    {
        Array<OneD, NekDouble> errNorms(3);
        errNorms = ErrorExtraPoints(field);
        error    = errNorms[2];
    }
    return error;
}
 
/**
 * Compute the error in the L2-, Linf- and H1-norm for a larger number of
 * quadrature points.
 * @param   field              The field to compare.
 * @returns                    Error in the L2-, Linf- and H1-norm.
 */
Array<OneD, NekDouble> EquationSystem::ErrorExtraPoints(unsigned int field)
{
    int NumModes = GetNumExpModes();
    Array<OneD, NekDouble> errNorms(3); // Expects 3 norms (L2, Linf, H1)
 
    const LibUtilities::PointsKey PkeyT1(m_NumQuadPointsError,
                                         LibUtilities::eGaussLobattoLegendre);
    const LibUtilities::PointsKey PkeyT2(m_NumQuadPointsError,
                                         LibUtilities::eGaussRadauMAlpha1Beta0);
    const LibUtilities::PointsKey PkeyQ1(m_NumQuadPointsError,
                                         LibUtilities::eGaussLobattoLegendre);
    const LibUtilities::PointsKey PkeyQ2(m_NumQuadPointsError,
                                         LibUtilities::eGaussLobattoLegendre);
    const LibUtilities::BasisKey BkeyT1(LibUtilities::eModified_A, NumModes,
                                        PkeyT1);
    const LibUtilities::BasisKey BkeyT2(LibUtilities::eModified_B, NumModes,
                                        PkeyT2);
    const LibUtilities::BasisKey BkeyQ1(LibUtilities::eModified_A, NumModes,
                                        PkeyQ1);
    const LibUtilities::BasisKey BkeyQ2(LibUtilities::eModified_A, NumModes,
                                        PkeyQ2);
 
    LibUtilities::BasisKeyVector Tkeys, Qkeys;
 
    // make a copy of the ExpansionInfoMap
    SpatialDomains::ExpansionInfoMap NewExpInfo = m_graph->GetExpansionInfo();
    SpatialDomains::ExpansionInfoMapShPtr ExpInfo =
        MemoryManager<SpatialDomains::ExpansionInfoMap>::AllocateSharedPtr(
            NewExpInfo);
 
    // reset new graph with new keys
    Tkeys.push_back(BkeyT1);
    Tkeys.push_back(BkeyT2);
    m_graph->ResetExpansionInfoToBasisKey(ExpInfo, LibUtilities::eTriangle,
                                          Tkeys);
    Qkeys.push_back(BkeyQ1);
    Qkeys.push_back(BkeyQ2);
    m_graph->ResetExpansionInfoToBasisKey(ExpInfo, LibUtilities::eQuadrilateral,
                                          Qkeys);
 
    MultiRegions::ExpListSharedPtr ErrorExplist =
        MemoryManager<MultiRegions::ExpList>::AllocateSharedPtr(m_session,
                                                                NewExpInfo);
 
    int ErrorCoordim = ErrorExplist->GetCoordim(0);
    int ErrorNq      = ErrorExplist->GetTotPoints();
 
    Array<OneD, NekDouble> ErrorXc0(ErrorNq, 0.0);
    Array<OneD, NekDouble> ErrorXc1(ErrorNq, 0.0);
    Array<OneD, NekDouble> ErrorXc2(ErrorNq, 0.0);
 
    switch (ErrorCoordim)
    {
        case 1:
            ErrorExplist->GetCoords(ErrorXc0);
            break;
        case 2:
            ErrorExplist->GetCoords(ErrorXc0, ErrorXc1);
            break;
        case 3:
            ErrorExplist->GetCoords(ErrorXc0, ErrorXc1, ErrorXc2);
            break;
    }
    LibUtilities::EquationSharedPtr exSol =
        m_session->GetFunction("ExactSolution", field);
 
    // Evaluate the exact solution
    Array<OneD, NekDouble> ErrorSol(ErrorNq);
 
    exSol->Evaluate(ErrorXc0, ErrorXc1, ErrorXc2, m_time, ErrorSol);
 
    // Calcualte spectral/hp approximation on the quadrature points
    // of this new expansion basis
    ErrorExplist->BwdTrans(m_fields[field]->GetCoeffs(),
                           ErrorExplist->UpdatePhys());
 
    errNorms[0] = ErrorExplist->L2(ErrorExplist->GetPhys(), ErrorSol);
    errNorms[1] = ErrorExplist->Linf(ErrorExplist->GetPhys(), ErrorSol);
    errNorms[2] = ErrorExplist->H1(ErrorExplist->GetPhys(), ErrorSol);
 
    return errNorms;
}
 
/**
 * Set the physical fields based on a restart file, or a function
 * describing the initial condition given in the session.
 * @param  initialtime           Time at which to evaluate the function.
 * @param  dumpInitialConditions Write the initial condition to file?
 */
void EquationSystem::v_SetInitialConditions(
    [[maybe_unused]] NekDouble initialtime, bool dumpInitialConditions,
    const int domain)
{
    if (m_session->GetComm()->GetRank() == 0)
    {
        cout << "Initial Conditions:" << endl;
    }
 
    if (m_session->DefinesFunction("InitialConditions"))
    {
        GetFunction("InitialConditions")
            ->Evaluate(m_session->GetVariables(), m_fields, m_time, domain);
        // Enforce C0 Continutiy of initial condiiton
        if ((m_projectionType == MultiRegions::eGalerkin) ||
            (m_projectionType == MultiRegions::eMixed_CG_Discontinuous))
        {
            for (int i = 0; i < m_fields.size(); ++i)
            {
                m_fields[i]->LocalToGlobal();
                m_fields[i]->GlobalToLocal();
                m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),
                                      m_fields[i]->UpdatePhys());
            }
        }
 
        if (m_session->GetComm()->GetRank() == 0)
        {
            for (int i = 0; i < m_fields.size(); ++i)
            {
                std::string varName = m_session->GetVariable(i);
                cout << "  - Field " << varName << ": "
                     << GetFunction("InitialConditions")
                            ->Describe(varName, domain)
                     << endl;
            }
        }
    }
    else
    {
        int nq = m_fields[0]->GetNpoints();
        for (int i = 0; i < m_fields.size(); i++)
        {
            Vmath::Zero(nq, m_fields[i]->UpdatePhys(), 1);
            m_fields[i]->SetPhysState(true);
            Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(),
                        1);
            if (m_session->GetComm()->GetRank() == 0)
            {
                cout << "  - Field " << m_session->GetVariable(i)
                     << ": 0 (default)" << endl;
            }
        }
    }
 
    if (dumpInitialConditions && m_checksteps && m_nchk == 0 &&
        !m_comm->IsParallelInTime())
    {
        Checkpoint_Output(m_nchk);
    }
    else if (dumpInitialConditions && m_nchk == 0 && m_comm->IsParallelInTime())
    {
        std::string newdir = m_sessionName + ".pit";
        if (!fs::is_directory(newdir))
        {
            fs::create_directory(newdir);
        }
        if (m_comm->GetTimeComm()->GetRank() == 0)
        {
            WriteFld(newdir + "/" + m_sessionName + "_0.fld");
        }
    }
    ++m_nchk;
}
 
/**
 *
 */
void EquationSystem::v_EvaluateExactSolution(unsigned int field,
                                             Array<OneD, NekDouble> &outfield,
                                             const NekDouble time)
{
    ASSERTL0(outfield.size() == m_fields[field]->GetNpoints(),
             "ExactSolution array size mismatch.");
    Vmath::Zero(outfield.size(), outfield, 1);
    if (m_session->DefinesFunction("ExactSolution"))
    {
        GetFunction("ExactSolution")
            ->Evaluate(m_session->GetVariable(field), outfield, time);
    }
}
 
/**
 * By default, nothing needs initialising at the EquationSystem level.
 */
void EquationSystem::v_DoInitialise([[maybe_unused]] bool dumpInitialConditions)
{
}
 
/**
 *
 */
void EquationSystem::v_DoSolve()
{
}
 
/**
 * Virtual function to define if operator in DoSolve is
 * negated with regard to the strong form. This is currently
 * only used in Arnoldi solves. Default is false.
 */
bool EquationSystem::v_NegatedOp(void)
{
    return false;
}
 
/**
 *
 */
void EquationSystem::v_TransCoeffToPhys()
{
}
 
/**
 *
 */
void EquationSystem::v_TransPhysToCoeff()
{
}
 
/// Virtual function for generating summary information.
void EquationSystem::v_GenerateSummary(SummaryList &l)
{
    SessionSummary(l);
}
 
/**
 * Write the field data to file. The file is named according to the
 * session name with the extension .fld appended.
 */
void EquationSystem::v_Output(void)
{
    if (!m_comm->IsParallelInTime())
    {
        // Serial-in-time
        WriteFld(m_sessionName + ".fld");
    }
    else
    {
        // Parallel-in-time
        std::string newdir = m_sessionName + ".pit";
        if (!fs::is_directory(newdir))
        {
            fs::create_directory(newdir);
        }
        WriteFld(newdir + "/" + m_sessionName + "_" +
                 std::to_string(m_windowPIT * m_comm->GetTimeComm()->GetSize() +
                                m_comm->GetTimeComm()->GetRank() + 1) +
                 ".fld");
    }
}
 
/**
 * Zero the physical fields.
 */
void EquationSystem::ZeroPhysFields(void)
{
    for (int i = 0; i < m_fields.size(); i++)
    {
        Vmath::Zero(m_fields[i]->GetNpoints(), m_fields[i]->UpdatePhys(), 1);
    }
}
 
/**
 * FwdTrans the m_fields members
 */
void EquationSystem::FwdTransFields(void)
{
    for (int i = 0; i < m_fields.size(); i++)
    {
        m_fields[i]->FwdTrans(m_fields[i]->GetPhys(),
                              m_fields[i]->UpdateCoeffs());
        m_fields[i]->SetPhysState(false);
    }
}
 
/**
 * Write the n-th checkpoint file.
 * @param   n   The index of the checkpoint file.
 */
void EquationSystem::Checkpoint_Output(const int n)
{
    if (!m_comm->IsParallelInTime())
    {
        // Serial-in-time
        std::string outname = m_sessionName + "_" + std::to_string(n);
        WriteFld(outname + ".chk");
    }
    else
    {
        // Parallel-in-time
        auto loc         = m_sessionName.find("_xml/");
        auto sessionName = m_sessionName.substr(0, loc);
        std::string paradir =
            sessionName + "_" + std::to_string(m_iterPIT) + ".pit";
        if (!fs::is_directory(paradir))
        {
            fs::create_directory(paradir);
        }
        std::string outname = paradir + "/" + sessionName + "_timeLevel" +
                              std::to_string(m_session->GetTimeLevel()) + "_" +
                              std::to_string(n);
        WriteFld(outname + ".chk");
    }
}
 
/**
 * Write the n-th checkpoint file.
 * @param   n   The index of the checkpoint file.
 */
void EquationSystem::Checkpoint_Output(
    const int n, MultiRegions::ExpListSharedPtr &field,
    std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
    std::vector<std::string> &variables)
{
    if (!m_comm->IsParallelInTime())
    {
        // Serial-in-time
        std::string outname = m_sessionName + "_" + std::to_string(n);
        WriteFld(outname, field, fieldcoeffs, variables);
    }
    else
    {
        ASSERTL0(false, "Not Implemented for Parallel-in-Time");
    }
}
 
/**
 * Write the n-th base flow into a .chk file
 * @param   n   The index of the base flow file.
 */
void EquationSystem::Checkpoint_BaseFlow(const int n)
{
    std::string outname = m_sessionName + "_BaseFlow_" + std::to_string(n);
 
    WriteFld(outname + ".chk");
}
 
/**
 * Writes the field data to a file with the given filename.
 * @param   outname     Filename to write to.
 */
void EquationSystem::WriteFld(const std::string &outname)
{
    std::vector<Array<OneD, NekDouble>> fieldcoeffs(m_fields.size());
    std::vector<std::string> variables(m_fields.size());
 
    for (int i = 0; i < m_fields.size(); ++i)
    {
        if (m_fields[i]->GetNcoeffs() == m_fields[0]->GetNcoeffs())
        {
            fieldcoeffs[i] = m_fields[i]->UpdateCoeffs();
        }
        else
        {
            fieldcoeffs[i] = Array<OneD, NekDouble>(m_fields[0]->GetNcoeffs());
            m_fields[0]->ExtractCoeffsToCoeffs(
                m_fields[i], m_fields[i]->GetCoeffs(), fieldcoeffs[i]);
        }
        variables[i] = m_boundaryConditions->GetVariable(i);
    }
 
    ExtraFldOutput(fieldcoeffs, variables);
 
    WriteFld(outname, m_fields[0], fieldcoeffs, variables);
}
 
/**
 * Writes the field data to a file with the given filename.
 * @param   outname         Filename to write to.
 * @param   field           ExpList on which data is based.
 * @param   fieldcoeffs     An array of array of expansion coefficients.
 * @param   variables       An array of variable names.
 */
void EquationSystem::WriteFld(const std::string &outname,
                              MultiRegions::ExpListSharedPtr &field,
                              std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
                              std::vector<std::string> &variables)
{
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
        field->GetFieldDefinitions();
    std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
 
    // Copy Data into FieldData and set variable
    for (int j = 0; j < fieldcoeffs.size(); ++j)
    {
        for (int i = 0; i < FieldDef.size(); ++i)
        {
            // Could do a search here to find correct variable
            FieldDef[i]->m_fields.push_back(variables[j]);
            field->AppendFieldData(FieldDef[i], FieldData[i], fieldcoeffs[j]);
        }
    }
 
    // Update time in field info if required
    if (m_fieldMetaDataMap.find("Time") != m_fieldMetaDataMap.end())
    {
        m_fieldMetaDataMap["Time"] = boost::lexical_cast<std::string>(m_time);
    }
 
    // Update step in field info if required
    if (m_fieldMetaDataMap.find("ChkFileNum") != m_fieldMetaDataMap.end())
    {
        m_fieldMetaDataMap["ChkFileNum"] = std::to_string(m_nchk);
    }
 
    // If necessary, add mapping information to metadata
    //      and output mapping coordinates
    Array<OneD, MultiRegions::ExpListSharedPtr> fields(1);
    fields[0] = field;
    GlobalMapping::MappingSharedPtr mapping =
        GlobalMapping::Mapping::Load(m_session, fields);
    LibUtilities::FieldMetaDataMap fieldMetaDataMap(m_fieldMetaDataMap);
    mapping->Output(fieldMetaDataMap, outname);
 
    // If necessary, add informaton for moving frame reference to metadata
    // X, Y, Z translational displacements
    // Theta_x, Theta_y, Theta_z angular displacements
    // U, V, W translational velocity
    // Omega_x, Omega_y, Omega_z angular velocity
    // A_x, A_y, A_z translational acceleration
    // DOmega_x, DOmega_y, DOmega_z angular acceleration
    // X0, Y0, Z0 pivot point
    for (size_t i = 0;
         i < m_strFrameData.size() && i < m_movingFrameData.size(); ++i)
    {
        fieldMetaDataMap[m_strFrameData[i]] =
            boost::lexical_cast<std::string>(m_movingFrameData[i]);
    }
 
    m_fld->Write(outname, FieldDef, FieldData, fieldMetaDataMap,
                 m_session->GetBackups());
}
 
/**
 * Import field from infile and load into \a m_fields. This routine will
 * also perform a \a BwdTrans to ensure data is in both the physical and
 * coefficient storage.
 * @param   infile  Filename to read.
 */
void EquationSystem::ImportFld(
    const std::string &infile,
    Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
{
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
    std::vector<std::vector<NekDouble>> FieldData;
    LibUtilities::FieldIOSharedPtr field_fld =
        LibUtilities::FieldIO::CreateForFile(m_session, infile);
    field_fld->Import(infile, FieldDef, FieldData);
 
    // Copy FieldData into m_fields
    for (int j = 0; j < pFields.size(); ++j)
    {
        Vmath::Zero(pFields[j]->GetNcoeffs(), pFields[j]->UpdateCoeffs(), 1);
 
        for (int i = 0; i < FieldDef.size(); ++i)
        {
            ASSERTL1(FieldDef[i]->m_fields[j] == m_session->GetVariable(j),
                     std::string("Order of ") + infile +
                         std::string(" data and that defined in "
                                     "m_boundaryconditions differs"));
 
            pFields[j]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                            FieldDef[i]->m_fields[j],
                                            pFields[j]->UpdateCoeffs());
        }
        pFields[j]->BwdTrans(pFields[j]->GetCoeffs(), pFields[j]->UpdatePhys());
    }
}
 
/**
 * Import field from infile and load into \a m_fields. This routine will
 * also perform a \a BwdTrans to ensure data is in both the physical and
 * coefficient storage.
 * @param   infile  Filename to read.
 * If optionan \a ndomains is specified it assumes we loop over nodmains
 * for each nvariables.
 */
void EquationSystem::ImportFldToMultiDomains(
    const std::string &infile,
    Array<OneD, MultiRegions::ExpListSharedPtr> &pFields, const int ndomains)
{
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
    std::vector<std::vector<NekDouble>> FieldData;
 
    LibUtilities::Import(infile, FieldDef, FieldData);
 
    int nvariables = GetNvariables();
 
    ASSERTL0(
        ndomains * nvariables == pFields.size(),
        "Number of fields does not match the number of variables and domains");
 
    // Copy FieldData into m_fields
    for (int j = 0; j < ndomains; ++j)
    {
        for (int i = 0; i < nvariables; ++i)
        {
            Vmath::Zero(pFields[j * nvariables + i]->GetNcoeffs(),
                        pFields[j * nvariables + i]->UpdateCoeffs(), 1);
 
            for (int n = 0; n < FieldDef.size(); ++n)
            {
                ASSERTL1(FieldDef[n]->m_fields[i] == m_session->GetVariable(i),
                         std::string("Order of ") + infile +
                             std::string(" data and that defined in "
                                         "m_boundaryconditions differs"));
 
                pFields[j * nvariables + i]->ExtractDataToCoeffs(
                    FieldDef[n], FieldData[n], FieldDef[n]->m_fields[i],
                    pFields[j * nvariables + i]->UpdateCoeffs());
            }
            pFields[j * nvariables + i]->BwdTrans(
                pFields[j * nvariables + i]->GetCoeffs(),
                pFields[j * nvariables + i]->UpdatePhys());
        }
    }
}
 
/**
 * Import field from infile and load into \a pField. This routine will
 * also perform a \a BwdTrans to ensure data is in both the physical and
 * coefficient storage.
 */
void EquationSystem::ImportFld(const std::string &infile,
                               MultiRegions::ExpListSharedPtr &pField,
                               std::string &pFieldName)
{
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
    std::vector<std::vector<NekDouble>> FieldData;
 
    LibUtilities::FieldIOSharedPtr field_fld =
        LibUtilities::FieldIO::CreateForFile(m_session, infile);
    field_fld->Import(infile, FieldDef, FieldData);
    int idx = -1;
 
    Vmath::Zero(pField->GetNcoeffs(), pField->UpdateCoeffs(), 1);
 
    for (int i = 0; i < FieldDef.size(); ++i)
    {
        // find the index of the required field in the file.
        for (int j = 0; j < FieldData.size(); ++j)
        {
            if (FieldDef[i]->m_fields[j] == pFieldName)
            {
                idx = j;
            }
        }
        ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");
 
        pField->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                    FieldDef[i]->m_fields[idx],
                                    pField->UpdateCoeffs());
    }
    pField->BwdTrans(pField->GetCoeffs(), pField->UpdatePhys());
}
 
/**
 * Import field from infile and load into the array \a coeffs.
 *
 * @param infile   Filename to read.
 * @param fieldStr an array of string identifying fields to be imported
 * @param coeffs   array of array of coefficients to store imported data
 */
void EquationSystem::ImportFld(const std::string &infile,
                               std::vector<std::string> &fieldStr,
                               Array<OneD, Array<OneD, NekDouble>> &coeffs)
{
 
    ASSERTL0(fieldStr.size() <= coeffs.size(),
             "length of fieldstr should be the same as coeffs");
 
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
    std::vector<std::vector<NekDouble>> FieldData;
 
    LibUtilities::FieldIOSharedPtr field_fld =
        LibUtilities::FieldIO::CreateForFile(m_session, infile);
    field_fld->Import(infile, FieldDef, FieldData);
 
    // Copy FieldData into m_fields
    for (int j = 0; j < fieldStr.size(); ++j)
    {
        Vmath::Zero(coeffs[j].size(), coeffs[j], 1);
        for (int i = 0; i < FieldDef.size(); ++i)
        {
            m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
                                             fieldStr[j], coeffs[j]);
        }
    }
}
 
/**
 * Write out a summary of the session data.
 * @param   out         Output stream to write data to.
 */
void EquationSystem::SessionSummary(SummaryList &s)
{
    if (m_session->DefinesSolverInfo("EQTYPE"))
    {
        AddSummaryItem(s, "EquationType", m_session->GetSolverInfo("EQTYPE"));
    }
    AddSummaryItem(s, "Session Name", m_sessionName);
    AddSummaryItem(s, "Spatial Dim.", m_spacedim);
    AddSummaryItem(s, "Max SEM Exp. Order",
                   m_fields[0]->EvalBasisNumModesMax());
 
    if (m_session->GetComm()->GetSize() > 1)
    {
        AddSummaryItem(s, "Num. Processes", m_session->GetComm()->GetSize());
    }
 
    if (m_HomogeneousType == eHomogeneous1D)
    {
        AddSummaryItem(s, "Quasi-3D", "Homogeneous in z-direction");
        AddSummaryItem(s, "Expansion Dim.", m_expdim + 1);
        AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
        AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
        AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
        if (m_halfMode)
        {
            AddSummaryItem(s, "ModeType", "Half Mode");
        }
        else if (m_singleMode)
        {
            AddSummaryItem(s, "ModeType", "Single Mode");
        }
        else if (m_multipleModes)
        {
            AddSummaryItem(s, "ModeType", "Multiple Modes");
        }
    }
    else if (m_HomogeneousType == eHomogeneous2D)
    {
        AddSummaryItem(s, "Quasi-3D", "Homogeneous in yz-plane");
        AddSummaryItem(s, "Expansion Dim.", m_expdim + 2);
        AddSummaryItem(s, "Num. Hom. Modes (y)", m_npointsY);
        AddSummaryItem(s, "Num. Hom. Modes (z)", m_npointsZ);
        AddSummaryItem(s, "Hom. length (LY)", m_LhomY);
        AddSummaryItem(s, "Hom. length (LZ)", m_LhomZ);
        AddSummaryItem(s, "FFT Type", m_useFFT ? "FFTW" : "MVM");
    }
    else
    {
        AddSummaryItem(s, "Expansion Dim.", m_expdim);
    }
 
    if (m_session->DefinesSolverInfo("UpwindType"))
    {
        AddSummaryItem(s, "Riemann Solver",
                       m_session->GetSolverInfo("UpwindType"));
    }
 
    if (m_session->DefinesSolverInfo("AdvectionType"))
    {
        std::string AdvectionType;
        AdvectionType = m_session->GetSolverInfo("AdvectionType");
        AddSummaryItem(
            s, "Advection Type",
            GetAdvectionFactory().GetClassDescription(AdvectionType));
    }
 
    if (m_session->DefinesSolverInfo("DiffusionType"))
    {
        std::string DiffusionType;
        DiffusionType = m_session->GetSolverInfo("DiffusionType");
        AddSummaryItem(
            s, "Diffusion Type",
            GetDiffusionFactory().GetClassDescription(DiffusionType));
    }
 
    if (m_projectionType == MultiRegions::eGalerkin)
    {
        AddSummaryItem(s, "Projection Type", "Continuous Galerkin");
    }
    else if (m_projectionType == MultiRegions::eDiscontinuous)
    {
        AddSummaryItem(s, "Projection Type", "Discontinuous Galerkin");
    }
    else if (m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
    {
        AddSummaryItem(s, "Projection Type",
                       "Mixed Continuous Galerkin and Discontinuous");
    }
}
 
/**
 *
 */
Array<OneD, bool> EquationSystem::v_GetSystemSingularChecks()
{
    return Array<OneD, bool>(m_session->GetVariables().size(), false);
}
 
/**
 *
 */
MultiRegions::ExpListSharedPtr EquationSystem::v_GetPressure()
{
    ASSERTL0(false, "This function is not valid for the Base class");
    MultiRegions::ExpListSharedPtr null;
    return null;
}
 
/**
 *
 */
void EquationSystem::v_ExtraFldOutput(
    [[maybe_unused]] std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
    [[maybe_unused]] std::vector<std::string> &variables)
{
}
 
} // namespace Nektar::SolverUtils