Skip to content
This repository has been archived by the owner on Mar 8, 2021. It is now read-only.

Commit

Permalink
Browse files Browse the repository at this point in the history
NaN debugging code for polynomial pushing.
  • Loading branch information
Andreas Kloeckner committed Mar 28, 2010
1 parent 57f20eb commit 8fec889
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 4 deletions.
51 changes: 47 additions & 4 deletions src/cpp/push_monomial.hpp
Expand Up @@ -110,6 +110,21 @@ namespace pyrticle
);
}

template <class VecType>
void debug(
const particle_number pn,
const mesh_data::mesh_data::element_number en,
const VecType &data) const
{
const unsigned basis_size = m_ldis.m_basis.size();
dyn_vector intp_coeffs(subrange(m_interpolation_coefficients,
pn*basis_size, (pn+1)*basis_size));
dyn_vector ldata(subrange(data,
en*basis_size, (en+1)*basis_size));
std::cout << "INTP" << intp_coeffs << std::endl;
std::cout << "DATA" << ldata << std::endl;
}

const double operator()(
const particle_number pn,
const mesh_data::mesh_data::element_number en,
Expand Down Expand Up @@ -249,14 +264,42 @@ namespace pyrticle

const double charge = ps.charges[pn];

bounded_vector el_force(3);
el_force[0] = charge*e[0];
el_force[1] = charge*e[1];
el_force[2] = charge*e[2];
bounded_vector el_force(charge*e);

const bounded_vector v = subrange(velocities, v_pstart, v_pend);
bounded_vector mag_force = cross(v, charge*b);

#if 0
// code for debugging NaNs
if (isnan_any(el_force) || isnan_any(mag_force))
{
const unsigned x_pstart = ps.xdim()*pn;
const unsigned x_pend = ps.xdim()*(pn+1);

const bounded_vector x = subrange(ps.positions, x_pstart, x_pend);

if (isnan_any(el_force))
{
std::cout << "EL FORCE HAD NAN" << std::endl;
interp.debug(pn, in_el, ex);
interp.debug(pn, in_el, ey);
interp.debug(pn, in_el, ez);
}

if (isnan_any(mag_force))
{
std::cout << "MAG FORCE HAD NAN" << std::endl;
interp.debug(pn, in_el, bx);
interp.debug(pn, in_el, by);
interp.debug(pn, in_el, bz);
}

std::cout
<< el_force << " pn " << pn << " at " << x << " in el " << in_el
<< std::endl;
}
#endif

// truncate forces to dimensions_velocity entries
subrange(result, v_pstart, v_pend) = subrange(
el_force + mag_force, 0, ps.vdim());
Expand Down
12 changes: 12 additions & 0 deletions src/cpp/tools.hpp
Expand Up @@ -324,6 +324,18 @@ namespace pyrticle



template <class VectorType>
bool isnan_any(VectorType const &vec)
{
BOOST_FOREACH(typename VectorType::value_type value, vec)
if (isnan(value))
return true;
return false;
}




template <class T>
class stats_gatherer
{
Expand Down

0 comments on commit 8fec889

Please sign in to comment.