Commit 37d2ad71 authored by David Moxey's avatar David Moxey

Update regression tests

parent 798924b4
......@@ -153,17 +153,20 @@ namespace Nektar
* @brief Set up the Stokes solution used to impose constant flowrate
* through a boundary.
*
* This routine solves a Stokes equation using a unit forcing that is
* determined by the user in order to obtain a correction field applied at
* the end of each timestep to impose a constant volumetric flow rate
* through a user-defined surface.
* This routine solves a Stokes equation using a unit forcing direction,
* specified by the user to be in the desired flow direction. This field can
* then be used to correct the end of each timestep to impose a constant
* volumetric flow rate through a user-defined boundary.
*
* There are essentially three modes of operation:
* There are three modes of operation:
*
* - Standard two-dimensional or three-dimensional simulations
* - Standard two-dimensional or three-dimensional simulations (e.g. pipes
* or channels)
* - 3DH1D simulations where the forcing is not in the homogeneous
* direction.
* - 3DH1D simulations where the forcing is in the homogeneous direction.
* direction (e.g. channel flow, where the y-direction of the 2D mesh
* is perpendicular to the wall);
* - 3DH1D simulations where the forcing is in the homogeneous direction
* (e.g. pipe flow in the z-direction).
*
* In the first two cases, the user should define:
* - the `Flowrate` parameter, which dictates the volumetric flux through
......@@ -173,11 +176,12 @@ namespace Nektar
* - define a `FlowrateForce` function with components `ForceX`, `ForceY`
* and `ForceZ` that defines a unit forcing in the appropriate direction.
*
* In the latter case, the user should define the `Flowrate` and the
* `FlowrateForce` function. The reference area is taken to be the
* homogeneous plane.
* In the latter case, the user should define only the `Flowrate`; the
* reference area is taken to be the homogeneous plane and the force is
* assumed to be the unit z-vector \f$ \hat{e}_z \f$.
*
* We then solve one timestep of the Stokes problem
* This routine solves a single timestep of the Stokes problem
* (premultiplied by the backwards difference coefficient):
*
* \f[ \frac{\partial\mathbf{u}}{\partial t} = -\nabla p +
* \nu\nabla^2\mathbf{u} + \mathbf{f} \f]
......@@ -188,12 +192,12 @@ namespace Nektar
*
* \f[ \alpha = \frac{\overline{Q} - Q(\mathbf{u^n})}{Q(\mathbf{u}_s)} \f]
*
* where \f$ Q(\cdot)\f$ is implemented in
* VelocityCorrectionScheme::MeasureFlowrate that calculates the volumetric
* flux. For more details, see chapter 3.2 of the thesis of D. Moxey
* (University of Warwick, 2011).
* where \f$ Q(\cdot)\f$ is the volumetric flux through the appropriate
* surface or line, which is implemented in
* VelocityCorrectionScheme::MeasureFlowrate. For more details, see chapter
* 3.2 of the thesis of D. Moxey (University of Warwick, 2011).
*/
void VelocityCorrectionScheme::SetupFlowrate()
void VelocityCorrectionScheme::SetupFlowrate(NekDouble aii_dt)
{
m_flowrateBndID = -1;
m_flowrateArea = 0.0;
......@@ -202,12 +206,21 @@ namespace Nektar
m_fields[0]->GetBndConditions();
std::string forces[] = { "X", "Y", "Z" };
Array<OneD, NekDouble> flowrateForce(m_spacedim);
Array<OneD, NekDouble> flowrateForce(m_spacedim, 0.0);
// Set up flowrate forces.
bool defined = true;
for (int i = 0; i < m_spacedim; ++i)
{
std::string varName = std::string("Force") + forces[i];
ASSERTL0(m_session->DefinesFunction("FlowrateForce", varName),
defined = m_session->DefinesFunction("FlowrateForce", varName);
if (!defined && m_HomogeneousType == eHomogeneous1D)
{
break;
}
ASSERTL0(defined,
"A 'FlowrateForce' function must defined with components "
"[ForceX, ...] to define direction of flowrate forcing");
......@@ -216,6 +229,14 @@ namespace Nektar
flowrateForce[i] = ffunc->Evaluate();
}
// For 3DH1D simulations, if force isn't defined then assume in
// z-direction.
if (!defined)
{
flowrateForce[2] = 1.0;
}
// Find the boundary condition that is tagged as the flowrate boundary.
for (int i = 0; i < bcs.num_elements(); ++i)
{
if (boost::iequals(bcs[i]->GetUserDefined(), "Flowrate"))
......@@ -231,13 +252,16 @@ namespace Nektar
"One boundary region must be marked using the 'Flowrate' "
"user-defined type to monitor the volumetric flowrate.");
// Extract an appropriate expansion list to represents the boundary.
if (m_flowrateBndID >= 0)
{
// For a boundary, extract the boundary itself.
m_flowrateBnd = m_fields[0]->GetBndCondExpansions()[m_flowrateBndID];
}
else if (m_HomogeneousType == eHomogeneous1D)
{
// Get z IDs to find zero-th plane
// For 3DH1D simulations with no force specified, find the mean
// (0th) plane.
Array<OneD, unsigned int> zIDs = m_fields[0]->GetZIDs();
int tmpId = -1;
......@@ -259,6 +283,11 @@ namespace Nektar
}
}
// At this point, some processors may not have m_flowrateBnd set if they
// don't contain the appropriate boundary. To calculate the area, we
// integrate 1.0 over the boundary (which has been set up with the
// appropriate subcommunicator to avoid deadlock), and then communicate
// this to the other processors with an AllReduce.
if (m_flowrateBnd)
{
Array<OneD, NekDouble> inArea(m_flowrateBnd->GetNpoints(), 1.0);
......@@ -266,6 +295,9 @@ namespace Nektar
}
m_comm->AllReduce(m_flowrateArea, LibUtilities::ReduceMax);
// Set up some storage for the Stokes solution (to be stored in
// m_flowrateStokes) and its initial condition (inTmp), which holds the
// unit forcing.
int nqTot = m_fields[0]->GetNpoints();
Array<OneD, Array<OneD, NekDouble> > inTmp(m_spacedim);
m_flowrateStokes = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
......@@ -273,8 +305,9 @@ namespace Nektar
for (int i = 0; i < m_spacedim; ++i)
{
inTmp[i] = Array<OneD, NekDouble>(
nqTot, flowrateForce[i] * m_timestep);
nqTot, flowrateForce[i] * aii_dt);
m_flowrateStokes[i] = Array<OneD, NekDouble>(nqTot, 0.0);
if (m_HomogeneousType == eHomogeneous1D)
{
Array<OneD, NekDouble> inTmp2(nqTot);
......@@ -282,17 +315,28 @@ namespace Nektar
m_fields[i]->SetWaveSpace(true);
inTmp[i] = inTmp2;
}
Vmath::Zero(m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
Vmath::Zero(
m_fields[i]->GetNcoeffs(), m_fields[i]->UpdateCoeffs(), 1);
}
// Create temporary extrapolation object to avoid issues with
// m_extrapolation for HOPBCs using higher order timestepping schemes.
ExtrapolateSharedPtr tmpExtrap = m_extrapolation;
m_extrapolation = GetExtrapolateFactory().CreateInstance(
"Standard", m_session, m_fields, m_pressure, m_velocity,
m_advObject);
// Finally, calculate the solution and the flux of the Stokes
// solution. We set m_greenFlux to maximum numeric limit, which signals
// to SolveUnsteadyStokesSystem that we don't need to apply a flowrate
// force.
m_greenFlux = numeric_limits<NekDouble>::max();
SolveUnsteadyStokesSystem(inTmp, m_flowrateStokes, 0.0, m_timestep);
m_flowrateAiidt = aii_dt;
SolveUnsteadyStokesSystem(inTmp, m_flowrateStokes, 0.0, aii_dt);
m_greenFlux = MeasureFlowrate(m_flowrateStokes);
// Reset extrapolation.
SetUpExtrapolation();
// Open field
// If the user specified IO_FlowSteps, open a handle to store output.
if (m_comm->GetRank() == 0 && m_flowrateSteps)
{
std::string filename = m_session->GetSessionName();
......@@ -303,6 +347,8 @@ namespace Nektar
<< "# -------------------------------------------"
<< endl;
}
m_extrapolation = tmpExtrap;
}
/**
......@@ -324,6 +370,8 @@ namespace Nektar
if (m_flowrateBnd && m_flowrateBndID >= 0)
{
// If we're an actual boundary, calculate the vector flux through
// the boundary.
Array<OneD, Array<OneD, NekDouble> > boundary(m_spacedim);
for (int i = 0; i < m_spacedim; ++i)
......@@ -337,7 +385,7 @@ namespace Nektar
}
else if (m_flowrateBnd)
{
// Homogeneous case
// 3DH1D case: compute flux through the zero-th (mean) plane.
flowrate = m_flowrateBnd->Integral(inarray[2]);
// Now communicate this with other planes
......@@ -352,6 +400,7 @@ namespace Nektar
}
else
{
// Remaining processors that do not contain boundary.
m_comm->AllReduce(flowrate, LibUtilities::ReduceSum);
}
......@@ -466,11 +515,7 @@ namespace Nektar
m_F[i] = Array< OneD, NekDouble> (m_fields[0]->GetTotPoints(), 0.0);
}
// Set up flowrate before m_fields are initialised.
if (m_flowrate > 0.0)
{
SetupFlowrate();
}
m_flowrateAiidt = 0.0;
AdvectionSystem::v_DoInitialise();
......@@ -582,6 +627,13 @@ namespace Nektar
const NekDouble time,
const NekDouble aii_Dt)
{
// Set up flowrate if we're starting for the first time or the value of
// aii_Dt has changed.
if (m_flowrate > 0.0 && (aii_Dt != m_flowrateAiidt))
{
SetupFlowrate(aii_Dt);
}
int physTot = m_fields[0]->GetTotPoints();
// Substep the pressure boundary condition if using substepping
......
......@@ -151,8 +151,10 @@ namespace Nektar
std::ofstream m_flowrateStream;
/// Interval at which to record flowrate data
int m_flowrateSteps;
/// Value of aii_dt used to compute Stokes flowrate solution.
NekDouble m_flowrateAiidt;
void SetupFlowrate();
void SetupFlowrate(NekDouble aii_dt);
NekDouble MeasureFlowrate(
const Array<OneD, Array<OneD, NekDouble> > &inarray);
......
......@@ -9,16 +9,16 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">8.7003e-16</value>
<value variable="v" tolerance="1e-12">1.18471e-15</value>
<value variable="w" tolerance="1e-12">8.94953e-07</value>
<value variable="p" tolerance="1e-8">1.91856e-15</value>
<value variable="u" tolerance="1e-12">0</value>
<value variable="v" tolerance="1e-12">0</value>
<value variable="w" tolerance="1e-12">6.57092e-15</value>
<value variable="p" tolerance="1e-8">0</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">5.70037e-15</value>
<value variable="v" tolerance="1e-12">5.63826e-15</value>
<value variable="w" tolerance="1e-12">2.77108e-06</value>
<value variable="p" tolerance="1e-8">9.82442e-15</value>
<value variable="u" tolerance="1e-12">2.93546e-18</value>
<value variable="v" tolerance="1e-12">2.93958e-18</value>
<value variable="w" tolerance="1e-12">2.37588e-14</value>
<value variable="p" tolerance="1e-8">2.74618e-17</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<NEKTAR>
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="6" FIELDS="u,v,w,p" TYPE="MODIFIED" />
<E COMPOSITE="C[0]" NUMMODES="4" FIELDS="u,v,w,p" TYPE="MODIFIED" />
</EXPANSIONS>
<CONDITIONS>
<SOLVERINFO>
......@@ -14,11 +14,11 @@
</SOLVERINFO>
<PARAMETERS>
<P> TimeStep = 0.001 </P>
<P> TimeStep = 0.01 </P>
<P> NumSteps = 100 </P>
<P> IO_CheckSteps = 1000 </P>
<P> IO_InfoSteps = 1000 </P>
<P> Kinvis = 1/200 </P>
<P> Kinvis = 1 </P>
<P> HomModesZ = 8 </P>
<P> LZ = 1.0 </P>
<P> Flowrate = 1.0 </P>
......@@ -58,16 +58,10 @@
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="FlowrateForce">
<E VAR="ForceX" VALUE="0" />
<E VAR="ForceY" VALUE="0" />
<E VAR="ForceZ" VALUE="1" />
</FUNCTION>
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="0" />
<E VAR="v" VALUE="0" />
<E VAR="w" VALUE="6*y*(1-y)" />
<E VAR="w" VALUE="0" />
<E VAR="p" VALUE="0" />
</FUNCTION>
......
......@@ -9,16 +9,16 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">1.54851e-15</value>
<value variable="v" tolerance="1e-12">1.92039e-15</value>
<value variable="w" tolerance="1e-12">8.94953e-07</value>
<value variable="p" tolerance="1e-8">3.23051e-15</value>
<value variable="u" tolerance="1e-12">0</value>
<value variable="v" tolerance="1e-12">0</value>
<value variable="w" tolerance="1e-12">7.58667e-15</value>
<value variable="p" tolerance="1e-8">0</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">9.74693e-15</value>
<value variable="v" tolerance="1e-12">7.87791e-15</value>
<value variable="w" tolerance="1e-12">2.77108e-06</value>
<value variable="p" tolerance="1e-8">1.64691e-14</value>
<value variable="u" tolerance="1e-12">1.76842e-18</value>
<value variable="v" tolerance="1e-12">1.77358e-18</value>
<value variable="w" tolerance="1e-12">3.01981e-14</value>
<value variable="p" tolerance="1e-8">2.08802e-17</value>
</metric>
</metrics>
</test>
......@@ -9,14 +9,14 @@
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">7.66614e-11</value>
<value variable="v" tolerance="1e-12">4.52553e-12</value>
<value variable="p" tolerance="1e-8">7.0151e-09</value>
<value variable="u" tolerance="1e-12">9.00299e-13</value>
<value variable="v" tolerance="1e-12">1.72538e-12</value>
<value variable="p" tolerance="1e-8">5.92475e-10</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">1.1156e-10</value>
<value variable="v" tolerance="1e-12">1.17928e-11</value>
<value variable="p" tolerance="1e-8">8.32526e-09</value>
<value variable="u" tolerance="1e-12">2.27274e-12</value>
<value variable="v" tolerance="1e-12">3.14404e-12</value>
<value variable="p" tolerance="1e-8">9.17508e-10</value>
</metric>
</metrics>
</test>
......@@ -10,15 +10,15 @@
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0</value>
<value variable="v" tolerance="1e-12">0</value>
<value variable="w" tolerance="1e-12">1.13277e-15</value>
<value variable="p" tolerance="1e-8">4.56093e-14</value>
<value variable="v" tolerance="1e-12">5.94828e-16</value>
<value variable="w" tolerance="1e-12">1.34946e-15</value>
<value variable="p" tolerance="1e-8">9.43679e-14</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">1.48905e-16</value>
<value variable="v" tolerance="1e-12">3.12187e-16</value>
<value variable="w" tolerance="1e-12">1.90958e-14</value>
<value variable="p" tolerance="1e-8">5.11546e-14</value>
<value variable="u" tolerance="1e-12">3.76927e-16</value>
<value variable="v" tolerance="1e-12">1.70832e-15</value>
<value variable="w" tolerance="1e-12">1.93179e-14</value>
<value variable="p" tolerance="1e-8">1.49283e-13</value>
</metric>
</metrics>
</test>
......@@ -154,7 +154,7 @@
<I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
<I PROPERTY="AdvectionForm" VALUE="Convective" />
<I PROPERTY="Projection" VALUE="Galerkin" />
<I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
<I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder3" />
</SOLVERINFO>
<PARAMETERS>
......
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