Commit 79863088 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'fix/linearised-wavespace' into 'master'

Fix/linearised wavespace

This MR fixes the linearised stability solver when using the homogeneous 1D code.

See merge request !454
parents 5b67aa99 5d601d93
......@@ -389,6 +389,7 @@ void LinearisedAdvection::v_Advect(
//3D
case 3:
{
grad1 = Array<OneD, NekDouble> (nPointsTot);
grad2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_u1 = Array<OneD, NekDouble> (nPointsTot);
......@@ -399,9 +400,16 @@ void LinearisedAdvection::v_Advect(
grad_base_v2 = Array<OneD, NekDouble> (nPointsTot);
grad_base_w2 = Array<OneD, NekDouble> (nPointsTot);
fields[0]->PhysDeriv(m_baseflow[0], grad_base_u0, grad_base_u1,grad_base_u2);
fields[0]->PhysDeriv(m_baseflow[1], grad_base_v0, grad_base_v1,grad_base_v2);
fields[0]->PhysDeriv(m_baseflow[2], grad_base_w0, grad_base_w1, grad_base_w2);
// Differentiate base flow in physical space
bool oldwavespace = fields[0]->GetWaveSpace();
fields[0]->SetWaveSpace(false);
fields[0]->PhysDeriv(m_baseflow[0],
grad_base_u0, grad_base_u1, grad_base_u2);
fields[0]->PhysDeriv(m_baseflow[1],
grad_base_v0, grad_base_v1, grad_base_v2);
fields[0]->PhysDeriv(m_baseflow[2],
grad_base_w0, grad_base_w1, grad_base_w2);
fields[0]->SetWaveSpace(oldwavespace);
//HalfMode has W(x,y,t)=0
if(m_HalfMode || m_SingleMode)
......@@ -415,6 +423,12 @@ void LinearisedAdvection::v_Advect(
}
fields[0]->PhysDeriv(inarray[n], grad0, grad1, grad2);
if (oldwavespace)
{
fields[0]->HomogeneousBwdTrans(grad0, grad0);
fields[0]->HomogeneousBwdTrans(grad1, grad1);
fields[0]->HomogeneousBwdTrans(grad2, grad2);
}
switch (n)
{
......@@ -447,15 +461,18 @@ void LinearisedAdvection::v_Advect(
else
{
//Evaluate U du'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
outarray[n], 1);
//Evaluate U du'/dx+ V du'/dy
Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
outarray[n], 1, outarray[n], 1);
//Evaluate (U du'/dx+ V du'/dy)+u' dU/dx
Vmath::Vvtvp(nPointsTot,grad_base_u0,1,advVel[0],1,outarray[n],1,outarray[n],1);
//Evaluate (U du'/dx+ V du'/dy +u' dU/dx)+v' dU/dy
Vmath::Vvtvp(nPointsTot,grad_base_u1,1,advVel[1],1,outarray[n],1,outarray[n],1);
//Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy) + W du'/dz
Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
outarray[n], 1, outarray[n], 1);
//Evaluate (U du'/dx+ V du'/dy +u' dU/dx +v' dU/dy + W du'/dz)+ w' dU/dz
Vmath::Vvtvp(nPointsTot,grad_base_u2,1,advVel[2],1,outarray[n],1,outarray[n],1);
}
......@@ -486,15 +503,18 @@ void LinearisedAdvection::v_Advect(
else
{
//Evaluate U dv'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
outarray[n], 1);
//Evaluate U dv'/dx+ V dv'/dy
Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
outarray[n], 1, outarray[n], 1);
//Evaluate (U dv'/dx+ V dv'/dy)+u' dV/dx
Vmath::Vvtvp(nPointsTot,grad_base_v0,1,advVel[0],1,outarray[n],1,outarray[n],1);
//Evaluate (U du'/dx+ V du'/dy +u' dV/dx)+v' dV/dy
Vmath::Vvtvp(nPointsTot,grad_base_v1,1,advVel[1],1,outarray[n],1,outarray[n],1);
//Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy) + W du'/dz
Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
outarray[n], 1, outarray[n], 1);
//Evaluate (U du'/dx+ V dv'/dy +u' dV/dx +v' dV/dy + W dv'/dz)+ w' dV/dz
Vmath::Vvtvp(nPointsTot,grad_base_v2,1,advVel[2],1,outarray[n],1,outarray[n],1);
}
......@@ -526,21 +546,30 @@ void LinearisedAdvection::v_Advect(
else
{
//Evaluate U dw'/dx
Vmath::Vmul (nPointsTot,grad0,1,m_baseflow[0],1,outarray[n],1);
Vmath::Vmul (nPointsTot, grad0, 1, m_baseflow[0], 1,
outarray[n], 1);
//Evaluate U dw'/dx+ V dw'/dx
Vmath::Vvtvp(nPointsTot,grad1,1,m_baseflow[1],1,outarray[n],1,outarray[n],1);
Vmath::Vvtvp(nPointsTot, grad1, 1, m_baseflow[1], 1,
outarray[n], 1, outarray[n], 1);
//Evaluate (U dw'/dx+ V dw'/dx)+u' dW/dx
Vmath::Vvtvp(nPointsTot,grad_base_w0,1,advVel[0],1,outarray[n],1,outarray[n],1);
//Evaluate (U dw'/dx+ V dw'/dx +w' dW/dx)+v' dW/dy
Vmath::Vvtvp(nPointsTot,grad_base_w1,1,advVel[1],1,outarray[n],1,outarray[n],1);
//Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy) + W dw'/dz
Vmath::Vvtvp(nPointsTot,grad2,1,m_baseflow[2],1,outarray[n],1,outarray[n],1);
Vmath::Vvtvp(nPointsTot, grad2, 1, m_baseflow[2], 1,
outarray[n], 1, outarray[n], 1);
//Evaluate (U dw'/dx+ V dw'/dx +u' dW/dx +v' dW/dy + W dw'/dz)+ w' dW/dz
Vmath::Vvtvp(nPointsTot,grad_base_w2,1,advVel[2],1,outarray[n],1,outarray[n],1);
}
break;
}
if (oldwavespace)
{
fields[0]->HomogeneousFwdTrans(outarray[n],outarray[n]);
}
break;
}
default:
ASSERTL0(false,"dimension unknown");
}
......@@ -672,8 +701,10 @@ void LinearisedAdvection::ImportFldBase(std::string pInfile,
}
else
{
bool oldwavespace = pFields[j]->GetWaveSpace();
pFields[j]->SetWaveSpace(false);
pFields[j]->BwdTrans(tmp_coeff, m_baseflow[j]);
pFields[j]->SetWaveSpace(oldwavespace);
}
}
......
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