PulseWavePropagation.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File PulseWavePropagation.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: Pulse Wave Propagation solve routines based on the weak
// formulation (1):
//
///////////////////////////////////////////////////////////////////////////////
 
#include <iostream>
 
#include <LibUtilities/BasicUtils/Timer.h>
#include <PulseWaveSolver/EquationSystems/PulseWavePropagation.h>
 
using namespace std;
 
namespace Nektar
{
 
string PulseWavePropagation::className =
    GetEquationSystemFactory().RegisterCreatorFunction(
        "PulseWavePropagation", PulseWavePropagation::create,
        "Pulse Wave Propagation equation.");
/**
 *  @class PulseWavePropagation
 *
 *  Set up the routines based on the weak formulation from
 *  "Computational Modelling of 1D blood flow with variable
 *  mechanical properties" by S. J. Sherwin et al. The weak
 *  formulation (1) reads:
 *  \f$ \sum_{e=1}^{N_{el}} \left[ \left( \frac{\partial \mathbf{U}^{\delta}
 * }{\partial t} , \mathbf{\psi}^{\delta} \right)_{\Omega_e} - \left(
 * \frac{\partial \mathbf{F(\mathbf{U})}^{\delta} }
 *    {\partial x}, \mathbf{\psi}^{\delta}  \right)_{\Omega_e} + \left[
 * \mathbf{\psi}^{\delta} \cdot \{ \mathbf{F}^u -
 * \mathbf{F}(\mathbf{U}^{\delta}) \} \right]_{x_e^l}^{x_eû} \right] = 0 \f$
 */
PulseWavePropagation::PulseWavePropagation(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const SpatialDomains::MeshGraphSharedPtr &pGraph)
    : PulseWaveSystem(pSession, pGraph)
{
}
 
void PulseWavePropagation::v_InitObject(bool DeclareField)
{
    // Will set up an array of vessels/fields in PulseWaveSystem::v_InitObject
    // so set DeclareField to false so that the fields are not set up in
    // EquationSystem unnecessarily. Note the number of fields in Equation
    // system is related to the number of variables. The number of vessels is
    // therefore held in PulwWaveSystem.
    PulseWaveSystem::v_InitObject(false);
 
    if (m_session->DefinesSolverInfo("PressureArea"))
    {
        m_pressureArea = GetPressureAreaFactory().CreateInstance(
            m_session->GetSolverInfo("PressureArea"), m_vessels, m_session);
    }
    else
    {
        m_pressureArea = GetPressureAreaFactory().CreateInstance(
            "Beta", m_vessels, m_session);
    }
 
    if (m_explicitAdvection)
    {
        m_ode.DefineOdeRhs(&PulseWavePropagation::DoOdeRhs, this);
        m_ode.DefineProjection(&PulseWavePropagation::DoOdeProjection, this);
    }
    else
    {
        ASSERTL0(false, "Implicit Pulse Wave Propagation not set up.");
    }
 
    // Create advection object
    string advName;
    string riemName;
    switch (m_upwindTypePulse)
    {
        case eUpwindPulse:
        {
            advName  = "WeakDG";
            riemName = "UpwindPulse";
        }
        break;
        default:
        {
            ASSERTL0(false, "populate switch statement for upwind flux");
        }
        break;
    }
    m_advObject =
        SolverUtils::GetAdvectionFactory().CreateInstance(advName, advName);
    m_advObject->SetFluxVector(&PulseWavePropagation::GetFluxVector, this);
    m_riemannSolver = SolverUtils::GetRiemannSolverFactory().CreateInstance(
        riemName, m_session);
    m_riemannSolver->SetScalar("A0", &PulseWavePropagation::GetA0, this);
    m_riemannSolver->SetScalar("beta", &PulseWavePropagation::GetBeta, this);
    m_riemannSolver->SetScalar("alpha", &PulseWavePropagation::GetAlpha, this);
    m_riemannSolver->SetScalar("N", &PulseWavePropagation::GetN, this);
    m_riemannSolver->SetParam("rho", &PulseWavePropagation::GetRho, this);
    m_riemannSolver->SetParam("domains", &PulseWavePropagation::GetDomains,
                              this);
 
    m_advObject->SetRiemannSolver(m_riemannSolver);
    m_advObject->InitObject(m_session, m_fields);
}
 
PulseWavePropagation::~PulseWavePropagation()
{
}
 
/**
 *  Computes the right hand side of (1). The RHS is everything
 *  except the term that contains the time derivative
 *  \f$\frac{\partial \mathbf{U}}{\partial t}\f$. In case of a
 *  Discontinuous Galerkin projection, m_advObject->Advect
 *  will be called
 *
 */
void PulseWavePropagation::DoOdeRhs(
    const Array<OneD, const Array<OneD, NekDouble>> &inarray,
    Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
{
    int i;
 
    Array<OneD, Array<OneD, NekDouble>> physarray(m_nVariables);
 
    // Dummy array for WeakDG advection
    Array<OneD, Array<OneD, NekDouble>> advVel(m_spacedim);
 
    // Output array for advection
    Array<OneD, Array<OneD, NekDouble>> out(m_nVariables);
 
    int cnt = 0;
 
    // Set up Inflow and Outflow boundary conditions.
    SetPulseWaveBoundaryConditions(inarray, outarray, time);
 
    // Set up any interface conditions and write into boundary condition
    EnforceInterfaceConditions(inarray);
 
    // do advection evaluation in all domains
    for (int omega = 0; omega < m_nDomains; ++omega)
    {
        LibUtilities::Timer timer;
        m_currentDomain = omega;
        int nq          = m_vessels[omega * m_nVariables]->GetTotPoints();
 
        timer.Start();
        for (i = 0; i < m_nVariables; ++i)
        {
            physarray[i] = inarray[i] + cnt;
            out[i]       = outarray[i] + cnt;
        }
 
        for (i = 0; i < m_nVariables; ++i)
        {
            m_fields[i] = m_vessels[omega * m_nVariables + i];
        }
 
        m_advObject->Advect(m_nVariables, m_fields, advVel, physarray, out,
                            time);
        for (i = 0; i < m_nVariables; ++i)
        {
            Vmath::Neg(nq, out[i], 1);
        }
        timer.Stop();
        timer.AccumulateRegion("PulseWavePropagation:_DoOdeRHS", 1);
        cnt += nq;
    }
}
 
void PulseWavePropagation::DoOdeProjection(
    const Array<OneD, const Array<OneD, NekDouble>> &inarray,
    Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
{
    // Just copy over array
    for (int i = 0; i < m_nVariables; ++i)
    {
        Vmath::Vcopy(inarray[i].size(), inarray[i], 1, outarray[i], 1);
    }
}
 
/**
 *  Does the projection between ... space and the ... space. Also checks for
 *Q-inflow boundary conditions at the inflow of the current arterial segment and
 *applies the Q-inflow if specified
 */
void PulseWavePropagation::SetPulseWaveBoundaryConditions(
    const Array<OneD, const Array<OneD, NekDouble>> &inarray,
    Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
 
{
    int omega;
 
    Array<OneD, MultiRegions::ExpListSharedPtr> vessel(2);
 
    int offset = 0;
 
    // This will be moved to the RCR boundary condition once factory is setup
    if (time == 0)
    {
        m_Boundary = Array<OneD, PulseWaveBoundarySharedPtr>(2 * m_nDomains);
 
        for (omega = 0; omega < m_nDomains; ++omega)
        {
            vessel[0] = m_vessels[2 * omega];
            vessel[1] = m_vessels[2 * omega + 1];
 
            for (int j = 0; j < 2; ++j)
            {
                std::string BCType;
 
                if (j < vessel[0]->GetBndConditions().size())
                {
                    BCType = vessel[0]->GetBndConditions()[j]->GetUserDefined();
                }
 
                // if no condition given define it to be NoUserDefined
                if (BCType.empty() || BCType == "Interface")
                {
                    BCType = "NoUserDefined";
                }
 
                m_Boundary[2 * omega + j] = GetBoundaryFactory().CreateInstance(
                    BCType, m_vessels, m_session, m_pressureArea);
 
                // turn on time dependent BCs
                if (BCType == "Q-inflow")
                {
                    vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
                }
                else if (BCType == "A-inflow")
                {
                    vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
                }
                else if (BCType == "U-inflow")
                {
                    vessel[1]->GetBndConditions()[j]->SetIsTimeDependent(true);
                }
                else if (BCType == "RCR-terminal")
                {
                    vessel[0]->GetBndConditions()[j]->SetIsTimeDependent(true);
                }
            }
        }
    }
 
    SetBoundaryConditions(time);
 
    // Loop over all vessels and set boundary conditions
    LibUtilities::Timer timer;
    for (omega = 0; omega < m_nDomains; ++omega)
    {
        timer.Start();
        for (int n = 0; n < 2; ++n)
        {
            m_Boundary[2 * omega + n]->DoBoundary(
                inarray, m_A_0, m_beta, m_alpha, time, omega, offset, n);
        }
 
        offset += m_vessels[2 * omega]->GetTotPoints();
        timer.Stop();
        timer.AccumulateRegion("PulseWavePropagation:_SetBCs", 1);
    }
}
 
/**
 *  Calculates the second term of the weak form (1): \f$
 *  \left( \frac{\partial \mathbf{F(\mathbf{U})}^{\delta}
 *  }{\partial x}, \mathbf{\psi}^{\delta} \right)_{\Omega_e}
 *  \f$
 *  The variables of the system are $\mathbf{U} = [A,u]^T$
 *  physfield[0] = A        physfield[1] = u
 *  flux[0] = F[0] = A*u    flux[1] = F[1] = u^2/2 + p/rho
 */
void PulseWavePropagation::GetFluxVector(
    const Array<OneD, Array<OneD, NekDouble>> &physfield,
    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &flux)
{
    int nq = m_vessels[m_currentDomain * m_nVariables]->GetTotPoints();
    NekDouble domain   = m_currentDomain;
    m_pressure[domain] = Array<OneD, NekDouble>(nq);
    Array<OneD, NekDouble> dAUdx(nq);
    NekDouble viscoelasticGradient = 0.0;
 
    LibUtilities::Timer timer;
 
    for (int j = 0; j < nq; ++j)
    {
        timer.Start();
        flux[0][0][j] = physfield[0][j] * physfield[1][j];
        timer.Stop();
        timer.AccumulateRegion("PulseWavePropagation:GetFluxVector-flux", 3);
    }
 
    // d/dx of AU, for the viscoelastic tube law and extra fields
    m_fields[0]->PhysDeriv(flux[0][0], dAUdx);
 
    for (int j = 0; j < nq; ++j)
    {
        if ((j == 0) || (j == nq - 1))
        {
            viscoelasticGradient = dAUdx[j];
        }
        else
        {
            viscoelasticGradient = (dAUdx[j] + dAUdx[j + 1]) / 2;
        }
 
        m_pressureArea->GetPressure(m_pressure[domain][j], m_beta[domain][j],
                                    physfield[0][j], m_A_0[domain][j],
                                    viscoelasticGradient, m_gamma[domain][j],
                                    m_alpha[domain][j]);
 
        flux[1][0][j] = physfield[1][j] * physfield[1][j] / 2 +
                        m_pressure[domain][j] / m_rho;
    }
 
    m_session->MatchSolverInfo("OutputExtraFields", "True", extraFields, true);
 
    if (extraFields)
    {
        /*
        Calculates wave speed and characteristic variables.
 
        Ideally this should be moved to PulseWaveSystem, but it's easiest to
        implement here.
        */
        int counter = 0;
 
        m_PWV[domain] = Array<OneD, NekDouble>(nq);
        m_W1[domain]  = Array<OneD, NekDouble>(nq);
        m_W2[domain]  = Array<OneD, NekDouble>(nq);
 
        for (int j = 0; j < nq; ++j)
        {
            m_pressureArea->GetC(m_PWV[domain][j], m_beta[domain][j],
                                 physfield[0][counter + j], m_A_0[domain][j],
                                 m_alpha[domain][j]);
            m_pressureArea->GetW1(m_W1[domain][j], physfield[1][counter + j],
                                  m_beta[domain][j], physfield[0][counter + j],
                                  m_A_0[domain][j], m_alpha[domain][j]);
            m_pressureArea->GetW2(m_W2[domain][j], physfield[1][counter + j],
                                  m_beta[domain][j], physfield[0][counter + j],
                                  m_A_0[domain][j], m_alpha[domain][j]);
        }
 
        counter += nq;
    }
}
 
Array<OneD, NekDouble> &PulseWavePropagation::GetA0()
{
    return m_A_0_trace[m_currentDomain];
}
 
Array<OneD, NekDouble> &PulseWavePropagation::GetBeta()
{
    return m_beta_trace[m_currentDomain];
}
 
Array<OneD, NekDouble> &PulseWavePropagation::GetAlpha()
{
    return m_alpha_trace[m_currentDomain];
}
 
Array<OneD, NekDouble> &PulseWavePropagation::GetN()
{
    return m_trace_fwd_normal[m_currentDomain];
}
 
NekDouble PulseWavePropagation::GetRho()
{
    return m_rho;
}
 
NekDouble PulseWavePropagation::GetDomains()
{
    return m_nDomains;
}
 
/**
 *  Print summary routine, calls virtual routine reimplemented in
 *  UnsteadySystem
 */
void PulseWavePropagation::v_GenerateSummary(SolverUtils::SummaryList &s)
{
    PulseWaveSystem::v_GenerateSummary(s);
}
 
} // namespace Nektar