Skip to content

GitLab

Commit 10455dda authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'master' into feature/nodal-cleanup

parents c88edcae 6b52a85e
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
View file @ 10455dda
......@@ -3,9 +3,11 @@ Changelog
v4.4.0
------
**Library:**
**Library**:
- Add support for variable polynomial order for 3D simulations with continuous
Galerkin discretisation (!604)
- Bump version of gsmpi to suppress autotuning output unless `--verbose` is
specified (!652)
- Add support for variable polynomial order with periodic boundary conditions
(!658)
- Statistics are now printed for lowest level of multi-level static condensation
......@@ -15,8 +17,12 @@ v4.4.0
- New FieldUtils library allows support for most `FieldConvert` post-processing
operations during simulation using a new filter (!589)
- Adjust CMake dependencies to reduce compile time (!671)
- Homogeneous1D dealiasing improvements (!622)
- Rework nodal utilities to support nodal prismatic elements (!660)
**ADRSolver:**
- Add a projection equation system for C^0 projections (!675)
**APESolver:**
- Use a continuous basefield projection and revert to constant c formulation (!664)
- Added ability to compute CFL number (!664)
......@@ -24,17 +30,25 @@ v4.4.0
**IncNavierStokesSolver:**
- Add ability to simulate additional scalar fields (!624)
- Improve performance when using homogeneous dealiasing (!622)
**NekMesh:**
- Modify curve module to allow for spline input (!628)
- New module for inserting an alternate high-order surface into the
the working mesh (!669)
- Add STL surface writer module (!668)
- New module for inserting an alternate high-order surface into the working
mesh (!669)
**FieldConvert:**
- Move all modules to a new library, FieldUtils, to support post-processing
during simulations (!589)
- Add module to stretch homogeneous direction (!609)
v4.3.4
------
**Library:**
- Fix performance issue with `v_ExtractDataToCoeffs` for post-processing of large
simulations (!672)
v4.3.3
------
**Library**:
......
cmake/ThirdPartyMPI.cmake
View file @ 10455dda
......@@ -51,23 +51,23 @@ IF( NEKTAR_USE_MPI )
IF (THIRDPARTY_BUILD_GSMPI)
EXTERNALPROJECT_ADD(
gsmpi-1.2
URL ${TPURL}/gsmpi-1.2.tar.bz2
URL_MD5 35901be16791bfdeafa9c4d0e06d189b
gsmpi-1.2.1
URL ${TPURL}/gsmpi-1.2.1.tar.bz2
URL_MD5 18dcb4cd1dcc7876173465c404b1142d
STAMP_DIR ${TPBUILD}/stamp
DOWNLOAD_DIR ${TPSRC}
SOURCE_DIR ${TPSRC}/gsmpi-1.2
BINARY_DIR ${TPBUILD}/gsmpi-1.2
TMP_DIR ${TPBUILD}/gsmpi-1.2-tmp
SOURCE_DIR ${TPSRC}/gsmpi-1.2.1
BINARY_DIR ${TPBUILD}/gsmpi-1.2.1
TMP_DIR ${TPBUILD}/gsmpi-1.2.1-tmp
INSTALL_DIR ${TPDIST}
CONFIGURE_COMMAND
CONFIGURE_COMMAND
${CMAKE_COMMAND}
-G ${CMAKE_GENERATOR}
-DCMAKE_C_COMPILER:FILEPATH=${CMAKE_C_COMPILER}
-DCMAKE_CXX_COMPILER:FILEPATH=${CMAKE_CXX_COMPILER}
-DCMAKE_BUILD_TYPE:STRING=Debug
-DCMAKE_INSTALL_PREFIX:PATH=${TPDIST}
${TPSRC}/gsmpi-1.2
${TPSRC}/gsmpi-1.2.1
)
SET(GSMPI_LIBRARY gsmpi CACHE FILEPATH
"GSMPI path" FORCE)
......@@ -78,7 +78,7 @@ IF( NEKTAR_USE_MPI )
MESSAGE(STATUS "Build GSMPI: ${TPDIST}/lib/lib${GSMPI_LIBRARY}.a")
MESSAGE(STATUS "Build XXT: ${TPDIST}/lib/lib${XXT_LIBRARY}.a")
ELSE (THIRDPARTY_BUILD_GSMPI)
ADD_CUSTOM_TARGET(gsmpi-1.2 ALL)
ADD_CUSTOM_TARGET(gsmpi-1.2.1 ALL)
INCLUDE (FindGSMPI)
INCLUDE (FindXXT)
ENDIF (THIRDPARTY_BUILD_GSMPI)
......
docs/user-guide/solvers/adr.tex
View file @ 10455dda
......@@ -24,6 +24,8 @@ see the table below.
\textbf{Equation to solve} & \textbf{EquationType} & \textbf{Dimensions} &
\textbf{Projections} \\
\midrule
$u = f$ &
\inltt{Projection} & All & Continuous/Discontinuous \\
$\nabla^2 u = 0$ &
\inltt{Laplace} & All & Continuous/Discontinuous \\
$\nabla^2 u = f$ &
......
library/LibUtilities/CMakeLists.txt
View file @ 10455dda
......@@ -419,7 +419,7 @@ IF( NEKTAR_USE_MPI )
SET_TARGET_PROPERTIES(LibUtilities
PROPERTIES COMPILE_FLAGS "${THE_COMPILE_FLAGS} ${MPI_COMPILE_FLAGS}")
ENDIF()
ADD_DEPENDENCIES(LibUtilities gsmpi-1.2)
ADD_DEPENDENCIES(LibUtilities gsmpi-1.2.1)
ENDIF( NEKTAR_USE_MPI )
# Lapack and Blas
......
library/LibUtilities/Communication/GsLib.hpp
View file @ 10455dda
......@@ -148,8 +148,10 @@ namespace Gs
* communication.
* @return GSLib data structure containing mapping information.
*/
static inline gs_data* Init ( const Nektar::Array<OneD, long> pId,
const LibUtilities::CommSharedPtr& pComm)
static inline gs_data* Init(
const Nektar::Array<OneD, long> pId,
const LibUtilities::CommSharedPtr &pComm,
bool verbose = true)
{
#ifdef NEKTAR_USE_MPI
if (pComm->GetSize() == 1)
......@@ -162,8 +164,8 @@ namespace Gs
MPI_Comm_dup(vCommMpi->GetComm(), &vComm.c);
vComm.id = vCommMpi->GetRank();
vComm.np = vCommMpi->GetSize();
gs_data* result = nektar_gs_setup(pId.get(),pId.num_elements(),
&vComm, 0, gs_auto, 1);
gs_data* result = nektar_gs_setup(pId.get(),pId.num_elements(),
&vComm, 0, gs_auto, (int)verbose);
MPI_Comm_free(&vComm.c);
return result;
#else
......
library/MultiRegions/AssemblyMap/AssemblyMapCG.cpp
View file @ 10455dda
......@@ -1112,12 +1112,14 @@ namespace Nektar
int unique_edges = foundEdges.size();
int unique_faces = foundFaces.size();
bool verbose = m_session->DefinesCmdLineArgument("verbose");
// Now construct temporary GS objects. These will be used to
// populate the arrays tmp3 and tmp4 with the multiplicity of
// the vertices and edges respectively to identify those
// vertices and edges which are located on partition boundary.
Array<OneD, long> vertArray(unique_verts, &procVerts[0]);
Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm);
Gs::gs_data *tmp1 = Gs::Init(vertArray, vComm, verbose);
Array<OneD, NekDouble> tmp4(unique_verts, 1.0);
Array<OneD, NekDouble> tmp5(unique_edges, 1.0);
Array<OneD, NekDouble> tmp6(unique_faces, 1.0);
......@@ -1127,7 +1129,7 @@ namespace Nektar
if (unique_edges > 0)
{
Array<OneD, long> edgeArray(unique_edges, &procEdges[0]);
Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm);
Gs::gs_data *tmp2 = Gs::Init(edgeArray, vComm, verbose);
Gs::Gather(tmp5, Gs::gs_add, tmp2);
Gs::Finalise(tmp2);
}
......@@ -1135,7 +1137,7 @@ namespace Nektar
if (unique_faces > 0)
{
Array<OneD, long> faceArray(unique_faces, &procFaces[0]);
Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm);
Gs::gs_data *tmp3 = Gs::Init(faceArray, vComm, verbose);
Gs::Gather(tmp6, Gs::gs_add, tmp3);
Gs::Finalise(tmp3);
}
......@@ -1324,6 +1326,8 @@ namespace Nektar
const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
bool verbose = m_session->DefinesCmdLineArgument("verbose");
m_signChange = false;
// Stores vertex, edge and face reordered vertices.
......@@ -1434,7 +1438,7 @@ namespace Nektar
edgeId[i] = dofIt->first + 1;
edgeDof[i] = (NekDouble) dofIt->second;
}
Gs::gs_data *tmp = Gs::Init(edgeId, vComm);
Gs::gs_data *tmp = Gs::Init(edgeId, vComm, verbose);
Gs::Gather(edgeDof, Gs::gs_min, tmp);
Gs::Finalise(tmp);
for (i=0; i < dofs[1].size(); i++)
......@@ -1465,7 +1469,7 @@ namespace Nektar
faceP[i] = (NekDouble) dofIt->second;
faceQ[i] = (NekDouble) dofIt2->second;
}
Gs::gs_data *tmp2 = Gs::Init(faceId, vComm);
Gs::gs_data *tmp2 = Gs::Init(faceId, vComm, verbose);
Gs::Gather(faceP, Gs::gs_min, tmp2);
Gs::Gather(faceQ, Gs::gs_min, tmp2);
Gs::Finalise(tmp2);
......@@ -2188,6 +2192,7 @@ namespace Nektar
const LocalRegions::ExpansionVector &locExpVector = *(locExp.GetExp());
LibUtilities::CommSharedPtr vCommRow = m_comm->GetRowComm();
const bool verbose = locExp.GetSession()->DefinesCmdLineArgument("verbose");
m_globalToUniversalMap = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
m_globalToUniversalMapUnique = Nektar::Array<OneD, int>(m_numGlobalCoeffs, -1);
......@@ -2355,8 +2360,8 @@ namespace Nektar
tmp[i] = m_globalToUniversalMap[i];
}
m_gsh = Gs::Init(tmp, vCommRow);
m_bndGsh = Gs::Init(tmp2, vCommRow);
m_gsh = Gs::Init(tmp, vCommRow, verbose);
m_bndGsh = Gs::Init(tmp2, vCommRow, verbose);
Gs::Unique(tmp, vCommRow);
for (unsigned int i = 0; i < m_numGlobalCoeffs; ++i)
{
......@@ -2386,6 +2391,7 @@ namespace Nektar
const boost::shared_ptr<LocalRegions::ExpansionVector> exp
= locexp.GetExp();
int nelmts = exp->size();
const bool verbose = locexp.GetSession()->DefinesCmdLineArgument("verbose");
// Get Default Map and turn off any searched values.
returnval = MemoryManager<AssemblyMapCG>
......@@ -2482,7 +2488,7 @@ namespace Nektar
{
tmp[i] = returnval->m_globalToUniversalMap[i];
}
returnval->m_gsh = Gs::Init(tmp, vCommRow);
returnval->m_gsh = Gs::Init(tmp, vCommRow, verbose);
Gs::Unique(tmp, vCommRow);
for (unsigned int i = 0; i < nglocoeffs; ++i)
{
......
library/MultiRegions/ExpList.cpp
View file @ 10455dda
......@@ -2179,17 +2179,19 @@ namespace Nektar
ASSERTL0(i != fielddef->m_fields.size(),
"Field (" + field + ") not found in file.");
// Determine mapping from element ids to location in expansion list
map<int, int> elmtToExpId;
// Loop in reverse order so that in case where using a Homogeneous
// expansion it sets geometry ids to first part of m_exp
// list. Otherwise will set to second (complex) expansion
for(i = (*m_exp).size()-1; i >= 0; --i)
if (m_elmtToExpId.size() == 0)
{
elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
// Loop in reverse order so that in case where using a
// Homogeneous expansion it sets geometry ids to first part of
// m_exp list. Otherwise will set to second (complex) expansion
for(i = (*m_exp).size()-1; i >= 0; --i)
{
m_elmtToExpId[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
}
}
boost::unordered_map<int, int>::iterator eIt;
for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
{
// Reset modes_offset in the case where all expansions of
......@@ -2203,14 +2205,16 @@ namespace Nektar
fielddef->m_numModes, modes_offset);
const int elmtId = fielddef->m_elementIDs[i];
if (elmtToExpId.count(elmtId) == 0)
eIt = m_elmtToExpId.find(elmtId);
if (eIt == m_elmtToExpId.end())
{
offset += datalen;
modes_offset += (*m_exp)[0]->GetNumBases();
continue;
}
expId = elmtToExpId[elmtId];
expId = eIt->second;
if (datalen == (*m_exp)[expId]->GetNcoeffs())
{
......@@ -2220,8 +2224,8 @@ namespace Nektar
else
{
(*m_exp)[expId]->ExtractDataToCoeffs(
&fielddata[offset], fielddef->m_numModes,
modes_offset, &coeffs[m_coeff_offset[expId]]);
&fielddata[offset], fielddef->m_numModes,
modes_offset, &coeffs[m_coeff_offset[expId]]);
}
offset += datalen;
......@@ -2447,7 +2451,20 @@ namespace Nektar
"This method is not defined or valid for this class type");
}
void ExpList::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,const Array<OneD, NekDouble> &inarray2,Array<OneD, NekDouble> &outarray,CoeffState coeffstate)
void ExpList::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,
const Array<OneD, NekDouble> &inarray2,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
}
void ExpList::v_DealiasedDotProd(
const Array<OneD, Array<OneD, NekDouble> > &inarray1,
const Array<OneD, Array<OneD, NekDouble> > &inarray2,
Array<OneD, Array<OneD, NekDouble> > &outarray,
CoeffState coeffstate)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
......
library/MultiRegions/ExpList.h
View file @ 10455dda
......@@ -348,7 +348,13 @@ namespace Nektar
const Array<OneD, NekDouble> &inarray2,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal);
inline void DealiasedDotProd(
const Array<OneD, Array<OneD, NekDouble> > &inarray1,
const Array<OneD, Array<OneD, NekDouble> > &inarray2,
Array<OneD, Array<OneD, NekDouble> > &outarray,
CoeffState coeffstate = eLocal);
inline void GetBCValues(
Array<OneD, NekDouble> &BndVals,
const Array<OneD, NekDouble> &TotField,
......@@ -875,7 +881,7 @@ namespace Nektar
}
/// Returns the session object
boost::shared_ptr<LibUtilities::SessionReader> GetSession()
boost::shared_ptr<LibUtilities::SessionReader> GetSession() const
{
return m_session;
}
......@@ -1020,6 +1026,9 @@ namespace Nektar
// or not
bool m_WaveSpace;
/// Mapping from geometry ID of element to index inside #m_exp
boost::unordered_map<int, int> m_elmtToExpId;
/// This function assembles the block diagonal matrix of local
/// matrices of the type \a mtype.
const DNekScalBlkMatSharedPtr GenBlockMatrix(
......@@ -1254,7 +1263,13 @@ namespace Nektar
const Array<OneD, NekDouble> &inarray2,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate = eLocal);
virtual void v_DealiasedDotProd(
const Array<OneD, Array<OneD, NekDouble> > &inarray1,
const Array<OneD, Array<OneD, NekDouble> > &inarray2,
Array<OneD, Array<OneD, NekDouble> > &outarray,
CoeffState coeffstate = eLocal);
virtual void v_GetBCValues(
Array<OneD, NekDouble> &BndVals,
const Array<OneD, NekDouble> &TotField,
......@@ -1777,7 +1792,19 @@ namespace Nektar
{
v_DealiasedProd(inarray1,inarray2,outarray,coeffstate);
}
/**
*
*/
inline void ExpList::DealiasedDotProd(
const Array<OneD, Array<OneD, NekDouble> > &inarray1,
const Array<OneD, Array<OneD, NekDouble> > &inarray2,
Array<OneD, Array<OneD, NekDouble> > &outarray,
CoeffState coeffstate)
{
v_DealiasedDotProd(inarray1,inarray2,outarray,coeffstate);
}
/**
*
*/
......
library/MultiRegions/ExpListHomogeneous1D.cpp
View file @ 10455dda
......@@ -88,8 +88,9 @@ namespace Nektar
{
if(m_useFFT)
{
NekDouble size = 1.5*m_homogeneousBasis->GetNumPoints();
m_padsize = int(size);
int size = m_homogeneousBasis->GetNumPoints() +
m_homogeneousBasis->GetNumPoints() / 2;
m_padsize = size + (size % 2);
m_FFT_deal = LibUtilities::GetNektarFFTFactory()
.CreateInstance("NekFFTW", m_padsize);
}
......@@ -168,26 +169,33 @@ namespace Nektar
/**
* Dealiasing routine
*
* @param inarray1 First term of the product
* @param inarray2 Second term of the product
* @param outarray Dealiased product
*/
void ExpListHomogeneous1D::v_DealiasedProd(const Array<OneD, NekDouble> &inarray1,
const Array<OneD, NekDouble> &inarray2,
Array<OneD, NekDouble> &outarray,
CoeffState coeffstate)
{
// inarray1 = first term of the product in full physical space
// inarray2 = second term of the product in full physical space
// dealiased product stored in outarray
int num_dofs = inarray1.num_elements();
int N = m_homogeneousBasis->GetNumPoints();
Array<OneD, NekDouble> V1(num_dofs);
Array<OneD, NekDouble> V2(num_dofs);
Array<OneD, NekDouble> V1V2(num_dofs);
HomogeneousFwdTrans(inarray1,V1,coeffstate);
HomogeneousFwdTrans(inarray2,V2,coeffstate);
if(m_WaveSpace)
{
V1 = inarray1;
V2 = inarray2;
}
else
{
HomogeneousFwdTrans(inarray1,V1,coeffstate);
HomogeneousFwdTrans(inarray2,V2,coeffstate);
}
int num_points_per_plane = num_dofs/m_planes.num_elements();
int num_proc;
......@@ -245,11 +253,189 @@ namespace Nektar
&(ShufV1V2[i*N]), 1);
}
m_transposition->Transpose(ShufV1V2, V1V2, false,
// Moving the results to the output
if (m_WaveSpace)
{
m_transposition->Transpose(ShufV1V2, outarray, false,
LibUtilities::eZtoXY);
}
else
{
m_transposition->Transpose(ShufV1V2, V1V2, false,
LibUtilities::eZtoXY);
HomogeneousBwdTrans(V1V2, outarray, coeffstate);
}
}
/**
* Dealiasing routine for dot product
*
* @param inarray1 First term of the product with dimension ndim
* (e.g. U)
* @param inarray2 Second term of the product with dimension ndim*nvec
* (e.g. grad U)
* @param outarray Dealiased product with dimension nvec
*/
void ExpListHomogeneous1D::v_DealiasedDotProd(
const Array<OneD, Array<OneD, NekDouble> > &inarray1,
const Array<OneD, Array<OneD, NekDouble> > &inarray2,
Array<OneD, Array<OneD, NekDouble> > &outarray,
CoeffState coeffstate)
{
int ndim = inarray1.num_elements();
ASSERTL1( inarray2.num_elements() % ndim == 0,
"Wrong dimensions for DealiasedDotProd.");
int nvec = inarray2.num_elements() / ndim;
int num_dofs = inarray1[0].num_elements();
int N = m_homogeneousBasis->GetNumPoints();
int num_points_per_plane = num_dofs/m_planes.num_elements();
int num_proc;
if(!m_session->DefinesSolverInfo("HomoStrip"))
{
num_proc = m_comm->GetColumnComm()->GetSize();
}
else
{
num_proc = m_StripZcomm->GetSize();
}
int num_dfts_per_proc = num_points_per_plane / num_proc
+ (num_points_per_plane % num_proc > 0);
// Get inputs in Fourier space
Array<OneD, Array<OneD, NekDouble> > V1(ndim);
Array<OneD, Array<OneD, NekDouble> > V2(ndim*nvec);
if(m_WaveSpace)
{
for (int i = 0; i < ndim; i++)
{
V1[i] = inarray1[i];
}
for (int i = 0; i < ndim*nvec; i++)
{
V2[i] = inarray2[i];
}
}
else
{
for (int i = 0; i < ndim; i++)
{
V1[i] = Array<OneD, NekDouble> (num_dofs);
HomogeneousFwdTrans(inarray1[i],V1[i],coeffstate);
}
for (int i = 0; i < ndim*nvec; i++)
{
V2[i] = Array<OneD, NekDouble> (num_dofs);
HomogeneousFwdTrans(inarray2[i],V2[i],coeffstate);
}
}
// Allocate variables for ffts
Array<OneD, Array<OneD, NekDouble> > ShufV1(ndim);
Array<OneD, NekDouble> ShufV1_PAD_coef(m_padsize,0.0);
Array<OneD, Array<OneD, NekDouble> > ShufV1_PAD_phys(ndim);
for (int i = 0; i < ndim; i++)
{
ShufV1[i] = Array<OneD, NekDouble>
(num_dfts_per_proc*N,0.0);
ShufV1_PAD_phys[i] = Array<OneD, NekDouble>
(m_padsize,0.0);
}
// Moving the results in physical space for the output
HomogeneousBwdTrans(V1V2, outarray, coeffstate);
Array<OneD, Array<OneD, NekDouble> > ShufV2(ndim*nvec);
Array<OneD, NekDouble> ShufV2_PAD_coef(m_padsize,0.0);
Array<OneD, Array<OneD, NekDouble> > ShufV2_PAD_phys(ndim*nvec);
for (int i = 0; i < ndim*nvec; i++)
{
ShufV2[i] = Array<OneD, NekDouble>
(num_dfts_per_proc*N,0.0);
ShufV2_PAD_phys[i] = Array<OneD, NekDouble>
(m_padsize,0.0);
}
Array<OneD, Array<OneD, NekDouble> > ShufV1V2(nvec);
Array<OneD, NekDouble> ShufV1V2_PAD_coef(m_padsize,0.0);
Array<OneD, NekDouble> ShufV1V2_PAD_phys(m_padsize,0.0);
for (int i = 0; i < nvec; i++)
{
ShufV1V2[i] = Array<OneD, NekDouble>
(num_dfts_per_proc*N,0.0);
}
// Transpose inputs
for (int i = 0; i < ndim; i++)
{
m_transposition->Transpose(V1[i], ShufV1[i], false,
LibUtilities::eXYtoZ);
}
for (int i = 0; i < ndim*nvec; i++)
{
m_transposition->Transpose(V2[i], ShufV2[i], false,
LibUtilities::eXYtoZ);
}
// Looping on the pencils
for(int i = 0 ; i < num_dfts_per_proc ; i++)
{
for (int j = 0; j < ndim; j++)
{
// Copying the i-th pencil pf lenght N into a bigger
// pencil of lenght 1.5N We are in Fourier space
Vmath::Vcopy(N, &(ShufV1[j][i*N]), 1,
&(ShufV1_PAD_coef[0]), 1);
// Moving to physical space using the padded system
m_FFT_deal->FFTBwdTrans(ShufV1_PAD_coef, ShufV1_PAD_phys[j]);
}
for (int j = 0; j < ndim*nvec; j++)
{
Vmath::Vcopy(N, &(ShufV2[j][i*N]), 1,
&(ShufV2_PAD_coef[0]), 1);
m_FFT_deal->FFTBwdTrans(ShufV2_PAD_coef, ShufV2_PAD_phys[j]);
}
// Performing the vectors multiplication in physical space on
// the padded system
for (int j = 0; j < nvec; j++)
{
Vmath::Zero(m_padsize, ShufV1V2_PAD_phys, 1);
for (int k = 0; k < ndim; k++)
{
Vmath::Vvtvp(m_padsize, ShufV1_PAD_phys[k], 1,
ShufV2_PAD_phys[j*ndim+k], 1,
ShufV1V2_PAD_phys, 1,
ShufV1V2_PAD_phys, 1);
}
// Moving back the result (V1*V2)_phys in Fourier space,
// padded system