Nektar::SolverUtils::DriverAdaptive Class Reference

Base class for the adaptive polynomial order driver. More...

#include <DriverAdaptive.h>

Inheritance diagram for Nektar::SolverUtils::DriverAdaptive:

Collaboration diagram for Nektar::SolverUtils::DriverAdaptive:

Static Public Member Functions
static DriverSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession)
	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	DriverAdaptive (const LibUtilities::SessionReaderSharedPtr pSession)
	Constructor. More...
virtual SOLVER_UTILS_EXPORT	~DriverAdaptive ()
	Destructor. More...
virtual SOLVER_UTILS_EXPORT void	v_InitObject (ostream &out=cout)
	Second-stage initialisation. More...
virtual SOLVER_UTILS_EXPORT void	v_Execute (ostream &out=cout)
	Virtual function for solve implementation. More...
SOLVER_UTILS_EXPORT void	ReplaceExpansion (Array< OneD, MultiRegions::ExpListSharedPtr > &fields, map< int, int > deltaP)
	Update EXPANSIONS tag inside XML schema to reflect new polynomial order distribution. More...
Protected Member Functions inherited from Nektar::SolverUtils::Driver
	Driver (const LibUtilities::SessionReaderSharedPtr pSession)
	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)

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< DriverAdaptive >

Additional Inherited Members
Public Member Functions inherited from Nektar::SolverUtils::Driver
virtual	~Driver ()
	Destructor. More...
SOLVER_UTILS_EXPORT void	InitObject (ostream &out=cout)
	Initialise Object. More...
SOLVER_UTILS_EXPORT void	Execute (ostream &out=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)
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...
Array< OneD, EquationSystemSharedPtr >	m_equ
	Equation system to solve. More...
int	m_nequ
	number of equations More...
enum EvolutionOperatorType	m_EvolutionOperator
	Evolution Operator. More...

Detailed Description

Base class for the adaptive polynomial order driver.

Definition at line 48 of file DriverAdaptive.h.

Constructor & Destructor Documentation

Nektar::SolverUtils::DriverAdaptive::DriverAdaptive ( const LibUtilities::SessionReaderSharedPtr  pSession )

protected

Constructor.

Definition at line 56 of file DriverAdaptive.cpp.

    : Driver(pSession)
{
}

Nektar::SolverUtils::DriverAdaptive::~DriverAdaptive ( )

protectedvirtual

Destructor.

Definition at line 65 of file DriverAdaptive.cpp.

{
}

Member Function Documentation

static DriverSharedPtr Nektar::SolverUtils::DriverAdaptive::create ( const LibUtilities::SessionReaderSharedPtr &  pSession )

inlinestatic

Creates an instance of this class.

Definition at line 54 of file DriverAdaptive.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

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

void Nektar::SolverUtils::DriverAdaptive::ReplaceExpansion ( Array< OneD, MultiRegions::ExpListSharedPtr > &  fields, map< int, int >  deltaP  )

protected

Update EXPANSIONS tag inside XML schema to reflect new polynomial order distribution.

Parameters

fields	Input fields.
deltaP	Map of polynomial order expansions

Definition at line 378 of file DriverAdaptive.cpp.

References Nektar::LibUtilities::BasisTypeMap, Nektar::ParseUtils::GenerateSeqString(), Nektar::SolverUtils::Driver::m_session, and Nektar::LibUtilities::ShapeTypeMap.

Referenced by v_Execute().

{
    int nExp, nDim;

    // Get field definitions
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs =
        fields[0]->GetFieldDefinitions();

    int expDim = 2;

    if (fielddefs[0]->m_numHomogeneousDir == 1)
    {
        nDim = 3;
    }
    else
    {
        nDim = 2;
    }

    // Add variables to field definition
    for (int i = 0; i < fielddefs.size(); ++i)
    {
        for (int j = 0; j < fields.num_elements(); ++j)
        {
            fielddefs[i]->m_fields.push_back(m_session->GetVariable(j));
        }
    }

    // Get tinyxml objects
    TiXmlElement *exp_tag = m_session->GetElement("NEKTAR/EXPANSIONS");

    // Clear current expansions
    exp_tag->Clear();

    // Write new expansion information
    for (int f = 0; f < fielddefs.size(); ++f)
    {
        nExp = fielddefs[f]->m_elementIDs.size();

        //---------------------------------------------
        // Write ELEMENTS
        TiXmlElement *elemTag = new TiXmlElement("ELEMENTS");
        exp_tag->LinkEndChild(elemTag);

        // Write FIELDS
        std::string fieldsString;
        {
            std::stringstream fieldsStringStream;
            bool first = true;
            for (std::vector<int>::size_type i = 0;
                 i < fielddefs[f]->m_fields.size(); i++)
            {
                if (!first)
                {
                    fieldsStringStream << ",";
                }
                fieldsStringStream << fielddefs[f]->m_fields[i];
                first = false;
            }
            fieldsString = fieldsStringStream.str();
        }
        elemTag->SetAttribute("FIELDS", fieldsString);

        // Write SHAPE
        std::string shapeString;
        {
            std::stringstream shapeStringStream;
            shapeStringStream
                << LibUtilities::ShapeTypeMap[fielddefs[f]->m_shapeType];
            if (fielddefs[f]->m_numHomogeneousDir == 1)
            {
                shapeStringStream << "-HomogenousExp1D";
            }
            else if (fielddefs[f]->m_numHomogeneousDir == 2)
            {
                shapeStringStream << "-HomogenousExp2D";
            }

            shapeString = shapeStringStream.str();
        }
        elemTag->SetAttribute("SHAPE", shapeString);

        // Write BASIS
        std::string basisString;
        {
            std::stringstream basisStringStream;
            bool first = true;
            for (std::vector<LibUtilities::BasisType>::size_type i = 0;
                 i < fielddefs[f]->m_basis.size(); i++)
            {
                if (!first)
                {
                    basisStringStream << ",";
                }
                basisStringStream
                    << LibUtilities::BasisTypeMap[fielddefs[f]->m_basis[i]];
                first = false;
            }
            basisString = basisStringStream.str();
        }
        elemTag->SetAttribute("BASIS", basisString);

        // Write homogeneuous length details
        if (fielddefs[f]->m_numHomogeneousDir)
        {
            std::string homoLenString;
            {
                std::stringstream homoLenStringStream;
                bool first = true;
                for (int i = 0; i < fielddefs[f]->m_numHomogeneousDir; ++i)
                {
                    if (!first)
                    {
                        homoLenStringStream << ",";
                    }
                    homoLenStringStream
                        << fielddefs[f]->m_homogeneousLengths[i];
                    first = false;
                }
                homoLenString = homoLenStringStream.str();
            }
            elemTag->SetAttribute("HOMOGENEOUSLENGTHS", homoLenString);
        }

        // Write homogeneuous planes/lines details
        if (fielddefs[f]->m_numHomogeneousDir)
        {
            if (fielddefs[f]->m_homogeneousYIDs.size() > 0)
            {
                std::string homoYIDsString;
                {
                    std::stringstream homoYIDsStringStream;
                    bool first = true;
                    for (int i = 0; i < fielddefs[f]->m_homogeneousYIDs.size();
                         i++)
                    {
                        if (!first)
                        {
                            homoYIDsStringStream << ",";
                        }
                        homoYIDsStringStream
                            << fielddefs[f]->m_homogeneousYIDs[i];
                        first = false;
                    }
                    homoYIDsString = homoYIDsStringStream.str();
                }
                elemTag->SetAttribute("HOMOGENEOUSYIDS", homoYIDsString);
            }

            if (fielddefs[f]->m_homogeneousZIDs.size() > 0)
            {
                std::string homoZIDsString;
                {
                    std::stringstream homoZIDsStringStream;
                    bool first = true;
                    for (int i = 0; i < fielddefs[f]->m_homogeneousZIDs.size();
                         i++)
                    {
                        if (!first)
                        {
                            homoZIDsStringStream << ",";
                        }
                        homoZIDsStringStream
                            << fielddefs[f]->m_homogeneousZIDs[i];
                        first = false;
                    }
                    homoZIDsString = homoZIDsStringStream.str();
                }
                elemTag->SetAttribute("HOMOGENEOUSZIDS", homoZIDsString);
            }
        }

        // Write NUMMODESPERDIR
        std::string numModesString;
        {
            std::stringstream numModesStringStream;

            numModesStringStream << "MIXORDER:";
            bool first = true;
            int eId;
            Array<OneD, int> order(nDim, 0);
            for (int n = 0; n < nExp; n++)
            {
                eId = fielddefs[f]->m_elementIDs[n];

                for (int i = 0; i < expDim; i++)
                {
                    order[i] = deltaP[eId];
                }

                for (int i = 0; i < nDim; i++)
                {
                    if (fielddefs[f]->m_uniOrder)
                    {
                        order[i] += fielddefs[f]->m_numModes[i];
                    }
                    else
                    {
                        order[i] += fielddefs[f]->m_numModes[n * nDim + i];
                    }

                    if (!first)
                    {
                        numModesStringStream << ",";
                    }

                    numModesStringStream << order[i];
                    first = false;
                }
            }
            numModesString = numModesStringStream.str();
        }
        elemTag->SetAttribute("NUMMODESPERDIR", numModesString);

        // Write IDs. Should ideally look at ways of compressing this stream if
        // just sequential.
        elemTag->SetAttribute(
            "ID", ParseUtils::GenerateSeqString(fielddefs[f]->m_elementIDs));
    }
}

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

protectedvirtual

Virtual function for solve implementation.

Implements Nektar::SolverUtils::Driver.

Definition at line 77 of file DriverAdaptive.cpp.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::eOrtho_A, Nektar::LibUtilities::eOrtho_B, Nektar::LibUtilities::eQuadrilateral, Nektar::LibUtilities::eTriangle, Nektar::GlobalMapping::Mapping::Load(), Nektar::SolverUtils::Driver::m_comm, Nektar::SolverUtils::Driver::m_equ, Nektar::SolverUtils::Driver::m_session, Nektar::GlobalMapping::MappingSharedPtr, Nektar::LibUtilities::ReduceMax, ReplaceExpansion(), Nektar::SolverUtils::Driver::v_InitObject(), Vmath::Vmul(), Vmath::Vsub(), and Vmath::Vsum().

{
    time_t starttime, endtime;
    NekDouble CPUtime;

    m_equ[0]->PrintSummary(out);

    // First run using original order
    time(&starttime);
    m_equ[0]->DoInitialise();

    // Obtain initial time in case a restart was used
    NekDouble startTime = m_equ[0]->GetFinalTime();
    m_equ[0]->DoSolve();

    // Load session parameters and solver info
    bool isHomogeneous1D;
    int nRuns, minP, maxP, sensorVar;
    NekDouble lowerTol, upperTol;
    m_session->LoadParameter  ("NumRuns",                nRuns,     1);
    m_session->LoadParameter  ("AdaptiveMinModes",       minP,      4);
    m_session->LoadParameter  ("AdaptiveMaxModes",       maxP,      12);
    m_session->LoadParameter  ("AdaptiveLowerTolerance", lowerTol,  1e-8);
    m_session->LoadParameter  ("AdaptiveUpperTolerance", upperTol,  1e-6);
    m_session->LoadParameter  ("AdaptiveSensorVariable", sensorVar, 0);
    m_session->MatchSolverInfo("Homogeneous", "1D", isHomogeneous1D, false);

    // Get number of elements and planes
    int nExp, nPlanes;
    if (isHomogeneous1D)
    {
        nExp    = m_equ[0]->UpdateFields()[0]->GetPlane(0)->GetExpSize();
        nPlanes = m_equ[0]->UpdateFields()[0]->GetZIDs().num_elements();
    }
    else
    {
        nExp    = m_equ[0]->UpdateFields()[0]->GetExpSize();
        nPlanes = 1;
    }

    int       nFields  = m_equ[0]->UpdateFields().num_elements();
    int       numSteps = m_session->GetParameter("NumSteps");
    NekDouble period   = m_session->GetParameter("TimeStep") * numSteps;

    Array<OneD, NekDouble> coeffs, phys, physReduced, tmpArray;

    // Get mapping
    GlobalMapping::MappingSharedPtr mapping;
    mapping = GlobalMapping::Mapping::Load(m_session,
                                           m_equ[0]->UpdateFields());

    // Adaptive loop
    LocalRegions::ExpansionSharedPtr Exp;
    for (int i = 1; i < nRuns; i++)
    {
        // Get field expansions
        Array<OneD, MultiRegions::ExpListSharedPtr> fields =
            m_equ[0]->UpdateFields();

        // Determine the change to be applied in the order
        map<int, int> deltaP;
        int offset = 0;
        for (int n = 0; n < nExp; n++)
        {
            offset = fields[sensorVar]->GetPhys_Offset(n);
            Exp    = fields[sensorVar]->GetExp(n);

            int P  = Exp->GetBasis(0)->GetNumModes();
            int Q  = Exp->GetBasis(1)->GetNumModes();
            int qa = Exp->GetBasis(0)->GetNumPoints();
            int qb = Exp->GetBasis(1)->GetNumPoints();

            // Declare orthogonal basis.
            LibUtilities::PointsKey pa(qa, Exp->GetBasis(0)->GetPointsType());
            LibUtilities::PointsKey pb(qb, Exp->GetBasis(1)->GetPointsType());

            StdRegions::StdExpansionSharedPtr OrthoExp;
            switch (Exp->GetGeom()->GetShapeType())
            {
                case LibUtilities::eQuadrilateral:
                {
                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P - 1,
                                              pa);
                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, Q - 1,
                                              pb);
                    OrthoExp = MemoryManager<
                        StdRegions::StdQuadExp>::AllocateSharedPtr(Ba, Bb);
                    break;
                }
                case LibUtilities::eTriangle:
                {
                    LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, P - 1,
                                              pa);
                    LibUtilities::BasisKey Bb(LibUtilities::eOrtho_B, Q - 1,
                                              pb);
                    OrthoExp =
                        MemoryManager<StdRegions::StdTriExp>::AllocateSharedPtr(
                            Ba, Bb);
                    break;
                }
                default:
                    ASSERTL0(false, "Shape not supported.");
                    break;
            }

            int nq = OrthoExp->GetTotPoints();

            NekDouble error = 0;
            NekDouble fac   = 0;
            NekDouble tmp   = 0;

            coeffs      = Array<OneD, NekDouble>(OrthoExp->GetNcoeffs());
            physReduced = Array<OneD, NekDouble>(OrthoExp->GetTotPoints());
            tmpArray    = Array<OneD, NekDouble>(OrthoExp->GetTotPoints(), 0.0);

            // Refinement based only on one variable
            for (int plane = 0; plane < nPlanes; plane++)
            {
                if (isHomogeneous1D)
                {
                    phys =
                        fields[sensorVar]->GetPlane(plane)->GetPhys() + offset;
                }
                else
                {
                    phys = fields[sensorVar]->GetPhys() + offset;
                }

                // Project solution to lower order
                OrthoExp->FwdTrans(phys, coeffs);
                OrthoExp->BwdTrans(coeffs, physReduced);

                // Calculate error =||phys-physReduced||^2 / ||phys||^2
                Vmath::Vsub(nq, phys,     1, physReduced, 1, tmpArray, 1);
                Vmath::Vmul(nq, tmpArray, 1, tmpArray,    1, tmpArray, 1);
                tmp = Exp->Integral(tmpArray);

                Vmath::Vmul(nq, phys, 1, phys, 1, tmpArray, 1);
                fac = Exp->Integral(tmpArray);

                tmp = abs(tmp / fac);

                if (tmp != tmp)
                {
                    ASSERTL0(false,
                             "Adaptive procedure encountered NaN value.");
                }

                // Get maximum value along planes
                error = (tmp > error) ? tmp : error;
            }

            // Combine planes from different processes
            m_session->GetComm()->GetColumnComm()->AllReduce(
                error, LibUtilities::ReduceMax);

            // Override tolerances if function is defined
            if (m_session->DefinesFunction("AdaptiveLowerTolerance"))
            {
                int nq = Exp->GetTotPoints();

                // Obtain points from the the element
                Array<OneD, NekDouble> xc0, xc1, xc2;
                xc0 = Array<OneD, NekDouble>(nq, 0.0);
                xc1 = Array<OneD, NekDouble>(nq, 0.0);
                xc2 = Array<OneD, NekDouble>(nq, 0.0);
                Exp->GetCoords(xc0, xc1, xc2);

                // Evaluate function from session file
                Array<OneD, NekDouble> tolerance(nq, 0.0);
                LibUtilities::EquationSharedPtr ffunc =
                    m_session->GetFunction("AdaptiveLowerTolerance", 0);
                ffunc->Evaluate(xc0, xc1, xc2, tolerance);
                lowerTol = Vmath::Vsum(nq, tolerance, 1) / nq;
            }

            if (m_session->DefinesFunction("AdaptiveUpperTolerance"))
            {
                int nq = Exp->GetTotPoints();

                // Obtain points from the the element
                Array<OneD, NekDouble> xc0, xc1, xc2;
                xc0 = Array<OneD, NekDouble>(nq, 0.0);
                xc1 = Array<OneD, NekDouble>(nq, 0.0);
                xc2 = Array<OneD, NekDouble>(nq, 0.0);
                Exp->GetCoords(xc0, xc1, xc2);

                // Evaluate function from session file
                Array<OneD, NekDouble> tolerance(nq, 0.0);
                LibUtilities::EquationSharedPtr ffunc =
                    m_session->GetFunction("AdaptiveUpperTolerance", 0);
                ffunc->Evaluate(xc0, xc1, xc2, tolerance);
                upperTol = Vmath::Vsum(nq, tolerance, 1) / nq;
            }

            // Determine what to do with the polynomial order
            if ((error > upperTol) && (P < maxP))
            {
                deltaP[Exp->GetGeom()->GetGlobalID()] = 1;
            }
            else if ((error < lowerTol) && P > minP)
            {
                deltaP[Exp->GetGeom()->GetGlobalID()] = -1;
            }
            else
            {
                deltaP[Exp->GetGeom()->GetGlobalID()] = 0;
            }
        }

        // Write new expansion section to the session reader
        ReplaceExpansion(fields, deltaP);

        // Reset GlobalLinSys Manager to avoid using too much memory
        //
        // @todo This could be made better by replacing individual matrices
        //       within the linear system.
        LibUtilities::NekManager<MultiRegions::GlobalLinSysKey,
                                 MultiRegions::GlobalLinSys>::
            ClearManager(std::string("GlobalLinSys"));

        int chkNumber = m_equ[0]->GetCheckpointNumber();
        int chkSteps  = m_equ[0]->GetCheckpointSteps();

        // Initialise driver again
        Driver::v_InitObject(out);

        // Update mapping (must be before m_equ[0]->DoInitialise();)
        mapping->ReplaceField(m_equ[0]->UpdateFields());

        // Set chkSteps to zero to avoid writing initial condition
        m_equ[0]->SetCheckpointSteps(0);

        // Initialise equation
        m_equ[0]->DoInitialise();
        m_equ[0]->SetInitialStep(i * numSteps);
        m_equ[0]->SetSteps(i * numSteps + numSteps);
        m_equ[0]->SetTime(startTime + i * period);
        m_equ[0]->SetBoundaryConditions(startTime + i * period);
        m_equ[0]->SetCheckpointNumber(chkNumber);
        m_equ[0]->SetCheckpointSteps(chkSteps);

        // Project solution to new expansion
        for (int n = 0; n < nFields; n++)
        {
            m_equ[0]->UpdateFields()[n]->ExtractCoeffsToCoeffs(
                fields[n], fields[n]->GetCoeffs(),
                m_equ[0]->UpdateFields()[n]->UpdateCoeffs());
            m_equ[0]->UpdateFields()[n]->BwdTrans_IterPerExp(
                m_equ[0]->UpdateFields()[n]->GetCoeffs(),
                m_equ[0]->UpdateFields()[n]->UpdatePhys());
        }

        // Solve equation
        m_equ[0]->DoSolve();
    }

    time(&endtime);

    m_equ[0]->Output();

    if (m_comm->GetRank() == 0)
    {
        CPUtime = difftime(endtime, starttime);
        cout << "-------------------------------------------" << endl;
        cout << "Total Computation Time = " << CPUtime << "s" << endl;
        cout << "-------------------------------------------" << endl;
    }

    // Evaluate and output computation time and solution accuracy.
    // The specific format of the error output is essential for the
    // regression tests to work.

    // Evaluate L2 Error
    for (int i = 0; i < m_equ[0]->GetNvariables(); ++i)
    {
        Array<OneD, NekDouble> exactsoln(m_equ[0]->GetTotPoints(), 0.0);

        // Evaluate "ExactSolution" function, or zero array
        m_equ[0]->EvaluateExactSolution(i, exactsoln, m_equ[0]->GetFinalTime());

        NekDouble vL2Error   = m_equ[0]->L2Error  (i, exactsoln);
        NekDouble vLinfError = m_equ[0]->LinfError(i, exactsoln);

        if (m_comm->GetRank() == 0)
        {
            out << "L 2 error (variable " << m_equ[0]->GetVariable(i)
                << ") : " << vL2Error << endl;
            out << "L inf error (variable " << m_equ[0]->GetVariable(i)
                << ") : " << vLinfError << endl;
        }
    }
}

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

protectedvirtual

Second-stage initialisation.

Reimplemented from Nektar::SolverUtils::Driver.

Definition at line 72 of file DriverAdaptive.cpp.

References Nektar::SolverUtils::Driver::v_InitObject().

{
    Driver::v_InitObject(out);
}

Friends And Related Function Documentation

friend class MemoryManager< DriverAdaptive >

friend

Definition at line 51 of file DriverAdaptive.h.

Member Data Documentation

string Nektar::SolverUtils::DriverAdaptive::className

static

Initial value:

= GetDriverFactory().RegisterCreatorFunction(
    "Adaptive", DriverAdaptive::create)

Name of the class.

Definition at line 64 of file DriverAdaptive.h.

string Nektar::SolverUtils::DriverAdaptive::driverLookupId

staticprotected

Initial value:

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

Definition at line 84 of file DriverAdaptive.h.