FilterBodyFittedVelocity.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File FilterBodyFittedVelocity.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: Read the body-fitted coordinate system from the file output by
//              FieldConvert module bodyFittedVelocity and compute local max/
//              min/original body-fitted velocity compinents.
//
///////////////////////////////////////////////////////////////////////////////
 
#include <CompressibleFlowSolver/Misc/VariableConverter.h>
#include <SolverUtils/Filters/FilterBodyFittedVelocity.h>
 
#include "FieldUtils/ProcessModules/ProcessBodyFittedVelocity.h"
 
#include <MultiRegions/ContField.h>
#include <MultiRegions/ExpList.h>
 
using std::cout;
using std::endl;
 
namespace Nektar::SolverUtils
{
std::string FilterBodyFittedVelocity::className =
    GetFilterFactory().RegisterCreatorFunction(
        "BodyFittedVelocity", FilterBodyFittedVelocity::create);
 
FilterBodyFittedVelocity::FilterBodyFittedVelocity(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const std::weak_ptr<EquationSystem> &pEquation, const ParamMap &pParams)
    : FilterFieldConvert(pSession, pEquation, pParams)
{
    // Load sampling frequency
    auto it = pParams.find("SampleFrequency");
    if (it == pParams.end())
    {
        m_sampleFrequency = 1;
    }
    else
    {
        LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
        m_sampleFrequency = round(equ.Evaluate());
    }
 
    // Load flag for filter type
    it = pParams.find("OriginalOrMaxOrMin");
    if (it == pParams.end())
    {
        m_filterType = eOriginal;
    }
    else
    {
        std::string sOptionType = it->second.c_str();
        if (boost::iequals(sOptionType, "maximum") ||
            boost::iequals(sOptionType, "max"))
        {
            m_filterType = eMax;
        }
        else if (boost::iequals(sOptionType, "minumun") ||
                 boost::iequals(sOptionType, "min"))
        {
            m_filterType = eMin;
        }
        else if (boost::iequals(sOptionType, "original"))
        {
            m_filterType = eOriginal;
        }
        else
        {
            WARNINGL0(false, "Detailed filter type is not found, use original");
            m_filterType = eOriginal;
        }
    }
 
    // Load bodyFittedCoordinate from a file generated by FieldConvert -m
    // bodyFittedVelocity
    it = pParams.find("BodyFittedCoordinateFile");
    ASSERTL0(it->second.length() > 0,
             "Missing parameter 'BodyFittedCoordinateFile'.");
 
    if (it->second.find_last_of('.') != std::string::npos)
    {
 
        m_bodyFittedCooriateFile = it->second;
    }
    else
    {
        std::stringstream outname;
        outname << it->second << ".fld";
        m_bodyFittedCooriateFile = outname.str();
    }
 
    m_initialized = m_restartFile != "";
}
 
FilterBodyFittedVelocity::~FilterBodyFittedVelocity()
{
}
 
void FilterBodyFittedVelocity::v_Initialise(
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
    const NekDouble &time)
{
    // Initialise output arrays
    FilterFieldConvert::v_Initialise(pFields, time);
 
    // Generate the body-fitted coordinate system and distance field
    // bfcsDir[i][j][k]
    //   i - dir vec:   0-main tangential, 1-normal, (2-minor tangential)
    //   j - component: 0-x, 1-y, (2-z)
    //   k - point id
    const int npoints = pFields[0]->GetTotPoints();
    m_bfcsDir = Array<OneD, Array<OneD, Array<OneD, NekDouble>>>(m_spaceDim);
    for (int i = 0; i < m_spaceDim; ++i)
    {
        m_bfcsDir[i] = Array<OneD, Array<OneD, NekDouble>>(m_spaceDim);
        for (int j = 0; j < m_spaceDim; ++j)
        {
            m_bfcsDir[i][j] = Array<OneD, NekDouble>(npoints);
        }
    }
 
    // Initialize the body fitted coordinate in the file
    // m_bfsVars holds the var names for the coordinate system and will be
    // used to find corresponding field to load from the input file.
    // Ref: FilterFieldConvert.cpp (v_Initialise) & ProcessInterpField.cpp
    // (Process)
    m_bfsVars.resize(m_spaceDim * m_spaceDim); // 4 or 9
    if (m_spaceDim == 2)
    {
        m_bfsVars[0] = "bfc_dir0_x";
        m_bfsVars[1] = "bfc_dir0_y";
        m_bfsVars[2] = "bfc_dir1_x";
        m_bfsVars[3] = "bfc_dir1_y";
    }
    else
    {
        m_bfsVars[0] = "bfc_dir0_x";
        m_bfsVars[1] = "bfc_dir0_y";
        m_bfsVars[2] = "bfc_dir0_z";
        m_bfsVars[3] = "bfc_dir1_x";
        m_bfsVars[4] = "bfc_dir1_y";
        m_bfsVars[5] = "bfc_dir1_z";
        m_bfsVars[6] = "bfc_dir2_x";
        m_bfsVars[7] = "bfc_dir2_y";
        m_bfsVars[8] = "bfc_dir2_z";
    }
 
    // Load file
    std::vector<LibUtilities::FieldDefinitionsSharedPtr> fieldDef;
    std::vector<std::vector<NekDouble>> fieldData;
    LibUtilities::FieldMetaDataMap fieldMetaData;
    LibUtilities::FieldIOSharedPtr fld = LibUtilities::FieldIO::CreateForFile(
        m_session, m_bodyFittedCooriateFile);
    fld->Import(m_bodyFittedCooriateFile, fieldDef, fieldData, fieldMetaData);
 
    // Get a tmp array to hold the coeffs from the input file
    Array<OneD, NekDouble> tmp_coeffs;
    tmp_coeffs = Array<OneD, NekDouble>(pFields[0]->GetNcoeffs(), 0.0);
 
    // Extract data and convert to physical space
    for (int k = 0; k < m_bfsVars.size(); ++k)
    {
        // In the fieldDef, for example fieldDef[0] is the quad and
        // fieldDef[1] is the tri region
        // Each of fieldDef[i] contains all the fields we expected in an fld
        for (int i = 0; i < fieldData.size(); ++i)
        {
            pFields[0]->ExtractDataToCoeffs(fieldDef[i], fieldData[i],
                                            m_bfsVars[k], tmp_coeffs);
        }
 
        // Bwd transorm the data and put them into m_bfcsDir[][]
        int ii, jj;
        jj = k % m_spaceDim;
        ii = k / m_spaceDim;
 
        pFields[0]->BwdTrans(tmp_coeffs, m_bfcsDir[ii][jj]);
    }
 
    // Allocate storage
    m_curFieldsVels_Car.resize(m_spaceDim);
    m_curFieldsVels.resize(m_spaceDim); // m_spaceDim = m_nAddFields-1
    m_outFieldsVels.resize(m_spaceDim);
 
    for (int n = 0; n < m_spaceDim; ++n)
    {
        m_curFieldsVels_Car[n] = Array<OneD, NekDouble>(npoints, 0.0);
        m_curFieldsVels[n]     = Array<OneD, NekDouble>(npoints, 0.0);
        m_outFieldsVels[n]     = Array<OneD, NekDouble>(npoints, 0.0);
    }
 
    // Allocate storage for rho,p,T
    if (m_problemType == eCompressible)
    {
        m_curFieldsThermalVars.resize(3);
        m_outFieldsThermalVars.resize(3);
 
        for (int n = 0; n < 3; ++n)
        {
            m_curFieldsThermalVars[n] = Array<OneD, NekDouble>(npoints, 0.0);
            m_outFieldsThermalVars[n] = Array<OneD, NekDouble>(npoints, 0.0);
        }
    }
 
    // m_outFields contains the initialization values
    // 2d: rho,rhou,rhov,E,u,v,p,T,s,a,Mach,Sensor,distanceToWall,u_bfc,v_bfc
    if (m_initialized)
    {
        for (int n = 0; n < m_spaceDim; ++n)
        {
            pFields[0]->BwdTrans(m_outFields[m_nVars + 1 + n],
                                 m_outFieldsVels[n]);
            if (pFields[0]->GetWaveSpace())
            {
                pFields[0]->HomogeneousBwdTrans(npoints, m_outFieldsVels[n],
                                                m_outFieldsVels[n]);
            }
        }
 
        if (m_problemType == eCompressible)
        {
            // For cfs, initialize thermal vars
            // shift_thermal for [rho, p, T] in m_outFields
            int shift_thermal[3] = {0, m_spaceDim * 2 + 2, m_spaceDim * 2 + 3};
 
            for (int n = 0; n < 3; ++n)
            {
                pFields[0]->BwdTrans(m_outFields[shift_thermal[n]],
                                     m_outFieldsThermalVars[n]);
                if (pFields[0]->GetWaveSpace())
                {
                    pFields[0]->HomogeneousBwdTrans(npoints,
                                                    m_outFieldsThermalVars[n],
                                                    m_outFieldsThermalVars[n]);
                }
            }
        }
    }
    else
    {
        ASSERTL0(false, "Restart file is expectd. It can be generated by \
                         FieldConvert in serial. For example: \
                         FieldConvert -m bodyFittedVelocity:bnd=0:bfcOut=1 \
                         mesh.xml session.xml baseflow.fld bfvFile.fld");
    }
}
 
void FilterBodyFittedVelocity::v_FillVariablesName(
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields)
{
    // Fill in the var names using the routine in base class
    FilterFieldConvert::v_FillVariablesName(pFields);
 
    // Check type of problem
    std::vector<std::string> varNames =
        pFields[0]->GetSession()->GetVariables();
    if (boost::iequals(varNames[0], "u") && boost::iequals(varNames[1], "v"))
    {
        m_problemType = eIncompressible;
        m_spaceDim    = pFields.size() - 1;
    }
    else if (boost::iequals(varNames[0], "rho") &&
             boost::iequals(varNames[1], "rhou"))
    {
        m_problemType = eCompressible;
        m_spaceDim    = pFields.size() - 2;
    }
    else
    {
        WARNINGL0(false, "Problem type is not claear. Please check");
        m_problemType = eOthers;
        m_spaceDim    = 2;
    }
 
    m_nVars      = m_variables.size(); // Include extra vars
    m_nFields    = pFields.size();     // Conservative vars for cfs
    m_nAddFields = 3;
 
    // Fill in the body-fitted velocity names
    // bfc represents body-fitted coordinate
    m_variables.push_back("distanceToWall");
    m_variables.push_back("u_bfc");
    m_variables.push_back("v_bfc");
    if (m_spaceDim == 3)
    {
        m_variables.push_back("w_bfc");
        ++m_nAddFields;
    }
    else if (m_spaceDim != 2)
    {
        ASSERTL0(false, "Unsupported dimension");
    }
}
 
// m_curFieldsVels_Car is the velocity in Cartesian coordinate system in current
//                     step
// m_curFieldsVels     is the velocity in body-fittd coordinate system in
//                     current step
// m_outFieldsVels     is the velocity in body-fittd coordinate
//                     sysetm for output
// m_outFields         is the coefficients for output
//
// Before execusion, the fields are expected to be initialized so that
// m_outFields and m_outFieldsVels are not empty.
//
// *** Note ***
// In order to get the amplitude field, take the difference of
// max (perturbed) and max (perturbation-free) for a better result.
// The difference of max (perturbed) and original (perturbation-free)
// contains noises for unknown reasons at the moment.
void FilterBodyFittedVelocity::v_ProcessSample(
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
    std::vector<Array<OneD, NekDouble>> &fieldcoeffs,
    [[maybe_unused]] const NekDouble &time)
{
    // Use shift_vel to get the current u,v,w in Cartesian coordinate
    // cfs_2D:
    // rho,rhou,rhov,E,u,v,p,T,s,a,Mach,Sensor,distanceToWall,u_bfc,v_bfc
    // inc_2D: u,v,p
    int shift_vel = (m_problemType == eCompressible) ? (m_spaceDim + 2) : 0;
 
    // Get current u,v,w in Cartesian coordinate
    for (int n = 0; n < m_spaceDim; ++n)
    {
        pFields[0]->BwdTrans(fieldcoeffs[shift_vel + n],
                             m_curFieldsVels_Car[n]);
        if (pFields[0]->GetWaveSpace())
        {
            pFields[0]->HomogeneousBwdTrans(pFields[0]->GetNpoints(),
                                            m_curFieldsVels_Car[n],
                                            m_curFieldsVels_Car[n]);
        }
    }
 
    // Project the velocity into the body-fitted coordinate system
    int npoints = pFields[0]->GetTotPoints();
    Array<OneD, NekDouble> vel_tmp(npoints);
 
    for (int i = 0; i < m_spaceDim; ++i) // loop for bfc velocity
    {
        Vmath::Zero(npoints, m_curFieldsVels[i], 1);
 
        for (int j = 0; j < m_spaceDim; ++j)
        {
            Vmath::Vmul(npoints, m_curFieldsVels_Car[j], 1, m_bfcsDir[i][j], 1,
                        vel_tmp, 1);
            Vmath::Vadd(npoints, vel_tmp, 1, m_curFieldsVels[i], 1,
                        m_curFieldsVels[i], 1);
        }
    }
 
    // Get max/min/original for the u/v/w_bfc field
    for (int n = 0; n < m_spaceDim; ++n)
    {
        size_t length = m_outFieldsVels[n].size();
 
        if (m_filterType == eMax)
        {
            // Compute max
            for (int i = 0; i < length; ++i)
            {
                if (m_curFieldsVels[n][i] > m_outFieldsVels[n][i])
                {
                    m_outFieldsVels[n][i] = m_curFieldsVels[n][i];
                }
            }
        }
        else if (m_filterType == eMin)
        {
            // Compute min
            for (int i = 0; i < length; ++i)
            {
                if (m_curFieldsVels[n][i] < m_outFieldsVels[n][i])
                {
                    m_outFieldsVels[n][i] = m_curFieldsVels[n][i];
                }
            }
        }
        else
        {
            // Original field
            m_outFieldsVels[n] = m_curFieldsVels[n];
        }
    }
 
    // Forward transform and put into m_outFields
    for (int n = 0; n < m_spaceDim; ++n)
    {
        pFields[0]->FwdTransLocalElmt(m_outFieldsVels[n],
                                      m_outFields[m_nVars + 1 + n]);
    }
 
    // Process the thermal variables for compressible flows
    if (m_problemType == eCompressible)
    {
        // Get current rho,p,T fields
        // shift_thermal for [rho, p, T] in fieldcoeffs
        int shift_thermal[3] = {0, m_spaceDim * 2 + 2, m_spaceDim * 2 + 3};
 
        for (int n = 0; n < 3; ++n)
        {
            pFields[0]->BwdTrans(fieldcoeffs[shift_thermal[n]],
                                 m_curFieldsThermalVars[n]);
            if (pFields[0]->GetWaveSpace())
            {
                pFields[0]->HomogeneousBwdTrans(pFields[0]->GetNpoints(),
                                                m_curFieldsThermalVars[n],
                                                m_curFieldsThermalVars[n]);
            }
        }
 
        // Get max/min/original for the rho/p/T fields
        for (int n = 0; n < 3; ++n)
        {
            size_t length = m_outFieldsThermalVars[n].size();
 
            if (m_filterType == eMax)
            {
                // Compute max
                for (int i = 0; i < length; ++i)
                {
                    if (m_curFieldsThermalVars[n][i] >
                        m_outFieldsThermalVars[n][i])
                    {
                        m_outFieldsThermalVars[n][i] =
                            m_curFieldsThermalVars[n][i];
                    }
                }
            }
            else if (m_filterType == eMin)
            {
                // Compute min
                for (int i = 0; i < length; ++i)
                {
                    if (m_curFieldsThermalVars[n][i] <
                        m_outFieldsThermalVars[n][i])
                    {
                        m_outFieldsThermalVars[n][i] =
                            m_curFieldsThermalVars[n][i];
                    }
                }
            }
            else
            {
                // Original field
                m_outFieldsThermalVars[n] = m_curFieldsThermalVars[n];
            }
        }
 
        // Forward transform and put into m_outFields
        for (int n = 0; n < 3; ++n)
        {
            pFields[0]->FwdTransLocalElmt(m_outFieldsThermalVars[n],
                                          m_outFields[shift_thermal[n]]);
        }
    }
}
 
void FilterBodyFittedVelocity::v_PrepareOutput(
    [[maybe_unused]] const Array<OneD, const MultiRegions::ExpListSharedPtr>
        &pFields,
    [[maybe_unused]] const NekDouble &time)
{
    m_fieldMetaData["NumberOfFieldDumps"] = std::to_string(m_numSamples);
}
 
NekDouble FilterBodyFittedVelocity::v_GetScale()
{
    return 1.0;
}
 
} // namespace Nektar::SolverUtils