Skip to content

Added pairwise gradient and improved instantiations #477

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
brbass wants to merge 1 commit into
base: develop
Choose a base branch
from
+108 −47
Open

Added pairwise gradient and improved instantiations #477

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 11 additions & 0 deletions src/FieldOperations/FieldListFunctions.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@

namespace Spheral {

template<typename Dimension> class ConnectivityMap;
template<typename Dimension> class TableKernel;

template<typename Dimension, typename DataType> class FieldList;
Expand Down Expand Up @@ -44,6 +45,16 @@ gradient(const FieldList<Dimension, std::vector<DataType>>& fieldList,
const FieldList<Dimension, typename Dimension::SymTensor>& Hfield,
const TableKernel<Dimension>& kernel);

template<typename Dimension, typename DataType>
void
gradientPairs(FieldList<Dimension, typename MathTraits<Dimension, DataType>::GradientType>& result,
const FieldList<Dimension, DataType>& field,
const FieldList<Dimension, typename Dimension::Vector>& position,
const FieldList<Dimension, typename Dimension::Scalar>& weight,
const FieldList<Dimension, typename Dimension::SymTensor>& H,
const ConnectivityMap<Dimension>& conn,
const TableKernel<Dimension>& kernel);

// Calculate the divergence of the given FieldList.
template<typename Dimension, typename DataType>
FieldList<Dimension, typename MathTraits<Dimension, DataType>::DivergenceType>
Expand Down
55 changes: 55 additions & 0 deletions src/FieldOperations/gradient.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "Field/FieldList.hh"
#include "Field/Field.hh"
#include "Field/NodeIterators.hh"
#include "Neighbor/ConnectivityMap.hh"
#include "NodeList/NodeList.hh"
#include "Neighbor/Neighbor.hh"
#include "Kernel/TableKernel.hh"
Expand Down Expand Up @@ -277,6 +278,60 @@ gradient(const FieldList<Dimension, std::vector<DataType>>& fieldList,
return result;
}

//------------------------------------------------------------------------------
// Calculate the gradient of a FieldList.
//------------------------------------------------------------------------------
template<typename Dimension, typename DataType>
void
gradientPairs(FieldList<Dimension, typename MathTraits<Dimension, DataType>::GradientType>& result,
const FieldList<Dimension, DataType>& field,
const FieldList<Dimension, typename Dimension::Vector>& position,
const FieldList<Dimension, typename Dimension::Scalar>& weight,
const FieldList<Dimension, typename Dimension::SymTensor>& H,
const ConnectivityMap<Dimension>& conn,
const TableKernel<Dimension>& kernel) {
typedef typename MathTraits<Dimension, DataType>::GradientType GradientType;
result = GradientType::zero();

const auto& pairs = conn.nodePairList();
const auto npairs = pairs.size();

for (auto k = 0u; k < npairs; ++k) {
const auto ni = pairs[k].i_list;
const auto nj = pairs[k].j_list;
const auto i = pairs[k].i_node;
const auto j = pairs[k].j_node;

const auto& ri = position(ni, i);
const auto& rj = position(nj, j);
const auto& vi = weight(ni, i);
const auto& vj = weight(nj, j);
const auto& hi = H(ni, i);
const auto& hj = H(nj, j);
const auto& fi = field(ni, i);
const auto& fj = field(nj, j);

const auto rij = ri - rj;
const auto hdeti = hi.Determinant();
const auto hdetj = hj.Determinant();
const auto etai = hi * rij;
const auto etaj = hj * rij;
const auto etaMagi = etai.magnitude();
const auto etaMagj = etaj.magnitude();
const auto etaUniti = etai.unitVector();
const auto etaUnitj = etaj.unitVector();
const auto hetaUniti = hi * etaUniti;
const auto hetaUnitj = hj * etaUnitj;

const auto dwi = hetaUniti * kernel.gradValue(etaMagi, hdeti);
const auto dwj = hetaUnitj * kernel.gradValue(etaMagj, hdetj);
const auto dwij = 0.5 * (dwi + dwj);

result(ni, i) += vj * (fj - fi) * dwij;
result(nj, j) -= vi * (fi - fj) * dwij;
}
}

//------------------------------------------------------------------------------
// The limiter method.
//------------------------------------------------------------------------------
Expand Down
89 changes: 42 additions & 47 deletions src/FieldOperations/gradientInst.cc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,62 +6,57 @@
#include "Geometry/Dimension.hh"

namespace Spheral {
"""

//============================== gradient() ==============================
template
FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>::GradientType>
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& fieldList,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
const TableKernel< Dim< %(ndim)s > >& kernel);
for Value in ("Scalar", "Vector"):
text += """
template
FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>::GradientType>
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& fieldList,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
const TableKernel< Dim< %(ndim)s > >& kernel);
FieldList<%%(Dim)s, MathTraits<%%(Dim)s, %(Value)s>::GradientType>
gradient<%%(Dim)s, %(Value)s>(const FieldList<%%(Dim)s, %(Value)s>& fieldList,
const FieldList<%%(Dim)s, %%(Vector)s>& position,
const FieldList<%%(Dim)s, %%(Scalar)s>& weight,
const FieldList<%%(Dim)s, %%(Scalar)s>& mass,
const FieldList<%%(Dim)s, %%(Scalar)s>& density,
const FieldList<%%(Dim)s, %%(SymTensor)s>& Hfield,
const TableKernel< %%(Dim)s >& kernel);

template
FieldList<Dim< %(ndim)s >, std::vector<MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>::GradientType>>
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>(const FieldList<Dim< %(ndim)s >, std::vector<Dim< %(ndim)s >::Scalar>>& fieldList,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
const TableKernel< Dim< %(ndim)s > >& kernel);
template
FieldList<Dim< %(ndim)s >, std::vector<MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>::GradientType>>
gradient<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>(const FieldList<Dim< %(ndim)s >, std::vector<Dim< %(ndim)s >::Vector>>& fieldList,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& weight,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& mass,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& density,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
const TableKernel< Dim< %(ndim)s > >& kernel);
FieldList<%%(Dim)s, std::vector<MathTraits<%%(Dim)s, %(Value)s>::GradientType>>
gradient<%%(Dim)s, %(Value)s>(const FieldList<%%(Dim)s, std::vector<%(Value)s>>& fieldList,
const FieldList<%%(Dim)s, %%(Vector)s>& position,
const FieldList<%%(Dim)s, %%(Scalar)s>& weight,
const FieldList<%%(Dim)s, %%(Scalar)s>& mass,
const FieldList<%%(Dim)s, %%(Scalar)s>& density,
const FieldList<%%(Dim)s, %%(SymTensor)s>& Hfield,
const TableKernel< %%(Dim)s >& kernel);

template
void
gradientPairs<%%(Dim)s, %(Value)s>(FieldList<%%(Dim)s, MathTraits<%%(Dim)s, %(Value)s>::GradientType>& result,
const FieldList<%%(Dim)s, %(Value)s>& field,
const FieldList<%%(Dim)s, %%(Vector)s>& position,
const FieldList<%%(Dim)s, %%(Scalar)s>& weight,
const FieldList<%%(Dim)s, %%(SymTensor)s>& Hfield,
const ConnectivityMap<%%(Dim)s>& conn,
const TableKernel<%%(Dim)s>& kernel);
""" % {"Value" : "%(" + Value + ")s"}

text += """
//============================== limiter() ==============================
template
FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>
limiter<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>& fieldList,
const FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Scalar>::GradientType>& gradient,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
const TableKernel< Dim< %(ndim)s > >& kernel);
FieldList<%(Dim)s, %(SymTensor)s>
limiter<%(Dim)s, %(Scalar)s>(const FieldList<%(Dim)s, %(Scalar)s>& fieldList,
const FieldList<%(Dim)s, MathTraits<%(Dim)s, %(Scalar)s>::GradientType>& gradient,
const FieldList<%(Dim)s, %(Vector)s>& position,
const FieldList<%(Dim)s, %(SymTensor)s>& Hfield,
const TableKernel< %(Dim)s >& kernel);
template
FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>
limiter<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>(const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& fieldList,
const FieldList<Dim< %(ndim)s >, MathTraits<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>::GradientType>& gradient,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::Vector>& position,
const FieldList<Dim< %(ndim)s >, Dim< %(ndim)s >::SymTensor>& Hfield,
const TableKernel< Dim< %(ndim)s > >& kernel);
FieldList<%(Dim)s, %(SymTensor)s>
limiter<%(Dim)s, %(Vector)s>(const FieldList<%(Dim)s, %(Vector)s>& fieldList,
const FieldList<%(Dim)s, MathTraits<%(Dim)s, %(Vector)s>::GradientType>& gradient,
const FieldList<%(Dim)s, %(Vector)s>& position,
const FieldList<%(Dim)s, %(SymTensor)s>& Hfield,
const TableKernel< %(Dim)s >& kernel);

}
"""