quadrature comparison #255
Conversation
alexfikl
left a comment
alexfikl
left a comment
There was a problem hiding this comment.
This is very cool! Just ran it locally to see the errors and there's a pretty big mismatch :\
Left some comments in there, but there's no big issues, just nitpicks. I imagine ruff is going to complain about some things, so might be worth running it locally until the CI recovers.
That's great! Looks good to me now 😁 EDIT: Hm, no idea what those segmentation faults are about, but they don't seem related to this PR :\ |
Me either. What's weirder is that I tried to reproduce them in an interactive tmate session (#256), and no dice. Maybe conda-forge was temporarily in a weird state? |
|
It seems that jemalloc was at least involved in the problem. I've disabled it. |
| # VR quadrature | ||
| qbx_order = vr_target_order | ||
| sources_h, weights_nodes_h = quadrature(level, | ||
| vr_target_order, qbx_order=qbx_order) | ||
| vr_value = np.dot(wave(sources_h), weights_nodes_h) | ||
| vr_result.append(vr_value) | ||
| vr_nodes.append(len(sources_h[0])) | ||
| vr_err.append(np.abs(vr_value - ref)) | ||
|
|
||
| # Tensor quadrature | ||
| qbx_order = tensor_target_order | ||
| sources_h, weights_nodes_h = quadrature(level, | ||
| tensor_target_order, qbx_order=qbx_order, | ||
| group_cls=TensorProductElementGroup) | ||
| tensor_value = np.dot(wave(sources_h), weights_nodes_h) | ||
| tensor_result.append(tensor_value) | ||
| tensor_nodes.append(len(sources_h[0])) | ||
| tensor_err.append(np.abs(tensor_value - ref)) |
There was a problem hiding this comment.
This should be a loop (since it's doing the same thing twice).
| n = Normalize[{-5, 4, 1}]; | ||
| r = 1; | ||
| wave[θ_, φ_] := | ||
| Exp[I r (n[[1]] Sin[θ] Cos[φ] + | ||
| n[[2]] Sin[θ] Sin[φ] + | ||
| n[[3]] Cos[θ])]; | ||
|
|
||
| NIntegrate[ | ||
| wave[θ, φ] * r^2 Sin[θ], | ||
| {θ, 0, Pi}, | ||
| {φ, 0, 2 Pi}, | ||
| WorkingPrecision -> 21 | ||
| ] |
There was a problem hiding this comment.
Could you try to avoid using Mathematica? The main problem is that I can't reproduce the result, or tweak the computation in any way.
https://mpmath.org/doc/current/calculus/integration.html#mpmath.quad should work just as well.
| actx = PyOpenCLArrayContext(queue) | ||
|
|
||
|
|
||
| def quadrature(level, target_order, qbx_order, group_cls=SimplexElementGroup): |
There was a problem hiding this comment.
| def quadrature(level, target_order, qbx_order, group_cls=SimplexElementGroup): | |
| def quadrature(level, target_order, qbx_order, group_cls): |
| Order VR Exact_to Tensor Exact_to | ||
| ----------------------------------------- | ||
| 1 2 3 | ||
| 2 4 5 | ||
| 3 5 7 | ||
| 4 7 9 | ||
| 5 8 11 | ||
| 6 10 13 | ||
| 7 12 15 | ||
| 8 14 17 | ||
| 9 15 19 | ||
| 10 17 21 | ||
| 11 19 23 | ||
| 12 20 25 | ||
| 13 22 27 | ||
| 14 24 29 | ||
| 15 25 31 | ||
| 16 27 33 | ||
| 17 28 35 | ||
| 18 30 37 | ||
| 19 32 39 |
There was a problem hiding this comment.
It's not useful to incorporate this table, since it's not guaranteed to match the code on an ongoing basis. Instead, the code should generate the table.
| actx = PyOpenCLArrayContext(queue) | ||
|
|
||
|
|
||
| def quadrature(level, target_order, qbx_order, group_cls=SimplexElementGroup): |
There was a problem hiding this comment.
The way this is currently expressed is flawed: At no point do you explicitly specify which quadrature should be used. And what pytential uses internally is an implementation detail that could change at any time. If you're intending to measure differences between specific quadrature rules, your code should exert positive control over what rule is used, instead of relying on implementation coincidence. It may be easier to do this by using meshmode primitives.
| qbx = QBXLayerPotentialSource( | ||
| pre_density_discr, target_order, qbx_order, | ||
| fmm_order=False) | ||
| discr_stage = sym.QBX_SOURCE_STAGE1 |
There was a problem hiding this comment.
Using pytential's stage-1 is a potential confounder, since it may apply refinement to the incoming discretization.
|
|
||
| sources = density_discr.nodes() | ||
| weights_nodes = bind(places, | ||
| sym.weights_and_area_elements(ambient_dim=3, dim=2, dofdesc=dofdesc))(actx) |
There was a problem hiding this comment.
Please document how you convinced yourself that the error is not dominated by geometry derivatives.
3 participants
No description provided.