Hermite conditions

[Davide-Miotti]

here we go!

[github-actions]

Benchmark Results

	main	b6d005e...	main / b6d005e...
Directional	2.4 ± 0.073 ms	2.37 ± 0.072 ms	1.01
Directional (per point)	2.4 ± 0.076 ms	2.37 ± 0.078 ms	1.01
Gradient	7.61 ± 0.23 ms	7.37 ± 0.26 ms	1.03
MonomialBasis/dim=1/deg=0	0.0461 ± 0.012 μs	0.0405 ± 0.012 μs	1.14
MonomialBasis/dim=1/deg=1	0.0794 ± 0.012 μs	0.0768 ± 0.014 μs	1.03
MonomialBasis/dim=1/deg=2	0.0673 ± 0.018 μs	0.0856 ± 0.019 μs	0.786
MonomialBasis/dim=2/deg=0	25 ± 9.8 ns	0.0367 ± 0.011 μs	0.682
MonomialBasis/dim=2/deg=1	0.0407 ± 0.012 μs	0.0355 ± 0.012 μs	1.15
MonomialBasis/dim=2/deg=2	0.047 ± 0.012 μs	0.0487 ± 0.011 μs	0.965
MonomialBasis/dim=3/deg=0	0.0408 ± 0.012 μs	0.0513 ± 0.012 μs	0.794
MonomialBasis/dim=3/deg=1	0.0468 ± 0.012 μs	0.0363 ± 0.012 μs	1.29
MonomialBasis/dim=3/deg=2	0.0491 ± 0.012 μs	0.0432 ± 0.013 μs	1.14
Partial	2.42 ± 0.1 ms	2.38 ± 0.092 ms	1.01
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∂	10.3 ± 0.051 ns	10.1 ± 0.09 ns	1.01
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∂²	10.6 ± 0.071 ns	10.6 ± 0.09 ns	0.993
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∇	17.3 ± 0.14 ns	17.3 ± 0.061 ns	1
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∇²	18.2 ± 0.07 ns	18.2 ± 0.07 ns	1
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∂	10.3 ± 0.07 ns	10.1 ± 0.09 ns	1.01
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∂²	10.6 ± 0.06 ns	10.6 ± 0.061 ns	0.992
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∇	17.3 ± 0.14 ns	17.3 ± 0.061 ns	1
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∇²	18.2 ± 0.08 ns	18.3 ± 0.15 ns	0.996
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∂	10.3 ± 0.08 ns	10.2 ± 0.19 ns	1
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∂²	10.6 ± 0.05 ns	10.6 ± 0.12 ns	0.997
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∇	17.3 ± 0.06 ns	17.3 ± 0.061 ns	1
RBF/Gaussian, exp(-(ε*r)²)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∇²	18.3 ± 0.05 ns	18.3 ± 0.14 ns	1
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∂	6.8 ± 0.09 ns	6.74 ± 0.15 ns	1.01
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∂²	14 ± 0.08 ns	13.9 ± 0.09 ns	1.01
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∇	8.36 ± 0.1 ns	8.69 ± 0.05 ns	0.962
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 0/0/∇²	16.3 ± 0.071 ns	16.3 ± 0.22 ns	1
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∂	6.79 ± 0.15 ns	6.66 ± 0.089 ns	1.02
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∂²	14 ± 0.08 ns	13.9 ± 0.09 ns	1.01
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∇	8.36 ± 0.11 ns	8.69 ± 0.06 ns	0.963
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 1/1/∇²	16.3 ± 0.08 ns	16.3 ± 0.089 ns	0.997
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∂	6.76 ± 0.15 ns	6.72 ± 0.16 ns	1.01
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∂²	14 ± 0.081 ns	13.9 ± 0.09 ns	1.01
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∇	8.36 ± 0.11 ns	8.69 ± 0.051 ns	0.963
RBF/Inverse Multiquadrics, 1/sqrt((r*ε)²+1)
├─Shape factor: ε = 1
└─Polynomial augmentation: degree 2/2/∇²	16.3 ± 0.07 ns	16.3 ± 0.11 ns	1
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 0/0/∂	3.73 ± 0.01 ns	3.42 ± 0.001 ns	1.09
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 0/0/∂²	5.27 ± 0.01 ns	4.7 ± 0.01 ns	1.12
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 0/0/∇	5.49 ± 0.03 ns	5.43 ± 0.02 ns	1.01
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 0/0/∇²	7 ± 0.04 ns	3.11 ± 0 ns	2.25
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 1/1/∂	3.73 ± 0.01 ns	3.42 ± 0.001 ns	1.09
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 1/1/∂²	5.27 ± 0.01 ns	4.7 ± 0.01 ns	1.12
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 1/1/∇	5.49 ± 0.03 ns	5.45 ± 4.5 ns	1.01
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 1/1/∇²	7 ± 0.05 ns	3.11 ± 0 ns	2.25
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 2/2/∂	3.73 ± 0.01 ns	3.42 ± 0.01 ns	1.09
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 2/2/∂²	5.27 ± 0.01 ns	4.7 ± 0.01 ns	1.12
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 2/2/∇	5.49 ± 0.021 ns	5.43 ± 0.02 ns	1.01
RBF/Polyharmonic spline (r³)
└─Polynomial augmentation: degree 2/2/∇²	7 ± 0.06 ns	3.11 ± 0 ns	2.25
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 0/0/∂	4.27 ± 0.01 ns	4.27 ± 0.01 ns	1
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 0/0/∂²	5.61 ± 0.061 ns	5.8 ± 0.011 ns	0.967
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 0/0/∇	6.38 ± 0.039 ns	6.28 ± 0.03 ns	1.02
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 0/0/∇²	7.07 ± 0.34 ns	4.5 ± 0.01 ns	1.57
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 1/1/∂	4.27 ± 0.01 ns	4.27 ± 0.01 ns	1
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 1/1/∂²	5.59 ± 0.07 ns	5.8 ± 0.011 ns	0.964
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 1/1/∇	6.38 ± 0.03 ns	6.28 ± 0.039 ns	1.02
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 1/1/∇²	7.07 ± 0.34 ns	4.5 ± 0.01 ns	1.57
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 2/2/∂	4.27 ± 0.01 ns	4.27 ± 0.01 ns	1
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 2/2/∂²	5.6 ± 0.07 ns	5.8 ± 0.02 ns	0.966
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 2/2/∇	6.38 ± 0.039 ns	6.28 ± 0.04 ns	1.02
RBF/Polyharmonic spline (r¹)
└─Polynomial augmentation: degree 2/2/∇²	7.07 ± 0.34 ns	4.5 ± 0.01 ns	1.57
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 0/0/∂	4.96 ± 0.01 ns	4.96 ± 0.01 ns	1
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 0/0/∂²	4.96 ± 0.01 ns	4.96 ± 0.01 ns	1
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 0/0/∇	6.02 ± 0.11 ns	5.93 ± 0.12 ns	1.02
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 0/0/∇²	5 ± 0.09 ns	3.11 ± 0 ns	1.61
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 1/1/∂	4.96 ± 0.01 ns	4.96 ± 0.001 ns	1
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 1/1/∂²	4.96 ± 0.01 ns	4.96 ± 0.01 ns	1
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 1/1/∇	6.02 ± 0.11 ns	5.92 ± 0.11 ns	1.02
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 1/1/∇²	4.99 ± 0.042 ns	3.11 ± 0 ns	1.61
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 2/2/∂	4.96 ± 0.01 ns	4.96 ± 0.001 ns	1
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 2/2/∂²	4.96 ± 0.01 ns	4.96 ± 0.01 ns	1
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 2/2/∇	6.02 ± 0.11 ns	5.93 ± 0.12 ns	1.02
RBF/Polyharmonic spline (r⁵)
└─Polynomial augmentation: degree 2/2/∇²	4.98 ± 0.05 ns	3.11 ± 0 ns	1.6
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 0/0/∂	10.5 ± 0.011 ns	10.5 ± 0.07 ns	1
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 0/0/∂²	4.96 ± 0.001 ns	5.27 ± 0.01 ns	0.941
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 0/0/∇	12.4 ± 0.12 ns	12.3 ± 0.04 ns	1.01
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 0/0/∇²	5.47 ± 0.019 ns	8.18 ± 0.021 ns	0.668
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 1/1/∂	10.9 ± 0.089 ns	10.5 ± 0.02 ns	1.04
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 1/1/∂²	4.96 ± 0.001 ns	5.26 ± 0.01 ns	0.943
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 1/1/∇	12.5 ± 0.11 ns	12.3 ± 0.04 ns	1.01
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 1/1/∇²	5.47 ± 0.019 ns	8.18 ± 0.021 ns	0.668
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 2/2/∂	10.5 ± 0.02 ns	10.5 ± 0.02 ns	1
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 2/2/∂²	4.96 ± 0.001 ns	5.27 ± 0.01 ns	0.941
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 2/2/∇	13 ± 0.1 ns	12.3 ± 0.05 ns	1.05
RBF/Polyharmonic spline (r⁷)
└─Polynomial augmentation: degree 2/2/∇²	5.47 ± 0.019 ns	8.18 ± 0.041 ns	0.668
time_to_load	0.525 ± 0.0038 s	0.521 ± 0.0029 s	1.01

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[codecov]

Codecov Report

Attention: Patch coverage is 99.74874% with 1 line in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/solve_with_Hermite.jl	99.08%	1 Missing ⚠️

Files with missing lines	Coverage Δ
src/RadialBasisFunctions.jl	100.00% <ø> (ø)
src/StencilData.jl	100.00% <100.00%> (ø)
src/basis/polyharmonic_spline.jl	100.00% <100.00%> (ø)
src/operators/monomial/partial.jl	98.25% <100.00%> (+0.17%)	⬆️
src/solve.jl	98.76% <100.00%> (-0.02%)	⬇️
src/solve_utils.jl	100.00% <100.00%> (ø)
src/solve_with_Hermite.jl	99.08% <99.08%> (ø)

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Davide-Miotti]

improvements to be made inside StencilData.jl in the function _update_stencil() there we should use an LDL solver which does not allocate anything (I can implement it fairly quickly) and use it instead of the backslash solver "", in order to do that we might include a matrix of weights inside StencilData, which is preallocated at the moment the struct is instantiated

[kylebeggs]

Thanks Davide for this PR. This is going to be awesome. So right now if you want to use Hermite to build an operator, how would that look like? I think you'd have to use all the low-level functions that aren't meant for the public users (i.e. they are not exported from the package), right? if so, that's fine, we can merge this after we make the minor changes I commented on. Then we can expose this functionality to the public API in another PR. This is a great start to get the ball rolling.

[kylebeggs]

there we should use an LDL solver which does not allocate anything (I can implement it fairly quickly) and use it instead of the backslash solver ""

the matrices are Symmetric types, so \ does use LDL. see https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize which is what \ calls.

Some future work could be implementing custom solvers to exploit the structure (and save some memory because we don't care about saving/returning the lagrange multipliers). As you know this is not a new idea (matrices coming from KKT conditions) but that usually involves bigger matrices. I think if we use MMatrix (https://juliaarrays.github.io/StaticArrays.jl/stable/api/#Mutable-arrays:-MVector,-MMatrix-and-MArray). Its something I have thought about some, but haven't made much progress.

[Davide-Miotti]

the matrices are Symmetric types, so \ does use LDL. see https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize which is what \ calls.

thank you, I pushed a version where the weights are a field of StencilData

[Davide-Miotti]

Some future work could be implementing custom solvers to exploit the structure (and save some memory because we don't care about saving/returning the lagrange multipliers). As you know this is not a new idea (matrices coming from KKT conditions) but that usually involves bigger matrices. I think if we use MMatrix (https://juliaarrays.github.io/StaticArrays.jl/stable/api/#Mutable-arrays:-MVector,-MMatrix-and-MArray). Its something I have thought about some, but haven't made much progress.

Thank you, I'll read this, we should also think about using an AMG solver for the sparse system.

[Davide-Miotti]

I think you'd have to use all the low-level functions that aren't meant for the public users (i.e. they are not exported from the package), right?

Correct, I did not think about implementing a higher-level function yet, I wanted to have this "backend" part finished and tested before we do that. We might face some issue, like the fact that we need access to normals and boundary flags, in my view these should be aggregated along with point coordinates (and possibly other info) in a larger struct.

[kylebeggs]

Ok, then we can merge this since it is non-breaking and we can make further modifications and update the user facing code at a later time. Thanks Davide!! 🎉