Commit e73e7858 authored by Zhenguo Yan's avatar Zhenguo Yan

DiffusionIP symm TraceJump calculation

parent 9951a0a3
......@@ -991,7 +991,8 @@ namespace Nektar
}
void DisContField1D::v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight)
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp)
{
int cnt, n, e, npts;
// Fill boundary conditions into missing elements.
......@@ -1004,7 +1005,8 @@ namespace Nektar
{
id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
// Bwd[id] = m_bndCondExpansions[n]->GetPhys()[0]; //this is not getting the correct value?
weight[id] = m_BndCondBwdWeight[n];
weightave[id] = m_BndCondBwdWeight[n];
weightjmp[id] = 0.0;
cnt++;
}
else if (m_bndConditions[n]->GetBoundaryConditionType() ==
......@@ -1017,7 +1019,8 @@ namespace Nektar
"boundary condition");
id = m_trace->GetPhys_Offset(m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt));
// Bwd[id] = Fwd[id];
weight[id] = m_BndCondBwdWeight[n];
weightave[id] = m_BndCondBwdWeight[n];
weightjmp[id] = 0.0;
cnt++;
}
else if (m_bndConditions[n]->GetBoundaryConditionType() ==
......
......@@ -276,7 +276,8 @@ namespace Nektar
const Array<OneD, const NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd);
virtual void v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight);
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp);
virtual void v_GetFwdBwdTracePhys(
Array<OneD, NekDouble> &Fwd,
......
......@@ -1579,7 +1579,8 @@ namespace Nektar
* NOTE: periodic boundary is considered interior traces and is not treated here.
*/
void DisContField2D::v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight)
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp)
{
int cnt, n, e, npts;
......@@ -1598,7 +1599,8 @@ namespace Nektar
id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
id2 = m_trace->GetPhys_Offset(m_traceMap->
GetBndCondTraceToGlobalTraceMap(cnt+e));
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weight[id2],1);
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weightave[id2],1);
Vmath::Fill(npts,0.0, &weightjmp[id2],1);
}
......@@ -1615,7 +1617,8 @@ namespace Nektar
id1 = m_bndCondExpansions[n]->GetPhys_Offset(e);
id2 = m_trace->GetPhys_Offset(m_traceMap->
GetBndCondTraceToGlobalTraceMap(cnt+e));
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weight[id2],1);
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weightave[id2],1);
Vmath::Fill(npts,0.0, &weightjmp[id2],1);
}
cnt += e;
......
......@@ -321,7 +321,8 @@ namespace Nektar
const Array<OneD, const NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd);
virtual void v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight);
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp);
virtual void v_GetFwdBwdTracePhys(
Array<OneD, NekDouble> &Fwd,
......
......@@ -1874,7 +1874,8 @@ using namespace std;
* NOTE: periodic boundary is considered interior traces and is not treated here.
*/
void DisContField3D::v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight)
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp)
{
int cnt, n, e, npts, phys_offset;
// Fill boundary conditions into missing elements
......@@ -1896,7 +1897,8 @@ using namespace std;
// &(m_bndCondExpansions[n]->GetPhys())[id1], 1,
// &Bwd[id2], 1);
// Vmath::Vcopy(npts,&Fwd[id2],1,&Bwd[id2],1);
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weight[id2],1);
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weightave[id2],1);
Vmath::Fill(npts,0.0, &weightjmp[id2],1);
}
cnt += e;
......@@ -1920,7 +1922,8 @@ using namespace std;
//Neumann " "boundary condition");
// Vmath::Vcopy(npts,&Fwd[id2],1,&Bwd[id2],1);
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weight[id2],1);
Vmath::Fill(npts,m_BndCondBwdWeight[n], &weightave[id2],1);
Vmath::Fill(npts,0.0, &weightjmp[id2],1);
}
cnt += e;
......
......@@ -275,7 +275,8 @@ namespace Nektar
const Array<OneD, const NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd);
virtual void v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight);
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp);
virtual void v_GetFwdBwdTracePhys(
Array<OneD, NekDouble> &Fwd,
......
......@@ -2502,15 +2502,18 @@ namespace Nektar
}
}
void ExpList::GetBwdWeight(Array<OneD, NekDouble> &weight)
void ExpList::GetBwdWeight(
Array<OneD, NekDouble> &weightAver,
Array<OneD, NekDouble> &weightJump)
{
int nTracePts = weight.num_elements();
int nTracePts = weightAver.num_elements();
// average for interior traces
for(int i = 0; i < nTracePts; ++i)
{
weight[i] = 0.5;
weightAver[i] = 0.5;
weightJump[i] = 1.0;
}
FillBwdWITHBwdWeight(weight);
FillBwdWITHBwdWeight(weightAver,weightJump);
}
const Array<OneD,const std::shared_ptr<ExpList> >
......@@ -3242,7 +3245,8 @@ namespace Nektar
}
void ExpList::v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight)
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp)
{
ASSERTL0(false,"v_FillBwdWITHBwdWeight not defined");
}
......
......@@ -769,7 +769,9 @@ namespace Nektar
inline void GetElmtNormalLength(Array<OneD, NekDouble> &lengths);
void GetBwdWeight(Array<OneD, NekDouble> &weight);
void GetBwdWeight(
Array<OneD, NekDouble> &weightAver,
Array<OneD, NekDouble> &weightJump);
inline void AddTraceIntegral(
const Array<OneD, const NekDouble> &Fx,
......@@ -831,7 +833,8 @@ namespace Nektar
Array<OneD, NekDouble> &Bwd);
inline void FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight);
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp);
inline void PeriodicBwdCopy(
const Array<OneD, const NekDouble> &Fwd,
......@@ -1320,7 +1323,8 @@ namespace Nektar
Array<OneD, NekDouble> &Bwd);
virtual void v_FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &weight);
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp);
virtual void v_PeriodicBwdCopy(
const Array<OneD, const NekDouble> &Fwd,
......@@ -2453,9 +2457,10 @@ namespace Nektar
}
inline void ExpList::FillBwdWITHBwdWeight(
Array<OneD, NekDouble> &Bwd)
Array<OneD, NekDouble> &weightave,
Array<OneD, NekDouble> &weightjmp)
{
v_FillBwdWITHBwdWeight(Bwd);
v_FillBwdWITHBwdWeight(weightave,weightjmp);
}
inline void ExpList::PeriodicBwdCopy(
......
......@@ -100,13 +100,15 @@ namespace Nektar
Vmath::Sdiv(nTracePts,1.0,m_traceNormDirctnElmtLength,1,m_traceNormDirctnElmtLengthRecip,1);
// }
m_tracBwdWeight = Array<OneD, NekDouble> (nTracePts,0.0);
pFields[0]->GetBwdWeight(m_tracBwdWeight);
m_tracBwdWeightAver = Array<OneD, NekDouble> (nTracePts,0.0);
m_tracBwdWeightJump = Array<OneD, NekDouble> (nTracePts,0.0);
pFields[0]->GetBwdWeight(m_tracBwdWeightAver,m_tracBwdWeightJump);
Array<OneD, NekDouble> tmpBwdWeight(nTracePts,0.0);
Array<OneD, NekDouble> tmpBwdWeightJump(nTracePts,0.0);
for(int i =1; i<nVariable;i++)
{
pFields[i]->GetBwdWeight(tmpBwdWeight);
Vmath::Vsub(nTracePts,tmpBwdWeight,1,m_tracBwdWeight,1,tmpBwdWeight,1);
pFields[i]->GetBwdWeight(tmpBwdWeight,tmpBwdWeightJump);
Vmath::Vsub(nTracePts,tmpBwdWeight,1,m_tracBwdWeightAver,1,tmpBwdWeight,1);
Vmath::Vabs(nTracePts,tmpBwdWeight,1,tmpBwdWeight,1);
NekDouble norm = 0.0;
for(int j = 0; j<nTracePts; j++)
......@@ -565,10 +567,22 @@ namespace Nektar
// note: here the jump is 2.0*(aver-vFwd)
// because Viscous wall use a symmetry value as the Bwd, not the target one
Array<OneD, NekDouble> tmpF (npnts,0.0);
Array<OneD, NekDouble> tmpB (npnts,0.0);
Array<OneD, NekDouble> Fweight (npnts,2.0);
Array<OneD, NekDouble> Bweight;
Bweight = m_tracBwdWeightJump;
Vmath::Vsub(npnts,Fweight,1,Bweight,1,Fweight,1);
for (int i = 0; i < nConvectiveFields; ++i)
{
Vmath::Vsub(npnts,aver[i],1,vFwd[i],1,jump[i],1);
Vmath::Smul(npnts,2.0,jump[i],1,jump[i],1);
Vmath::Vsub(npnts,aver[i],1,vFwd[i],1,tmpF,1);
Vmath::Vsub(npnts,vBwd[i],1,aver[i],1,tmpB,1);
Vmath::Vmul(npnts,tmpF,1,Fweight,1,tmpF,1);
Vmath::Vmul(npnts,tmpB,1,Bweight,1,tmpB,1);
Vmath::Vadd(npnts,tmpF,1,tmpB,1,jump[i],1);
}
}
......@@ -586,7 +600,7 @@ namespace Nektar
Array<OneD, NekDouble> Fweight (npnts,1.0);
Array<OneD, NekDouble> Bweight;
Bweight = m_tracBwdWeight;
Bweight = m_tracBwdWeightAver;
Vmath::Vsub(npnts,Fweight,1,Bweight,1,Fweight,1);
......
......@@ -71,7 +71,8 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > m_traceNormals;
Array<OneD, Array<OneD, NekDouble> > m_traceAver;
Array<OneD, Array<OneD, NekDouble> > m_traceJump;
Array<OneD, NekDouble> m_tracBwdWeight;
Array<OneD, NekDouble> m_tracBwdWeightAver;
Array<OneD, NekDouble> m_tracBwdWeightJump;
Array<OneD, NekDouble> m_traceNormDirctnElmtLength;
Array<OneD, NekDouble> m_traceNormDirctnElmtLengthRecip;
LibUtilities::SessionReaderSharedPtr m_session;
......@@ -198,7 +199,7 @@ namespace Nektar
const Array<OneD, const Array<OneD, NekDouble> > &vBwd,
Array<OneD, Array<OneD, NekDouble> > &aver,
Array<OneD, Array<OneD, NekDouble> > &jump);
virtual const Array<OneD, const Array<OneD, NekDouble> > &v_GetTraceNormal()
{
return m_traceNormals;
......
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