Skip to content
Snippets Groups Projects
Commit 52a97236 authored by Andreas Klöckner's avatar Andreas Klöckner
Browse files

Fix direct-QBX for rscale

parent 5a15a2c9
No related branches found
No related tags found
1 merge request!46Scale expansion coefficients with expansion radius
Pipeline #
...@@ -61,7 +61,9 @@ def stringify_expn_index(i): ...@@ -61,7 +61,9 @@ def stringify_expn_index(i):
def expand(expansion_nr, sac, expansion, avec, bvec): def expand(expansion_nr, sac, expansion, avec, bvec):
coefficients = expansion.coefficients_from_source(avec, bvec) rscale = sym.Symbol("rscale")
coefficients = expansion.coefficients_from_source(avec, bvec, rscale)
assigned_coeffs = [ assigned_coeffs = [
sym.Symbol( sym.Symbol(
...@@ -71,7 +73,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec): ...@@ -71,7 +73,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec):
for i in expansion.get_coefficient_identifiers()] for i in expansion.get_coefficient_identifiers()]
return sac.assign_unique("expn%d_result" % expansion_nr, return sac.assign_unique("expn%d_result" % expansion_nr,
expansion.evaluate(assigned_coeffs, bvec)) expansion.evaluate(assigned_coeffs, bvec, rscale))
# {{{ layer potential computation # {{{ layer potential computation
...@@ -101,6 +103,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): ...@@ -101,6 +103,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper):
return """ return """
<> a[idim] = center[idim,itgt] - src[idim,isrc] {dup=idim} <> a[idim] = center[idim,itgt] - src[idim,isrc] {dup=idim}
<> b[idim] = tgt[idim,itgt] - center[idim,itgt] {dup=idim} <> b[idim] = tgt[idim,itgt] - center[idim,itgt] {dup=idim}
<> rscale = expansion_radii[itgt]
""" """
def get_src_tgt_arguments(self): def get_src_tgt_arguments(self):
...@@ -111,6 +114,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): ...@@ -111,6 +114,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper):
shape=(self.dim, "ntargets"), order="C"), shape=(self.dim, "ntargets"), order="C"),
lp.GlobalArg("center", None, lp.GlobalArg("center", None,
shape=(self.dim, "ntargets"), order="C"), shape=(self.dim, "ntargets"), order="C"),
lp.GlobalArg("expansion_radii", None, shape="ntargets"),
lp.ValueArg("nsources", None), lp.ValueArg("nsources", None),
lp.ValueArg("ntargets", None), lp.ValueArg("ntargets", None),
] ]
...@@ -242,7 +246,8 @@ class LayerPotential(LayerPotentialBase): ...@@ -242,7 +246,8 @@ class LayerPotential(LayerPotentialBase):
for iknl in range(len(self.expansions)) for iknl in range(len(self.expansions))
] ]
def __call__(self, queue, targets, sources, centers, strengths, **kwargs): def __call__(self, queue, targets, sources, centers, strengths, expansion_radii,
**kwargs):
""" """
:arg strengths: are required to have area elements and quadrature weights :arg strengths: are required to have area elements and quadrature weights
already multiplied in. already multiplied in.
...@@ -253,7 +258,8 @@ class LayerPotential(LayerPotentialBase): ...@@ -253,7 +258,8 @@ class LayerPotential(LayerPotentialBase):
for i, dens in enumerate(strengths): for i, dens in enumerate(strengths):
kwargs["strength_%d" % i] = dens kwargs["strength_%d" % i] = dens
return knl(queue, src=sources, tgt=targets, center=centers, **kwargs) return knl(queue, src=sources, tgt=targets, center=centers,
expansion_radii=expansion_radii, **kwargs)
# }}} # }}}
...@@ -283,7 +289,7 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase): ...@@ -283,7 +289,7 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase):
for iknl in range(len(self.expansions)) for iknl in range(len(self.expansions))
] ]
def __call__(self, queue, targets, sources, centers, **kwargs): def __call__(self, queue, targets, sources, centers, expansion_radii, **kwargs):
knl = self.get_optimized_kernel() knl = self.get_optimized_kernel()
return knl(queue, src=sources, tgt=targets, center=centers, return knl(queue, src=sources, tgt=targets, center=centers,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment