Add synthetic turbulence generation for the compressible solver

CHANGELOG.md

@@ -29,6 +29,9 @@ v5.7.0
 **ShallowWaterSolver**
 - Implement implicit time-discritization (!1784)
 
+**CompressibleSolver**
+- Add synthetic turbulence generator for the compressible solver (!1859)
+
 **NekMesh**
 - Added revolve module (!1825)
 

docs/user-guide/xml/xml-forcing.tex

@@ -40,10 +40,10 @@ This force type specifies the name of a body forcing function expressed in the \
 </FORCE>
 \end{lstlisting}
 
-\subsection{IncNSSyntheticTurbulence}
+\subsection{Synthetic turbulence generator}
 This force type allows the user to apply synthetic turbulence generation in the flow field. The Synthetic Eddy Method is implemented. The approach developed here is based on a source term formulation. This formulation allows the user to apply synthetic turbulence generation in any specific area of the domain, not only in the boundary condition as most methodologies do. So that, after defining a synthetic eddy region (box of eddies), the user can randomly release eddies inside this box which are going to be convected downstream and will produce turbulence depending on the flow conditions. Each eddy  that leaves the synthetic eddy region is reintroduced in the inlet plane of the box, so this mechanism re-energise the system, roughly speaking. 
 
-Below it is shown how to define the Synthetic Eddy Method for a fully three-dimensional Navier-Stokes simulation. Note that this definition is under the \inltt{FORCING} tag. Firstly, in the \inltt{TYPE} entry, we define the force type as \texttt{IncNSSyntheticTurbulence}. In the \inltt{BoxOfEddies} tag, under the \inltt{FORCE} tag, the center plane of the synthetic eddy region is defined. The coordinates of its center are given by \texttt{x0, y0, z0} and lengths of its sides are \texttt{lyref}  and \texttt{lzref} in the $y$- and $z$-directions, respectively. In the \inltt{Sigma} tag, we define the standard deviation (\texttt{sigma}) of the Gaussian function with zero mean, which is used to compute the stochastic signal. After that, the bulk velocity (\texttt{Ub}) of the flow must be provided in the \inltt{BulkVelocity} tag. 
+Below it is shown how to define the Synthetic Eddy Method for a fully three-dimensional Navier-Stokes simulation. Note that this definition is under the \inltt{FORCING} tag. Firstly, in the \inltt{TYPE} entry, we define the force type as \texttt{IncNSSyntheticTurbulence} for the incompressible solver and \texttt{CFSSyntheticTurbulence} for the compressible solver. In the \inltt{BoxOfEddies} tag, under the \inltt{FORCE} tag, the center plane of the synthetic eddy region is defined. The coordinates of its center are given by \texttt{x0, y0, z0} and lengths of its sides are \texttt{lyref}  and \texttt{lzref} in the $y$- and $z$-directions, respectively. Note that the length in the x-direction is defined in the characteristic length scale function (see below), so that \inltt{l00}  defines the value of $l_{x}$. In the \inltt{Sigma} tag, we define the standard deviation (\texttt{sigma}) of the Gaussian function with zero mean, which is used to compute the stochastic signal. After that, the bulk velocity (\texttt{Ub}) of the flow must be provided in the \inltt{BulkVelocity} tag. 
 
 \begin{lstlisting}[style=XMLStyle] 
 <FORCE TYPE="IncNSSyntheticTurbulence">
@@ -59,12 +59,12 @@ In order to define the Reynolds stresses (\inltt{ReynoldsStresses} tag) and the
 
 \begin{lstlisting}[style=XMLStyle] 
 <FUNCTION NAME="ReynoldsStresses">
-    <E VAR="r00" VALUE="1e-6" />
-    <E VAR="r10" VALUE="10*y" />
+    <E VAR="r00" VALUE="1e-3" /> 
+    <E VAR="r10" VALUE="10*y+y^2+5*y^3" />
     <E VAR="r20" VALUE="0.0"  />
-    <E VAR="r11" VALUE="1e-6" />
+    <E VAR="r11" VALUE="1e-3" />
     <E VAR="r21" VALUE="0.0"  />
-    <E VAR="r22" VALUE="1e-6" />
+    <E VAR="r22" VALUE="1e-3" />
  </FUNCTION>
 \end{lstlisting}
 
@@ -84,7 +84,7 @@ Also, in the Synthetic Eddy Method implemented here, an isotropic or anisotropic
 </FUNCTION>
 \end{lstlisting}
 
-Note that the Synthetic Eddy Method is only supported for a fully three-dimensional Incompressible Navier-Stokes simulation.
+Note that the synthetic turbulence generator is only supported for fully three-dimensional simulations.
 
 \subsection{MovingReferenceFrame}
 This force type allows the solution of incompressilbe Navier-Stokes in moving frame of reference. The moving frame is attached the to body and can have translational, rotational or both motions. Although the Navier-Stokes equations are solved in a moving reference frame, our formulation is based on the absolute velocity and pressure (in inertial frame). However, note that these absolute velocities and any other vector quantities are expressed using the coordinate basis of the moving frame. Further, note that if you are using the FilterAeroForces, the force vector $\left(F_x, F_y, F_z\right)$ is automatically converted and output in the inertial frame (ground reference frame).

library/SolverUtils/CMakeLists.txt

@@ -55,6 +55,7 @@ SET(SOLVER_UTILS_SOURCES
   Forcing/ForcingMovingReferenceFrame.cpp
   Forcing/ForcingNoise.cpp
   Forcing/ForcingProgrammatic.cpp
+  Forcing/ForcingSyntheticEddy.cpp
 )
 
 SET(SOLVER_UTILS_HEADERS
@@ -116,6 +117,7 @@ SET(SOLVER_UTILS_HEADERS
   Forcing/ForcingMovingReferenceFrame.h
   Forcing/ForcingNoise.h
   Forcing/ForcingProgrammatic.h
+  Forcing/ForcingSyntheticEddy.h
 )
 
 IF (NEKTAR_USE_ARPACK)

library/SolverUtils/Forcing/ForcingSyntheticEddy.h

+///////////////////////////////////////////////////////////////////////////////
+//
+// File:  ForcingSyntheticEddy.h
+//
+// 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: Derived base class - Synthetic turbulence generation.
+//              This code implements the Synthetic Eddy Method (SEM).
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NEKTAR_SOLVERUTILS_FORCINGSYNTHETICEDDY
+#define NEKTAR_SOLVERUTILS_FORCINGSYNTHETICEDDY
+
+#include <SolverUtils/Forcing/Forcing.h>
+#include <string>
+
+namespace Nektar::SolverUtils
+{
+
+class ForcingSyntheticEddy : virtual public SolverUtils::Forcing
+{
+public:
+    friend class MemoryManager<ForcingSyntheticEddy>;
+
+    /// Creates an instance of this class
+    SOLVER_UTILS_EXPORT static ForcingSharedPtr create(
+        const LibUtilities::SessionReaderSharedPtr &pSession,
+        const std::weak_ptr<EquationSystem> &pEquation,
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
+    {
+        ForcingSharedPtr p =
+            MemoryManager<ForcingSyntheticEddy>::AllocateSharedPtr(pSession,
+                                                                   pEquation);
+        p->InitObject(pFields, pNumForcingFields, pForce);
+        return p;
+    }
+
+    /// Name of the class
+    SOLVER_UTILS_EXPORT static std::string className;
+
+protected:
+    // Initial object
+    SOLVER_UTILS_EXPORT void v_InitObject(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const unsigned int &pNumForcingFields,
+        const TiXmlElement *pForce) override;
+    // Apply forcing term
+    SOLVER_UTILS_EXPORT void v_Apply(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const Array<OneD, Array<OneD, NekDouble>> &inarray,
+        Array<OneD, Array<OneD, NekDouble>> &outarray,
+        const NekDouble &time) override;
+    // Apply forcing term
+    SOLVER_UTILS_EXPORT void v_ApplyCoeff(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const Array<OneD, Array<OneD, NekDouble>> &inarray,
+        Array<OneD, Array<OneD, NekDouble>> &outarray,
+        const NekDouble &time) override;
+    // Set Cholesky decomposition of the Reynolds Stresses in the domain
+    SOLVER_UTILS_EXPORT void SetCholeskyReyStresses(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Set the Number of the Eddies in the box
+    SOLVER_UTILS_EXPORT void SetNumberOfEddies();
+    /// Set reference lengths
+    SOLVER_UTILS_EXPORT void ComputeRefLenghts();
+    /// Set Xi max
+    SOLVER_UTILS_EXPORT void ComputeXiMax();
+    /// Compute Random 1 or -1 value
+    SOLVER_UTILS_EXPORT Array<OneD, Array<OneD, int>> GenerateRandomOneOrMinusOne();
+    /// Compute Constant C
+    SOLVER_UTILS_EXPORT NekDouble ComputeConstantC(int row, int col);
+    /// Compute Gaussian
+    SOLVER_UTILS_EXPORT NekDouble ComputeGaussian(NekDouble coord,
+                                                  NekDouble xiMaxVal,
+                                                  NekDouble constC = 1.0);
+    /// Check if point is inside the box of eddies
+    SOLVER_UTILS_EXPORT bool InsideBoxOfEddies(NekDouble coord0,
+                                               NekDouble coord1,
+                                               NekDouble coord2);
+    /// Compute Stochastic Signal
+    SOLVER_UTILS_EXPORT Array<OneD, Array<OneD, NekDouble>> ComputeStochasticSignal(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Compute Velocity Fluctuation
+    SOLVER_UTILS_EXPORT Array<OneD, Array<OneD, NekDouble>>
+    ComputeVelocityFluctuation(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        Array<OneD, Array<OneD, NekDouble>> stochasticSignal);
+    /// Compute Characteristic Convective Turbulent Time
+    SOLVER_UTILS_EXPORT Array<OneD, Array<OneD, NekDouble>> ComputeCharConvTurbTime(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Compute Smoothing Factor
+    SOLVER_UTILS_EXPORT Array<OneD, Array<OneD, NekDouble>> ComputeSmoothingFactor(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        Array<OneD, Array<OneD, NekDouble>> convTurb);
+    /// Set box of eddies mask
+    SOLVER_UTILS_EXPORT void SetBoxOfEddiesMask(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Initialise forcing term for each eddy
+    SOLVER_UTILS_EXPORT void InitialiseForcingEddy(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Compute the initial position of the eddies inside the box
+    SOLVER_UTILS_EXPORT void ComputeInitialRandomLocationOfEddies();
+    /// Update positions of the eddies inside the box
+    SOLVER_UTILS_EXPORT void UpdateEddiesPositions();
+    /// Remove eddy from forcing term
+    SOLVER_UTILS_EXPORT void RemoveEddiesFromForcing(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Compute the initial location of the eddies for the test case
+    SOLVER_UTILS_EXPORT void ComputeInitialLocationTestCase();
+
+    // Members
+    // Expressions (functions) of the prescribed Reynolds stresses
+    std::map<int, LibUtilities::EquationSharedPtr> m_R;
+    /// Cholesky decomposition of the Reynolds Stresses
+    /// for all points in the mesh
+    Array<OneD, Array<OneD, NekDouble>> m_Cholesky;
+    /// Bulk velocity
+    NekDouble m_Ub;
+    /// Standard deviation
+    NekDouble m_sigma;
+    /// Inlet length in the y-direction and z-direction
+    Array<OneD, NekDouble> m_lyz;
+    /// Center of the inlet plane
+    Array<OneD, NekDouble> m_rc;
+    /// Number of eddies in the box
+    int m_N;
+    /// Characteristic lenght scales
+    Array<OneD, NekDouble> m_l;
+    /// Reference lenghts
+    Array<OneD, NekDouble> m_lref;
+    /// Xi max
+    Array<OneD, NekDouble> m_xiMax;
+    // XiMaxMin - Value form Giangaspero et al. 2022
+    NekDouble m_xiMaxMin = 0.4;
+    /// Space dimension
+    int m_spacedim;
+    /// Ration of specific heats
+    NekDouble m_gamma;
+    /// Box of eddies mask
+    Array<OneD, int> m_mask;
+    /// Eddy position
+    Array<OneD, Array<OneD, NekDouble>> m_eddyPos;
+    /// Check when the forcing should be applied
+    bool m_calcForcing{true};
+    /// Eddies that add to the domain
+    std::vector<unsigned int> m_eddiesIDForcing;
+    /// Current time
+    NekDouble m_currTime;
+    /// Check for test case
+    bool m_tCase;
+    /// Forcing for each eddy
+    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_ForcingEddy;
+
+    SOLVER_UTILS_EXPORT ForcingSyntheticEddy(
+        const LibUtilities::SessionReaderSharedPtr &pSession,
+        const std::weak_ptr<EquationSystem> &pEquation);
+    SOLVER_UTILS_EXPORT ~ForcingSyntheticEddy(void) override = default;
+};
+
+} // namespace Nektar::SolverUtils
+
+#endif

solvers/CompressibleFlowSolver/CMakeLists.txt

@@ -38,6 +38,7 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
         EquationSystems/NavierStokesCFEAxisym.cpp
         EquationSystems/NavierStokesImplicitCFE.cpp
         Forcing/ForcingAxiSymmetric.cpp
+        Forcing/ForcingCFSSyntheticEddy.cpp
         Forcing/ForcingQuasi1D.cpp
        )
 
@@ -95,6 +96,7 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
         EquationSystems/NavierStokesCFE.h
         EquationSystems/NavierStokesImplicitCFE.h
         Forcing/ForcingAxiSymmetric.h
+        Forcing/ForcingCFSSyntheticEddy.h
         Forcing/ForcingQuasi1D.h
         Misc/EquationOfState.h
         Misc/IdealGasEoS.h
@@ -151,6 +153,8 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
     ADD_NEKTAR_TEST(IsentropicVortex_WeakDG_HexDeformed)
     ADD_NEKTAR_TEST(RinglebFlow_P3)
     ADD_NEKTAR_TEST(RinglebFlow_P8 LENGTHY)
+    ADD_NEKTAR_TEST(ChanFlow3D_infTurbExpl)
+    ADD_NEKTAR_TEST(ChanFlow3D_infTurbImpl)
     #ADD_NEKTAR_TEST(Couette_WeakDG_LDG_MODIFIED)
     ADD_NEKTAR_TEST(Couette_WeakDG_IP_MODIFIED_IM)
     #ADD_NEKTAR_TEST(Couette_WeakDG_LDGNS_MODIFIED_IM)

solvers/CompressibleFlowSolver/Forcing/ForcingCFSSyntheticEddy.cpp

+///////////////////////////////////////////////////////////////////////////////
+//
+// File: ForcingCFSSyntheticEddy.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: Derived base class - Synthetic turbulence forcing for the
+//              Compressible solver.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#include <CompressibleFlowSolver/Forcing/ForcingCFSSyntheticEddy.h>
+
+namespace Nektar::SolverUtils
+{
+
+std::string ForcingCFSSyntheticEddy::className =
+    GetForcingFactory().RegisterCreatorFunction(
+        "CFSSyntheticTurbulence", ForcingCFSSyntheticEddy::create,
+        "Compressible Synthetic Eddy Turbulence Forcing (Generation)");
+
+ForcingCFSSyntheticEddy::ForcingCFSSyntheticEddy(
+    const LibUtilities::SessionReaderSharedPtr &pSession,
+    const std::weak_ptr<EquationSystem> &pEquation)
+    : Forcing(pSession, pEquation), ForcingSyntheticEddy(pSession, pEquation)
+{
+}
+
+/**
+ * @brief Apply forcing term if an eddy left the box of eddies and
+ *        update the eddies positions.
+ *
+ * @param pFields   Pointer to fields.
+ * @param inarray   Given fields. The fields are in in physical space.
+ * @param outarray  Calculated solution after forcing term being applied
+ *                  in physical space.
+ * @param time      time.
+ */
+void ForcingCFSSyntheticEddy::v_Apply(
+    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+    [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
+    Array<OneD, Array<OneD, NekDouble>> &outarray,
+    [[maybe_unused]] const NekDouble &time)
+{
+    // Number of Variables
+    int nVars = pFields.size();
+    // Total number of coefficients
+    unsigned int nqTot = pFields[0]->GetTotPoints();
+
+    // Only calculates the forcing in the first time step and when an eddy
+    // leaves the synthetic eddy region (box).
+    if (m_calcForcing)
+    {
+        CalculateForcing(pFields);
+        m_calcForcing = false;
+    }
+
+    for (size_t i = 0; i < nVars; ++i)
+    {
+        Vmath::Vadd(nqTot, m_Forcing[i], 1, outarray[i], 1, outarray[i], 1);
+    }
+
+    // Update eddies position inside the box.
+    UpdateEddiesPositions();
+}
+
+/**
+ * @brief Apply forcing term if an eddy left the box of eddies and
+ *        update the eddies positions.
+ *
+ * @param pFields   Pointer to fields.
+ * @param inarray   Given fields. The fields are in in physical space.
+ * @param outarray  Calculated solution after forcing term being applied
+ *                  in coeficient space.
+ * @param time      time.
+ */
+void ForcingCFSSyntheticEddy::v_ApplyCoeff(
+    const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
+    [[maybe_unused]] const Array<OneD, Array<OneD, NekDouble>> &inarray,
+    Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble &time)
+{
+    // Number of Variables
+    int nVars = fields.size();
+    // Total number of coefficients
+    unsigned int nCoeff = fields[0]->GetNcoeffs();
+
+    // Only calculates the forcing in the first time step and when an eddy
+    // leaves the box
+    if (m_calcForcing)
+    {
+        CalculateForcing(fields);
+        m_calcForcing = false;
+    }
+
+    Array<OneD, NekDouble> forcingCoeff(nCoeff);
+    for (size_t i = 0; i < nVars; ++i)
+    {
+        fields[i]->FwdTrans(m_Forcing[i], forcingCoeff);
+        Vmath::Vadd(nCoeff, forcingCoeff, 1, outarray[i], 1, outarray[i], 1);
+    }
+
+    // Check for update: implicit solver
+    if (m_currTime != time)
+    {
+        UpdateEddiesPositions();
+        m_currTime = time;
+    }
+}
+
+/**
+ * @brief Calculate forcing term.
+ *
+ * @param pFfields  Pointer to fields.
+ */
+void ForcingCFSSyntheticEddy::CalculateForcing(
+    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
+{
+    // Total number of quadrature points
+    int nqTot = pFields[0]->GetTotPoints();
+    // Number of Variables
+    int nVars = pFields.size();
+    // Define parameters
+    Array<OneD, NekDouble> rho(nqTot), temperature(nqTot), pressure(nqTot);
+    // Velocity fluctuation module
+    NekDouble velFlucMod;
+    // Variable converter
+    m_varConv = MemoryManager<VariableConverter>::AllocateSharedPtr(m_session,
+                                                                    m_spacedim);
+
+    // Compute Stochastic Signal
+    Array<OneD, Array<OneD, NekDouble>> stochasticSignal;
+    stochasticSignal = ComputeStochasticSignal(pFields);
+
+    // Compute rho and mach mean inside the box of eddies
+    std::pair<NekDouble, NekDouble> rhoMachMean;
+    rhoMachMean = ComputeRhoMachMean(pFields);
+
+    // Compute velocity fluctuation
+    Array<OneD, Array<OneD, NekDouble>> velFluc;
+    velFluc = ComputeVelocityFluctuation(pFields, stochasticSignal);
+
+    // Compute density fluctuation
+    Array<OneD, Array<OneD, NekDouble>> rhoFluc;
+    rhoFluc = ComputeDensityFluctuation(pFields, velFluc, rhoMachMean);
+
+    // Compute characteristic convective turbulent time
+    Array<OneD, Array<OneD, NekDouble>> convTurbTime;
+    convTurbTime = ComputeCharConvTurbTime(pFields);
+
+    // Compute smoothing factor
+    Array<OneD, Array<OneD, NekDouble>> smoothFac;
+    smoothFac = ComputeSmoothingFactor(pFields, convTurbTime);
+
+    // Get physical data
+    Array<OneD, Array<OneD, NekDouble>> physFields(nVars);
+    for (size_t i = 0; i < nVars; ++i)
+    {
+        physFields[i] = pFields[i]->GetPhys();
+    }
+
+    // Calculate parameters
+    m_varConv->GetTemperature(physFields, temperature);
+    m_varConv->GetPressure(physFields, pressure);
+    m_varConv->GetRhoFromPT(pressure, temperature, rho);
+
+    // Check if eddies left the box (reinjected). Note that the member
+    // m_eddiesIDForcing is populate with the eddies that left the box.
+    // If any left it is going to be empty.
+    if (!m_eddiesIDForcing.empty())
+    {
+        // Forcing must stop applying for eddies that left the box
+        RemoveEddiesFromForcing(pFields);
+
+        // Update Forcing term which are going to be applied until
+        // the eddy leave the leave the domain
+        // Select the eddies to apply the forcing. Superposition.
+        for (auto &n : m_eddiesIDForcing)
+        {
+            for (size_t i = 0; i < nqTot; ++i)
+            {
+                if (m_mask[i])
+                {
+                    // density term
+                    m_ForcingEddy[n][0][i] =
+                        (rhoFluc[n][i] * smoothFac[0][i]) / convTurbTime[0][i];
+                    // Update forcing
+                    m_Forcing[0][i] += m_ForcingEddy[n][0][i];
+                    //  velocity term
+                    for (size_t j = 0; j < m_spacedim; ++j)
+                    {
+                        m_ForcingEddy[n][j + 1][i] =
+                            (rho[i] * velFluc[n][j * nqTot + i] *
+                             smoothFac[j][i]) /
+                            convTurbTime[j][i];
+                        // Update forcing
+                        m_Forcing[j + 1][i] += m_ForcingEddy[n][j + 1][i];
+                    }
+                    // energy term
+                    velFlucMod =
+                        velFluc[n][i] * velFluc[n][i] +
+                        velFluc[n][nqTot + i] * velFluc[n][nqTot + i] +
+                        velFluc[n][2 * nqTot + i] * velFluc[n][2 * nqTot + i];
+
+                    m_ForcingEddy[n][nVars - 1][i] =
+                        (0.5 * rho[i] * velFlucMod * smoothFac[0][i]) /
+                        convTurbTime[0][i];
+                    // Update forcing
+                    m_Forcing[nVars - 1][i] += m_ForcingEddy[n][nVars - 1][i];
+                }
+            }
+        }
+        // delete eddies
+        m_eddiesIDForcing.erase(m_eddiesIDForcing.begin(),
+                                m_eddiesIDForcing.end());
+    }
+    else
+    {
+        NEKERROR(ErrorUtil::efatal, "Failed: Eddies ID vector is empty.");
+    }
+}
+
+/**
+ * @brief Compute density fluctuation for the source term
+ *
+ * @param pFfields  Pointer to fields.
+ */
+Array<OneD, Array<OneD, NekDouble>> ForcingCFSSyntheticEddy::
+    ComputeDensityFluctuation(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const Array<OneD, Array<OneD, NekDouble>> &velFluc,
+        std::pair<NekDouble, NekDouble> rhoMachMean)
+{
+    // Total number of quadrature points
+    int nqTot = pFields[0]->GetTotPoints();
+    // Density fluctuation
+    Array<OneD, Array<OneD, NekDouble>> rhoFluc(m_N);
+
+    for (auto &n : m_eddiesIDForcing)
+    {
+        rhoFluc[n] = Array<OneD, NekDouble>(nqTot, 0.0);
+
+        for (size_t i = 0; i < nqTot; ++i)
+        {
+            if (m_mask[i])
+            {
+                // only need velocity fluctation in the x-direction
+                rhoFluc[n][i] = rhoMachMean.first * (m_gamma - 1) *
+                                rhoMachMean.second * rhoMachMean.second *
+                                (velFluc[n][i] / m_Ub);
+            }
+        }
+    }
+
+    return rhoFluc;
+}
+
+/**
+ * @brief Calculate the density and mach number mean
+ *        inside the box of eddies.
+ *
+ * @param  pFfields  Pointer to fields.
+ * @return           Pair with density and mach nember means.
+ */
+std::pair<NekDouble, NekDouble> ForcingCFSSyntheticEddy::ComputeRhoMachMean(
+    const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields)
+{
+    // Total number of quadrature points
+    int nqTot = pFields[0]->GetTotPoints();
+    // Number of Variables
+    int nVars = pFields.size();
+    // <rho> and <mach>^{2}
+    NekDouble rhoMean = 0.0, machMean = 0.0;
+    // Counter for the mean
+    int count = 0;
+
+    // Get physical field
+    Array<OneD, Array<OneD, NekDouble>> physFields(nVars);
+    for (size_t i = 0; i < nVars; ++i)
+    {
+        physFields[i] = pFields[i]->GetPhys();
+    }
+
+    // Define parameters
+    Array<OneD, NekDouble> soundSpeed(nqTot), mach(nqTot), rho(nqTot),
+        temperature(nqTot), pressure(nqTot);
+    // Calculate parameters
+
+    m_varConv->GetTemperature(physFields, temperature);
+    m_varConv->GetPressure(physFields, pressure);
+    m_varConv->GetRhoFromPT(pressure, temperature, rho);
+    m_varConv->GetSoundSpeed(physFields, soundSpeed);
+    m_varConv->GetMach(physFields, soundSpeed, mach);
+
+    // Sum
+    for (size_t i = 0; i < nqTot; ++i)
+    {
+        if (m_mask[i])
+        {
+            rhoMean += rho[i];
+            machMean += mach[i];
+            count += 1;
+        }
+    }
+
+    // Density mean
+    rhoMean = rhoMean / count;
+    // Mach number mean
+    machMean = machMean / count;
+
+    return std::make_pair(rhoMean, machMean);
+}
+
+} // namespace Nektar::SolverUtils

solvers/CompressibleFlowSolver/Forcing/ForcingCFSSyntheticEddy.h

+//////////////////////////////////////////////////////////////////////////////
+//
+// File:  ForcingCFSSyntheticEddy.h
+//
+// 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: Derived base class - Synthetic turbulence forcing for the
+//              Compressible solver.
+//
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef NEKTAR_SOLVERUTILS_FORCINGCFSSYNTHETICEDDY
+#define NEKTAR_SOLVERUTILS_FORCINGCFSSYNTHETICEDDY
+
+#include <CompressibleFlowSolver/Misc/VariableConverter.h>
+#include <SolverUtils/Forcing/Forcing.h>
+#include <SolverUtils/Forcing/ForcingSyntheticEddy.h>
+#include <string>
+
+namespace Nektar::SolverUtils
+{
+
+class ForcingCFSSyntheticEddy : virtual public SolverUtils::Forcing,
+                                virtual public SolverUtils::ForcingSyntheticEddy
+{
+public:
+    friend class MemoryManager<ForcingCFSSyntheticEddy>;
+
+    /// Creates an instance of this class
+    static ForcingSharedPtr create(
+        const LibUtilities::SessionReaderSharedPtr &pSession,
+        const std::weak_ptr<EquationSystem> &pEquation,
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const unsigned int &pNumForcingFields, const TiXmlElement *pForce)
+    {
+        ForcingSharedPtr p =
+            MemoryManager<ForcingCFSSyntheticEddy>::AllocateSharedPtr(
+                pSession, pEquation);
+        p->InitObject(pFields, pNumForcingFields, pForce);
+        return p;
+    }
+
+    /// Name of the class
+    static std::string className;
+
+protected:
+    // Apply forcing term
+    void v_Apply(const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+                 const Array<OneD, Array<OneD, NekDouble>> &inarray,
+                 Array<OneD, Array<OneD, NekDouble>> &outarray,
+                 const NekDouble &time) override;
+    // Apply forcing term
+    void v_ApplyCoeff(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const Array<OneD, Array<OneD, NekDouble>> &inarray,
+        Array<OneD, Array<OneD, NekDouble>> &outarray,
+        const NekDouble &time) override;
+    /// Calculate Forcing
+    void CalculateForcing(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+    /// Compute Velocity Fluctuation
+    Array<OneD, Array<OneD, NekDouble>> ComputeDensityFluctuation(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
+        const Array<OneD, Array<OneD, NekDouble>> &velFLuc,
+        std::pair<NekDouble, NekDouble> rhoMachMean);
+    /// Compute rho and mach mean
+    std::pair<NekDouble, NekDouble> ComputeRhoMachMean(
+        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
+
+    /// Auxiliary object to convert variables
+    VariableConverterSharedPtr m_varConv;
+
+private:
+    ForcingCFSSyntheticEddy(
+        const LibUtilities::SessionReaderSharedPtr &pSession,
+        const std::weak_ptr<EquationSystem> &pEquation);
+    ~ForcingCFSSyntheticEddy(void) override = default;
+};
+
+} // namespace Nektar::SolverUtils
+
+#endif

solvers/CompressibleFlowSolver/Tests/ChanFlow3D_infTurbExpl.tst

+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>Fully 3D Compressible Channel Flow. Synthetic turbulence generation is introduced in the flow field using forcing tag. This test case uses the explicit solver. </description>
+    <executable>CompressibleFlowSolver</executable>
+    <parameters>ChanFlow3D_infTurbExpl.xml</parameters>
+    <files>
+        <file description="Session File">ChanFlow3D_infTurbExpl.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="rho"  tolerance="1e-7">0.0496262</value>
+            <value variable="rhou" tolerance="1e-7">0.0510315</value>
+            <value variable="rhov" tolerance="1e-8">0.00044049</value>
+            <value variable="rhow" tolerance="1e-9">0.000513142</value>
+            <value variable="E"    tolerance="1e-5">8.88794</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="rho"  tolerance="1e-7">0.0142084</value>
+            <value variable="rhou" tolerance="1e-7">0.0220761</value>
+            <value variable="rhov" tolerance="1e-8">0.00487243</value>
+            <value variable="rhow" tolerance="1e-8">0.00275159</value>
+            <value variable="E"    tolerance="1e-5">2.55629</value>
+        </metric>
+    </metrics>
+</test>
+
+

solvers/CompressibleFlowSolver/Tests/ChanFlow3D_infTurbImpl.tst

+<?xml version="1.0" encoding="utf-8"?>
+<test>
+    <description>Fully 3D Compressible Channel Flow. Synthetic turbulence generation is introduced in the flow field using forcing tag. This test cases uses the implicit solver. </description>
+    <executable>CompressibleFlowSolver</executable>
+    <parameters>ChanFlow3D_infTurbImpl.xml</parameters>
+    <files>
+        <file description="Session File">ChanFlow3D_infTurbImpl.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="rho"  tolerance="1e-7">0.0496262</value>
+            <value variable="rhou" tolerance="1e-7">0.0510306</value>
+            <value variable="rhov" tolerance="1e-8">0.000247923</value>
+            <value variable="rhow" tolerance="1e-8">0.000463459</value>
+            <value variable="E"    tolerance="1e-5">8.88793</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="rho"  tolerance="1e-7">0.0140256</value>
+            <value variable="rhou" tolerance="1e-7">0.0188351</value>
+            <value variable="rhov" tolerance="1e-8">0.00117766</value>
+            <value variable="rhow" tolerance="1e-8">0.00196196</value>
+            <value variable="E"    tolerance="1e-5"> 2.51336</value>
+        </metric>
+    </metrics>
+</test>
+
+

solvers/IncNavierStokesSolver/CMakeLists.txt

@@ -34,7 +34,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
         ./Filters/FilterReynoldsStresses.cpp
         ./Filters/FilterMovingBody.cpp
         ./Filters/FilterAeroForcesSPM.cpp
-	./Forcing/ForcingIncNSSyntheticEddy.cpp
+        ./Forcing/ForcingIncNSSyntheticEddy.cpp
         ./Forcing/ForcingMovingBody.cpp
         ./Forcing/ForcingStabilityCoupledLNS.cpp
         ./BoundaryConditions/IncBoundaryConditions.cpp
@@ -71,6 +71,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
             Filters/FilterAeroForcesSPM.h
             Filters/FilterMovingBody.h
             Filters/FilterReynoldsStresses.h
+            Forcing/ForcingIncNSSyntheticEddy.h
             Forcing/ForcingMovingBody.h
             Forcing/ForcingStabilityCoupledLNS.h
             BoundaryConditions/IncBoundaryConditions.h

solvers/IncNavierStokesSolver/Forcing/ForcingIncNSSyntheticEddy.h

@@ -28,7 +28,8 @@
 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 // DEALINGS IN THE SOFTWARE.
 //
-// Description: Derived base class - Synthetic turbulence forcing
+// Description: Derived base class - Synthetic turbulence forcing for the
+//              Incompressible solver
 //
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -36,12 +37,14 @@
 #define NEKTAR_SOLVERUTILS_FORCINGINCNSSYNTHETICEDDY
 
 #include <SolverUtils/Forcing/Forcing.h>
+#include <SolverUtils/Forcing/ForcingSyntheticEddy.h>
 #include <string>
 
 namespace Nektar::SolverUtils
 {
-
-class ForcingIncNSSyntheticEddy : public SolverUtils::Forcing
+class ForcingIncNSSyntheticEddy
+    : virtual public SolverUtils::Forcing,
+      virtual public SolverUtils::ForcingSyntheticEddy
 {
 public:
     friend class MemoryManager<ForcingIncNSSyntheticEddy>;
@@ -64,112 +67,14 @@ public:
     static std::string className;
 
 protected:
-    void v_InitObject(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
-        const unsigned int &pNumForcingFields,
-        const TiXmlElement *pForce) override;
-
+    // Apply forcing term
     void v_Apply(const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
                  const Array<OneD, Array<OneD, NekDouble>> &inarray,
                  Array<OneD, Array<OneD, NekDouble>> &outarray,
                  const NekDouble &time) override;
-
-    // Set Cholesky decomposition of the Reynolds Stresses in the domain
-    void SetCholeskyReyStresses(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-    /// Set the Number of the Eddies in the box
-    void SetNumberOfEddies();
-    /// Set reference lengths
-    void ComputeRefLenghts();
-    /// Set Xi max
-    void ComputeXiMax();
-    /// Compute Random 1 or -1 value
-    Array<OneD, Array<OneD, int>> GenerateRandomOneOrMinusOne();
-    /// Compute Constant C
-    NekDouble ComputeConstantC(int row, int col);
-    /// Compute Gaussian
-    NekDouble ComputeGaussian(NekDouble coord, NekDouble xiMaxVal,
-                              NekDouble constC = 1.0);
-    /// Check if point is inside the box of eddies
-    bool InsideBoxOfEddies(NekDouble coord0, NekDouble coord1,
-                           NekDouble coord2);
-    /// Compute Stochastic Signal
-    Array<OneD, Array<OneD, NekDouble>> ComputeStochasticSignal(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-    /// Compute Velocity Fluctuation
-    Array<OneD, Array<OneD, NekDouble>> ComputeVelocityFluctuation(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
-        Array<OneD, Array<OneD, NekDouble>> stochasticSignal);
-    /// Compute Characteristic Convective Turbulent Time
-    Array<OneD, Array<OneD, NekDouble>> ComputeCharConvTurbTime(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-    /// Compute Smoothing Factor
-    Array<OneD, Array<OneD, NekDouble>> ComputeSmoothingFactor(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
-        Array<OneD, Array<OneD, NekDouble>> convTurb);
-    /// Set box of eddies mask
-    void SetBoxOfEddiesMask(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-    /// Compute the initial position of the eddies inside the box
-    void ComputeInitialRandomLocationOfEddies();
-    /// Update positions of the eddies inside the box
-    void UpdateEddiesPositions();
     /// Calculate Forcing
     void CalculateForcing(
         const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-    /// Compute the initial location of the eddies for the test case
-    void ComputeInitialLocationTestCase();
-    /// Remove eddy from forcing term
-    void RemoveEddiesFromForcing(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-    /// Initialise forcing term for each eddy
-    void InitialiseForcingEddy(
-        const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields);
-
-    // Members
-    // Expressions (functions) of the prescribed Reynolds stresses
-    std::map<int, LibUtilities::EquationSharedPtr> m_R;
-    /// Cholesky decomposition of the Reynolds Stresses
-    /// for all points in the mesh
-    Array<OneD, Array<OneD, NekDouble>> m_Cholesky;
-    /// Bulk velocity
-    NekDouble m_Ub;
-    /// Standard deviation
-    NekDouble m_sigma;
-    /// Inlet length in the y-direction and z-direction
-    Array<OneD, NekDouble> m_lyz;
-    /// Center of the inlet plane
-    Array<OneD, NekDouble> m_rc;
-    /// Number of eddies in the box
-    int m_N;
-    /// Characteristic lenght scales
-    Array<OneD, NekDouble> m_l;
-    /// Reference lenghts
-    Array<OneD, NekDouble> m_lref;
-    /// Xi max
-    Array<OneD, NekDouble> m_xiMax;
-    // XiMaxMin - Value form Giangaspero et al. 2022
-    NekDouble m_xiMaxMin = 0.4;
-    /// Space dimension
-    int m_spacedim;
-    /// Ration of specific heats
-    NekDouble m_gamma;
-    /// Box of eddies mask
-    Array<OneD, int> m_mask;
-    /// Eddy position
-    Array<OneD, Array<OneD, NekDouble>> m_eddyPos;
-    /// Check when the forcing should be applied
-    bool m_calcForcing{true};
-    /// Eddies that add to the domain
-    std::vector<unsigned int> m_eddiesIDForcing;
-    /// Current time
-    NekDouble m_currTime;
-    /// Keep applying force during GMRES iteration
-    bool m_implicitForcing{false};
-    /// Check for test case
-    bool m_tCase;
-    /// Forcing for each eddy
-    Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_ForcingEddy;
 
 private:
     ForcingIncNSSyntheticEddy(

solvers/IncNavierStokesSolver/Tests/ChanFlow3D_infTurb.tst

@@ -8,16 +8,16 @@
     </files>
     <metrics>
         <metric type="L2" id="1">
-            <value variable="u" tolerance="1e-5">3.64539</value>
-            <value variable="v" tolerance="1e-7">0.0268234</value>
-            <value variable="w" tolerance="1e-6">0.030717</value>
-            <value variable="p" tolerance="1e-7">0.0224162</value>
+            <value variable="u" tolerance="1e-5">3.64585</value>
+            <value variable="v" tolerance="1e-7">0.0316134</value>
+            <value variable="w" tolerance="1e-7">0.0361575</value>
+            <value variable="p" tolerance="1e-7">0.0214355</value>
         </metric>
         <metric type="Linf" id="2">
-            <value variable="u" tolerance="1e-5">1.40919</value>
-            <value variable="v" tolerance="1e-6">0.145714</value>
-            <value variable="w" tolerance="1e-6">0.119218</value>
-            <value variable="p" tolerance="1e-7">0.0906175</value>
+            <value variable="u" tolerance="1e-5">1.46953</value>
+            <value variable="v" tolerance="1e-6">0.183361</value>
+            <value variable="w" tolerance="1e-6">0.160353</value>
+            <value variable="p" tolerance="1e-7">0.0869825</value>
         </metric>
     </metrics>
 </test>

solvers/IncNavierStokesSolver/Tests/ChanFlow3D_infTurb.xml

@@ -41,7 +41,7 @@
 
         <PARAMETERS>
             <P> TimeStep      = 0.01            </P>
-            <P> NumSteps      = 75              </P>
+            <P> NumSteps      = 85              </P>
             <P> IO_CheckSteps = 0               </P>
             <P> IO_InfoSteps  = 1               </P>
             <P> IO_CFLSteps   = 1               </P>
@@ -118,7 +118,7 @@
             <E VAR="r00"   VALUE="(+2194869.72*(1-abs(y))^14-15843749.2*(1-abs(y))^13+51288368.6*(1-abs(y))^12-98280270.7*(1-abs(y))^11+123921817.0*(1-abs(y))^10-108083933.0*(1-abs(y))^9+66700892.4*(1-abs(y))^8-29258748.2*(1-abs(y))^7+9017078.56*(1-abs(y))^6-1889548.13*(1-abs(y))^5+250183.446*(1-abs(y))^4-17029.8457*(1-abs(y))^3+1.26582088*(1-abs(y))^2+68.7492352*(1-abs(y))^1-0.0148027907*(1-abs(y))^0)^2*(tauw/rhom)"     />
             <E VAR="r11"   VALUE="(+115374.19*(1-abs(y))^14-835012.705*(1-abs(y))^13+2718672.87*(1-abs(y))^12-5262526.56*(1-abs(y))^11+6743594.93*(1-abs(y))^10-6028822.18*(1-abs(y))^9+3861137.61*(1-abs(y))^8-1790967.72*(1-abs(y))^7+601467.524*(1-abs(y))^6-144772.848*(1-abs(y))^5+24433.5085*(1-abs(y))^4-2749.69383*(1-abs(y))^3+170.501429*(1-abs(y))^2+1.19079597*(1-abs(y))^1-0.00240344727*(1-abs(y))^0)^2*(tauw/rhom)"     />
             <E VAR="r22"   VALUE="(-467724.204*(1-abs(y))^14+3358371.37*(1-abs(y))^13-10821834.7*(1-abs(y))^12+20672126.8*(1-abs(y))^11-26050692.2*(1-abs(y))^10+22809003.9*(1-abs(y))^9-14237511.6*(1-abs(y))^8+6400502.89*(1-abs(y))^7-2069981.91*(1-abs(y))^6+476610.822*(1-abs(y))^5-76731.1913*(1-abs(y))^4+8464.95414*(1-abs(y))^3-638.190408*(1-abs(y))^2+33.7641881*(1-abs(y))^1+0.00362081709*(1-abs(y))^0)^2*(tauw/rhom)"     />
-            <E VAR="r10"   VALUE="(y>0)*(-y*(+143887.186*(1-abs(y))^14-924421.506*(1-abs(y))^13+2558324.72*(1-abs(y))^12-3912224.86*(1-abs(y))^11+3423219.44*(1-abs(y))^10-1362528.94*(1-abs(y))^9-423777.124*(1-abs(y))^8+908513.536*(1-abs(y))^7-583322.471*(1-abs(y))^6+215251.323*(1-abs(y))^5-49187.1036*(1-abs(y))^4+6726.52324*(1-abs(y))^3-466.587078*(1-abs(y))^2+5.87033236*(1-abs(y))^1-0.0123273856*(1-abs(y))^0))*(tauw/rhom)  +  (y<0)*(y*(+143887.186*(1-abs(y))^14-924421.506*(1-abs(y))^13+2558324.72*(1-abs(y))^12-3912224.86*(1-abs(y))^11+3423219.44*(1-abs(y))^10-1362528.94*(1-abs(y))^9-423777.124*(1-abs(y))^8+908513.536*(1-abs(y))^7-583322.471*(1-abs(y))^6+215251.323*(1-abs(y))^5-49187.1036*(1-abs(y))^4+6726.52324*(1-abs(y))^3-466.587078*(1-abs(y))^2+5.87033236*(1-abs(y))^1-0.0123273856*(1-abs(y))^0))*(tauw/rhom)"     />
+            <E VAR="r10"   VALUE="(y>0)*(-1.0*(+143887.186*(1-abs(y))^14-924421.506*(1-abs(y))^13+2558324.72*(1-abs(y))^12-3912224.86*(1-abs(y))^11+3423219.44*(1-abs(y))^10-1362528.94*(1-abs(y))^9-423777.124*(1-abs(y))^8+908513.536*(1-abs(y))^7-583322.471*(1-abs(y))^6+215251.323*(1-abs(y))^5-49187.1036*(1-abs(y))^4+6726.52324*(1-abs(y))^3-466.587078*(1-abs(y))^2+5.87033236*(1-abs(y))^1-0.0123273856*(1-abs(y))^0))*(tauw/rhom)  +  (y<0)*(1.0*(+143887.186*(1-abs(y))^14-924421.506*(1-abs(y))^13+2558324.72*(1-abs(y))^12-3912224.86*(1-abs(y))^11+3423219.44*(1-abs(y))^10-1362528.94*(1-abs(y))^9-423777.124*(1-abs(y))^8+908513.536*(1-abs(y))^7-583322.471*(1-abs(y))^6+215251.323*(1-abs(y))^5-49187.1036*(1-abs(y))^4+6726.52324*(1-abs(y))^3-466.587078*(1-abs(y))^2+5.87033236*(1-abs(y))^1-0.0123273856*(1-abs(y))^0))*(tauw/rhom)"     />
             <E VAR="r20"   VALUE="0.0"     />
             <E VAR="r21"   VALUE="0.0"     />
         </FUNCTION>