Commit 965e1cb3 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo

Solved conflict

parents 0ede3856 f9450fcb
......@@ -23,6 +23,7 @@ v4.4.0
to boundary conditions (!615)
- Allow expansions to be loaded directly from field file (!617)
- New options for load balancing (DOF or BOUNDARY) in mesh partitioner (!617)
- Update Body/Field forces at each timestep (!665)
**ADRSolver:**
- Add a projection equation system for C^0 projections (!675)
......@@ -31,6 +32,7 @@ v4.4.0
- Use a continuous basefield projection and revert to constant c formulation (!664)
- Added ability to compute CFL number (!664)
- Output Sourceterm (!664)
- Use the Forcing framework to define source terms (!665)
**IncNavierStokesSolver:**
- Add ability to simulate additional scalar fields (!624)
......
......@@ -118,6 +118,16 @@ namespace Nektar
return vForceList;
}
Nektar::Array<OneD, Array<OneD, NekDouble> > &Forcing::UpdateForces()
{
return m_Forcing;
}
const Nektar::Array<OneD, const Array<OneD, NekDouble> > &Forcing::GetForces()
{
return m_Forcing;
}
void Forcing::EvaluateTimeFunction(
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
......
......@@ -90,6 +90,11 @@ namespace SolverUtils
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const unsigned int& pNumForcingFields = 0);
SOLVER_UTILS_EXPORT const Array<OneD, const Array<OneD, NekDouble> >
&GetForces();
SOLVER_UTILS_EXPORT Array<OneD, Array<OneD, NekDouble> > &UpdateForces();
protected:
/// Session reader
LibUtilities::SessionReaderSharedPtr m_session;
......
......@@ -65,7 +65,6 @@ namespace SolverUtils
const TiXmlElement* pForce)
{
m_NumVariable = pNumForcingFields;
int nq = pFields[0]->GetTotPoints();
const TiXmlElement* funcNameElmt = pForce->FirstChildElement("BODYFORCE");
if(!funcNameElmt)
......@@ -76,9 +75,15 @@ namespace SolverUtils
"specifying function name which prescribes body force.");
}
string funcName = funcNameElmt->GetText();
ASSERTL0(m_session->DefinesFunction(funcName),
"Function '" + funcName + "' not defined.");
m_funcName = funcNameElmt->GetText();
ASSERTL0(m_session->DefinesFunction(m_funcName),
"Function '" + m_funcName + "' not defined.");
bool singleMode, halfMode;
m_session->MatchSolverInfo("ModeType","SingleMode",singleMode,false);
m_session->MatchSolverInfo("ModeType","HalfMode", halfMode, false);
bool homogeneous = pFields[0]->GetExpType() == MultiRegions::e3DH1D;
m_transform = (singleMode || halfMode || homogeneous);
// Time function is optional
funcNameElmt = pForce->FirstChildElement("BODYFORCETIMEFCN");
......@@ -104,43 +109,48 @@ namespace SolverUtils
}
m_Forcing = Array<OneD, Array<OneD, NekDouble> > (m_NumVariable);
for (int i = 0; i < m_NumVariable; ++i)
{
m_Forcing[i] = Array<OneD, NekDouble> (pFields[0]->GetTotPoints(), 0.0);
}
Update(pFields, 0.0);
}
std::string s_FieldStr;
void ForcingBody::Update(
const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields,
const NekDouble &time)
{
for (int i = 0; i < m_NumVariable; ++i)
{
m_Forcing[i] = Array<OneD, NekDouble> (nq, 0.0);
s_FieldStr = m_session->GetVariable(i);
ASSERTL0(m_session->DefinesFunction(funcName, s_FieldStr),
std::string s_FieldStr = m_session->GetVariable(i);
ASSERTL0(m_session->DefinesFunction(m_funcName, s_FieldStr),
"Variable '" + s_FieldStr + "' not defined.");
EvaluateFunction(pFields, m_session, s_FieldStr,
m_Forcing[i], funcName);
m_Forcing[i], m_funcName, time);
}
bool singleMode, halfMode;
m_session->MatchSolverInfo("ModeType","SingleMode",singleMode,false);
m_session->MatchSolverInfo("ModeType","HalfMode", halfMode, false);
bool homogeneous = pFields[0]->GetExpType() == MultiRegions::e3DH1D;
// If singleMode or halfMode, transform the forcing term to be in
// physical space in the plane, but Fourier space in the homogeneous
// direction
if (singleMode || halfMode || homogeneous)
if (m_transform)
{
for (int i = 0; i < m_NumVariable; ++i)
{
pFields[0]->HomogeneousFwdTrans(m_Forcing[i],m_Forcing[i]);
pFields[0]->HomogeneousFwdTrans(m_Forcing[i], m_Forcing[i]);
}
}
}
void ForcingBody::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
if(m_hasTimeFcnScaling)
{
Array<OneD, NekDouble> TimeFcn(1);
......@@ -157,6 +167,8 @@ namespace SolverUtils
}
else
{
Update(fields, time);
for (int i = 0; i < m_NumVariable; i++)
{
Vmath::Vadd(outarray[i].num_elements(), outarray[i], 1,
......
......@@ -85,11 +85,18 @@ namespace SolverUtils
const NekDouble &time);
private:
std::string m_funcName;
bool m_hasTimeFcnScaling;
LibUtilities::EquationSharedPtr m_timeFcnEqn;
bool m_transform;
ForcingBody(const LibUtilities::SessionReaderSharedPtr& pSession);
virtual ~ForcingBody(void){};
bool m_hasTimeFcnScaling;
LibUtilities::EquationSharedPtr m_timeFcnEqn;
void Update(
const Array<OneD, MultiRegions::ExpListSharedPtr>& pFields,
const NekDouble &time);
};
}
......
......@@ -138,22 +138,8 @@ void APE::v_InitObject()
m_bf[i] = Array<OneD, NekDouble>(GetTotPoints());
}
EvaluateFunction(m_bfNames, m_bf, "Baseflow", m_time);
Array<OneD, NekDouble> tmpC(GetNcoeffs());
for (int i = 0; i < m_spacedim + 2; ++i)
{
m_bfField->IProductWRTBase(m_bf[i], tmpC);
m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
m_bfField->LocalToGlobal(tmpC, tmpC);
m_bfField->GlobalToLocal(tmpC, tmpC);
m_bfField->BwdTrans(tmpC, m_bf[i]);
}
// Initialize the sourceterm
m_sourceTerms = Array<OneD, NekDouble>(GetTotPoints(), 0.0);
EvaluateFunction("S", m_sourceTerms, "Source", m_time);
m_fields[0]->FwdTrans(m_sourceTerms, tmpC);
m_fields[0]->BwdTrans(tmpC, m_sourceTerms);
m_forcing = SolverUtils::Forcing::Load(m_session, m_fields, m_spacedim + 1);
// Do not forwards transform initial condition
m_homoInitialFwd = false;
......@@ -313,24 +299,24 @@ void APE::GetFluxVector(
/**
* @brief v_PostIntegrate
*/
bool APE::v_PostIntegrate(int step)
bool APE::v_PreIntegrate(int step)
{
if (m_cflsteps && !((step + 1) % m_cflsteps))
EvaluateFunction(m_bfNames, m_bf, "Baseflow", m_time);
Array<OneD, NekDouble> tmpC(GetNcoeffs());
std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
{
NekDouble cfl = GetCFLEstimate();
if (m_comm->GetRank() == 0)
for (int i = 0; i < (*x)->GetForces().num_elements(); ++i)
{
cout << "CFL: " << cfl << endl;
m_bfField->IProductWRTBase((*x)->GetForces()[i], tmpC);
m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
m_bfField->LocalToGlobal(tmpC, tmpC);
m_bfField->GlobalToLocal(tmpC, tmpC);
m_bfField->BwdTrans(tmpC, (*x)->UpdateForces()[i]);
}
}
EvaluateFunction("S", m_sourceTerms, "Source", m_time);
EvaluateFunction(m_bfNames, m_bf, "Baseflow", m_time);
Array<OneD, NekDouble> tmpC(GetNcoeffs());
m_fields[0]->FwdTrans(m_sourceTerms, tmpC);
m_fields[0]->BwdTrans(tmpC, m_sourceTerms);
for (int i = 0; i < m_spacedim + 2; ++i)
{
// ensure the field is C0-continuous
......@@ -341,6 +327,24 @@ bool APE::v_PostIntegrate(int step)
m_bfField->BwdTrans(tmpC, m_bf[i]);
}
return UnsteadySystem::v_PreIntegrate(step);
}
/**
* @brief v_PostIntegrate
*/
bool APE::v_PostIntegrate(int step)
{
if (m_cflsteps && !((step + 1) % m_cflsteps))
{
NekDouble cfl = GetCFLEstimate();
if (m_comm->GetRank() == 0)
{
cout << "CFL: " << cfl << endl;
}
}
return UnsteadySystem::v_PostIntegrate(step);
}
......@@ -365,7 +369,11 @@ void APE::DoOdeRhs(const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Vmath::Neg(nq, outarray[i], 1);
}
AddSource(outarray);
std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
{
(*x)->Apply(m_fields, outarray, outarray, m_time);
}
}
......@@ -493,15 +501,6 @@ void APE::WallBC(int bcRegion, int cnt,
}
/**
* @brief sourceterm for p' equation obtained from GetSource
*/
void APE::AddSource(Array< OneD, Array< OneD, NekDouble > > &outarray)
{
Vmath::Vadd(GetTotPoints(), m_sourceTerms, 1, outarray[0], 1, outarray[0], 1);
}
/**
* @brief Compute the advection velocity in the standard space
* for each element of the expansion.
......@@ -602,8 +601,6 @@ void APE::v_ExtraFldOutput(
std::vector<Array<OneD, NekDouble> > &fieldcoeffs,
std::vector<std::string> &variables)
{
const int nCoeffs = m_fields[0]->GetNcoeffs();
for (int i = 0; i < m_spacedim + 2; i++)
{
Array<OneD, NekDouble> tmpC(GetNcoeffs());
......@@ -613,16 +610,30 @@ void APE::v_ExtraFldOutput(
m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
m_bfField->LocalToGlobal(tmpC, tmpC);
m_bfField->GlobalToLocal(tmpC, tmpC);
m_bfField->BwdTrans(tmpC, m_bf[i]);
variables.push_back(m_bfNames[i]);
fieldcoeffs.push_back(tmpC);
}
variables.push_back("S");
Array<OneD, NekDouble> FwdS(nCoeffs);
m_fields[0]->FwdTrans(m_sourceTerms, FwdS);
fieldcoeffs.push_back(FwdS);
int f = 0;
std::vector<SolverUtils::ForcingSharedPtr>::const_iterator x;
for (x = m_forcing.begin(); x != m_forcing.end(); ++x)
{
for (int i = 0; i < (*x)->GetForces().num_elements(); ++i)
{
Array<OneD, NekDouble> tmpC(GetNcoeffs());
m_bfField->IProductWRTBase((*x)->GetForces()[i], tmpC);
m_bfField->MultiplyByElmtInvMass(tmpC, tmpC);
m_bfField->LocalToGlobal(tmpC, tmpC);
m_bfField->GlobalToLocal(tmpC, tmpC);
variables.push_back("F_" + boost::lexical_cast<string>(f) +
"_" + m_session->GetVariable(i));
fieldcoeffs.push_back(tmpC);
}
f++;
}
}
......
......@@ -39,6 +39,7 @@
#include <SolverUtils/UnsteadySystem.h>
#include <SolverUtils/Advection/Advection.h>
#include <SolverUtils/Forcing/Forcing.h>
#include <SolverUtils/RiemannSolvers/RiemannSolver.h>
using namespace Nektar::SolverUtils;
......@@ -72,12 +73,12 @@ class APE : public UnsteadySystem
protected:
SolverUtils::AdvectionSharedPtr m_advection;
std::vector<SolverUtils::ForcingSharedPtr> m_forcing;
SolverUtils::RiemannSolverSharedPtr m_riemannSolver;
Array<OneD, Array<OneD, NekDouble> > m_traceBasefield;
Array<OneD, Array<OneD, NekDouble> > m_vecLocs;
/// Isentropic coefficient, Ratio of specific heats (APE)
NekDouble m_gamma;
Array<OneD, NekDouble> m_sourceTerms;
Array<OneD, Array<OneD, NekDouble> > m_bf;
MultiRegions::ExpListSharedPtr m_bfField;
std::vector<std::string> m_bfNames;
......@@ -101,9 +102,9 @@ class APE : public UnsteadySystem
const Array<OneD, Array<OneD, NekDouble> > &physfield,
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > &flux);
virtual bool v_PostIntegrate(int step);
virtual bool v_PreIntegrate(int step);
void AddSource(Array< OneD, Array< OneD, NekDouble > >& outarray);
virtual bool v_PostIntegrate(int step);
void GetStdVelocity(Array< OneD, NekDouble >& stdV);
......
......@@ -806,10 +806,6 @@
<E VAR="p0" VALUE="Pinfinity" />
<E VAR="rho0" VALUE="Rho0" />
</FUNCTION>
<FUNCTION NAME="Source"> <!-- Incompressible base flow -->
<E VAR="S" VALUE="0" />
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0" />
......
......@@ -358,11 +358,7 @@
<E VAR="p0" VALUE="Pinfinity" />
<E VAR="rho0" VALUE="Rho0" />
</FUNCTION>
<FUNCTION NAME="Source"> <!-- Incompressible base flow -->
<E VAR="S" VALUE="0" />
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0" />
<E VAR="u" VALUE="0" />
......
......@@ -75,7 +75,7 @@
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="pmax*omega0*cos(omega0*t)*exp(-32*x^2)"/>
<E VAR="p" VALUE="pmax*omega0*cos(omega0*t)*exp(-32*x^2)"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
......@@ -86,4 +86,9 @@
<E VAR="u" VALUE="0"/>
</FUNCTION>
</CONDITIONS>
<FORCING>
<FORCE TYPE="Field">
<FIELDFORCE> Source <FIELDFORCE/>
</FORCE>
</FORCING>
</NEKTAR>
......@@ -79,9 +79,6 @@
<E VAR="p0" VALUE="Pinfinity"/>
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
......
......@@ -796,9 +796,6 @@
<E VAR="u0" VALUE="200" />
<E VAR="v0" VALUE="0" />
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
......
......@@ -8,14 +8,14 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="p" tolerance="1e-4">21.3387</value>
<value variable="u" tolerance="1e-7">0.0101445</value>
<value variable="v" tolerance="1e-7">0.0105278</value>
<value variable="p" tolerance="1e-4">6.77218</value>
<value variable="u" tolerance="1e-7">0.00160109</value>
<value variable="v" tolerance="1e-7">0.00163104</value>
</metric>
<metric type="Linf" id="2">
<value variable="p" tolerance="1e-4">91.6371</value>
<value variable="u" tolerance="1e-7">0.0383799</value>
<value variable="v" tolerance="1e-7">0.0381269</value>
<value variable="p" tolerance="1e-4">30.1206</value>
<value variable="u" tolerance="1e-7">0.00609606</value>
<value variable="v" tolerance="1e-7">0.00607659</value>
</metric>
</metrics>
</test>
......@@ -317,7 +317,7 @@
<I PROPERTY="UpwindType" VALUE="LaxFriedrichs"/>
</SOLVERINFO>
<PARAMETERS>
<P> TimeStep = 0.00001 </P>
<P> TimeStep = 1e-05 </P>
<P> NumSteps = 10 </P>
<P> FinTime = TimeStep*NumSteps </P>
<P> IO_CheckSteps = 10 </P>
......@@ -347,18 +347,25 @@
<E VAR="p0" VALUE="Pinfinity"/>
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
<E VAR="v" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="InitialConditions">
<E VAR="p" VALUE="100*exp(-32*(x^2+y^2))"/> <!-- Gaussian pulse located at the origin -->
<E VAR="p" VALUE="0"/> <!-- Gaussian pulse located at the origin -->
<E VAR="u" VALUE="0"/>
<E VAR="v" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="p" VALUE="100 * 2*PI*5E2 * cos(2*PI*5E2 * t) * exp(-32*(x^2+y^2))"/>
<E VAR="u" VALUE="0"/>
<E VAR="v" VALUE="0"/>
</FUNCTION>
</CONDITIONS>
<FORCING>
<FORCE TYPE="Field">
<FIELDFORCE> Source <FIELDFORCE/>
</FORCE>
</FORCING>
</NEKTAR>
......@@ -347,9 +347,6 @@
<E VAR="p0" VALUE="Pinfinity"/>
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p,u,v" VALUE="0"/>
</FUNCTION>
......
......@@ -353,9 +353,6 @@
<E VAR="p0" VALUE="Pinfinity"/>
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
......
......@@ -130,9 +130,6 @@
<E VAR="u0" VALUE="100" />
<E VAR="v0" VALUE="0" />
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
......
......@@ -1405,9 +1405,6 @@
<E VAR="p0" VALUE="Pinfinity"/>
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
......
......@@ -1398,9 +1398,6 @@
<E VAR="p0" VALUE="Pinfinity"/>
<E VAR="rho0" VALUE="Rho0"/>
</FUNCTION>
<FUNCTION NAME="Source">
<E VAR="S" VALUE="0"/>
</FUNCTION>
<FUNCTION NAME="ExactSolution"> <!-- Not really the exact solution -->
<E VAR="p" VALUE="0"/>
<E VAR="u" VALUE="0"/>
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment