Skip to content
GitLab
Commit cb59f55d authored by Matthias Messner's avatar Matthias Messner
Browse files

adjusted sign in P2P for force computation; interleaved potential and force in 

L2P; inserted mutual P2P (symmetry of matrix)
parent 3ce09f30
Hide whitespace changes
Inline Side-by-side
Src/Kernels/Chebyshev/FAbstractChebKernel.hpp
View file @ cb59f55d
......@@ -114,38 +114,80 @@ public:
{
// loop: target particles
typename ContainerClass::BasicIterator iTargets(*TargetParticles);
while (iTargets.hasNotFinished()) {
ParticleClass& Target = iTargets.data();
{ // loop: source particles (target leaf cell == source leaf cell)
typename ContainerClass::ConstBasicIterator iSources(*SourceParticles);
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// only if target and source are not identical
if (&Target != &Source) this->directInteraction(Target, Source);
// progress sources
iSources.gotoNext();
if (TargetParticles != SourceParticles) {
while (iTargets.hasNotFinished()) {
ParticleClass& Target = iTargets.data();
{ // loop: source particles (target leaf cell == source leaf cell)
typename ContainerClass::ConstBasicIterator iSources(*SourceParticles);
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// only if target and source are not identical
this->directInteraction(Target, Source);
// progress sources
iSources.gotoNext();
}
}
{ // loop: source particles (target leaf cell != source leaf cell)
for (unsigned int idx=0; idx<27; ++idx) {
if (NeighborSourceParticles[idx]) {
typename ContainerClass::ConstBasicIterator iSources(*NeighborSourceParticles[idx]);
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// target and source cannot be identical
this->directInteraction(Target, Source);
// progress sources
iSources.gotoNext();
}
}
}
}
// progress targets
iTargets.gotoNext();
}
{ // loop: source particles (target leaf cell != source leaf cell)
for (unsigned int idx=0; idx<27; ++idx) {
if (NeighborSourceParticles[idx]) {
typename ContainerClass::ConstBasicIterator iSources(*NeighborSourceParticles[idx]);
while (iSources.hasNotFinished()) {
const ParticleClass& Source = iSources.data();
// target and source cannot be identical
this->directInteraction(Target, Source);
// progress sources
iSources.gotoNext();
} else {
while (iTargets.hasNotFinished()) {
ParticleClass& Target = iTargets.data();
{ // loop: source particles (target leaf cell == source leaf cell)
typename ContainerClass::BasicIterator iSources = iTargets;
iSources.gotoNext();
while (iSources.hasNotFinished()) {
ParticleClass& Source = iSources.data();
// only if target and source are not identical
this->directInteractionMutual(Target, Source);
// progress sources
iSources.gotoNext();
}
}
{ // loop: source particles (target leaf cell != source leaf cell)
for (unsigned int idx=0; idx<=13; ++idx) {
if (NeighborSourceParticles[idx]) {
typename ContainerClass::BasicIterator iSources(*NeighborSourceParticles[idx]);
while (iSources.hasNotFinished()) {
ParticleClass& Source = iSources.data();
// target and source cannot be identical
this->directInteractionMutual(Target, Source);
// progress sources
iSources.gotoNext();
}
}
}
}
// progress targets
iTargets.gotoNext();
}
// progress targets
iTargets.gotoNext();
}
}
......@@ -164,6 +206,21 @@ private:
Target.incForces(force.getX(), force.getY(), force.getZ());
}
void directInteractionMutual(ParticleClass& Target, ParticleClass& Source) const
{
const FReal one_over_r = MatrixKernel->evaluate(Target.getPosition(), Source.getPosition());
const FReal wt = Target.getPhysicalValue();
const FReal ws = Source.getPhysicalValue();
// potential
Target.incPotential(one_over_r * ws);
Source.incPotential(one_over_r * wt);
// force
F3DPosition force(Source.getPosition() - Target.getPosition());
force *= ((ws*wt) * (one_over_r*one_over_r*one_over_r));
Target.incForces( force.getX(), force.getY(), force.getZ());
Source.incForces((-force.getX()), (-force.getY()), (-force.getZ()));
}
};
......
Src/Kernels/Chebyshev/FChebInterpolator.hpp
View file @ cb59f55d
......@@ -169,6 +169,16 @@ public:
const FReal *const localExpansion,
ContainerClass *const localParticles) const;
/**
* Local to particle operation: application of \f$S_\ell(x,\bar x_m)\f$ and
* \f$\nabla_x S_\ell(x,\bar x_m)\f$ (interpolation)
*/
template <class ContainerClass>
void applyL2PTotal(const F3DPosition& center,
const FReal width,
const FReal *const localExpansion,
ContainerClass *const localParticles) const;
void applyM2M(const unsigned int ChildIndex,
const FReal *const ChildExpansion,
......@@ -305,8 +315,11 @@ inline void FChebInterpolator<ORDER>::applyL2P(const F3DPosition& center,
S[0] *= FReal(2.); S[0] += FReal(1.);
S[1] *= FReal(2.); S[1] += FReal(1.);
S[2] *= FReal(2.); S[2] += FReal(1.);
targetValue += FReal(1.) / nnodes * S[0] * S[1] * S[2] * localExpansion[n];
targetValue += S[0] * S[1] * S[2] * localExpansion[n];
}
// scale
targetValue /= nnodes;
// set potential
iter.data().setPotential(targetValue);
// increment target iterator
......@@ -438,4 +451,136 @@ inline void FChebInterpolator<ORDER>::applyL2PGradient(const F3DPosition& center
}
/**
* Local to particle operation: application of \f$S_\ell(x,\bar x_m)\f$ and
* \f$\nabla_x S_\ell(x,\bar x_m)\f$ (interpolation)
*/
template <int ORDER>
template <class ContainerClass>
inline void FChebInterpolator<ORDER>::applyL2PTotal(const F3DPosition& center,
const FReal width,
const FReal *const localExpansion,
ContainerClass *const localParticles) const
{
// setup local to global mapping
const map_glob_loc map(center, width);
F3DPosition Jacobian;
map.computeJacobian(Jacobian);
const FReal jacobian[3] = {Jacobian.getX(), Jacobian.getY(), Jacobian.getZ()};
F3DPosition localPosition;
FReal T_of_x[ORDER][3];
FReal U_of_x[ORDER][3];
FReal P[3];
typename ContainerClass::BasicIterator iter(*localParticles);
while(iter.hasNotFinished()){
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
// evaluate chebyshev polynomials of source particle
// T_0(x_i) and T_1(x_i)
T_of_x[0][0] = FReal(1.); T_of_x[1][0] = localPosition.getX();
T_of_x[0][1] = FReal(1.); T_of_x[1][1] = localPosition.getY();
T_of_x[0][2] = FReal(1.); T_of_x[1][2] = localPosition.getZ();
// U_0(x_i) and U_1(x_i)
U_of_x[0][0] = FReal(1.); U_of_x[1][0] = localPosition.getX() * FReal(2.);
U_of_x[0][1] = FReal(1.); U_of_x[1][1] = localPosition.getY() * FReal(2.);
U_of_x[0][2] = FReal(1.); U_of_x[1][2] = localPosition.getZ() * FReal(2.);
for (unsigned int o=2; o<ORDER; ++o) {
// T_o(x_i)
T_of_x[o][0] = FReal(2.)*localPosition.getX()*T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = FReal(2.)*localPosition.getY()*T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = FReal(2.)*localPosition.getZ()*T_of_x[o-1][2] - T_of_x[o-2][2];
// U_o(x_i)
U_of_x[o][0] = FReal(2.)*localPosition.getX()*U_of_x[o-1][0] - U_of_x[o-2][0];
U_of_x[o][1] = FReal(2.)*localPosition.getY()*U_of_x[o-1][1] - U_of_x[o-2][1];
U_of_x[o][2] = FReal(2.)*localPosition.getZ()*U_of_x[o-1][2] - U_of_x[o-2][2];
}
// scale, because dT_o/dx = oU_{o-1}
for (unsigned int o=2; o<ORDER; ++o) {
U_of_x[o-1][0] *= FReal(o);
U_of_x[o-1][1] *= FReal(o);
U_of_x[o-1][2] *= FReal(o);
}
// apply P and increment forces
FReal potential = FReal(0.);
FReal forces[3] = {FReal(0.), FReal(0.), FReal(0.)};
for (unsigned int n=0; n<nnodes; ++n) {
// tensor indices of chebyshev nodes
const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};
// f0 component //////////////////////////////////////
P[0] = U_of_x[0][0] * T_of_roots[1][j[0]];
P[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
P[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
// for potential
FReal S0 = T_of_x[1][0] * T_of_roots[1][j[0]];
for (unsigned int o=2; o<ORDER; ++o) {
P[0] += U_of_x[o-1][0] * T_of_roots[o][j[0]];
P[1] += T_of_x[o ][1] * T_of_roots[o][j[1]];
P[2] += T_of_x[o ][2] * T_of_roots[o][j[2]];
// for potential
S0 += T_of_x[o][0] * T_of_roots[o][j[0]];
}
P[0] *= FReal(2.);
P[1] *= FReal(2.); P[1] += FReal(1.);
P[2] *= FReal(2.); P[2] += FReal(1.);
forces[0] += P[0] * P[1] * P[2] * localExpansion[n];
// for potential
S0 *= FReal(2.); S0 += FReal(1.);
potential += FReal(1.) / nnodes * S0 * P[1] * P[2] * localExpansion[n];
// f1 component //////////////////////////////////////
P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
P[1] = U_of_x[0][1] * T_of_roots[1][j[1]];
P[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
for (unsigned int o=2; o<ORDER; ++o) {
P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]];
P[1] += U_of_x[o-1][1] * T_of_roots[o][j[1]];
P[2] += T_of_x[o ][2] * T_of_roots[o][j[2]];
}
P[0] *= FReal(2.); P[0] += FReal(1.);
P[1] *= FReal(2.);
P[2] *= FReal(2.); P[2] += FReal(1.);
forces[1] += P[0] * P[1] * P[2] * localExpansion[n];
// f2 component //////////////////////////////////////
P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
P[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
P[2] = U_of_x[0][2] * T_of_roots[1][j[2]];
for (unsigned int o=2; o<ORDER; ++o) {
P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]];
P[1] += T_of_x[o ][1] * T_of_roots[o][j[1]];
P[2] += U_of_x[o-1][2] * T_of_roots[o][j[2]];
}
P[0] *= FReal(2.); P[0] += FReal(1.);
P[1] *= FReal(2.); P[1] += FReal(1.);
P[2] *= FReal(2.);
forces[2] += P[0] * P[1] * P[2] * localExpansion[n];
}
// scale forces
forces[0] *= jacobian[0] / nnodes;
forces[1] *= jacobian[1] / nnodes;
forces[2] *= jacobian[2] / nnodes;
// set computed potential
iter.data().incPotential(potential);
// set computed forces
iter.data().incForces(forces[0] * iter.data().getPhysicalValue(),
forces[1] * iter.data().getPhysicalValue(),
forces[2] * iter.data().getPhysicalValue());
// increment iterator
iter.gotoNext();
}
}
#endif
Src/Kernels/Chebyshev/FChebKernel.hpp
View file @ cb59f55d
......@@ -165,13 +165,20 @@ public:
M2LHandler->applyU(LeafCell->getLocal() + AbstractBaseClass::nnodes,
const_cast<CellClass*>(LeafCell)->getLocal());
// 2.a) apply Sx
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
// 2.b) apply Px (grad Sx)
//// 2.a) apply Sx
//AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// TargetParticles);
//// 2.b) apply Px (grad Sx)
//AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// TargetParticles);
// 2.c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
......
Src/Kernels/Chebyshev/FChebSymKernel.hpp
View file @ cb59f55d
......@@ -313,17 +313,24 @@ public:
void L2P(const CellClass* const LeafCell,
ContainerClass* const TargetParticles)
{
// a) apply Sx
const F3DPosition LeafCellCenter(getLeafCellCenter(LeafCell->getCoordinate()));
AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
// b) apply Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
// // a) apply Sx
// AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// TargetParticles);
// // b) apply Px (grad Sx)
// AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
// LeafCell->getLocal(),
// TargetParticles);
// c) apply Sx and Px (grad Sx)
AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter,
AbstractBaseClass::BoxWidthLeaf,
LeafCell->getLocal(),
TargetParticles);
}
......
Tests/Kernels/testCompareKernels.cpp
View file @ cb59f55d
......@@ -190,8 +190,8 @@ int main(int argc, char* argv[])
{ // begin Chebyshef kernel
// accuracy
const unsigned int ORDER = 5;
const FReal epsilon = FReal(1e-5);
const unsigned int ORDER = 7;
const FReal epsilon = FReal(1e-7);
FReal* p1; p1 = NULL;
FReal* f1; p1 = NULL;
......
Tests/Utils/testChebSxUCBSy.cpp
View file @ cb59f55d
......@@ -26,6 +26,7 @@
#include "../../Src/Utils/FTic.hpp"
#include "../../Src/Utils/FMath.hpp"
#include "../../Src/Utils/FParameters.hpp"
#include "../../Src/Containers/FVector.hpp"
......@@ -136,7 +137,7 @@ int main(int argc, char* argv[])
////////////////////////////////////////////////////////////////////
// approximative computation
const unsigned int ORDER = 5;
const FReal epsilon = FReal(atof(argv[1]));
const FReal epsilon = FParameters::getValue(argc, argv, "-eps", 1e-5);
const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
typedef FChebInterpolator<ORDER> InterpolatorClass;
InterpolatorClass S;
......@@ -181,7 +182,7 @@ int main(int argc, char* argv[])
*/
// Interpolate f_i = \sum_m^L S(x_i,\bar x_m) * F_m
S.applyL2P(cx, width, F, X.getTargets());
S.applyL2PTotal(cx, width, F, X.getTargets());
time.tac();
std::cout << "Done in " << time.elapsed() << "sec." << std::endl;
......
Supports Markdown
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