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
Floating point robustness fixes in grid deposition.
  • Loading branch information
Andreas Kloeckner committed Mar 28, 2010
1 parent 4f81a9e commit 4e16263
Showing 1 changed file with 22 additions and 7 deletions.
29 changes: 22 additions & 7 deletions pyrticle/deposition/grid.py
Expand Up @@ -175,7 +175,8 @@ def prepare_average_groups(self):
avg_group_finder[idx] = ag

for ag in avg_groups:
self.backend.average_groups.extend(ag)
self.backend.average_groups.extend(
[int(ag_i) for ag_i in ag])
self.backend.average_group_starts.append(
len(self.backend.average_groups))

Expand Down Expand Up @@ -226,11 +227,16 @@ def make_pointwise_interpolation_matrix(self, eog, eg, el, ldis, svd, scaled_vdm
inv_s = numpy.zeros((len(s),), dtype=float)
inv_s[nonzero_flags] = 1/s[nonzero_flags]

if not nonzero_flags.all():
from warnings import warn
warn("grid dep: taking pseudoinverse of poorly conditioned matrix--"
"this should have been caught earlier")

# compute the pseudoinverse of the structured Vandermonde matrix
inv_s_diag = numpy.zeros(
(node_count, point_count),
dtype=float)
inv_s_diag[:len(s),:len(s)] = numpy.diag(1/s)
inv_s_diag[:len(s),:len(s)] = numpy.diag(inv_s)

svdm_pinv = numpy.dot(numpy.dot(vt.T, inv_s_diag), u.T)

Expand All @@ -254,6 +260,11 @@ def make_pointwise_interpolation_matrix(self, eog, eg, el, ldis, svd, scaled_vdm
if self.filter is not None:
imat = numpy.dot(self.filter.get_filter_matrix(eg), imat)

assert not numpy.isnan(imat).any(), \
"encountered NaN in element %d's interpolation matrix" % el.id
assert not numpy.isinf(imat).any(), \
"encountered infinity in element %d's interpolation matrix" % el.id

eog.interpolation_matrix = numpy.asarray(imat, order="F")

if basis_subset is None:
Expand Down Expand Up @@ -585,7 +596,10 @@ def prepare_with_pointwise_projection_and_basis_reduction(self):
assert s[-1] == numpy.min(s)
assert s[0] == numpy.max(s)

if len(basis) > len(points) or s[0]/s[-1] > 10:
# Imagine that--s can have negative entries.
# (AK: I encountered one negative zero.)

if len(basis) > len(points) or numpy.abs(s[0]/s[-1]) > 10:
retry = True

# badly conditioned, kill a basis entry
Expand Down Expand Up @@ -620,8 +634,8 @@ def prepare_with_pointwise_projection_and_basis_reduction(self):
% el.id)

if ldis.node_count() > len(basis):
print "element %d: #nodes=%d, leftover modes=%d" % (
el.id, ldis.node_count(), len(basis),)
print "element %d: #nodes=%d, killed modes=%d" % (
el.id, ldis.node_count(), ldis.node_count()-len(basis),)

basis_len_vec[discr.find_el_range(el.id)] = len(basis)
el_condition_vec[discr.find_el_range(el.id)] = s[0]/s[-1]
Expand Down Expand Up @@ -794,8 +808,9 @@ def _deposit_densites(self, state, velocities, pslice):
state, velocities, pslice))

def _deposit_j(self, state, velocities, pslice):
return self.remap_grid_to_mesh(self.deposit_grid_j(
state, velocities, pslice))
grid_j = self.deposit_grid_j(state, velocities, pslice)

return self.remap_grid_to_mesh(grid_j)

def _deposit_rho(self, state, pslice):
return self.remap_grid_to_mesh(
Expand Down

0 comments on commit 4e16263

Please sign in to comment.