Mercurial > repos > public > sbplib_julia
changeset 1679:529b533a1dbb feature/sbp_operators/laplace_curvilinear
Merge feature/sbp_operators/laplace_curvilinear
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Sun, 30 Jun 2024 10:57:05 +0200 |
parents | de6300bd36cc (current diff) 13a7a4ff49e3 (diff) |
children | b30db2ea34ed |
files | Project.toml src/Grids/Grids.jl |
diffstat | 13 files changed, 85 insertions(+), 37 deletions(-) [+] |
line wrap: on
line diff
--- a/Project.toml Fri Jun 28 21:40:19 2024 +0200 +++ b/Project.toml Sun Jun 30 10:57:05 2024 +0200 @@ -10,12 +10,12 @@ TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac" [weakdeps] +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" [extensions] +SbplibPlotsExt = "Plots" SbplibMakieExt = "Makie" -SbplibPlotsExt = "Plots" [compat] julia = "1.5"
--- a/docs/src/grids_and_grid_functions.md Fri Jun 28 21:40:19 2024 +0200 +++ b/docs/src/grids_and_grid_functions.md Sun Jun 30 10:57:05 2024 +0200 @@ -7,6 +7,24 @@ ## Interface for grids All grids are expected to work as a grid function for the coordinate function, and thus implements Julia's Indexing- and Iteration-interfaces. Notably they are *not* abstract arrays because that inteface is too restrictive for the types of grids we wish to implement. + +## Plotting +Plotting of grids and grid functions is supported through a package extension with Makie.jl. + +For grids we have: +* `plot(::Grid{<:Any,2})` (same as `lines`) +* `lines(::Grid{<:Any,2})` +* `scatter(::Grid{<:Any,2})` + +For 1D grid functions we have: +* `plot(::Grid{<:Any,1}, ::AbstractVector)` (same as `lines`) +* `lines(::Grid{<:Any,1}, ::AbstractVector)` +* `scatter(::Grid{<:Any,1}, ::AbstractVector)` + +For 2D grid functions we have: +* `plot(::Grid{<:Any,2}, ::AbstractArray{<:Any,2})` (constructs a 2d mesh) +* `surface(::Grid{<:Any,2}, ::AbstractArray{<:Any,2})` + ## To write about <!-- # TODO: --> * Grid functions
--- a/ext/SbplibMakieExt.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/ext/SbplibMakieExt.jl Sun Jun 30 10:57:05 2024 +0200 @@ -8,9 +8,6 @@ function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2}) ps = map(Tuple, g)[:] values = gf[:] - - N = length(ps) - faces = Vector{NTuple{3,Int}}() n,m = size(g)
--- a/src/Grids/Grids.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/src/Grids/Grids.jl Sun Jun 30 10:57:05 2024 +0200 @@ -30,6 +30,30 @@ export ConcreteChart export parameterspace +export ParameterSpace +export HyperBox +export Simplex +export Interval +export Rectangle +export Box +export Triangle +export Tetrahedron + +export limits +export unitinterval +export unitsquare +export unitcube +export unithyperbox + +export verticies +export unittriangle +export unittetrahedron +export unitsimplex + +export Chart +export ConcreteChart +export parameterspace + # Grid export Grid export coordinate_size @@ -39,6 +63,7 @@ export boundary_indices export boundary_identifiers export boundary_grid +export min_spacing export coarsen export refine export eval_on
--- a/src/Grids/equidistant_grid.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/src/Grids/equidistant_grid.jl Sun Jun 30 10:57:05 2024 +0200 @@ -47,6 +47,8 @@ """ inverse_spacing(g::EquidistantGrid) = 1/step(g.points) +min_spacing(g::EquidistantGrid) = spacing(g) + boundary_identifiers(::EquidistantGrid) = (Lower(), Upper()) boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
--- a/src/Grids/grid.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/src/Grids/grid.jl Sun Jun 30 10:57:05 2024 +0200 @@ -82,6 +82,14 @@ # TODO: Implement `setindex!`? # TODO: Implement a more general ComponentView that can handle non-AbstractArrays. + +""" + min_spacing(g::Grid) + +The smallest distance between any pair of grid points in `g`. +""" +function min_spacing end + """ refine(g::Grid, r)
--- a/src/Grids/tensor_grid.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/src/Grids/tensor_grid.jl Sun Jun 30 10:57:05 2024 +0200 @@ -50,6 +50,12 @@ Base.size(g::TensorGrid, d) = size(g)[d] +function min_spacing(g::TensorGrid) + relevant_grids = filter(g->!isa(g,ZeroDimGrid),g.grids) + d = min_spacing.(relevant_grids) + return minimum(d) +end + refine(g::TensorGrid, r::Int) = mapreduce(g->refine(g,r), TensorGrid, g.grids) coarsen(g::TensorGrid, r::Int) = mapreduce(g->coarsen(g,r), TensorGrid, g.grids)
--- a/src/RegionIndices/RegionIndices.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/src/RegionIndices/RegionIndices.jl Sun Jun 30 10:57:05 2024 +0200 @@ -13,7 +13,7 @@ Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i) Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i) Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i) - Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed? + Index(t::Tuple{T, DataType}) where T<:Integer = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed? end export Index
--- a/test/Grids/equidistant_grid_test.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/test/Grids/equidistant_grid_test.jl Sun Jun 30 10:57:05 2024 +0200 @@ -57,6 +57,11 @@ @test inverse_spacing(EquidistantGrid(0:0.1:10)) == 10 end + @testset "min_spacing" begin + @test min_spacing(EquidistantGrid(0:10)) == 1 + @test min_spacing(EquidistantGrid(0:0.1:10)) == 0.1 + end + @testset "boundary_identifiers" begin g = EquidistantGrid(0:0.1:10) @test boundary_identifiers(g) == (Lower(), Upper())
--- a/test/Grids/tensor_grid_test.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/test/Grids/tensor_grid_test.jl Sun Jun 30 10:57:05 2024 +0200 @@ -138,6 +138,15 @@ @test axes(g) == (1:11,1:6) end + @testset "min_spacing" begin + g₁ = EquidistantGrid(range(0,1,length=11)) + g₂ = EquidistantGrid(range(2,3,length=6)) + g₃ = ZeroDimGrid(@SVector[1,2]) + + @test min_spacing(TensorGrid(g₁, g₂)) == 1/10 + @test min_spacing(TensorGrid(g₂, g₃)) == 1/5 + end + @testset "refine" begin g1(n) = EquidistantGrid(range(0,1,length=n)) g2(n) = EquidistantGrid(range(2,3,length=n))
--- a/test/LazyTensors/lazy_tensor_operations_test.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/test/LazyTensors/lazy_tensor_operations_test.jl Sun Jun 30 10:57:05 2024 +0200 @@ -4,13 +4,13 @@ using Tullio -struct DummyMapping{T,R,D} <: LazyTensor{T,R,D} end +struct TransposableDummyMapping{T,R,D} <: LazyTensor{T,R,D} end -LazyTensors.apply(m::DummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply -LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose +LazyTensors.apply(m::TransposableDummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply +LazyTensors.apply_transpose(m::TransposableDummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose -LazyTensors.range_size(m::DummyMapping) = :range_size -LazyTensors.domain_size(m::DummyMapping) = :domain_size +LazyTensors.range_size(m::TransposableDummyMapping) = :range_size +LazyTensors.domain_size(m::TransposableDummyMapping) = :domain_size struct SizeDoublingMapping{T,R,D} <: LazyTensor{T,R,D} @@ -24,7 +24,7 @@ @testset "Mapping transpose" begin - m = DummyMapping{Float64,2,3}() + m = TransposableDummyMapping{Float64,2,3}() @test m' isa LazyTensor{Float64, 3,2} @test m'' == m @test apply(m',zeros(Float64,(0,0)), 0, 0, 0) == :apply_transpose
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl Sun Jun 30 10:57:05 2024 +0200 @@ -13,20 +13,9 @@ using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils -""" - monomial(x,k) - -Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``. -Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)`` -""" -function monomial(x,k) - if k < 0 - return zero(x) - end - x^k/factorial(k) -end @testset "undivided_skewed04" begin + monomial(x,k) = k < 0 ? zero(x) : x^k/factorial(k) g = equidistant_grid(0., 11., 20) D,Dᵀ = undivided_skewed04(g, 1)
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Fri Jun 28 21:40:19 2024 +0200 +++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Sun Jun 30 10:57:05 2024 +0200 @@ -7,19 +7,6 @@ using Sbplib.SbpOperators: closure_size, Stencil, VolumeOperator -""" - monomial(x,k) - -Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``. -Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)`` -""" -function monomial(x,k) - if k < 0 - return zero(x) - end - x^k/factorial(k) -end - @testset "first_derivative" begin @testset "Constructors" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) @@ -39,6 +26,8 @@ @testset "Accuracy conditions" begin N = 20 g = equidistant_grid(0//1, 2//1, N) + + monomial(x,k) = k < 0 ? zero(x) : x^k/factorial(k) @testset for order ∈ [2,4] stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) D₁ = first_derivative(g, stencil_set)