Commit 4627afbf authored by Spencer Sherwin's avatar Spencer Sherwin

Merge branch 'fix/Homo1D_ReadBCFromFile' into 'master'

Fix/homo1D read bc from file

See merge request !1069
parents 823a9a84 2023ef4e
Pipeline #1164 passed with stages
in 7 minutes and 48 seconds
......@@ -1214,13 +1214,13 @@ void Mapping::v_UpdateBCs( const NekDouble time)
if( BndConds[n]->GetUserDefined() =="" ||
BndConds[n]->GetUserDefined() =="MovingBody")
{
BndExp[n]->FwdTrans_BndConstrained(BndExp[n]->GetPhys(),
BndExp[n]->UpdateCoeffs());
if (m_fields[i]->GetExpType() == MultiRegions::e3DH1D)
{
BndExp[n]->HomogeneousFwdTrans(BndExp[n]->GetCoeffs(),
BndExp[n]->UpdateCoeffs());
BndExp[n]->SetWaveSpace(false);
}
BndExp[n]->FwdTrans_BndConstrained(
BndExp[n]->GetPhys(), BndExp[n]->UpdateCoeffs());
}
}
}
......
......@@ -705,6 +705,12 @@ namespace Nektar
// Now fill in all other Dirichlet coefficients.
for(int i = 0; i < m_bndCondExpansions.num_elements(); ++i)
{
if (m_bndConditions[i]->GetBoundaryConditionType() ==
SpatialDomains::ePeriodic)
{
continue;
}
Array<OneD, NekDouble>& coeffs = m_bndCondExpansions[i]->UpdateCoeffs();
for(int j = 0; j < (m_bndCondExpansions[i])->GetNcoeffs(); ++j)
......@@ -734,7 +740,13 @@ namespace Nektar
Array<OneD, NekDouble>& coeffs = m_bndCondExpansions[nreg]->UpdateCoeffs();
for(int j = 0; j < nreg; ++j)
{
{
if (m_bndConditions[j]->GetBoundaryConditionType() ==
SpatialDomains::ePeriodic)
{
continue;
}
bndcnt += m_bndCondExpansions[j]->GetNcoeffs();
}
......
......@@ -217,42 +217,7 @@ namespace Nektar
m_useFFT, false,
PlanesBndCondExp, comm);
}
EvaluateBoundaryConditions(0.0, variable);
}
void DisContField3DHomogeneous1D::EvaluateBoundaryConditions(
const NekDouble time,
const std::string varName)
{
int n;
const Array<OneD, const NekDouble> z = m_homogeneousBasis->GetZ();
Array<OneD, NekDouble> local_z(m_planes.num_elements());
for (n = 0; n < m_planes.num_elements(); n++)
{
local_z[n] = z[m_transposition->GetPlaneID(n)];
}
for (n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->EvaluateBoundaryConditions(
time, varName, 0.5*m_lhom*(1.0+local_z[n]));
}
// Fourier transform coefficient space boundary values
// This will only be undertaken for time dependent
// boundary conditions unless time == 0.0 which is the
// case when the method is called from the constructor.
for (n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
if (time == 0.0 ||
m_bndConditions[n]->IsTimeDependent() )
{
m_bndCondExpansions[n]->HomogeneousFwdTrans(
m_bndCondExpansions[n]->GetCoeffs(),
m_bndCondExpansions[n]->UpdateCoeffs());
}
}
v_EvaluateBoundaryConditions(0.0, variable);
}
void DisContField3DHomogeneous1D::v_HelmSolve(
......@@ -315,9 +280,127 @@ namespace Nektar
const NekDouble x2_in,
const NekDouble x3_in)
{
boost::ignore_unused(x2_in, x3_in);
boost::ignore_unused(x2_in,x3_in);
int i;
int npoints;
int nbnd = m_bndCondExpansions.num_elements();
MultiRegions::ExpListSharedPtr locExpList;
for (i = 0; i < nbnd; ++i)
{
if (time == 0.0 || m_bndConditions[i]->IsTimeDependent())
{
locExpList = m_bndCondExpansions[i];
npoints = locExpList->GetNpoints();
Array<OneD, NekDouble> x0(npoints, 0.0);
Array<OneD, NekDouble> x1(npoints, 0.0);
Array<OneD, NekDouble> x2(npoints, 0.0);
Array<OneD, NekDouble> valuesFile(npoints, 1.0), valuesExp(npoints, 1.0);
locExpList->GetCoords(x0, x1, x2);
if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eDirichlet)
{
SpatialDomains::DirichletBCShPtr bcPtr = std::static_pointer_cast<
SpatialDomains::DirichletBoundaryCondition>(
m_bndConditions[i]);
std::string filebcs = bcPtr->m_filename;
std::string exprbcs = bcPtr->m_expr;
if (filebcs != "")
{
ExtractFileBCs(filebcs, bcPtr->GetComm(), varName, locExpList);
valuesFile = locExpList->GetPhys();
}
if (exprbcs != "")
{
LibUtilities::Equation condition =
std::static_pointer_cast<SpatialDomains::
DirichletBoundaryCondition >(
m_bndConditions[i])->m_dirichletCondition;
condition.Evaluate(x0, x1, x2, time, valuesExp);
}
EvaluateBoundaryConditions(time, varName);
Vmath::Vmul(npoints, valuesExp, 1, valuesFile, 1,
locExpList->UpdatePhys(), 1);
// set wave space to false since have set up phys values
locExpList->SetWaveSpace(false);
locExpList->FwdTrans_BndConstrained(
locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eNeumann)
{
SpatialDomains::NeumannBCShPtr bcPtr = std::static_pointer_cast<
SpatialDomains::NeumannBoundaryCondition>(
m_bndConditions[i]);
std::string filebcs = bcPtr->m_filename;
if (filebcs != "")
{
ExtractFileBCs(filebcs, bcPtr->GetComm(), varName, locExpList);
}
else
{
LibUtilities::Equation condition = std::
static_pointer_cast<SpatialDomains::
NeumannBoundaryCondition>(
m_bndConditions[i])->m_neumannCondition;
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
}
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::eRobin)
{
SpatialDomains::RobinBCShPtr bcPtr = std::static_pointer_cast<
SpatialDomains::RobinBoundaryCondition>(
m_bndConditions[i]);
std::string filebcs = bcPtr->m_filename;
if (filebcs != "")
{
ExtractFileBCs(filebcs, bcPtr->GetComm(), varName, locExpList);
}
else
{
LibUtilities::Equation condition = std::
static_pointer_cast<SpatialDomains::
RobinBoundaryCondition>(
m_bndConditions[i])->m_robinFunction;
condition.Evaluate(x0, x1, x2, time,
locExpList->UpdatePhys());
}
locExpList->IProductWRTBase(locExpList->GetPhys(),
locExpList->UpdateCoeffs());
}
else if (m_bndConditions[i]->GetBoundaryConditionType()
== SpatialDomains::ePeriodic)
{
continue;
}
else
{
ASSERTL0(false, "This type of BC not implemented yet");
}
}
}
}
std::shared_ptr<ExpList> &DisContField3DHomogeneous1D::
......
......@@ -801,7 +801,7 @@ namespace Nektar
}
void ExpList::FwdTrans_BndConstrained(
void ExpList::v_FwdTrans_BndConstrained(
const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray)
{
......
......@@ -1297,6 +1297,10 @@ namespace Nektar
const Array<OneD,const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
virtual void v_FwdTrans_BndConstrained(
const Array<OneD,const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray);
virtual void v_SmoothField(Array<OneD,NekDouble> &field);
virtual void v_IProductWRTBase(
......@@ -1746,6 +1750,17 @@ namespace Nektar
v_FwdTrans_IterPerExp(inarray,outarray);
}
/**
*
*/
inline void ExpList::FwdTrans_BndConstrained (
const Array<OneD, const NekDouble> &inarray,
Array<OneD,NekDouble> &outarray)
{
v_FwdTrans_BndConstrained(inarray,outarray);
}
/**
*
*/
......
......@@ -483,7 +483,9 @@ namespace Nektar
/**
* Forward transform element by element
*/
void ExpListHomogeneous1D::v_FwdTrans_IterPerExp(const Array<OneD, const NekDouble> &inarray, Array<OneD, NekDouble> &outarray)
void ExpListHomogeneous1D::v_FwdTrans_IterPerExp(const Array<OneD,
const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int cnt = 0, cnt1 = 0;
Array<OneD, NekDouble> tmparray;
......@@ -500,6 +502,29 @@ namespace Nektar
HomogeneousFwdTrans(outarray,outarray);
}
}
/**
* Forward transform element by element with boundaries constrained
*/
void ExpListHomogeneous1D::v_FwdTrans_BndConstrained(const Array<OneD,
const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int cnt = 0, cnt1 = 0;
Array<OneD, NekDouble> tmparray;
//spectral element FwdTrans plane by plane
for(int n = 0; n < m_planes.num_elements(); ++n)
{
m_planes[n]->FwdTrans_BndConstrained(inarray+cnt, tmparray = outarray + cnt1);
cnt += m_planes[n]->GetTotPoints();
cnt1 += m_planes[n]->GetNcoeffs();
}
if(!m_WaveSpace)
{
HomogeneousFwdTrans(outarray,outarray);
}
}
/**
* Backward transform
......
......@@ -190,7 +190,10 @@ namespace Nektar
virtual void v_FwdTrans_IterPerExp(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
virtual void v_FwdTrans_BndConstrained(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
virtual void v_BwdTrans(const Array<OneD,const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate);
......
......@@ -43,6 +43,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
ADD_NEKTAR_TEST(2DFlow_lineforcing_bcfromfile)
ADD_NEKTAR_TEST(BercovierEngelman)
ADD_NEKTAR_TEST(ChanFlow2D_bcsfromfiles)
ADD_NEKTAR_TEST(ChanFlow_3DH1D_bcsfromfiles)
ADD_NEKTAR_TEST(ChanFlow_3DH1D_MVM LENGTHY)
ADD_NEKTAR_TEST(ChanFlow_3DH2D_MVM LENGTHY)
ADD_NEKTAR_TEST(ChanFlow_LinNS_m8)
......
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>Channel Flow 3DH1D P=2 Boundary Conditions from files</description>
<executable>IncNavierStokesSolver</executable>
<parameters>ChanFlow_3DH1D_bcsfromfiles.xml</parameters>
<files>
<file description="Session File">ChanFlow_3DH1D_bcsfromfiles.xml</file>
<file description="Session File">ChanFlow_3DH1D_bcsfromfiles_walls.bc</file>
<file description="Session File">ChanFlow_3DH1D_bcsfromfiles_inflow.bc</file>
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-6">1.40307e-16</value>
<value variable="v" tolerance="1e-6">0</value>
<value variable="w" tolerance="1e-6">0</value>
<value variable="p" tolerance="1e-6">6.68056e-15</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-6">1.249e-15</value>
<value variable="v" tolerance="1e-6">1.89879e-16</value>
<value variable="w" tolerance="1e-6">9.51393e-18</value>
<value variable="p" tolerance="1e-6">3.68594e-14</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<NEKTAR>
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="3" FIELDS="u,v,w,p" TYPE="MODIFIED" />
</EXPANSIONS>
<CONDITIONS>
<SOLVERINFO>
<I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme"/>
<I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes"/>
<I PROPERTY="AdvectionForm" VALUE="Convective"/>
<I PROPERTY="Projection" VALUE="Galerkin"/>
<I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder2"/>
<I PROPERTY="HOMOGENEOUS" VALUE="1D"/>
</SOLVERINFO>
<PARAMETERS>
<P> TimeStep = 0.001 </P>
<P> NumSteps = 1000 </P>
<P> IO_CheckSteps = 1000 </P>
<P> IO_InfoSteps = 1000 </P>
<P> Kinvis = 1 </P>
<P> HomModesZ = 10 </P>
<P> LZ = 1.0 </P>
</PARAMETERS>
<VARIABLES>
<V ID="0"> u </V>
<V ID="1"> v </V>
<V ID="2"> w </V>
<V ID="3"> p </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
<B ID="2"> C[3] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" FILE="ChanFlow_3DH1D_bcsfromfiles_walls.bc" /> // u = 0
<D VAR="v" FILE="ChanFlow_3DH1D_bcsfromfiles_walls.bc" /> // v = 0
<D VAR="w" FILE="ChanFlow_3DH1D_bcsfromfiles_walls.bc" /> // w = 0
<N VAR="p" USERDEFINEDTYPE="H" VALUE="0" /> // High Order Pressure BC
</REGION>
<REGION REF="1">
<D VAR="u" FILE="ChanFlow_3DH1D_bcsfromfiles_inflow.bc" /> // u = y*(1-y)
<D VAR="v" FILE="ChanFlow_3DH1D_bcsfromfiles_inflow.bc" /> // v = 0
<D VAR="w" FILE="ChanFlow_3DH1D_bcsfromfiles_inflow.bc" /> // w = 0
<N VAR="p" USERDEFINEDTYPE="H" VALUE="0" /> // High Order Pressure BC
</REGION>
<REGION REF="2">
<N VAR="u" VALUE="0" />
<N VAR="v" VALUE="0" />
<N VAR="w" VALUE="0" />
<D VAR="p" VALUE="0" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="0" />
<E VAR="v" VALUE="0" />
<E VAR="w" VALUE="0" />
<E VAR="p" VALUE="0" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="y*(1-y)" />
<E VAR="v" VALUE="0" />
<E VAR="w" VALUE="0" />
<E VAR="p" VALUE="-2*Kinvis*(x-1)" />
</FUNCTION>
</CONDITIONS>
<GEOMETRY DIM="2" SPACE="2">
<VERTEX>
<!-- Always must have four values per entry. -->
<V ID="0"> 0.0 0.0 0.0 </V>
<V ID="1"> 0.5 0.0 0.0 </V>
<V ID="2"> 1.0 0.0 0.0 </V>
<V ID="3"> 0.0 0.5 0.0 </V>
<V ID="4"> 0.5 0.5 0.0 </V>
<V ID="5"> 1.0 0.5 0.0 </V>
<V ID="6"> 0.0 1.0 0.0 </V>
<V ID="7"> 0.5 1.0 0.0 </V>
<V ID="8"> 1.0 1.0 0.0 </V>
</VERTEX>
<EDGE>
<E ID="0"> 0 1 </E>
<E ID="1"> 1 2 </E>
<E ID="2"> 0 3 </E>
<E ID="3"> 1 4 </E>
<E ID="4"> 2 5 </E>
<E ID="5"> 3 4 </E>
<E ID="6"> 4 5 </E>
<E ID="7"> 3 6 </E>
<E ID="8"> 4 7 </E>
<E ID="9"> 5 8 </E>
<E ID="10"> 6 7 </E>
<E ID="11"> 7 8 </E>
</EDGE>
<!-- Q - quads, T - triangles, S - segments, E - tet, P - pyramid, R - prism, H - hex -->
<!-- Only certain element types are appropriate for the given dimension (dim on mesh) -->
<!-- Can also use faces to define 3-D elements. Specify with F[1] for face 1, for example. -->
<ELEMENT>
<Q ID="0"> 0 3 5 2 </Q>
<Q ID="1"> 1 4 6 3 </Q>
<Q ID="2"> 5 8 10 7 </Q>
<Q ID="3"> 6 9 11 8 </Q>
</ELEMENT>
<COMPOSITE>
<C ID="0"> Q[0-3] </C>
<C ID="1"> E[0,1,10,11] </C> // Walls
<C ID="2"> E[2,7] </C> // Inflow
<C ID="3"> E[4,9] </C> // Outflow
</COMPOSITE>
<DOMAIN> C[0] </DOMAIN>
</GEOMETRY>
</NEKTAR>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>
<Metadata>
<Provenance>
<GitBranch>refs/heads/fix/Homo1D_ReadBCFromFile</GitBranch>
<GitSHA1>fda88fc8f3eb4071a6485fb79bce8e4ff7aa1225</GitSHA1>
<Hostname>dyn1134-150.insecure.ic.ac.uk</Hostname>
<NektarVersion>4.5.0</NektarVersion>
<Timestamp>09-May-2018 13:08:09</Timestamp>
</Provenance>
<ChkFileNum>2</ChkFileNum>
<Kinvis>1</Kinvis>
<SessionName0>ChanFlow_3DH1D_bcsfromfiles.xml</SessionName0>
<Time>1.0000000000000007</Time>
<TimeStep>0.001</TimeStep>
</Metadata>
<ELEMENTS FIELDS="u,v,w,p" SHAPE="Segment-HomogenousExp1D" BASIS="Modified_A,Fourier" HOMOGENEOUSLENGTHS="1" HOMOGENEOUSZIDS="0,1,2,3,4,5,6,7,8,9" NUMMODESPERDIR="UNIORDER:3,10" ID="2,7" COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYEAGF+xRaUIgwQaVhoGCPag0DOTYoNKEgMMeVJpc91TsQaUR/mNB4+M3l4EhoFRS9en1JCgf4b8zV+/+OP68EMpH+E/Y6dnHjUvSMfzrILyxxd4yFsN/UzhmHj8+yQ/Dvw6p5beYzdQw/New6rSG5n4zNP8xMFwHu6eUyHAbBaNgFAxGYAChHAqhdIoh65QjETx7cam/VqB+uiesyyb3l+97hUcLbS7eWelrd6LBhp+FOfGn7aw9n5+FWqa/3LgncTsbw+57E/bMuC+6w2HWRJuj7FK/7RbvtNnDtcgn++FqG/UiZn+pXKc91ute7vFoDNujG/d2rYBPpg1vlL4Q08rIPfb9jvrb1vXvOb3sVPODpc17/C5ypDN9YrERdBHnub1G3uZ96qsdjzks9xSLqS1MvJ5uc7HoWIbS4dU2zy7O3CC7e4VNz/32aGuFlj3Wj5p3/rI6uudreKpBmdTBPTB/SkDprsmbir3dhWxx+Rfmz9K2mHXqy3tsBD50ajrO7raB+dOeN4NH4vTsPTO4JrxfLj9jD8yfn9YfL7VQm2DDuzZJdc3JpTYwf86Ikz7dcKZ4j24r24/pnxbD/aly1OlR7vHEPXOTRG2EZaftgfnzyaafv7i3GNpMPb9TMkurEO7P6IXTN+heCLZpnHeFf2dwO9yfKwR/xpVbNu+Jdlwv4Miwfw8A7pH6bQAA</ELEMENTS>
</NEKTAR>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>
<Metadata>
<Provenance>
<GitBranch>refs/heads/fix/Homo1D_ReadBCFromFile</GitBranch>
<GitSHA1>fda88fc8f3eb4071a6485fb79bce8e4ff7aa1225</GitSHA1>
<Hostname>dyn1134-150.insecure.ic.ac.uk</Hostname>
<NektarVersion>4.5.0</NektarVersion>
<Timestamp>09-May-2018 13:07:54</Timestamp>
</Provenance>
<ChkFileNum>2</ChkFileNum>
<Kinvis>1</Kinvis>
<SessionName0>ChanFlow_3DH1D_bcsfromfiles.xml</SessionName0>
<Time>1.0000000000000007</Time>
<TimeStep>0.001</TimeStep>
</Metadata>
<ELEMENTS FIELDS="u,v,w,p" SHAPE="Segment-HomogenousExp1D" BASIS="Modified_A,Fourier" HOMOGENEOUSLENGTHS="1" HOMOGENEOUSZIDS="0,1,2,3,4,5,6,7,8,9" NUMMODESPERDIR="UNIORDER:3,10" ID="0-1,10-11" COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYBgFo2AUjIJRMApGwSgYBaNgFBALDCCUgyWY+mD/qzhZbNt/9b241F8rUD/dE9Zl47XG7Vbi63CbOdEHqyLCyvbwszAn/rSdtWch05uVHIfS9xyatfzt1dhWmxn3RXc4zJpo81jH+XdVMNOeRr4i5ac3JuxRL2L2l8p12pO6s9fzj6fknsWhZ3snRCXa8EbpCzGtjNzz+nXTz9cyH3ffmd3DIfotxcbvIkc60ycWm9y3q1PaAh5aW0isuLvJXGlPsZjawsTr6Tad589tkbz33Hr7PeE7/RNL9vTcb4+2VmjZc73hak3TBf09pdte8gcwTbaB+RPmH00P48ZeOT5bXP6F+RPGd2U0VNx4JRXuT5h4+j7B7CX7S+H+hIk37lscYdSrYAPzJ0z8a9+22PYQTRuYP2HiTnf3bdEMerQb5k+YuHHHHr980S+7Yf6Eif/8f4its0zSBuZPmHhBJe+UuFvaNhLQ+DWH+ltUJ+wI2xchnPFb2hazTn15j81lnjOh4ZoRNl1SMaU5+jV77HkzeCROz96z4GXWzd7X6Xv0I4OeCHh22nxaf7zUQm2CTb9h6kS2h+J7Nh37r5zysGfPjDjp0w1nivdcUBJZ26Qsvafz/9zTc2sbbFSOOj3KPZ64p/2sinVd9IfddjNsV80+l2PzZNPPX9xbDG1ctn67Lr/3gXXNfyt90WaHPdELp2/QvRBso5Waa/ntwOvdageefXIKi9+zQvBnXLll8564gBeznB7p70mKmzn9v9QkG5g/Yf5pcM1dLKEvgDN+Yf6E8VOVtgkLK6fB/QkPf5aJ8tpfSuH+hMfX5NVqD6/r28D8CRNn6mac9PO1nA3MnzDxW7znW89KM+6B+RMmrlxp6rJX6OlumD9h4vl3luhnvhWzgfkTJl71rTY68r2uDQACU0+9</ELEMENTS>
</NEKTAR>
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