Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • mosaic/software/parallel-computing/openfpm/openfpm_numerics
  • argupta/openfpm_numerics
2 results
Show changes
Commits on Source (1)
Hide whitespace changes
Inline Side-by-side
src/DCPSE/Dcpse.cuh
View file @ 5de43ebb
......@@ -4,11 +4,10 @@
#ifndef OPENFPM_PDATA_DCPSE_CUH
#define OPENFPM_PDATA_DCPSE_CUH
#if defined(__NVCC__) && defined(HAVE_EIGEN)
#if defined(__NVCC__)
#include "Vector/vector_dist.hpp"
#include "MonomialBasis.hpp"
#include "DMatrix/EMatrix.hpp"
#include "SupportBuilder.hpp"
#include "SupportBuilder.cuh"
#include "Support.hpp"
......@@ -103,11 +102,7 @@ public:
opt(opt)
{
particles.ghost_get_subset();
if (supportSizeFactor < 1)
initializeAdaptive(particles, convergenceOrder, rCut);
else
initializeStaticSize(particles, convergenceOrder, rCut, supportSizeFactor);
initializeStaticSize(particles, convergenceOrder, rCut, supportSizeFactor);
}
Dcpse_gpu(vector_type &particles,
......@@ -130,11 +125,7 @@ public:
isSharedSupport(true)
{
particles.ghost_get_subset();
if (supportSizeFactor < 1)
initializeAdaptive(particles, convergenceOrder, rCut);
else
initializeStaticSize(particles, convergenceOrder, rCut, supportSizeFactor);
initializeStaticSize(particles, convergenceOrder, rCut, supportSizeFactor);
}
template<unsigned int prp>
......@@ -576,79 +567,6 @@ public:
}
private:
template <typename U>
void initializeAdaptive(vector_type &particles,
unsigned int convergenceOrder,
U rCut) {
// Still need to be tested
#ifdef SE_CLASS1
this->update_ctr=particles.getMapCtr();
#endif
if (!isSharedSupport) {
subsetKeyPid.resize(particles.size_local_orig());
supportRefs.resize(particles.size_local());
}
localEps.resize(particles.size_local());
localEpsInvPow.resize(particles.size_local());
kerOffsets.resize(particles.size_local()+1);
const T condVTOL = 1e2;
if (!isSharedSupport) {
SupportBuilder<vector_type,vector_type>
supportBuilder(particles, particles, differentialSignature, rCut, differentialOrder == 0);
unsigned int requiredSupportSize = monomialBasis.size();
// need to resize supportKeys1D to yet unknown supportKeysTotalN
// add() takes too long
openfpm::vector<openfpm::vector<size_t>> tempSupportKeys(supportRefs.size());
auto it = particles.getDomainIterator();
while (it.isNext()) {
auto key_o = particles.getOriginKey(it.get());
subsetKeyPid.get(key_o.getKey()) = it.get().getKey();
Support support = supportBuilder.getSupport(it, requiredSupportSize, opt);
supportRefs.get(key_o.getKey()) = key_o.getKey();
tempSupportKeys.get(key_o.getKey()) = support.getKeys();
kerOffsets.get(key_o.getKey()) = supportKeysTotalN;
if (maxSupportSize < support.size())
maxSupportSize = support.size();
supportKeysTotalN += support.size();
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
// Vandermonde matrix computation
Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
vandermonde(support, monomialBasis, particles, particles,HOverEpsilon);
vandermonde.getMatrix(V);
T condV = conditionNumber(V, condVTOL);
T eps = vandermonde.getEps();
if (condV > condVTOL) {
requiredSupportSize *= 2;
std::cout << "INFO: Increasing, requiredSupportSize = " << requiredSupportSize << std::endl; // debug
continue;
} else requiredSupportSize = monomialBasis.size();
++it;
}
kerOffsets.get(supportRefs.size()) = supportKeysTotalN;
supportKeys1D.resize(supportKeysTotalN);
size_t offset = 0;
for (size_t i = 0; i < tempSupportKeys.size(); ++i)
for (size_t j = 0; j < tempSupportKeys.get(i).size(); ++j, ++offset)
supportKeys1D.get(offset) = tempSupportKeys.get(i).get(j);
}
kerOffsets.hostToDevice(); supportKeys1D.hostToDevice();
assembleLocalMatrices_t(rCut);
}
template <typename U>
void initializeStaticSize(vector_type &particles,
unsigned int convergenceOrder,
......
src/DCPSE/Dcpse.hpp
View file @ 5de43ebb
......@@ -184,14 +184,7 @@ public:
opt(opt)
{
particles.ghost_get_subset(); // This communicates which ghost particles to be excluded from support
if (supportSizeFactor < 1)
{
initializeAdaptive(particles, particles, convergenceOrder, rCut);
}
else
{
initializeStaticSize(particles, particles, convergenceOrder, rCut, supportSizeFactor);
}
initializeStaticSize(particles, particles, convergenceOrder, rCut, supportSizeFactor);
}
//Surface DCPSE Constructor
......@@ -259,11 +252,9 @@ public:
isSharedLocalSupport(true)
{
particles.ghost_get_subset();
if (supportSizeFactor < 1)
initializeAdaptive(particles, particles, convergenceOrder, rCut);
else
initializeStaticSize(particles, particles, convergenceOrder, rCut, supportSizeFactor);
initializeStaticSize(particles, particles, convergenceOrder, rCut, supportSizeFactor);
}
Dcpse(vector_type &particlesFrom,vector_type2 &particlesTo,
Point<dim, unsigned int> differentialSignature,
unsigned int convergenceOrder,
......@@ -277,10 +268,7 @@ public:
opt(opt)
{
particlesFrom.ghost_get_subset();
if (supportSizeFactor < 1)
initializeAdaptive(particlesFrom,particlesTo,convergenceOrder, rCut);
else
initializeStaticSize(particlesFrom,particlesTo,convergenceOrder, rCut, supportSizeFactor);
initializeStaticSize(particlesFrom,particlesTo,convergenceOrder, rCut, supportSizeFactor);
}
......@@ -762,95 +750,6 @@ public:
}
private:
void initializeAdaptive(vector_type &particlesFrom,
vector_type2 &particlesTo,
unsigned int convergenceOrder,
T rCut) {
SupportBuilder<vector_type,vector_type2>
supportBuilder(particlesFrom, particlesTo, differentialSignature, rCut, differentialOrder == 0);
unsigned int requiredSupportSize = monomialBasis.size();
if (!isSharedLocalSupport)
localSupports.resize(particlesTo.size_local_orig());
localEps.resize(particlesTo.size_local_orig());
localEpsInvPow.resize(particlesTo.size_local_orig());
kerOffsets.resize(particlesTo.size_local_orig());
kerOffsets.fill(-1);
auto it = particlesTo.getDomainIterator();
while (it.isNext()) {
const T condVTOL = 1e2;
auto key_o = particlesTo.getOriginKey(it.get());
if (!isSharedLocalSupport)
localSupports.get(key_o.getKey()) = supportBuilder.getSupport(it, requiredSupportSize,opt);
Support& support = localSupports.get(key_o.getKey());
// Get the points in the support of the DCPSE kernel and store the support for reuse
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> V(support.size(), monomialBasis.size());
// Vandermonde matrix computation
Vandermonde<dim, T, EMatrix<T, Eigen::Dynamic, Eigen::Dynamic>>
vandermonde(support, monomialBasis,particlesFrom, particlesTo,HOverEpsilon);
vandermonde.getMatrix(V);
T eps = vandermonde.getEps();
if (!isSharedLocalSupport) {
T condV = conditionNumber(V, condVTOL);
if (condV > condVTOL) {
requiredSupportSize *= 2;
std::cout
<< "INFO: Increasing, requiredSupportSize = " << requiredSupportSize
<< std::endl; // debug
continue;
} else
requiredSupportSize = monomialBasis.size();
}
localSupports.get(key_o.getKey()) = support;
localEps.get(key_o.getKey()) = eps;
localEpsInvPow.get(key_o.getKey()) = 1.0 / openfpm::math::intpowlog(eps,differentialOrder);
// Compute the diagonal matrix E
DcpseDiagonalScalingMatrix<dim> diagonalScalingMatrix(monomialBasis);
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> E(support.size(), support.size());
diagonalScalingMatrix.buildMatrix(E, support, eps, particlesFrom, particlesTo);
// Compute intermediate matrix B
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> B = E * V;
// Compute matrix A
EMatrix<T, Eigen::Dynamic, Eigen::Dynamic> A = B.transpose() * B;
// Compute RHS vector b
DcpseRhs<dim> rhs(monomialBasis, differentialSignature);
EMatrix<T, Eigen::Dynamic, 1> b(monomialBasis.size(), 1);
rhs.template getVector<T>(b);
// Get the vector where to store the coefficients...
EMatrix<T, Eigen::Dynamic, 1> a(monomialBasis.size(), 1);
// ...solve the linear system...
a = A.colPivHouseholderQr().solve(b);
// ...and store the solution for later reuse
kerOffsets.get(key_o.getKey()) = calcKernels.size();
Point<dim, T> xp = particlesTo.getPosOrig(key_o);
const auto& support_keys = support.getKeys();
size_t N = support_keys.size();
for (size_t i = 0; i < N; ++i)
{
const auto& xqK = support_keys.get(i);
Point<dim, T> xq = particlesFrom.getPosOrig(xqK);
Point<dim, T> normalizedArg = (xp - xq) / eps;
calcKernels.add(computeKernel(normalizedArg, a));
}
//
++it;
}
}
void initializeStaticSize(vector_type &particlesFrom,vector_type2 &particlesTo,
unsigned int convergenceOrder,
T rCut,
......