Commit 037e5873 authored by Ale's avatar Ale
Browse files

new AeroForces into wavyforce branch

parent 2b50ea76
......@@ -65,9 +65,9 @@ namespace Nektar
m_outputFile = pParams.find("OutputFile")->second;
}
if (!(m_outputFile.length() >= 4
&& m_outputFile.substr(m_outputFile.length() - 4) == ".frc"))
&& m_outputFile.substr(m_outputFile.length() - 4) == ".fce"))
{
m_outputFile += ".frc";
m_outputFile += ".fce";
}
if (pParams.find("OutputFrequency") == pParams.end())
......@@ -108,12 +108,6 @@ namespace Nektar
"Missing parameter 'Boundary'.");
m_BoundaryString = pParams.find("Boundary")->second;
}
m_wavygeometry = false;
if(m_session->DefinesFunction("WavyGeometry"))
{
m_wavygeometry = true;
}
}
......@@ -175,27 +169,15 @@ namespace Nektar
}
}
// Fields are in physical space
for (int i = 0; i < pFields.num_elements(); ++i)
{
if (m_isHomogeneous1D)
{
pFields[i]->SetWaveSpace(false);
}
pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
pFields[i]->UpdatePhys());
pFields[i]->SetPhysState(true);
pFields[i]->PutPhysInToElmtExp();
}
LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
if (vComm->GetRank() == 0)
{
// Open output stream
m_outputStream.open(m_outputFile.c_str());
m_outputStream << "# Time, (drag-press , drag-visc) "
"Total Drag, (lift-press ,lift-visc) Total Lift:" ;
m_outputStream << "# Time\t Fx (press)\t Fx (visc)\t Fx (tot)\t"
" Fy (press)\t Fy (visc)\t Fy (tot)\t"
" Fz (press)\t Fz (visc)\t Fz (tot)\t";
m_outputStream << endl;
}
......@@ -241,24 +223,32 @@ namespace Nektar
Array<OneD, NekDouble> values;
LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
NekDouble D,L,D_p,D_t,L_p,L_t;
D=0.0;
L=0.0;
D_p=0.0;
D_t=0.0;
L_p=0.0;
L_t=0.0;
NekDouble Fx,Fy,Fz,Fxp,Fxv,Fyp,Fyv,Fzp,Fzv;
Fxp = 0.0; // x-component of the force due to pressure difference
Fxv = 0.0; // x-component of the force due to viscous stress
Fx = 0.0; // x-component of the force (total) Fx = Fxp + Fxv (Drag)
Fyp = 0.0; // y-component of the force due to pressure difference
Fyv = 0.0; // y-component of the force due to viscous stress
Fy = 0.0; // y-component of the force (total) Fy = Fyp + Fyv (Lift)
Fzp = 0.0; // z-component of the force due to pressure difference
Fzv = 0.0; // z-component of the force due to viscous stress
Fz = 0.0; // z-component of the force (total) Fz = Fzp + Fzv (Side)
NekDouble rho = (m_session->DefinesParameter("rho"))
? (m_session->GetParameter("rho"))
: 1;
NekDouble mu = rho*m_session->GetParameter("Kinvis");
for(int i = 0; i < pFields.num_elements(); ++i)
{
pFields[i]->SetWaveSpace(false);
pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
pFields[i]->UpdatePhys());
pFields[i]->SetPhysState(true);
pFields[i]->PutPhysInToElmtExp();
}
// Homogeneous 1D case Compute forces on all WALL boundaries
......@@ -267,7 +257,8 @@ namespace Nektar
{
if(vComm->GetColumnComm()->GetRank() == 0)
{
pFields[0]->GetPlane(0)->GetBoundaryToElmtMap(BoundarytoElmtID,BoundarytoTraceID);
pFields[0]->GetPlane(0)->GetBoundaryToElmtMap(
BoundarytoElmtID,BoundarytoTraceID);
BndExp = pFields[0]->GetPlane(0)->GetBndCondExpansions();
StdRegions::StdExpansion1DSharedPtr bc;
......@@ -283,7 +274,7 @@ namespace Nektar
elmt = pFields[0]->GetPlane(0)->GetExp(elmtid);
nq = elmt->GetTotPoints();
offset = pFields[0]->GetPlane(0)->GetPhys_Offset(elmtid);
// Initialise local arrays for the velocity
// gradients size of total number of quadrature
// points for each element (hence local).
......@@ -339,9 +330,9 @@ namespace Nektar
for(int j = 0; j < dim; ++j)
{
elmt->GetEdgePhysVals(boundary,bc,gradU[j],
fgradU[j]);
fgradU[j]);
elmt->GetEdgePhysVals(boundary,bc,gradV[j],
fgradV[j]);
fgradV[j]);
}
//normals of the element
......@@ -398,19 +389,16 @@ namespace Nektar
// boundaries
Vmath::Vvtvp(nbc,Pb,1,normals[0],1,
drag_p,1,drag_p, 1);
drag_p,1,drag_p, 1);
Vmath::Vvtvp(nbc,Pb,1,normals[1],1,
lift_p,1,lift_p,1);
lift_p,1,lift_p,1);
//integration over the boundary
D_t += bc->Integral(drag_t);
L_t += bc->Integral(lift_t);
D_p += bc->Integral(drag_p);
L_p += bc->Integral(lift_p);
Fxv += bc->Integral(drag_t);
Fyv += bc->Integral(lift_t);
D=D_t+D_p;
L=L_t+L_p;
Fxp += bc->Integral(drag_p);
Fyp += bc->Integral(lift_p);
}
}
else
......@@ -419,7 +407,15 @@ namespace Nektar
}
}
}
}
for(int i = 0; i < pFields.num_elements(); ++i)
{
pFields[i]->SetWaveSpace(true);
pFields[i]->BwdTrans(pFields[i]->GetCoeffs(),
pFields[i]->UpdatePhys());
pFields[i]->SetPhysState(false);
}
}
//3D WALL case
else if(dim==3 && !m_isHomogeneous1D)
{
......@@ -429,7 +425,7 @@ namespace Nektar
pFields[0]->GetBoundaryToElmtMap(BoundarytoElmtID,
BoundarytoTraceID);
BndExp = pFields[0]->GetBndCondExpansions();
LocalRegions:: Expansion2DSharedPtr bc;
LocalRegions::Expansion2DSharedPtr bc;
// loop over the types of boundary conditions
for(cnt = n = 0; n < BndExp.num_elements(); ++n)
......@@ -488,8 +484,10 @@ namespace Nektar
Array<OneD, NekDouble> drag_t(nbc);
Array<OneD, NekDouble> lift_t(nbc);
Array<OneD, NekDouble> side_t(nbc);
Array<OneD, NekDouble> drag_p(nbc);
Array<OneD, NekDouble> lift_p(nbc);
Array<OneD, NekDouble> side_p(nbc);
Array<OneD, NekDouble> temp(nbc);
Array<OneD, NekDouble> temp2(nbc);
......@@ -567,9 +565,30 @@ namespace Nektar
Vmath::Vmul(nbc,fgradV[1],1,normals[1],1,temp2,1);
Vmath::Smul(nbc,-0.5,fgradV[1],1,fgradV[1],1);
Vmath::Vadd(nbc,temp2,1,lift_t,1,lift_t,1);
Vmath::Smul(nbc,mu,lift_t,1,lift_t,1);
//zero temporary storage vector
Vmath::Zero(nbc,temp,0);
Vmath::Zero(nbc,temp2,0);
//b) SIDE TERMS
//-rho*kinvis*
// (2*dv/dy*nx+(du/dy+dv/dx)*nx+(dv/dz+dw/dy)*nz)
Vmath::Vadd(nbc,fgradV[2],1,fgradW[1],1,temp,1);
Vmath::Neg(nbc,temp,1);
Vmath::Vmul(nbc,temp,1,normals[1],1,temp,1);
Vmath::Vadd(nbc,fgradU[2],1,fgradW[1],1,side_t,1);
Vmath::Neg(nbc,side_t,1);
Vmath::Vmul(nbc,side_t,1,normals[0],1,side_t,1);
Vmath::Smul(nbc,-2.0,fgradW[2],1,fgradW[2],1);
Vmath::Vmul(nbc,fgradW[2],1,normals[2],1,temp2,1);
Vmath::Smul(nbc,-0.5,fgradW[2],1,fgradW[2],1);
Vmath::Vadd(nbc,temp2,1,side_t,1,side_t,1);
Vmath::Smul(nbc,mu,side_t,1,side_t,1);
// Compute normal tractive forces on all WALL
......@@ -578,16 +597,17 @@ namespace Nektar
drag_p,1,drag_p,1);
Vmath::Vvtvp(nbc,Pb,1,normals[1],1,
lift_p,1,lift_p,1);
Vmath::Vvtvp(nbc,Pb,1,normals[2],1,
side_p,1,side_p,1);
//integration over the boundary
D_t += bc->Expansion::Integral(drag_t);
L_t += bc->Expansion::Integral(lift_t);
Fxv += bc->Expansion::Integral(drag_t);
Fyv += bc->Expansion::Integral(lift_t);
Fzv += bc->Expansion::Integral(side_t);
D_p += bc->Expansion::Integral(drag_p);
L_p += bc->Expansion::Integral(lift_p);
D_t +=D_t+D_p;
L=L_t+L_p;
Fxp += bc->Expansion::Integral(drag_p);
Fyp += bc->Expansion::Integral(lift_p);
Fzp += bc->Expansion::Integral(side_p);
}
}
else
......@@ -632,8 +652,8 @@ namespace Nektar
elmt->PhysDeriv(U,gradU[0],gradU[1]);
elmt->PhysDeriv(V,gradV[0],gradV[1]);
bc = boost::dynamic_pointer_cast<StdRegions
::StdExpansion1D> (BndExp[n]->GetExp(i));
bc = boost::dynamic_pointer_cast<LocalRegions
::Expansion1D> (BndExp[n]->GetExp(i));
int nbc = bc->GetTotPoints();
Array<OneD, NekDouble> Pb(nbc);
......@@ -688,14 +708,11 @@ namespace Nektar
Vmath::Vvtvp(nbc,Pb,1,normals[1],1,
lift_p,1,lift_p,1);
D_p=D_p+bc->Integral(drag_p);
L_p=L_p+bc->Integral(lift_p);
D_t=D_t+bc->Integral(drag_t);
L_t=L_t+bc->Integral(lift_t);
Fxp += bc->Integral(drag_p);
Fyp += bc->Integral(lift_p);
D=D_p+D_t;
L=L_p+L_t;
Fxv += bc->Integral(drag_t);
Fyp += bc->Integral(lift_t);
}
}
else
......@@ -707,32 +724,44 @@ namespace Nektar
}
vComm->AllReduce(D_p, LibUtilities::ReduceSum);
vComm->AllReduce(D_t, LibUtilities::ReduceSum);
vComm->AllReduce(D, LibUtilities::ReduceSum);
vComm->AllReduce(L_p, LibUtilities::ReduceSum);
vComm->AllReduce(L_t, LibUtilities::ReduceSum);
vComm->AllReduce(L, LibUtilities::ReduceSum);
vComm->AllReduce(Fxp, LibUtilities::ReduceSum);
vComm->AllReduce(Fxv, LibUtilities::ReduceSum);
Fx = Fxp + Fxv;
vComm->AllReduce(Fyp, LibUtilities::ReduceSum);
vComm->AllReduce(Fyv, LibUtilities::ReduceSum);
Fy = Fyp + Fyv;
vComm->AllReduce(Fzp, LibUtilities::ReduceSum);
vComm->AllReduce(Fzv, LibUtilities::ReduceSum);
Fz = Fzp + Fzv;
if (vComm->GetRank() == 0)
{
m_outputStream.width(8);
m_outputStream << setprecision(6) << time;
m_outputStream.width(25);
m_outputStream << setprecision(8) << D_p;
m_outputStream << setprecision(8) << Fxp;
m_outputStream.width(25);
m_outputStream << setprecision(8) << D_t;
m_outputStream << setprecision(8) << Fxv;
m_outputStream.width(25);
m_outputStream << setprecision(16) << D;
m_outputStream << setprecision(16) << Fx;
m_outputStream.width(25);
m_outputStream << setprecision(8) << L_p;
m_outputStream << setprecision(8) << Fyp;
m_outputStream.width(25);
m_outputStream << setprecision(8) << Fyv;
m_outputStream.width(25);
m_outputStream << setprecision(8) << Fy;
m_outputStream.width(25);
m_outputStream << setprecision(8) << Fzp;
m_outputStream.width(25);
m_outputStream << setprecision(8) << L_t;
m_outputStream << setprecision(8) << Fzv;
m_outputStream.width(25);
m_outputStream << setprecision(8) << L;
m_outputStream << setprecision(8) << Fz;
m_outputStream << endl;
}
......
......@@ -91,7 +91,6 @@ namespace Nektar
std::string m_BoundaryString;
/// number of planes for homogeneous1D expansion
int m_nplanes;
bool m_wavygeometry;
};
}
}
......
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