Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use hcubature_buffer for decomposed geometries #157

Closed
wants to merge 15 commits into from
Closed

Use hcubature_buffer for decomposed geometries #157

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
2f9d0d0
Use hcubature with buffer for Rope
mikeingold Dec 19, 2024
c0bff58
Add type param
mikeingold Dec 19, 2024
179da14
Fix typo
mikeingold Dec 19, 2024
1727b79
Use a sample integrand to inform buffer
mikeingold Dec 19, 2024
abdb0ce
Update name
mikeingold Dec 19, 2024
bbb6fd6
Update accessor to generator
mikeingold Dec 19, 2024
b37cff0
Also benchmark HCubature rules
mikeingold Dec 19, 2024
a04b3ba
Generate geometry names
mikeingold Dec 19, 2024
b09776c
Drop field name
mikeingold Dec 19, 2024
795be4c
Try without explicit rule names
mikeingold Dec 19, 2024
2cd7b0c
Update rule name
mikeingold Dec 19, 2024
0f9f110
Use interp
mikeingold Dec 19, 2024
ca1c9ab
Change dict structures
mikeingold Dec 19, 2024
53e86a8
Update naming
mikeingold Dec 19, 2024
2e1038d
Implement for CylinderSurface
mikeingold Dec 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 18 additions & 15 deletions benchmark/benchmarks.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,8 @@ end
spec = (
f = p -> norm(to(p)),
f_exp = p::Point -> exp(-norm(to(p))^2 / u"m^2"),
g = (
geometries = (
bezier = BezierCurve([Point(t, sin(t), 0) for t in range(-pi, pi, length = 361)]),
line = Line(Point(0, 0, 0), Point(1, 1, 1)),
plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)),
ray = Ray(Point(0, 0, 0), Vec(0, 0, 1)),
rope = Rope([Point(t, t, t) for t in 1:32]...),
triangle = Triangle(Point(1, 0, 0), Point(0, 1, 0), Point(0, 0, 1)),
tetrahedron = let
Expand All @@ -56,19 +53,25 @@ spec = (
Tetrahedron(a, b, c, a + ẑ)
end
),
rule_gl = GaussLegendre(100),
rule_gk = GaussKronrod(),
rule_hc = HAdaptiveCubature()
geometries_exp = (
line = Line(Point(0, 0, 0), Point(1, 1, 1)),
plane = Plane(Point(0, 0, 0), Vec(0, 0, 1)),
ray = Ray(Point(0, 0, 0), Vec(0, 0, 1))
),
rules = (
( name = "GaussLegendre", rule = GaussLegendre(100) ),
( name = "HAdaptiveCubature", rule = HAdaptiveCubature() )
)
)

SUITE["Specializations/Scalar GaussLegendre"] = let s = BenchmarkGroup()
s["BezierCurve"] = @benchmarkable integral($spec.f, $spec.g.bezier, $spec.rule_gl)
s["Line"] = @benchmarkable integral($spec.f_exp, $spec.g.line, $spec.rule_gl)
s["Plane"] = @benchmarkable integral($spec.f_exp, $spec.g.plane, $spec.rule_gl)
s["Ray"] = @benchmarkable integral($spec.f_exp, $spec.g.ray, $spec.rule_gl)
s["Rope"] = @benchmarkable integral($spec.f, $spec.g.rope, $spec.rule_gl)
s["Triangle"] = @benchmarkable integral($spec.f, $spec.g.triangle, $spec.rule_gl)
s["Tetrahedron"] = @benchmarkable integral($spec.f, $spec.g.tetrahedron, $spec.rule_gl)
SUITE["Specializations"] = let s = BenchmarkGroup()
s[]
for r in spec.rules, geometry in spec.geometries
s[geometry.name, r.name] = @benchmarkable integral($spec.f, geometry, r.rule)
end
for r in spec.rules, geometry in spec.geometries_exp
s[geometry.name, r.name] = @benchmarkable integral($spec.f_exp, geometry, r.rule)
end
s
end

Expand Down
24 changes: 24 additions & 0 deletions src/specializations/Rope.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,3 +35,27 @@ function integral(
# Convert the Rope into Segments, sum the integrals of those
return sum(segment -> integral(f, segment, rule; kwargs...), Meshes.segments(rope))
end

# Use HCubature.hcubature_buffer to reduce allocations
function integral(
f,
rope::Meshes.Rope,
rule::HAdaptiveCubature;
FP::Type{T} = Float64,
kwargs...
Comment on lines +41 to +45
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[JuliaFormatter] reported by reviewdog 🐶

Suggested change
f,
rope::Meshes.Rope,
rule::HAdaptiveCubature;
FP::Type{T} = Float64,
kwargs...
f,
rope::Meshes.Rope,
rule::HAdaptiveCubature;
FP::Type{T} = Float64,
kwargs...

) where {T <: AbstractFloat}
# Geometry information
N = Meshes.paramdim(rope)
segments = Meshes.segments(rope)

# Use a sample integrand to develop and append a buffer to the given rule
sample = first(segments)
integrand(ts) = f(sample(ts...)) * differential(sample, ts)
uintegrand(ts) = Unitful.ustrip.(integrand(ts))
buffer = HCubature.hcubature_buffer(uintegrand, _zeros(FP, N), _ones(FP, N))
rule = HAdaptiveCubature(rule.kwargs..., buffer = buffer)

# Convert the Rope into Segments, sum the integrals of those
_subintegral(seg) = _integral(f, seg, rule; FP = FP, kwargs...)
return sum(_subintegral, segments)
end
Loading