SmoothedProfileMethod.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: SmoothedProfileMethod.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).
//
// License for the specific language governing rights and limitations under
// 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: Smoothed Profile Method for the Incompressible
// Navier Stokes equations
///////////////////////////////////////////////////////////////////////////////
 
#include <IncNavierStokesSolver/EquationSystems/SmoothedProfileMethod.h>
#include <IncNavierStokesSolver/Filters/FilterAeroForcesSPM.h>
#include <MultiRegions/ContField.h>
#include <MultiRegions/ContField3DHomogeneous1D.h>
#include <MultiRegions/ContField3DHomogeneous2D.h>
 
using namespace std;
 
namespace Nektar
{
using namespace MultiRegions;
 
string SmoothedProfileMethod::className =
    SolverUtils::GetEquationSystemFactory().RegisterCreatorFunction(
        "SmoothedProfileMethod", SmoothedProfileMethod::create);
 
string SmoothedProfileMethod::solverTypeLookupId =
    LibUtilities::SessionReader::RegisterEnumValue(
        "SolverType", "SmoothedProfileMethod", eSmoothedProfileMethod);
 
/**
 * @brief Construct a new Smoothed Profile Method object
 *
 * @param pSession
 * @param pGraph
 */
SmoothedProfileMethod::SmoothedProfileMethod(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const SpatialDomains::MeshGraphSharedPtr &pGraph)
    : UnsteadySystem(pSession, pGraph),
      VelocityCorrectionScheme(pSession, pGraph)
{
}
 
/**
 * @brief Destroy the Smoothed Profile Method object
 *
 */
SmoothedProfileMethod::~SmoothedProfileMethod(void)
{
}
 
void SmoothedProfileMethod::v_InitObject(bool DeclareField)
{
    VelocityCorrectionScheme::v_InitObject(DeclareField);
 
    // Update implicit time-intregration class operators
    m_ode.DefineImplicitSolve(
        &SmoothedProfileMethod::v_SolveUnsteadyStokesSystem, this);
 
    // Number of dims as number of velocity vectors
    size_t nvel = m_velocity.size();
 
    // Initialization of correction pressure and shape function
    switch (nvel)
    {
        case 1:
            if (m_projectionType == eGalerkin)
            {
                SetUpExpansions<ContField>(nvel);
            }
            else if (m_projectionType == eDiscontinuous)
            {
                SetUpExpansions<DisContField>(nvel);
            }
            break;
 
        case 2:
            if (m_projectionType == eGalerkin)
            {
                SetUpExpansions<ContField>(nvel);
            }
            else if (m_projectionType == eDiscontinuous)
            {
                SetUpExpansions<DisContField>(nvel);
            }
            break;
 
        case 3:
            if (m_projectionType == eGalerkin)
            {
                if (m_HomogeneousType == EquationSystem::eNotHomogeneous)
                {
                    SetUpExpansions<ContField>(nvel);
                }
                else if (m_HomogeneousType == EquationSystem::eHomogeneous1D)
                {
                    SetUpExpansions<ContField3DHomogeneous1D>(nvel);
                }
                else if (m_HomogeneousType == EquationSystem::eHomogeneous2D ||
                         m_HomogeneousType == EquationSystem::eHomogeneous3D)
                {
                    SetUpExpansions<ContField3DHomogeneous2D>(nvel);
                }
            }
            else if (m_projectionType == eDiscontinuous)
            {
                if (m_HomogeneousType == EquationSystem::eNotHomogeneous)
                {
                    SetUpExpansions<DisContField>(nvel);
                }
                else if (m_HomogeneousType == EquationSystem::eHomogeneous1D)
                {
                    SetUpExpansions<DisContField3DHomogeneous1D>(nvel);
                }
                else if (m_HomogeneousType == EquationSystem::eHomogeneous2D ||
                         m_HomogeneousType == EquationSystem::eHomogeneous3D)
                {
                    SetUpExpansions<DisContField3DHomogeneous2D>(nvel);
                }
            }
            break;
    }
 
    // Read 'm_phi' and its velocity
    ASSERTL0(m_session->DefinesFunction("ShapeFunction"),
             "ShapeFunction must be defined in the session file.")
    ReadPhi();
 
    // Allocate the vector 'm_up'
    size_t physTot = m_phi->GetTotPoints();
    m_velName.push_back("Up");
    if (nvel > 1)
    {
        m_velName.push_back("Vp");
    }
    if (nvel == 3)
    {
        m_velName.push_back("Wp");
    }
 
    m_up = Array<OneD, Array<OneD, NekDouble>>(nvel);
    for (size_t i = 0; i < nvel; ++i)
    {
        m_up[i] = Array<OneD, NekDouble>(physTot, 0.0);
    }
 
    // Make sure that m_phi and m_up are defined
    UpdatePhiUp(0.0);
 
    // Get the time integration scheme.
    LibUtilities::TimeIntScheme timeInt;
    if (m_session->DefinesTimeIntScheme())
    {
        timeInt = m_session->GetTimeIntScheme();
    }
    else
    {
        timeInt.method = m_session->GetSolverInfo("TimeIntegrationMethod");
        timeInt.order  = timeInt.method.back() - '0';
 
        // Remove everything past the IMEX.
        timeInt.method = timeInt.method.substr(0, 4);
    }
 
    // Select 'm_gamma0' depending on IMEX order
    ASSERTL0(
        boost::iequals(timeInt.method, "IMEX") && 1 <= timeInt.order &&
            timeInt.order <= 4,
        "The TimeIntegrationMethod scheme must be IMEX with order '1' to '4'.")
 
    switch (timeInt.order)
    {
        case 1:
            m_gamma0 = 1.0;
            break;
 
        case 2:
            m_gamma0 = 3.0 / 2.0;
            break;
 
        case 3:
            m_gamma0 = 11.0 / 6.0;
            break;
 
        case 4:
            m_gamma0 = 25.0 / 12.0;
            break;
    }
 
    // Check if the aeroforces filter is active, negative if inactive
    m_forcesFilter = -1;
    for (size_t i = 0; i < m_session->GetFilters().size(); ++i)
    {
        if (boost::iequals(m_session->GetFilters()[i].first, "AeroForcesSPM"))
        {
            m_forcesFilter = i;
            break;
        }
    }
}
 
/**
 * @brief Generates the summary of the current simulation
 *
 * @param s
 */
void SmoothedProfileMethod::v_GenerateSummary(SolverUtils::SummaryList &s)
{
    VelocityCorrectionScheme::v_GenerateSummary(s);
    SolverUtils::AddSummaryItem(s, "IB formulation",
                                "Smoothed Profile Method (SPM)");
}
 
/**
 * @brief Linear terms due to pressure and visosity are calculated here.
 * After solving the velocity filed without taking into account the
 * immersed boundaries, a new correction is applied through the force
 * \f$f_s\f$:
 *
 * \f[ \mathbf{f_s} = \frac{\Phi^{n+1}(\mathbf{u_p}-\mathbf{u^*})}
 * {\Delta t} \f]
 *
 * @param inarray
 * @param outarray
 * @param time
 * @param a_iixDt
 */
void SmoothedProfileMethod::v_SolveUnsteadyStokesSystem(
    const Array<OneD, const Array<OneD, NekDouble>> &inarray,
    Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time,
    const NekDouble a_iixDt)
{
    VelocityCorrectionScheme::v_SolveUnsteadyStokesSystem(inarray, outarray,
                                                          time, a_iixDt);
 
    size_t physTot = m_pressureP->GetNpoints();
 
    /* SPM correction of velocity */
    // Update 'm_phi' and 'm_up' if needed (evaluated at next time step)
    UpdatePhiUp(time + a_iixDt);
    // Update calculation of IB forcing 'm_fs'
    UpdateForcing(outarray, a_iixDt);
    // Estimate forces only if requested
    if (m_forcesFilter >= 0)
    {
        static_pointer_cast<FilterAeroForcesSPM>(
            m_filters[m_forcesFilter].second)
            ->CalculateForces(outarray, m_upPrev, m_phi, time, a_iixDt);
    }
    // Set BC conditions for pressure p_p
    SetUpCorrectionPressure(outarray, m_F);
    // Solve Poisson equation for pressure p_p
    SolveCorrectionPressure(m_F[0]);
    // Solve velocity in the next step with IB
    SolveCorrectedVelocity(m_F, outarray, a_iixDt);
 
    // Add pressures to get final value
    Vmath::Vadd(physTot, m_pressure->GetPhys(), 1, m_pressureP->GetPhys(), 1,
                m_pressure->UpdatePhys(), 1);
    m_pressure->FwdTrans(m_pressure->GetPhys(), m_pressure->UpdateCoeffs());
 
    // Add presure to outflow bc if using convective like BCs
    m_extrapolation->AddPressureToOutflowBCs(m_kinvis);
}
 
/**
 * @brief Sets the forcing term of the equation for the correction pressure
 * \f$p_p\f$:
 *
 * \f[ \nabla\cdot\mathbf{f_s} \f]
 *
 * @param fields
 * @param Forcing
 */
void SmoothedProfileMethod::SetUpCorrectionPressure(
    const Array<OneD, const Array<OneD, NekDouble>> &fields,
    Array<OneD, Array<OneD, NekDouble>> &Forcing)
{
    boost::ignore_unused(fields);
 
    size_t physTot = m_fs[0]->GetNpoints();
    size_t nvel    = m_velocity.size();
 
    // Set boundary conditions
    SetCorrectionPressureBCs();
 
    // Divergence of 'fs'
    m_fields[m_velocity[0]]->PhysDeriv(eX, m_fs[0]->GetPhys(), Forcing[0]);
 
    // Using 'Forcing[1]' as storage
    for (size_t i = 1; i < nvel; ++i)
    {
        size_t ind = m_velocity[i];
        m_fields[ind]->PhysDeriv(DirCartesianMap[i], m_fs[i]->GetPhys(),
                                 Forcing[1]);
        Vmath::Vadd(physTot, Forcing[1], 1, Forcing[0], 1, Forcing[0], 1);
    }
}
 
/**
 * @brief Solves the Poisson equation for the correction pressure
 * \f$p_p\f$:
 *
 * \f[ \nabla^2 p_p = \nabla\cdot\mathbf{f_s} \f]
 *
 * @param Forcing
 */
void SmoothedProfileMethod::SolveCorrectionPressure(
    const Array<OneD, NekDouble> &Forcing)
{
    StdRegions::ConstFactorMap factors;
    // Factor 'lambda=0' in Helmholtz equation to get the Poisson form
    factors[StdRegions::eFactorLambda] = 0.0;
 
    // Solve the Poisson equation
    m_pressureP->HelmSolve(Forcing, m_pressureP->UpdateCoeffs(), factors);
 
    // Update node values from coefficients
    m_pressureP->BwdTrans(m_pressureP->GetCoeffs(), m_pressureP->UpdatePhys());
}
 
/**
 * @brief Corrects the velocity field so that the IBs are taken into
 * account. Solves the explicit equation:
 *
 * \f[ \frac{\gamma_0(\mathbf{u_p}^{n+1} - \mathbf{u}^*)}{\Delta t} =
 *     \mathbf{f_s} - \nabla p_p \f]
 *
 * @param Forcing
 * @param fields
 * @param dt
 */
void SmoothedProfileMethod::SolveCorrectedVelocity(
    Array<OneD, Array<OneD, NekDouble>> &Forcing,
    Array<OneD, Array<OneD, NekDouble>> &fields, NekDouble dt)
{
    size_t physTot = m_phi->GetNpoints();
 
    // Gradient of p_p
    size_t nvel = m_velocity.size();
    if (nvel == 2)
    {
        m_pressureP->PhysDeriv(m_pressureP->GetPhys(), Forcing[0], Forcing[1]);
    }
    else
    {
        m_pressureP->PhysDeriv(m_pressureP->GetPhys(), Forcing[0], Forcing[1],
                               Forcing[2]);
    }
 
    // Velocity correction
    for (size_t i = 0; i < nvel; ++i)
    {
        // Adding -(1-m_phi)*grad(p_p) instead of -grad(p_p) reduces the
        // flux through the walls, but the flow is not incompressible
        if (m_session->DefinesSolverInfo("ForceBoundary") &&
            boost::iequals(m_session->GetSolverInfo("ForceBoundary"), "True"))
        {
            Vmath::Vvtvm(physTot, m_phi->GetPhys(), 1, Forcing[i], 1,
                         Forcing[i], 1, Forcing[i], 1);
            Vmath::Vadd(physTot, m_fs[i]->GetPhys(), 1, Forcing[i], 1,
                        Forcing[i], 1);
        }
        else
        {
            Vmath::Vsub(physTot, m_fs[i]->GetPhys(), 1, Forcing[i], 1,
                        Forcing[i], 1);
        }
        Blas::Daxpy(physTot, dt / m_gamma0, Forcing[i], 1, fields[i], 1);
    }
}
 
/**
 * @brief Updates the BCs for boundaries with Neumann or periodic BCs in
 * the pressure:
 *
 * \f[ \frac{\partial p_p}{\partial\mathbf{n}} =
 *     \mathbf{f_s}\cdot\mathbf{n} \f]
 */
void SmoothedProfileMethod::SetCorrectionPressureBCs()
{
    size_t nvel = m_velocity.size();
    Array<OneD, ExpListSharedPtr> BndExp;
    Array<OneD, SpatialDomains::BoundaryConditionShPtr> BndCond;
 
    // Get the BC expansions
    BndExp  = m_pressureP->GetBndCondExpansions();
    BndCond = m_pressureP->GetBndConditions();
 
    // For each boundary...
    for (size_t b = 0; b < BndExp.size(); ++b)
    {
        // Only for BCs based on the derivative
        if (BndCond[b]->GetBoundaryConditionType() ==
                SpatialDomains::eNeumann ||
            BndCond[b]->GetBoundaryConditionType() == SpatialDomains::ePeriodic)
        {
            // Calculate f_s values
            Array<OneD, Array<OneD, NekDouble>> f_s(nvel);
            for (size_t i = 0; i < nvel; ++i)
            {
                f_s[i] = m_fs[0]->GetBndCondExpansions()[b]->GetPhys();
            }
 
            // BC is f_s * n
            BndExp[b]->NormVectorIProductWRTBase(f_s, BndExp[b]->UpdatePhys());
        }
    }
}
 
/**
 * @brief Calculates the values of the shape function
 *
 * @param t
 */
void SmoothedProfileMethod::UpdatePhiUp(NekDouble t)
{
    // Initialise 'm_up' and 'm_phi' during first step
    if (t <= 0.0)
    {
        if (!m_filePhi)
        {
            // Update 'm_phi' only if it was provided as a function
            m_phiEvaluator->Evaluate("Phi", m_phi->UpdatePhys(), t);
        }
 
        // Initialize both variables for the first step
        m_phiEvaluator->Evaluate(m_velName, m_up, t);
 
        // Initialise 'm_upPrev' in all cases
        m_upPrev = m_up;
    }
    // If timedependent 'm_phi'
    // Phi functions from files are not timedependent
    else if (m_timeDependentPhi)
    {
        m_phiEvaluator->Evaluate("Phi", m_phi->UpdatePhys(), t);
 
        // And if velocities are timedependent as well
        if (m_timeDependentUp)
        {
            // Store previous value of u_p during simulation
            m_upPrev = m_up;
            m_phiEvaluator->Evaluate(m_velName, m_up, t);
        }
    }
}
 
/**
 * @brief For a body with a velocity \f$\mathbf{u_p}\f$, the force
 * \f$\mathbf{f_s}\f$ applied to the fluid ensures that the IBC are met:
 *
 * \f[ \mathbf{f_s} = \frac{\Phi^{n+1}\left(\mathbf{u_p}^{n+1} -
 * \mathbf{u^*}\right)}{\Delta t} \f]
 *
 * @param fields
 * @param dt
 * @param f_s
 */
void SmoothedProfileMethod::UpdateForcing(
    const Array<OneD, const Array<OneD, NekDouble>> &fields, NekDouble dt)
{
    size_t nvel = m_velocity.size();
    size_t nq   = m_phi->GetNpoints();
 
    for (size_t i = 0; i < nvel; ++i)
    {
        // In homogeneous cases, switch out of wave space
        Array<OneD, NekDouble> tmpField(nq);
        size_t ind = m_velocity[i];
 
        if (m_HomogeneousType != EquationSystem::eNotHomogeneous &&
            m_fields[ind]->GetWaveSpace())
        {
            m_fields[ind]->HomogeneousBwdTrans(nq, fields[i], tmpField);
            m_fs[i]->HomogeneousBwdTrans(nq, m_fs[i]->GetPhys(),
                                         m_fs[i]->UpdatePhys());
        }
        else
        {
            tmpField = fields[i];
        }
 
        Vmath::Vsub(nq, m_up[i], 1, tmpField, 1, m_fs[i]->UpdatePhys(), 1);
        Vmath::Vmul(nq, m_phi->GetPhys(), 1, m_fs[i]->GetPhys(), 1,
                    m_fs[i]->UpdatePhys(), 1);
        Vmath::Smul(nq, m_gamma0 / dt, m_fs[i]->GetPhys(), 1,
                    m_fs[i]->UpdatePhys(), 1);
 
        // And go back to wave space if the 'if' was executed
        if (m_HomogeneousType != EquationSystem::eNotHomogeneous &&
            m_fields[ind]->GetWaveSpace())
        {
            m_fs[i]->HomogeneousFwdTrans(nq, m_fs[i]->GetPhys(),
                                         m_fs[i]->UpdatePhys());
        }
    }
}
 
/**
 * @brief True if the function is timedependent, false otherwise
 *
 * @param name
 * @param type
 * @param attribute
 * @return string
 */
bool SmoothedProfileMethod::GetVarTimeDependence(string funcName,
                                                 string elemName)
{
    // Get the handler of the function
    TiXmlElement *function = GetFunctionHdl(funcName);
 
    // Go to the first element
    TiXmlElement *functionDef = function->FirstChildElement();
    ASSERTL0(functionDef, "At least one element must be defined in " + funcName)
 
    // And search the element with name 'elemName'
    string varName = functionDef->Attribute("VAR");
    while (functionDef && !boost::iequals(varName, elemName))
    {
        functionDef = functionDef->NextSiblingElement();
        varName     = functionDef->Attribute("VAR");
    }
 
    ASSERTL0(functionDef,
             "Variable " + elemName + " must be defined in " + funcName + ".");
 
    // And return the value of USERDEFINEDTYPE
    string attr;
    int err     = functionDef->QueryStringAttribute("USERDEFINEDTYPE", &attr);
    bool output = boost::iequals(attr, "TimeDependent");
 
    ASSERTL0((err == TIXML_NO_ATTRIBUTE) || (err == TIXML_SUCCESS && output),
             "USERDEFINEDTYPE in " + elemName +
                 " must be TimeDependent if defined");
 
    return output;
}
 
/**
 * @brief Returns a handle to the requested function. Returns NULL if it
 * does not exist
 *
 * @param functionName
 * @return TiXmlElement*
 */
TiXmlElement *SmoothedProfileMethod::GetFunctionHdl(string functionName)
{
    // Get the handler of first function block
    TiXmlElement *conds    = m_session->GetElement("Nektar/Conditions");
    TiXmlElement *function = conds->FirstChildElement("FUNCTION");
 
    // Loop over functions until the block 'name' is found
    string functionType = function->Attribute("NAME");
    while (function && !boost::iequals(functionType, functionName))
    {
        function     = function->NextSiblingElement("FUNCTION");
        functionType = function->Attribute("NAME");
    }
 
    return function;
}
 
void SmoothedProfileMethod::ReadPhi()
{
    // Function evaluator for Phi and Up
    m_phiEvaluator = GetFunction("ShapeFunction");
 
    TiXmlElement *function = GetFunctionHdl("ShapeFunction");
    TiXmlElement *child    = function->FirstChildElement();
    m_filePhi              = false;
 
    // If defined by using a file
    if (boost::iequals(child->ValueStr(), "F"))
    {
        // Get name of STL file
        string fileName;
        int status = child->QueryStringAttribute("FILE", &fileName);
        ASSERTL0(status == TIXML_SUCCESS,
                 "An FLD file with the values "
                 "of the phi function has to be supplied.")
        ASSERTL0(boost::iequals(fileName.substr(fileName.length() - 4), ".fld"),
                 "A valid FLD file must be supplied in the "
                 "'ShapeFunction' field.")
 
        // Get phi values from XML file (after "FieldConvert" the STL file)
        // First, load the data
        std::vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
        std::vector<std::vector<NekDouble>> fieldData;
        LibUtilities::FieldMetaDataMap fieldMetaData;
        LibUtilities::FieldIOSharedPtr phiFile =
            LibUtilities::FieldIO::CreateForFile(m_session, fileName);
        phiFile->Import(fileName, fieldDef, fieldData, fieldMetaData);
 
        // Only Phi field should be defined in the file
        ASSERTL0(fieldData.size() == 1, "Only one field (phi) must be "
                                        "defined in the FLD file.")
 
        // Extract Phi field to output
        string tmp("phi");
        m_phi->ExtractDataToCoeffs(fieldDef[0], fieldData[0], tmp,
                                   m_phi->UpdateCoeffs());
        m_phi->BwdTrans(m_phi->GetCoeffs(), m_phi->UpdatePhys());
        m_filePhi          = true;
        m_timeDependentPhi = false;
        m_timeDependentUp  = false;
    }
    else
    {
        // Check if Phi is timedependent
        m_timeDependentPhi = GetVarTimeDependence("ShapeFunction", "Phi");
 
        // If so, check if its velocity changes as well
        m_timeDependentUp = GetVarTimeDependence("ShapeFunction", "Up");
        switch (m_velocity.size())
        {
            case 2:
                m_timeDependentUp |=
                    GetVarTimeDependence("ShapeFunction", "Vp");
                break;
            case 3:
                m_timeDependentUp |=
                    GetVarTimeDependence("ShapeFunction", "Vp");
                m_timeDependentUp |=
                    GetVarTimeDependence("ShapeFunction", "Wp");
                break;
        }
    }
}
 
} // namespace Nektar