Open DanielVandH opened 2 months ago
Perhaps as a simple first step it would be useful to just run these scripts as parts of the CI to just see if they run completely, e.g. in runtests.jl
dir = readdir("./examples")
filter!(file -> endswith(file, ".jl"), dir)
for example_path in dir
script = joinpath("./examples", example_path)
mod = @eval module $(gensym()) end
Base.include(mod, script) # make sure each script is self-contained
end
unless you have something else you usually do instead for this.
That’s a great idea. Can you add a PR?
Note that’s using the old ApproxFun version. We can probably change it to Laplacian(axes(S,1))
we should also add laplacian(S) which would match the diff(S) in 1D
PR made. Yeah, I was trying to convert it to the newer version, and got to here before making this issue:
using MultivariateOrthogonalPolynomials, SpecialFunctions
a, b, c = 1, 1, 1
S = WeightedTriangle(a, b, c)
xy = axes(S, 1)
x, y = first.(xy), last.(xy)
∆ = Laplacian(xy)
f = @. 1.0 + erf(5.0(1.0 - 10.0((x - 1/2)^2 + (y-1/2)^2)))
N = 500
M = sparse(∆[Block.(1:N), Block.(1:N)])
StackOverflowError:
Stacktrace:
[1] unblock(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ BlockArrays C:\Users\User\.julia\packages\BlockArrays\K0HfQ\src\views.jl:10
[2] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:81
[3] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:85
[4] _getindex(::Type{Tuple{StaticArraysCore.SVector{2, Float64}}}, A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\abstractquasiarray.jl:372
[5] getindex(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}})
@ QuasiArrays C:\Users\User\.julia\packages\QuasiArrays\ZCn0F\src\abstractquasiarray.jl:367
--- the last 5 lines are repeated 15995 more times ---
[79981] unblock(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ BlockArrays C:\Users\User\.julia\packages\BlockArrays\K0HfQ\src\views.jl:10
[79982] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, inds::Tuple{QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:81
[79983] to_indices(A::QuasiArrays.Inclusion{StaticArraysCore.SVector{2, Float64}, DomainSets.EuclideanUnitSimplex{2, Float64, :closed}}, I::Tuple{BlockArrays.BlockRange{1, Tuple{UnitRange{Int64}}}})
@ ContinuumArrays C:\Users\User\.julia\packages\ContinuumArrays\Uyiyf\src\ContinuumArrays.jl:85
What's the correct code for reproducing Example 1 in the paper, i.e. $u{xx} + u{yy} = f(x, y)$ for $(x, y) \in T$ with $\left.u\right|_{\partial T} = 0$, where $f(x, y) = 1 + \mathrm{erf}(5(1 - 10((x - 1/2)^2 + (y-1/2)^2)))$?
Hmm I just checked and I realised that the code for partial derivatives for WeightedTriangle
isn't implemented yet. So actually the first task is to get that working. That is, we want to first add a function like:
@simplify function *(Dy::PartialDerivative{2}, P::WeightedTriangle)
a,b,c = P.P.a,P.P.b,P.P.c
k = mortar(Base.OneTo.(oneto(∞)))
T = promote_type(eltype(Dy), eltype(P))
M = _BandedBlockBandedMatrix((k .+ convert(T, b+c))', axes(k,1), (-1,1), (-1,1)) # TODO: replace with correct formula
WeightedTriangle(a,b-1,c-1) * M # TODO: deal with the case where b or c are 0
end
The examples here fail (below is Lines 1-12 from that file)
I believe in this case the fix is to do
S = WeightedTriangle(1, 1, 1)
? This error appears later in the script as well (might be some others too -Laplacian(S)
also errors -, but I haven't run it that far). I could try and fix it up later once I've looked into it all a bit more.