Commit 2484daf5 authored by Daniele de Grazia's avatar Daniele de Grazia
Browse files

Change in InitObject

parent 9d39cc34
......@@ -136,10 +136,68 @@ namespace Nektar
m_planeAdv->SetRiemannSolver(m_riemann);
m_planeAdv->SetFluxVectorVec(m_fluxVector);
fluxvector = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
fluxvector[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
}
fluxvector_homo =
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
(num_planes);
outarray_homo = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
fields_plane =
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
(num_planes);
inarray_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
outarray_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
advVel_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
for (i = 0; i < num_planes; ++i)
{
fields_plane[i] = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
inarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
outarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
fluxvector_homo[i] =
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(nConvectiveFields);
advVel_plane[i] = Array<OneD, Array<OneD, NekDouble> >(3);
for (j = 0; j < nConvectiveFields; j ++)
{
inarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_homo[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
fluxvector_homo[i][j] = Array<OneD, Array<OneD, NekDouble> >
(3);
for (int k = 0; k < 3 ; ++k)
{
fluxvector_homo[i][j][k] =
Array<OneD, NekDouble>(nPointsTot_plane, 0.0);
advVel_plane[i][k] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
}
}
}
}
/**
* @brief Compute the advection term at each time-step using the Flux
* Reconstruction approach (FR) looping on the planes.
......@@ -159,48 +217,13 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int i, j, k;
int nVel = advVel.num_elements();
Array<OneD, Array<OneD, NekDouble> > fluxvector(nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
fluxvector[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
}
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
fluxvector_homo(num_planes);
Array<OneD, Array<OneD, NekDouble> >
outarray_homo(nConvectiveFields);
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
fields_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
inarray_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
outarray_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
advVel_plane(num_planes);
for (i = 0; i < num_planes; ++i)
{
fields_plane[i] = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
inarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
outarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
fluxvector_homo[i] =
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(nConvectiveFields);
advVel_plane[i] = Array<OneD, Array<OneD, NekDouble> >(nVel);
for (j = 0; j < nConvectiveFields; j ++)
{
fields_plane[i][j]= fields[j]->GetPlane(i);
inarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
Vmath::Vcopy(nPointsTot_plane,
&inarray[j][i * nPointsTot_plane], 1,
......@@ -209,10 +232,9 @@ namespace Nektar
for (j = 0; j < nVel; j ++)
{
advVel_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
if (advVel[j].num_elements() != 0 )
{ Vmath::Vcopy(nPointsTot_plane,
{
Vmath::Vcopy(nPointsTot_plane,
&advVel[j][i * nPointsTot_plane], 1,
&advVel_plane[i][j][0], 1);
}
......@@ -224,15 +246,6 @@ namespace Nektar
for (j = 0; j < nConvectiveFields; j ++)
{
fluxvector_homo[i][j] =
Array<OneD, Array<OneD, NekDouble> >(nVel);
for (int k = 0; k < nVel ; ++k)
{
fluxvector_homo[i][j][k] =
Array<OneD, NekDouble>(nPointsTot_plane, 0.0);
}
Vmath::Vcopy(nPointsTot_plane,
&outarray_plane[i][j][0], 1,
&outarray[j][i * nPointsTot_plane], 1);
......@@ -250,8 +263,6 @@ namespace Nektar
for (i = 0; i < nConvectiveFields; ++i)
{
outarray_homo[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
fields[0]->PhysDeriv(2, fluxvector[i], outarray_homo[i]);
Vmath::Vadd(nPointsTot,
......
......@@ -65,8 +65,18 @@ namespace Nektar
int nPointsTot_plane;
int nCoeffs_plane;
int num_planes;
int i, j, k;
Array<OneD, unsigned int> planes;
Array<OneD, unsigned int> planes;
Array<OneD, Array<OneD, NekDouble> > fluxvector;
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
fluxvector_homo;
Array<OneD, Array<OneD, NekDouble> > outarray_homo;
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
fields_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > inarray_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > outarray_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > advVel_plane;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
......
......@@ -88,7 +88,44 @@ namespace Nektar
nPointsTot_plane = nPointsTot/num_planes;
nCoeffs_plane = nCoeffs/num_planes;
fields_plane =
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
(num_planes);
inarray_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
outarray_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
for (i = 0; i < num_planes; ++i)
{
fields_plane[i] = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
inarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
outarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
inarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
}
}
flux = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
flux_homo =
Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
outarray_z =
Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
for (j = 0; j < nConvectiveFields; j++)
{
flux[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
flux_homo[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
outarray_z[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
}
}
/**
......@@ -102,31 +139,11 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int i, j, k;
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
fields_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
inarray_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
outarray_plane(num_planes);
for (i = 0; i < num_planes; ++i)
{
fields_plane[i] = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
inarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
outarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
fields_plane[i][j]= fields[j]->GetPlane(i);
inarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
Vmath::Vcopy(nPointsTot_plane,
&inarray[j][i * nPointsTot_plane], 1,
......@@ -146,19 +163,11 @@ namespace Nektar
}
}
Array<OneD, Array<OneD, NekDouble> > flux(nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > flux_homo(nConvectiveFields);
Array<OneD, Array<OneD, NekDouble> > outarray_z(nConvectiveFields);
NekDouble beta;
int Homolen = fields[0]->GetHomoLen();
for (j = 0; j < nConvectiveFields; j++)
{
flux[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
flux_homo[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
outarray_z[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
// Transform flux in Fourier space
fields[0]->HomogeneousFwdTrans(inarray[j], flux[j]);
}
......
......@@ -64,8 +64,17 @@ namespace Nektar
int num_planes;
int nPointsTot_plane;
int nCoeffs_plane;
int i, j, k;
Array<OneD, unsigned int> planes;
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
fields_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > inarray_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > outarray_plane;
Array<OneD, Array<OneD, NekDouble> > flux;
Array<OneD, Array<OneD, NekDouble> > flux_homo;
Array<OneD, Array<OneD, NekDouble> > outarray_z;
Array<OneD, unsigned int> planes;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
......
......@@ -73,11 +73,12 @@ namespace Nektar
Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
{
int nConvectiveFields = pFields.num_elements();
nVariables = nConvectiveFields - 1;
Array<OneD, MultiRegions::ExpListSharedPtr>
pFields_plane0(nConvectiveFields);
for (int i = 0; i < nConvectiveFields; ++i)
for (i = 0; i < nConvectiveFields; ++i)
{
pFields_plane0[i] = pFields[i]->GetPlane(0);
}
......@@ -93,95 +94,68 @@ namespace Nektar
nPointsTot_plane = nPointsTot/num_planes;
nCoeffs_plane = nCoeffs/num_planes;
}
/**
* @brief Calculate WeakDG Diffusion for the Navier-Stokes (NS) equations
* using an LDG interface flux and the the flux in the third direction.
*
* The equations that need a diffusion operator are those related
* with the velocities and with the energy.
*/
void DiffusionLDGNS3DHomogeneous1D ::v_Diffuse(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
int i, j, k;
int nVariables = inarray.num_elements();
DiffusionLDGNSSharedPtr diffLDGNS = boost::dynamic_pointer_cast<
DiffusionLDGNS>(m_planeDiff);
Array<OneD, Array<OneD, NekDouble> > fluxvector(nConvectiveFields);
diffLDGNS =
boost::dynamic_pointer_cast<DiffusionLDGNS>(m_planeDiff);
fluxvector =
Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
fluxvector[j] = Array<OneD, NekDouble>(nPointsTot);
}
Array<OneD, Array<OneD, NekDouble> > derivatives01_homo(nVariables);
derivatives01_homo =
Array<OneD, Array<OneD, NekDouble> > (nVariables);
for (i = 0; i < nVariables; ++i)
{
derivatives01_homo[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
fields[0]->PhysDeriv(2, inarray[i], derivatives01_homo[i]);
}
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
fields_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
inarray_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
derivatives01_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
outarray_plane(num_planes);
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
fluxvector_homo(num_planes);
m_planeDiff->SetRiemannSolver(m_riemann);
m_planeDiff->SetFluxVectorVecNS(m_fluxVectorNS);
fields_plane =
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
(num_planes);
inarray_plane =
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (num_planes);
outarray_plane =
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (num_planes);
derivatives01_plane =
Array<OneD, Array<OneD, Array<OneD, NekDouble> > > (num_planes);
fluxvector_homo =
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
(num_planes);
for (i = 0; i < num_planes; ++i)
{
fields_plane[i] = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
(nConvectiveFields);
inarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nVariables);
(nVariables);
derivatives01_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nVariables);
(nVariables);
outarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
(nConvectiveFields);
fluxvector_homo[i] =
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >(spaceDim);
for (j = 0; j < nConvectiveFields; j ++)
{
fields_plane[i][j]= fields[j]->GetPlane(i);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
}
for (j = 0; j < nVariables; j ++)
{
inarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
Vmath::Vcopy(nPointsTot_plane,
&inarray[j][i * nPointsTot_plane], 1,
&inarray_plane[i][j][0], 1);
(nPointsTot_plane, 0.0);
derivatives01_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
Vmath::Vcopy(nPointsTot_plane,
&derivatives01_homo[j][i * nPointsTot_plane], 1,
&derivatives01_plane[i][j][0], 1);
(nPointsTot_plane, 0.0);
}
for (j = 0; j < spaceDim ; ++j)
{
fluxvector_homo[i][j] =
Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
fluxvector_homo[i][j] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
for (int k = 0; k < nConvectiveFields ; ++k)
{
......@@ -189,6 +163,54 @@ namespace Nektar
Array<OneD, NekDouble>(nPointsTot_plane, 0.0);
}
}
}
outarray_homo = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
outarray_homo[j] = Array<OneD, NekDouble>(nPointsTot, 0.0);
}
m_planeDiff->SetRiemannSolver(m_riemann);
m_planeDiff->SetFluxVectorVecNS(m_fluxVectorNS);
}
/**
* @brief Calculate WeakDG Diffusion for the Navier-Stokes (NS) equations
* using an LDG interface flux and the the flux in the third direction.
*
* The equations that need a diffusion operator are those related
* with the velocities and with the energy.
*/
void DiffusionLDGNS3DHomogeneous1D ::v_Diffuse(
const int nConvectiveFields,
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
for (i = 0; i < nVariables; ++i)
{
fields[0]->PhysDeriv(2, inarray[i], derivatives01_homo[i]);
}
for (i = 0; i < num_planes; ++i)
{
for (j = 0; j < nConvectiveFields; j ++)
{
fields_plane[i][j]= fields[j]->GetPlane(i);
}
for (j = 0; j < nVariables; j ++)
{
Vmath::Vcopy(nPointsTot_plane,
&inarray[j][i * nPointsTot_plane], 1,
&inarray_plane[i][j][0], 1);
Vmath::Vcopy(nPointsTot_plane,
&derivatives01_homo[j][i * nPointsTot_plane], 1,
&derivatives01_plane[i][j][0], 1);
}
// Set the first order derivatives in 3rd direction
diffLDGNS->SetHomoDerivs(derivatives01_plane[i]);
......@@ -215,12 +237,8 @@ namespace Nektar
}
}
Array<OneD, Array<OneD, NekDouble> > outarray_homo(nConvectiveFields);
for (i = 0; i < nConvectiveFields; ++i)
{
outarray_homo[i] = Array<OneD, NekDouble>(nPointsTot, 0.0);
fields[0]->PhysDeriv(2, fluxvector[i], outarray_homo[i]);
Vmath::Vadd(nPointsTot,
......
......@@ -64,8 +64,25 @@ namespace Nektar
int num_planes;
int nPointsTot_plane;
int nCoeffs_plane;
int nVariables;
int i, j, k;
Array<OneD, unsigned int> planes;
Array<OneD, Array<OneD, NekDouble> > fluxvector;
Array<OneD, Array<OneD, NekDouble> > derivatives01_homo;
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
fields_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
inarray_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
derivatives01_plane;
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
outarray_plane;
Array<OneD, Array<OneD, Array<OneD, Array<OneD, NekDouble> > > >
fluxvector_homo;
Array<OneD, Array<OneD, NekDouble> > outarray_homo;
DiffusionLDGNSSharedPtr diffLDGNS;
virtual void v_InitObject(
LibUtilities::SessionReaderSharedPtr pSession,
......
......@@ -152,6 +152,44 @@ namespace Nektar
nPointsTot_plane = nPointsTot/num_planes;
nCoeffs_plane = nCoeffs/num_planes;
fields_plane =
Array <OneD, Array<OneD, MultiRegions::ExpListSharedPtr> >
(num_planes);
inarray_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
outarray_plane = Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
(num_planes);
for (i = 0; i < num_planes; ++i)
{
fields_plane[i] = Array<OneD, MultiRegions::ExpListSharedPtr>
(nConvectiveFields);
inarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
outarray_plane[i] = Array<OneD, Array<OneD, NekDouble> >
(nConvectiveFields);
for (j = 0; j < nConvectiveFields; j ++)
{
inarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
outarray_plane[i][j] = Array<OneD, NekDouble>
(nPointsTot_plane, 0.0);
}
}
flux = Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
flux_homo =
Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);
outarray_z =
Array<OneD, Array<OneD, NekDouble> > (nConvectiveFields);