From 34f6d4d62bb24b98de7d0379a576f85022f412a2 Mon Sep 17 00:00:00 2001
From: Ankang Gao <ankanggao@ustc.edu.cn>
Date: Fri, 16 Aug 2024 13:54:54 +0000
Subject: [PATCH] Fix initial conditions in MRF when set-start-time is used
---
CHANGELOG.md | 2 +-
.../Forcing/ForcingMovingReferenceFrame.cpp | 152 ++++++++++--------
.../Forcing/ForcingMovingReferenceFrame.h | 23 ++-
.../BoundaryConditions/IncBaseCondition.cpp | 21 ++-
.../BoundaryConditions/IncBaseCondition.h | 2 +
.../BoundaryConditions/MovingFrameWall.cpp | 14 +-
.../BoundaryConditions/StaticWall.cpp | 8 +-
.../BoundaryConditions/StaticWall.h | 2 -
.../BoundaryConditions/TransMovingWall.cpp | 18 +--
.../Tests/FreeFallCyl.tst | 2 +-
10 files changed, 134 insertions(+), 110 deletions(-)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 13a78b4fde..f27ca90d1d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,7 +10,7 @@ v5.7.0
- Fix AdaptiveSFD for MPI (!1821)
**IncNavierStokesSolver**
-- Fix initial and boundary conditions in the moving reference frame (!1692)
+- Fix initial and boundary conditions in the moving reference frame (!1692, !1820)
- Fix memory-leak for the Mixed_CG_Discontinuous projection when initializing the traceMep (!1806)
- Add synthetic turbulence generation for the incompressible solver (!1664)
diff --git a/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.cpp b/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.cpp
index 42aff113d6..c9d5acc44e 100644
--- a/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.cpp
+++ b/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.cpp
@@ -107,18 +107,18 @@ void ForcingMovingReferenceFrame::v_InitObject(
m_hasPlane0 = pFields[0]->GetZIDs()[0] == 0;
}
// initialize variables
- m_velxyz = Array<OneD, NekDouble>(3, 0.0);
- m_omegaxyz = Array<OneD, NekDouble>(3, 0.0);
- m_extForceXYZ = Array<OneD, NekDouble>(6, 0.0);
- m_hasVel = Array<OneD, bool>(3, false);
- m_hasOmega = Array<OneD, bool>(3, false);
- m_travelWave = Array<OneD, NekDouble>(3, 0.);
- m_pivotPoint = Array<OneD, NekDouble>(3, 0.0);
- m_index = 0;
- m_currentTime = -1.;
+ m_velxyz = Array<OneD, NekDouble>(3, 0.0);
+ m_omegaxyz = Array<OneD, NekDouble>(3, 0.0);
+ m_body.extForceXYZ = Array<OneD, NekDouble>(6, 0.0);
+ m_hasVel = Array<OneD, bool>(3, false);
+ m_hasOmega = Array<OneD, bool>(3, false);
+ m_travelWave = Array<OneD, NekDouble>(3, 0.);
+ m_pivotPoint = Array<OneD, NekDouble>(3, 0.0);
+ m_index = 0;
+ m_currentTime = -1.;
LoadParameters(pForce);
InitBodySolver(pForce);
- if (m_circularCylinder || m_expdim == 1)
+ if (m_body.isCircular || m_expdim == 1)
{
for (int i = 0; i < 3; ++i)
{
@@ -162,7 +162,7 @@ void ForcingMovingReferenceFrame::v_InitObject(
auto FluidEq = std::dynamic_pointer_cast<FluidInterface>(equ);
FluidEq->SetMovingFramePivot(m_pivotPoint);
// initialize aeroforce filter
- if (m_hasFreeMotion)
+ if (m_body.hasFreeMotion)
{
InitialiseFilter(m_session, pFields, pForce);
}
@@ -193,11 +193,11 @@ void ForcingMovingReferenceFrame::LoadParameters(const TiXmlElement *pForce)
pForce->FirstChildElement("CIRCULARCYLINDER");
if (mssgTagSpecial)
{
- m_circularCylinder = true;
+ m_body.isCircular = true;
}
else
{
- m_circularCylinder = false;
+ m_body.isCircular = false;
}
// load frame velocity
funcNameElmt = pForce->FirstChildElement("FRAMEVELOCITY");
@@ -284,7 +284,7 @@ void ForcingMovingReferenceFrame::LoadParameters(const TiXmlElement *pForce)
std::string var = forceVar[i];
if (m_session->DefinesFunction(FuncName, var))
{
- m_extForceFunction[i] =
+ m_body.extForceFunction[i] =
m_session->GetFunction(FuncName, var);
}
}
@@ -294,7 +294,7 @@ void ForcingMovingReferenceFrame::LoadParameters(const TiXmlElement *pForce)
std::string var = momentVar[i];
if (m_session->DefinesFunction(FuncName, var))
{
- m_extForceFunction[i + 3] =
+ m_body.extForceFunction[i + 3] =
m_session->GetFunction(FuncName, var);
}
}
@@ -385,18 +385,18 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
const TiXmlElement *mssgTag;
std::string mssgStr;
// allocate memory and initialise
- m_bodyVel = Array<OneD, Array<OneD, NekDouble>>(3);
+ m_body.vel = Array<OneD, Array<OneD, NekDouble>>(3);
for (size_t i = 0; i < 3; ++i)
{
- m_bodyVel[i] = Array<OneD, NekDouble>(NumDof, 0.);
+ m_body.vel[i] = Array<OneD, NekDouble>(NumDof, 0.);
}
SetInitialConditions();
- UpdateFluidInterface(m_bodyVel, 1);
+ UpdateFluidInterface(m_body.vel, 1);
// read free motion DoFs
- m_DirDoFs.clear();
+ m_body.dirDoFs.clear();
for (int i = 0; i < NumDof; ++i)
{
- m_DirDoFs.insert(i);
+ m_body.dirDoFs.insert(i);
}
mssgTag = pForce->FirstChildElement("MOTIONPRESCRIBED");
if (mssgTag)
@@ -411,7 +411,7 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
{
if (EvaluateExpression(values[i]) == 0)
{
- m_DirDoFs.erase(i);
+ m_body.dirDoFs.erase(i);
if (i < m_spacedim)
{
m_hasVel[i] = true;
@@ -424,11 +424,12 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
}
}
}
- m_hasFreeMotion = m_DirDoFs.size() < NumDof;
+ m_body.hasFreeMotion = m_body.dirDoFs.size() < NumDof;
// read mass matrix
- Array<OneD, NekDouble> M(NumDof * NumDof, 0.);
- mssgTag = pForce->FirstChildElement("MASS");
- ASSERTL0(m_DirDoFs.size() == NumDof || mssgTag, "Mass matrix is required.");
+ m_body.M = Array<OneD, NekDouble>(NumDof * NumDof, 0.);
+ mssgTag = pForce->FirstChildElement("MASS");
+ ASSERTL0(m_body.dirDoFs.size() == NumDof || mssgTag,
+ "Mass matrix is required.");
if (mssgTag)
{
std::vector<std::string> values;
@@ -442,14 +443,14 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
{
for (int j = 0; j < NumDof; ++j)
{
- M[count] = EvaluateExpression(values[count]);
+ m_body.M[count] = EvaluateExpression(values[count]);
++count;
}
}
}
// read damping matrix
- Array<OneD, NekDouble> C(NumDof * NumDof, 0.);
- mssgTag = pForce->FirstChildElement("DAMPING");
+ m_body.C = Array<OneD, NekDouble>(NumDof * NumDof, 0.);
+ mssgTag = pForce->FirstChildElement("DAMPING");
if (mssgTag)
{
std::vector<std::string> values;
@@ -463,14 +464,14 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
{
for (int j = 0; j < NumDof; ++j)
{
- C[count] = EvaluateExpression(values[count]);
+ m_body.C[count] = EvaluateExpression(values[count]);
++count;
}
}
}
// read rigidity matrix
- Array<OneD, NekDouble> K(NumDof * NumDof, 0.);
- mssgTag = pForce->FirstChildElement("RIGIDITY");
+ m_body.K = Array<OneD, NekDouble>(NumDof * NumDof, 0.);
+ mssgTag = pForce->FirstChildElement("RIGIDITY");
if (mssgTag)
{
std::vector<std::string> values;
@@ -484,7 +485,7 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
{
for (int j = 0; j < NumDof; ++j)
{
- K[count] = EvaluateExpression(values[count]);
+ m_body.K[count] = EvaluateExpression(values[count]);
++count;
}
}
@@ -501,7 +502,8 @@ void ForcingMovingReferenceFrame::InitBodySolver(const TiXmlElement *pForce)
{
gamma = m_session->GetParameter("NewmarkGamma");
}
- m_bodySolver.SetNewmarkBeta(beta, gamma, m_timestep, M, C, K, m_DirDoFs);
+ m_bodySolver.SetNewmarkBeta(beta, gamma, m_timestep, m_body.M, m_body.C,
+ m_body.K, m_body.dirDoFs);
}
void ForcingMovingReferenceFrame::UpdatePrescribed(
@@ -537,7 +539,7 @@ void ForcingMovingReferenceFrame::UpdatePrescribed(
it.second->Evaluate(0., 0., 0., time);
}
}
- for (auto i : m_DirDoFs)
+ for (auto i : m_body.dirDoFs)
{
if (Dirs.find(i) == Dirs.end())
{
@@ -571,46 +573,49 @@ void ForcingMovingReferenceFrame::UpdateFrameVelocity(
{
// compute the velocites whoes functions are provided in inertial frame
Array<OneD, NekDouble> forcebody(6, 0.); // fluid force
- if (m_hasFreeMotion)
+ if (m_body.hasFreeMotion)
{
- for (auto it : m_extForceFunction)
+ for (auto it : m_body.extForceFunction)
{
- m_extForceXYZ[it.first] = it.second->Evaluate(0., 0., 0., time);
+ m_body.extForceXYZ[it.first] =
+ it.second->Evaluate(0., 0., 0., time);
}
- m_aeroforceFilter->GetForces(pFields, NullNekDouble1DArray, time);
+ m_body.aeroforceFilter->GetForces(pFields, NullNekDouble1DArray,
+ time);
auto equ = m_equ.lock();
ASSERTL0(equ, "Weak pointer to the equation system is expired");
auto FluidEq = std::dynamic_pointer_cast<FluidInterface>(equ);
FluidEq->GetAeroForce(forcebody);
}
- SolveBodyMotion(m_bodyVel, forcebody, Dirs);
+ SolveBodyMotion(m_body.vel, forcebody, Dirs);
}
if (m_isRoot && m_index % m_outputFrequency == 0)
{
m_outputStream << boost::format("%25.19e") % time << " ";
- for (size_t i = 0; i < m_bodyVel[0].size(); ++i)
+ for (size_t i = 0; i < m_body.vel[0].size(); ++i)
{
- m_outputStream << boost::format("%25.19e") % m_bodyVel[0][i] << " "
- << boost::format("%25.19e") % m_bodyVel[1][i] << " "
- << boost::format("%25.19e") % m_bodyVel[2][i] << " ";
+ m_outputStream << boost::format("%25.19e") % m_body.vel[0][i] << " "
+ << boost::format("%25.19e") % m_body.vel[1][i] << " "
+ << boost::format("%25.19e") % m_body.vel[2][i]
+ << " ";
}
m_outputStream << endl;
}
// extract values and transform to body frame
for (int i = 0; i < m_spacedim; ++i)
{
- m_velxyz[i] = m_bodyVel[1][i];
+ m_velxyz[i] = m_body.vel[1][i];
}
if (m_hasRotation)
{
- m_omegaxyz[2] = m_bodyVel[1][m_spacedim];
+ m_omegaxyz[2] = m_body.vel[1][m_spacedim];
Array<OneD, NekDouble> angle(3, 0.);
- angle[2] = m_bodyVel[0][m_spacedim];
+ angle[2] = m_body.vel[0][m_spacedim];
m_frame.SetAngle(angle);
m_frame.IneritalToBody(3, m_velxyz, m_velxyz);
}
// set displacements at the current time step
- UpdateFluidInterface(m_bodyVel, 0);
+ UpdateFluidInterface(m_body.vel, 0);
++m_index;
}
@@ -624,7 +629,7 @@ void ForcingMovingReferenceFrame::UpdateFluidInterface(
{
disp[i] = bodyVel[0][i];
}
- if (!m_circularCylinder)
+ if (!m_body.isCircular)
{
disp[5] = bodyVel[0][m_spacedim];
}
@@ -655,7 +660,7 @@ void ForcingMovingReferenceFrame::SolveBodyMotion(
Array<OneD, Array<OneD, NekDouble>> &bodyVel,
const Array<OneD, NekDouble> &forcebody, std::map<int, NekDouble> &Dirs)
{
- if (!m_hasFreeMotion)
+ if (!m_body.hasFreeMotion)
{
m_bodySolver.SolvePrescribed(bodyVel, Dirs);
}
@@ -677,9 +682,9 @@ void ForcingMovingReferenceFrame::SolveBodyMotion(
}
for (int i = 0; i < m_spacedim; ++i)
{
- force[i] += m_extForceXYZ[i];
+ force[i] += m_body.extForceXYZ[i];
}
- force[m_spacedim] = force[5] + m_extForceXYZ[5];
+ force[m_spacedim] = force[5] + m_body.extForceXYZ[5];
m_bodySolver.SolveFree(bodyVel, force);
}
else
@@ -705,9 +710,9 @@ void ForcingMovingReferenceFrame::SolveBodyMotion(
m_frame.BodyToInerital(3, forcebody + 3, tmp = force + 3);
for (int i = 0; i < m_spacedim; ++i)
{
- force[i] += m_extForceXYZ[i];
+ force[i] += m_body.extForceXYZ[i];
}
- force[m_spacedim] = force[5] + m_extForceXYZ[5];
+ force[m_spacedim] = force[5] + m_body.extForceXYZ[5];
m_bodySolver.Solve(tmpbodyVel, force, Dirs);
angle[2] = tmpbodyVel[0][m_spacedim];
}
@@ -750,11 +755,11 @@ void ForcingMovingReferenceFrame::UpdateBoundaryConditions(NekDouble time)
// compute the velocities whose functions are provided in inertial frame
std::map<int, NekDouble> Dirs;
UpdatePrescribed(time, Dirs);
- Array<OneD, Array<OneD, NekDouble>> bodyVel(m_bodyVel.size());
- for (size_t i = 0; i < m_bodyVel.size(); ++i)
+ Array<OneD, Array<OneD, NekDouble>> bodyVel(m_body.vel.size());
+ for (size_t i = 0; i < m_body.vel.size(); ++i)
{
- bodyVel[i] = Array<OneD, NekDouble>(m_bodyVel[i].size());
- Vmath::Vcopy(m_bodyVel[i].size(), m_bodyVel[i], 1, bodyVel[i], 1);
+ bodyVel[i] = Array<OneD, NekDouble>(m_body.vel[i].size());
+ Vmath::Vcopy(m_body.vel[i].size(), m_body.vel[i], 1, bodyVel[i], 1);
}
m_bodySolver.SolvePrescribed(bodyVel, Dirs);
// set displacements at the next time step
@@ -931,18 +936,23 @@ void ForcingMovingReferenceFrame::SetInitialConditions()
}
}
}
+ if (m_session->DefinesCmdLineArgument("set-start-time"))
+ {
+ time = std::stod(
+ m_session->GetCmdLineArgument<std::string>("set-start-time"));
+ }
if (fileData.size() == strFrameData.size())
{
- int NumDofm1 = m_bodyVel[0].size() - 1;
+ int NumDofm1 = m_body.vel[0].size() - 1;
for (int i = 0; i < m_spacedim; ++i)
{
- m_bodyVel[0][i] = fileData[strFrameData[i]];
- m_bodyVel[1][i] = fileData[strFrameData[i + 6]];
- m_bodyVel[2][i] = fileData[strFrameData[i + 12]];
+ m_body.vel[0][i] = fileData[strFrameData[i]];
+ m_body.vel[1][i] = fileData[strFrameData[i + 6]];
+ m_body.vel[2][i] = fileData[strFrameData[i + 12]];
}
- m_bodyVel[0][NumDofm1] = fileData[strFrameData[5]];
- m_bodyVel[1][NumDofm1] = fileData[strFrameData[11]];
- m_bodyVel[2][NumDofm1] = fileData[strFrameData[17]];
+ m_body.vel[0][NumDofm1] = fileData[strFrameData[5]];
+ m_body.vel[1][NumDofm1] = fileData[strFrameData[11]];
+ m_body.vel[2][NumDofm1] = fileData[strFrameData[17]];
}
std::map<int, NekDouble> Dirs;
UpdatePrescribed(time, Dirs);
@@ -954,19 +964,19 @@ void ForcingMovingReferenceFrame::SetInitialConditions(
{
for (auto it : Dirs)
{
- int NumDof = m_bodyVel[0].size();
+ int NumDof = m_body.vel[0].size();
int NumDof2 = NumDof << 1;
- if (it.first < m_bodyVel[0].size())
+ if (it.first < m_body.vel[0].size())
{
- m_bodyVel[1][it.first] = it.second;
+ m_body.vel[1][it.first] = it.second;
}
else if (it.first < NumDof2)
{
- m_bodyVel[0][it.first - NumDof] = it.second;
+ m_body.vel[0][it.first - NumDof] = it.second;
}
else
{
- m_bodyVel[2][it.first - NumDof2] = it.second;
+ m_body.vel[2][it.first - NumDof2] = it.second;
}
}
}
@@ -989,10 +999,10 @@ void ForcingMovingReferenceFrame::InitialiseFilter(
std::replace(pstr.begin(), pstr.end(), ',', ' ');
vParams["MomentPoint"] = pstr;
}
- m_aeroforceFilter = MemoryManager<FilterAeroForces>::AllocateSharedPtr(
+ m_body.aeroforceFilter = MemoryManager<FilterAeroForces>::AllocateSharedPtr(
pSession, m_equ.lock(), vParams);
- m_aeroforceFilter->Initialise(pFields, 0.0);
+ m_body.aeroforceFilter->Initialise(pFields, 0.0);
}
void Newmark_BetaSolver::SetNewmarkBeta(NekDouble beta, NekDouble gamma,
diff --git a/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.h b/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.h
index 6737b2a343..db7eda00f0 100644
--- a/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.h
+++ b/library/SolverUtils/Forcing/ForcingMovingReferenceFrame.h
@@ -153,7 +153,6 @@ private:
// prescribed functions in the session file
std::map<int, LibUtilities::EquationSharedPtr> m_frameVelFunction;
- std::map<int, LibUtilities::EquationSharedPtr> m_extForceFunction;
std::ofstream m_outputStream;
bool m_isRoot;
@@ -172,8 +171,6 @@ private:
Array<OneD, NekDouble> m_omegaxyz;
// coordinate vector
Array<OneD, Array<OneD, NekDouble>> m_coords;
- // externel force
- Array<OneD, NekDouble> m_extForceXYZ;
NekDouble m_currentTime;
NekDouble m_timestep;
@@ -181,17 +178,27 @@ private:
bool m_isH1d;
bool m_hasPlane0;
bool m_isH2d;
- bool m_hasFreeMotion;
NekInt m_spacedim;
NekInt m_expdim;
unsigned int m_index;
unsigned int m_outputFrequency;
+ struct
+ {
+ Array<OneD, Array<OneD, NekDouble>> vel;
+ bool hasFreeMotion;
+ std::set<int> dirDoFs;
+ bool isCircular;
+ // fluid force filter
+ FilterAeroForcesSharedPtr aeroforceFilter;
+ // externel force
+ std::map<int, LibUtilities::EquationSharedPtr> extForceFunction;
+ Array<OneD, NekDouble> extForceXYZ;
+ Array<OneD, NekDouble> M;
+ Array<OneD, NekDouble> C;
+ Array<OneD, NekDouble> K;
+ } m_body;
Newmark_BetaSolver m_bodySolver;
- Array<OneD, Array<OneD, NekDouble>> m_bodyVel;
- std::set<int> m_DirDoFs;
- bool m_circularCylinder;
- FilterAeroForcesSharedPtr m_aeroforceFilter;
ForcingMovingReferenceFrame(
const LibUtilities::SessionReaderSharedPtr &pSession,
diff --git a/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.cpp b/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.cpp
index a46841b134..94b479663b 100644
--- a/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.cpp
+++ b/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.cpp
@@ -73,6 +73,7 @@ void IncBaseCondition::v_Initialise(
if (pSession->DefinesParameter("ExtrapolateOrder"))
{
m_intSteps = std::round(pSession->GetParameter("ExtrapolateOrder"));
+ m_intSteps = std::min(3, std::max(0, m_intSteps));
}
else if (pSession->DefinesSolverInfo("TimeIntegrationMethod"))
{
@@ -92,7 +93,7 @@ void IncBaseCondition::v_Initialise(
}
else
{
- m_intSteps = 1;
+ m_intSteps = 0;
}
}
@@ -150,6 +151,10 @@ void IncBaseCondition::InitialiseCoords(
void IncBaseCondition::ExtrapolateArray(
const int numCalls, Array<OneD, Array<OneD, Array<OneD, NekDouble>>> &array)
{
+ if (m_intSteps == 1)
+ {
+ return;
+ }
int nint = std::min(numCalls, m_intSteps);
int nlevels = array.size();
int dim = array[0].size();
@@ -220,8 +225,8 @@ void IncBaseCondition::AddVisPressureBCs(
Array<OneD, Array<OneD, NekDouble>> &N,
std::map<std::string, NekDouble> ¶ms)
{
- if (params.find("Kinvis") == params.end() || params["Kinvis"] == 0. ||
- fields.size() == 0)
+ if (m_intSteps == 0 || params.find("Kinvis") == params.end() ||
+ params["Kinvis"] <= 0. || fields.size() == 0)
{
return;
}
@@ -247,8 +252,14 @@ void IncBaseCondition::AddVisPressureBCs(
Array<OneD, NekDouble> temp(m_npoints);
for (int i = 0; i < m_bnddim; i++)
{
- m_field->ExtractElmtToBndPhys(m_nbnd, Q[i], temp);
- Vmath::Svtvp(m_npoints, -kinvis, temp, 1, N[i], 1, N[i], 1);
+ m_field->ExtractElmtToBndPhys(m_nbnd, Q[i],
+ m_viscous[m_intSteps - 1][i]);
+ }
+ ExtrapolateArray(m_numCalls, m_viscous);
+ for (int i = 0; i < m_bnddim; i++)
+ {
+ Vmath::Svtvp(m_npoints, -kinvis, m_viscous[m_intSteps - 1][i], 1, N[i],
+ 1, N[i], 1);
}
}
diff --git a/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.h b/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.h
index 121a4582bc..9b272a1fd8 100644
--- a/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.h
+++ b/solvers/IncNavierStokesSolver/BoundaryConditions/IncBaseCondition.h
@@ -130,6 +130,8 @@ protected:
Array<OneD, Array<OneD, NekDouble>> m_coords;
MultiRegions::ExpListSharedPtr m_bndElmtExps;
MultiRegions::ExpListSharedPtr m_field;
+ Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_viscous;
+ int m_pressure;
static NekDouble StifflyStable_Betaq_Coeffs[3][3];
static NekDouble StifflyStable_Alpha_Coeffs[3][3];
diff --git a/solvers/IncNavierStokesSolver/BoundaryConditions/MovingFrameWall.cpp b/solvers/IncNavierStokesSolver/BoundaryConditions/MovingFrameWall.cpp
index 9b9acea4ce..0fba3e7ce2 100644
--- a/solvers/IncNavierStokesSolver/BoundaryConditions/MovingFrameWall.cpp
+++ b/solvers/IncNavierStokesSolver/BoundaryConditions/MovingFrameWall.cpp
@@ -114,19 +114,15 @@ void MovingFrameWall::v_Update(
if (m_hasPressure && fields.size() > 0)
{
++m_numCalls;
- for (int i = 0; i < m_bnddim; ++i)
- {
- Vmath::Fill(m_npoints, 0., m_viscous[m_intSteps - 1][i], 1);
- }
- AddVisPressureBCs(fields, m_viscous[m_intSteps - 1], params);
- ExtrapolateArray(m_numCalls, m_viscous);
- // copy extrapolated data
+
Array<OneD, Array<OneD, NekDouble>> rhs(m_bnddim);
for (int i = 0; i < m_bnddim; ++i)
{
- rhs[i] = Array<OneD, NekDouble>(m_npoints);
- Vmath::Vcopy(m_npoints, m_viscous[m_intSteps - 1][i], 1, rhs[i], 1);
+ rhs[i] = Array<OneD, NekDouble>(m_npoints, 0.);
}
+ // add viscous term
+ AddVisPressureBCs(fields, rhs, params);
+ // add DuDt term
AddRigidBodyAcc(rhs, params, nptsPlane0);
m_BndExp[m_pressure]->NormVectorIProductWRTBase(
rhs, m_BndExp[m_pressure]->UpdateCoeffs());
diff --git a/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.cpp b/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.cpp
index 730f6bd26c..f19bd6515e 100644
--- a/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.cpp
+++ b/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.cpp
@@ -94,14 +94,14 @@ void StaticWall::v_Update(
}
++m_numCalls;
// pressure
+ Array<OneD, Array<OneD, NekDouble>> rhs(m_bnddim);
for (int i = 0; i < m_bnddim; ++i)
{
- Vmath::Fill(m_npoints, 0., m_viscous[m_intSteps - 1][i], 1);
+ rhs[i] = Array<OneD, NekDouble>(m_npoints, 0.);
}
- AddVisPressureBCs(fields, m_viscous[m_intSteps - 1], params);
- ExtrapolateArray(m_numCalls, m_viscous);
+ AddVisPressureBCs(fields, rhs, params);
m_BndExp[m_pressure]->NormVectorIProductWRTBase(
- m_viscous[m_intSteps - 1], m_BndExp[m_pressure]->UpdateCoeffs());
+ rhs, m_BndExp[m_pressure]->UpdateCoeffs());
}
} // namespace Nektar
diff --git a/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.h b/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.h
index c697813730..96ab793a60 100644
--- a/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.h
+++ b/solvers/IncNavierStokesSolver/BoundaryConditions/StaticWall.h
@@ -76,8 +76,6 @@ protected:
const Array<OneD, const Array<OneD, NekDouble>> &Adv,
std::map<std::string, NekDouble> ¶ms) override;
- Array<OneD, Array<OneD, Array<OneD, NekDouble>>> m_viscous;
- int m_pressure;
StaticWall(const LibUtilities::SessionReaderSharedPtr pSession,
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
Array<OneD, SpatialDomains::BoundaryConditionShPtr> cond,
diff --git a/solvers/IncNavierStokesSolver/BoundaryConditions/TransMovingWall.cpp b/solvers/IncNavierStokesSolver/BoundaryConditions/TransMovingWall.cpp
index 4dac0b91c1..ea5ffa382f 100644
--- a/solvers/IncNavierStokesSolver/BoundaryConditions/TransMovingWall.cpp
+++ b/solvers/IncNavierStokesSolver/BoundaryConditions/TransMovingWall.cpp
@@ -106,13 +106,16 @@ void TransMovingWall::v_Update(
++m_numCalls;
int nptsPlane0 = 0;
SetNumPointsOnPlane0(nptsPlane0);
+ Array<OneD, Array<OneD, NekDouble>> rhs(m_bnddim);
for (int i = 0; i < m_bnddim; ++i)
{
- Vmath::Fill(m_npoints, 0., m_viscous[m_intSteps - 1][i], 1);
+ rhs[i] = Array<OneD, NekDouble>(m_npoints, 0.);
}
- AddVisPressureBCs(fields, m_viscous[m_intSteps - 1], params);
- ExtrapolateArray(m_numCalls, m_viscous);
- // Add dudt
+ // add viscous term
+ AddVisPressureBCs(fields, rhs, params);
+ // Add DuDt
+ std::map<std::string, NekDouble> transParams;
+ std::vector<std::string> accStr = {"A_x", "A_y", "A_z"};
NekDouble time = 0., dt2 = 2. * m_dt;
if (params.find("Time") != params.end())
{
@@ -120,11 +123,8 @@ void TransMovingWall::v_Update(
}
std::vector<NekDouble> times = {time - dt2, time - m_dt, time + m_dt,
time + dt2};
- Array<OneD, Array<OneD, NekDouble>> rhs(m_bnddim);
for (int i = 0; i < m_bnddim; ++i)
{
- rhs[i] = Array<OneD, NekDouble>(m_npoints);
- Vmath::Vcopy(m_npoints, m_viscous[m_intSteps - 1][i], 1, rhs[i], 1);
if (m_BndConds.find(i) != m_BndConds.end() && nptsPlane0)
{
NekDouble dudt = 0.;
@@ -136,10 +136,10 @@ void TransMovingWall::v_Update(
{
dudt += Fourth_Coeffs[j] * equ.Evaluate(0, 0., 0., times[j]);
}
- Vmath::Sadd(nptsPlane0, -dudt, m_viscous[m_intSteps - 1][i], 1,
- rhs[i], 1);
+ transParams[accStr[i]] = dudt;
}
}
+ AddRigidBodyAcc(rhs, transParams, nptsPlane0);
m_BndExp[m_pressure]->NormVectorIProductWRTBase(
rhs, m_BndExp[m_pressure]->UpdateCoeffs());
}
diff --git a/solvers/IncNavierStokesSolver/Tests/FreeFallCyl.tst b/solvers/IncNavierStokesSolver/Tests/FreeFallCyl.tst
index 19e7b5be22..5782d11b30 100644
--- a/solvers/IncNavierStokesSolver/Tests/FreeFallCyl.tst
+++ b/solvers/IncNavierStokesSolver/Tests/FreeFallCyl.tst
@@ -2,7 +2,7 @@
<test>
<description>Moving reference frame formualtion of a free falling circular cylinder</description>
<executable>IncNavierStokesSolver</executable>
- <parameters>FreeFallCyl.xml FreeFallCylc.xml</parameters>
+ <parameters>FreeFallCyl.xml FreeFallCylc.xml --set-start-time 0 --set-start-chknumber 0</parameters>
<files>
<file description="Mesh File">FreeFallCyl.xml</file>
<file description="Session File">FreeFallCylc.xml</file>
--
GitLab