Mercurial > repos > public > sbplib_julia
changeset 1328:f00a205ae347 refactor/grids
Merge default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 02 May 2023 20:14:39 +0200 |
| parents | a956fc8e79a6 (diff) 95cac1ee8476 (current diff) |
| children | e94ddef5e72f |
| files | |
| diffstat | 43 files changed, 1100 insertions(+), 804 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Mon May 01 11:37:09 2023 +0200 +++ b/Manifest.toml Tue May 02 20:14:39 2023 +0200 @@ -2,13 +2,13 @@ julia_version = "1.8.5" manifest_format = "2.0" -project_hash = "b024d6898b484706c36ee3b2a041918f3a9d2088" +project_hash = "a36735c53cfa4453f39635046eeaa47a4ea1231b" [[deps.Adapt]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c" +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.5.0" +version = "3.6.1" [[deps.ArrayInterface]] deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"] @@ -27,9 +27,9 @@ [[deps.Compat]] deps = ["Dates", "LinearAlgebra", "UUIDs"] -git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04" +git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.6.0" +version = "4.6.1" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -77,6 +77,12 @@ deps = ["SHA", "Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" @@ -96,9 +102,24 @@ [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c" +git-tree-sha1 = "08be5ee09a7632c32695d954a602df96a877bf0d" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.3" +version = "0.8.6" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "7756ce473bd10b67245bdebdc8d8670a85f6230b" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.18" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
--- a/Notes.md Mon May 01 11:37:09 2023 +0200 +++ b/Notes.md Tue May 02 20:14:39 2023 +0200 @@ -148,6 +148,7 @@ ## Identifiers for regions The identifiers (`Upper`, `Lower`, `Interior`) used for region indecies should probabily be included in the grid module. This allows new grid types to come with their own regions. +We implement this by refactoring RegionIndices to be agnostic to the region types and then moving the actual types to Grids. ## Regions and tensormappings - [ ] Use a trait to indicate if a LazyTensor uses indices with regions. @@ -199,79 +200,11 @@ Preferably dimensions and sizes should be checked when lazy objects are created, for example TensorApplication, TensorComposition and so on. If dimension checks decreases performance we can make them skippable later. ## Vector valued grid functions -Från slack konversation: - -Jonatan Werpers: -Med vektorvärda gridfunktioner vill vi ju fortfarande att grid funktionen ska vara till exempel AbstractArray{LitenVektor,2} -Och att man ska kunna göra allt man vill med LitenVektor -typ addera, jämföra osv -Och då borde points returnera AbstractArray{LitenVektor{Float,2},2} för ett 2d nät -Men det kanske bara ska vara Static arrays? - -Vidar Stiernström: -Ja, jag vet inte riktigt vad som är en rimlig representation -Du menar en vektor av static arrays då? - -Jonatan Werpers: -Ja, att LitenVektor är en StaticArray - -Vidar Stiernström: -Tuplar känns typ rätt inuitivt för att representera värdet i en punkt -men -det suger att man inte har + och - för dem - -Jonatan Werpers: -Ja precis - -Vidar Stiernström: -så kanske är bra med static arrays i detta fall - -Jonatan Werpers: -Man vill ju kunna köra en Operator rakt på och vara klar eller? - -Vidar Stiernström: -Har inte alls tänkt på hur det vi gör funkar mot vektorvärda funktioner -men känns som staticarrays är hur man vill göra det -tuplar är ju immutable också -blir jobbigt om man bara agerar på en komponent då - -Jonatan Werpers: -Hm… -Tål att tänkas på -Men det lär ju bli mer indirektion med mutables eller? -Hur fungerar det? -Det finns ju hur som helst både SVector och MVector i StaticArrays - -Vidar Stiernström: -När vi jobbat i c/c++ och kollat runt lite hur man brukar göra så lagrar man i princip alla sina obekanta i en lång vektor och så får man specificera i funktioerna vilken komponent man agerar på och till vilken man skriver -så man lagrar grejer enl: w = [u1, v1, u2, v2, …] i 1D. -Men alltså har ingen aning hur julia hanterar detta - -Jonatan Werpers: -Det vi är ute efter kanske är att en grid funcktion är en AbstractArray{T,2} where T<:NågotSomViKanRäknaMed -Och så får den typen var lite vad som helst. - -Vidar Stiernström: -Tror det kan vara farligt att ha nåt som är AbstractArray{LitenArray{NDof},Dim} -Jag gissar att det kompilatorn vill ha är en stor array med doubles - -Jonatan Werpers: -Och sen är det upp till den som använder grejerna att vara smart -Vill man vara trixig kan man väl då imlementera SuperHaxxorGridFunction <: AbstractArray{Array{…},2} som lagrar allt linjärt eller något sånt -Det kommer väl lösa sig när man börjar implementera vektorvärda saker -Euler nästa! -New -Vidar Stiernström: -Det vore skönt att inte behöva skriva såhär varje gång man testar mot en tupel :smile: @test [gp[i]...] ≈ [p[i]...] atol=5e-13 - -Jonatan Werpers: -https://github.com/JuliaArrays/ArraysOfArrays.jl -https://github.com/jw3126/Setfield.jl ### Test-applikationer -div och grad operationer +div- och grad-operationer -Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensor fält och `div()` på sista. +Enligt Wikipedia verkar det som att `∇⋅` agerar på första dimensionen av ett tensorfält och `div()` på sista. Om man generaliserar kanske `∇` i så fall bara lägger till en dimension i början. Kan vi implementera `⋅`(\cdot) så att de fungerar som man vill för både tensor-fält och tensor-operatorer? @@ -279,8 +212,8 @@ Är `∇` ett tensor-fält av tensor-operatorer? Vad är ett tensor-fält i vår kod? Är det en special-fall av en tensor-mapping? ### Grid-funktionen -Grid-funktionon har typen `AbstractArray{T,2} where T`. -`T` kan vara lite vad som helst, tillexemel en SVector eller Array, eller tuple. TensorOperatorerna bryr sig inte om exakt vad det är, mer än att typen måste stödja de operationer som operatorn använder. +Grid-funktioner har typen `AbstractArray{T,2} where T`. +`T` kan vara lite vad som helst, tillexemel en SVector eller Array, eller Tuple. Tensoroperatorerna bryr sig inte om exakt vad det är, mer än att typen måste stödja de operationer som operatorn använder. En nackdel kan vara hur man ska få ut gridfunktionen för tex andra komponenten. @@ -292,7 +225,7 @@ gf[2,3][2] # x̄[2] för en viss gridpunkt ``` -Note: Behöver bestämma om eval on skickar in `x̄` eller `x̄...` till `f`. Eller om man kan stödja båda. +Note: Behöver bestämma om `eval_on` skickar in `x̄` eller `x̄...` till `f`. Eller om man kan stödja båda. ### Tensor operatorer Vi kan ha tensor-operatorer som agerar på ett skalärt fält och ger ett vektorfält eller tensorfält. @@ -309,34 +242,8 @@ Kanske kan man implementera `⋅(tm::LazyTensor{R,D}, v::AbstractArray{T,D})` där T är en AbstractArray, tm på något sätt har komponenter, lika många som T har element. -### Ratade alternativ: - - -#### 2.AbstractArray{T,2+1} where T (NOPE!) -Blir inte den här. Bryter mot alla tankar om hur grid funktioner ska fungera. Om de tillåts ha en annan dimension än nätet blir allt hemskt. - -Man låter helt enkelt arrayen ha en extra dimension. En fördel är att man har en väldigt "native" typ. En nackdel kan vara att det eventuellt blir rörigt vilken dimension olika operatorer ska agera på. I värsta fall behöver vi "kroneckra in" de tillagda dimensionerna. Vektorfältets index kommer också att bli det första eftersom vi vill att de ska lagras kontinuerligt i minnet pga chachen. (Går kanske att lösa med en custom typ men då krånglar man till det för sig). En fördel skulle vara att man enkelt får ut olika komponenter. - -Syntax: -``` -gf = eval_on_grid(g,f) -gf[:,2,3] # Hela vektorn för en gridpunkt -gf[2,2,3] # Andra komponenten av vektor fältet i en punkt. -gf[2,:,:] # -``` - -### Evaluering av funktioner på nät -Hur ska man skriva funktioner som evalueras på nätet? `f(x,y) = ...` eller `f(x̄) = ...`? Eller båda? Kan eval_on_grid se skillnad eller får användaren specificera? - -``` -f(x,y) = [x^2, y^2] -f(x̄) = [x̄[1]^2, x̄[2]^2] -``` - -Påverkas detta av hur vi förväntar oss kunna skapa lata gridfunktioner? - ### Komponenter som gridfunktioner -En viktig operation för vektor fält är att kunna få ut komponenter som grid-funktioner. Detta behöver antagligen kunna ske lazy. +En viktig operation för vektorfält är att kunna få ut komponenter som grid-funktioner. Detta behöver antagligen kunna ske lazy. Det finns ett par olika lösningar: * Implementera en egen typ av view som tar hand om detta. Eller Accessors.jl? * Använda en LazyTensor @@ -354,18 +261,25 @@ @ourview gf[:,:][2] ``` -## Grids embedded in higher dimensions +### Prestanda-aspekter +[Vidar, Discord, 2023-03-03] +Typiskt sett finns det två sätt att representera vektorvärda gridfunktioner AbstractArray{T,Dim} där T är en vektor över komponenterna. Man skulle alltså i 1D ha +u = [ [u1[x1], u2[x1]] , [u1[x2], u2[x2]], ... [u1[xN], u2[xN]]]. Detta brukar kallas array of structs (AoS). Alternativet är struct of arrays (SoA), där man har alla gridpunkter för en given komponent u = [[u1[x1], u1[x2]],... u1[xN]], [u2[x1], u2[x2], ... u2[xN]]]. + +Personligen tycker jag att AoS känns som den mer naturliga representationen? Det skulle göra det enklarare att parallelisera en vektorvärd gridfunktion över gridpunkterna, och om man opererar på olika komponenter i samma funktion så är det också bra ur en minnesaccess-synpunkt då dessa kommer ligga nära vandra i minnet. Problemet är att AoS sabbar vektorisering på CPU då två gridpunkter i en komponent ligger långt bort från varandra. Efter lite eftersökningar (och efter att snackat lite med Ossian) så verkar det ändå som att AoS är dåligt på GPU, där man vill att trådar typiskt sett utföra samma operation på närliggande minne. + +Vad tänker du kring detta ur ett interface-perspektiv? Jag hittade paketet https://github.com/JuliaArrays/StructArrays.jl som verkar erbjuda AoS-interface men SoA-minneslayout så det kanske kan vara något vi kan använda? Inte native-stödd på samma sätt som SVector, men verkar iaf utvecklas aktivt. -For grids generated by asking for boundary grids for a regular grid, it would -make sense if these grids knew they were embedded in a higher dimension. They -would return coordinates in the full room. This would make sense when -drawing points for example, or when evaluating functions on the boundary. +[Efter telefonsamtal] För optimal prestanda behöver vi antagligen se till att man kan räkna ut varje komponent i en punkt individuellt. Detta så att man har frihet att till exempel låta den innersta loopen hålla komponentindexet konstant för att underlätta intruktionsvektorisering. + -Implementation of this is an issue that requires some thought. Adding an extra -"Embedded" type for each grid would make it easy to understand each type but -contribute to "type bloat". On the other hand adapting existing types to -handle embeddedness would complicate the now very simple grid types. Are there -other ways of doing the implentation? +[Vidare tankar] + * Det borde bara vara output-gridfunktionen som behöver special-indexeras? Det viktiga på inputsidan är att den är lagrad på rätt sätt i minnet. + * Det borde inte vara några problem att behålla det "optimala" interfacet (gf[1,1,1][2]) till gridfunktionerna. Om man verkligen behöver kan skapa parallella indexeringsmetoder som gör det man behöver, i.e, "deep indexing". + * Det är inte säkert att vi behöver göra något speciellt på outputsidan överhuvudtaget. Det känns inte orimligt att kompilatorn skulle kunna optimera bort den koden som räknar ut onödiga komponenter. + * Om vi behöver special-indexering kommer till exempel LazyTensorApplication att behöva implementera det. + * För att komma vidare med något mer avancerat behöver vi antagligen implementera några operatorer som ger och agerar på vektorvärda funktioner. Tex grad, elastiska operatorn, andra? + ## Performance measuring We should be measuring performance early. How does our effective cpu and memory bandwidth utilization compare to peak performance?
--- a/Project.toml Mon May 01 11:37:09 2023 +0200 +++ b/Project.toml Tue May 02 20:14:39 2023 +0200 @@ -4,6 +4,7 @@ version = "0.1.0" [deps] +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
--- a/docs/Manifest.toml Mon May 01 11:37:09 2023 +0200 +++ b/docs/Manifest.toml Tue May 02 20:14:39 2023 +0200 @@ -10,10 +10,10 @@ version = "0.0.1" [[deps.Adapt]] -deps = ["LinearAlgebra"] -git-tree-sha1 = "0310e08cb19f5da31d08341c6120c047598f5b9c" +deps = ["LinearAlgebra", "Requires"] +git-tree-sha1 = "cc37d689f599e8df4f464b2fa3870ff7db7492ef" uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -version = "3.5.0" +version = "3.6.1" [[deps.ArrayInterface]] deps = ["ArrayInterfaceCore", "Compat", "IfElse", "LinearAlgebra", "SnoopPrecompile", "Static"] @@ -35,9 +35,9 @@ [[deps.Compat]] deps = ["Dates", "LinearAlgebra", "UUIDs"] -git-tree-sha1 = "61fdd77467a5c3ad071ef8277ac6bd6af7dd4c04" +git-tree-sha1 = "7a60c856b9fa189eb34f5f8a6f6b5529b7942957" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "4.6.0" +version = "4.6.1" [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] @@ -119,9 +119,9 @@ [[deps.Parsers]] deps = ["Dates", "SnoopPrecompile"] -git-tree-sha1 = "151d91d63d8d6c1a5789ecb7de51547e00480f1b" +git-tree-sha1 = "478ac6c952fddd4399e71d4779797c538d0ff2bf" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "2.5.4" +version = "2.5.8" [[deps.Preferences]] deps = ["TOML"] @@ -141,12 +141,18 @@ deps = ["SHA", "Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +[[deps.Requires]] +deps = ["UUIDs"] +git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7" +uuid = "ae029012-a4dd-5104-9daa-d747884805df" +version = "1.3.0" + [[deps.SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" version = "0.7.0" [[deps.Sbplib]] -deps = ["TOML", "TiledIteration"] +deps = ["StaticArrays", "TOML", "TiledIteration"] path = ".." uuid = "5a373a26-915f-4769-bcab-bf03835de17b" version = "0.1.0" @@ -169,9 +175,24 @@ [[deps.Static]] deps = ["IfElse"] -git-tree-sha1 = "c35b107b61e7f34fa3f124026f2a9be97dea9e1c" +git-tree-sha1 = "08be5ee09a7632c32695d954a602df96a877bf0d" uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" -version = "0.8.3" +version = "0.8.6" + +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "7756ce473bd10b67245bdebdc8d8670a85f6230b" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.18" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + +[[deps.Statistics]] +deps = ["LinearAlgebra", "SparseArrays"] +uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/grid_refactor.md Tue May 02 20:14:39 2023 +0200 @@ -0,0 +1,125 @@ +# Grids refactor + +## Goals before merging +* A somewhat clear path towards multi-block grids and their grid functions. +* A somewhat clear path towards implementations of div() and rot() or the elastic operator (See Notes.md) + +## Change summary +* `EquidistantGrid` is now only a 1D thing. +* Higher dimensions are supported through `TensorGrid`. +* The old behavior of `EquidistantGrid` has been moved to the function `equidistant_grid`. +* Grids embedded in higher dimensions are now supported through tensor products with `ZeroDimGrid`s. +* Vector valued grid functions are now supported and the default element type is `SVector`. +* Grids are now expected to support Julia's indexing and iteration interface. +* `eval_on` can be called with both `f(x,y,...)` and `f(x̄)`. + + +## TODO +* Add benchmarks or allocation tests for eval_on and indexing grids. +* Document the expected behavior of grid functions +* Write down the thinking around Grid being an AbstractArray. Why it doesn't work. +* Write about the design choices in the docs. +* Merge and run benchmarks + +* Check all the docstring of all types that have been changed +* Clean out Notes.md of any solved issues +* Delete this document, move remaining notes to Notes.md + +## Remaining work for feature branches +* Multi-block grids +* Periodic grids +* Grids with modified boundary nodes +* Unstructured grids? + +## Frågor + +### Should we move utility functions to their own file? + +### Ska `Grid` vara en AbstractArray? +Efter som alla nät ska agera som en gridfunktion för koordinaterna måste man +svara på frågan hur vi hanterar generellla gridfunktioner samtidigt. + +Några saker att förhålla sig till: + - Multiblock nät? + - Unstructured? + - Triangular structured grids? + - Non-simply connected? + - CG/DG-nät + +Om alla nät är någon slags AbstractArray så kan tillexempel ett multiblock nät vara en AbstractArray{<:Grid, 1} och motsvarande gridfunktioner AbstractArray{<:AbstractArray,1}. +Där man alltså först indexerar för vilket när man vill åt och sedan indexerar nätet. Tex `mg[2][32,12]`. + +Ett ostrukturerat nät med till exempel trianglar skulle vi kunna se på samma sätt som ett multiblocknät. Antagligen har de två typerna av nät olika underliggande datastruktur anpassade efter ändamålet. + +Hur fungerar tankarna ovan om man har tex tensorprodukten av ett ostrukturerat nät och en ekvidistant nät? +```julia +m = Mesh2DTriangle(10) +e = EqudistantGrid(range(1:10) + +e[4] # fourth point + +m[3][5] # Fifth node in third triangle +m[3,5] # Fifth node in third triangle # Funkar bara om alla nät är samma, (stämmer inte i mb-fallet) + +g = TensorGrid(m, e) + +g[3,4][5] # ?? +g[3,4] # ?? + +g[3,5,4] # ?? + + + +``` + +Alla grids kanske inte är AbstractArray? Måste de vara rektangulära? Det blir svårt för strukturerade trianglar och NSC-griddar. Men de ska i allafall vara indexerbara? + +Man skulle kunna utesluta MultiblockGrid i tensorgrids + +CG-nät och DG-nät blir olika. +På CG-nät kanske man både vill indexera noder och trianglar beroende på vad man håller på med? + + +Om griddarna inte ska vara AbstractArray finns det många andra ställen som blir konstiga om de är AbstractArray. TensorApplication?! LazyArrays?! Är alla saker vi jobbar med egentligen mer generella object? Finns det något sätt att uttrycka koden så att man kan välja? + + +Det vi är ute efter är kanske att griddarna uppfyller Iteration och Indexing interfacen. + +#### Försök till slutsater + * Multiblock-nät indexeras i två nivåer tex `g[3][3,4]` + * Vi struntar i att implementera multiblock-nät som en del av ett tensorgrid till att börja med. + * En grid kan inte alltid vara en AbstractArray eftersom till exempel ett NCS eller strukturerad triangel inte har rätt form. + * Om vi har nod-indexerade ostrukturerade nät borde de fungera med TensorGrid. + * Griddar ska uppfylla Indexing och Iteration interfacen + +### Should Grid have function for the target manifold dimension? +Where would it be used? + In the constructor for TensorGrid + In eval on if we want to allow multiargument functions + Elsewhere? + +An alternative is to analyze T in Grid{T,D} to find the answer. (See combined_coordinate_vector_type in tensor_grid.jl) + +### Lazy version of map for our needs? +Could be used to + * evaulate functions on grids + * pick out components of grid functions + * More? + +Maybe this: +```julia +struct LazyMappedArray <: LazyArray + f::F + v::AT +end +``` + +Could allow us to remove eval_on. + +### Do we need functions like `getcomponent`? +Perhaps this can be more cleanly solved using map or a lazy version of map? +That approach would be more flexible and more general requiring few specialized functions. + +(see "Lazy version of map for our needs?" above) + +## Implement the tensor product operator for grids?
--- a/src/Grids/Grids.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/Grids/Grids.jl Tue May 02 20:14:39 2023 +0200 @@ -1,31 +1,46 @@ module Grids using Sbplib.RegionIndices +using Sbplib.LazyTensors +using StaticArrays # Grid export Grid -export dims -export points -export evalOn +export target_manifold_dim +export component_type + +export TensorGrid +export ZeroDimGrid + +export TensorGridBoundary + +export grid_id +export boundary_id + +export eval_on +export getcomponent # BoundaryIdentifier export BoundaryIdentifier -export CartesianBoundary -export dim -export region + # EquidistantGrid export EquidistantGrid export spacing export inverse_spacing -export restrict export boundary_identifiers export boundary_grid export refine export coarsen +export equidistant_grid +export CartesianBoundary + +abstract type BoundaryIdentifier end include("grid.jl") -include("boundary_identifier.jl") +include("tensor_grid.jl") include("equidistant_grid.jl") +include("zero_dim_grid.jl") + end # module
--- a/src/Grids/boundary_identifier.jl Mon May 01 11:37:09 2023 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ - -abstract type BoundaryIdentifier end - -struct CartesianBoundary{Dim, R<:Region} <: BoundaryIdentifier end -dim(::CartesianBoundary{Dim, R}) where {Dim, R} = Dim -region(::CartesianBoundary{Dim, R}) where {Dim, R} = R() \ No newline at end of file
--- a/src/Grids/equidistant_grid.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/Grids/equidistant_grid.jl Tue May 02 20:14:39 2023 +0200 @@ -1,71 +1,25 @@ - """ - EquidistantGrid{Dim,T<:Real} <: Grid - -`Dim`-dimensional equidistant grid with coordinates of type `T`. +# TODO: """ -struct EquidistantGrid{Dim,T<:Real} <: Grid - size::NTuple{Dim, Int} - limit_lower::NTuple{Dim, T} - limit_upper::NTuple{Dim, T} - - function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T} - if any(size .<= 0) - throw(DomainError("all components of size must be postive")) - end - if any(limit_upper.-limit_lower .<= 0) - throw(DomainError("all side lengths must be postive")) - end - return new{Dim,T}(size, limit_lower, limit_upper) - end +#TODO: Document recomendations for type of range. (LinRange is faster?) +struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1} + points::R end - -""" - EquidistantGrid(size, limit_lower, limit_upper) - -Construct an equidistant grid with corners at the coordinates `limit_lower` and -`limit_upper`. - -The length of the domain sides are given by the components of -`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and `limit_upper=(1,2)` the domain is defined -as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative. +Base.eachindex(g::EquidistantGrid) = eachindex(g.points) -The number of equidistantly spaced points in each coordinate direction are given -by the tuple `size`. -""" -function EquidistantGrid(size, limit_lower, limit_upper) - return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper) -end - - -""" - EquidistantGrid{T}() +# Indexing interface +Base.getindex(g::EquidistantGrid, i) = g.points[i] +Base.firstindex(g::EquidistantGrid) = firstindex(g.points) +Base.lastindex(g::EquidistantGrid) = lastindex(g.points) -Constructs a 0-dimensional grid. -""" -EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid - - -""" - EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) +# Iteration interface +Base.iterate(g::EquidistantGrid) = iterate(g.points) +Base.iterate(g::EquidistantGrid, state) = iterate(g.points, state) -Convenience constructor for 1D grids. -""" -function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T - return EquidistantGrid((size,),(limit_lower,),(limit_upper,)) -end - -Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T - -Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size) - -Base.size(g::EquidistantGrid) = g.size - -Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim - - - +Base.IteratorSize(::Type{<:EquidistantGrid}) = Base.HasShape{1}() +Base.length(g::EquidistantGrid) = length(g.points) +Base.size(g::EquidistantGrid) = size(g.points) """ @@ -73,7 +27,7 @@ The spacing between grid points. """ -spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1) +spacing(g::EquidistantGrid) = step(g.points) """ @@ -81,94 +35,27 @@ The reciprocal of the spacing between grid points. """ -inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid) +inverse_spacing(g::EquidistantGrid) = 1/step(g.points) -""" - points(grid::EquidistantGrid) - -The point of the grid as an array of tuples with the same dimension as the grid. -The points are stored as [(x1,y1), (x1,y2), … (x1,yn); - (x2,y1), (x2,y2), … (x2,yn); - ⋮ ⋮ ⋮ - (xm,y1), (xm,y2), … (xm,yn)] -""" -function points(grid::EquidistantGrid) - indices = Tuple.(CartesianIndices(grid.size)) - h = spacing(grid) - return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices) -end - - -""" - restrict(::EquidistantGrid, dim) - -Pick out given dimensions from the grid and return a grid for them. -""" -function restrict(grid::EquidistantGrid, dim) - size = grid.size[dim] - limit_lower = grid.limit_lower[dim] - limit_upper = grid.limit_upper[dim] - - return EquidistantGrid(size, limit_lower, limit_upper) -end +boundary_identifiers(::EquidistantGrid) = (Lower(), Upper()) +boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin]) +boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end]) """ - orthogonal_dims(grid::EquidistantGrid,dim) - -Returns the dimensions of grid orthogonal to that of dim. -""" -function orthogonal_dims(grid::EquidistantGrid, dim) - orth_dims = filter(i -> i != dim, dims(grid)) - if orth_dims == dims(grid) - throw(DomainError(string("dimension ",string(dim)," not matching grid"))) - end - return orth_dims -end - - -""" - boundary_identifiers(::EquidistantGrid) - -Returns a tuple containing the boundary identifiers for the grid, stored as - (CartesianBoundary(1,Lower), - CartesianBoundary(1,Upper), - CartesianBoundary(2,Lower), - ...) -""" -boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,) - - -""" - boundary_grid(grid::EquidistantGrid, id::CartesianBoundary) - -Creates the lower-dimensional restriciton of `grid` spanned by the dimensions -orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional -grid is a zero-dimensional grid. -""" -function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary) - orth_dims = orthogonal_dims(grid, dim(id)) - return restrict(grid, orth_dims) -end -boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}() - - -""" - refine(grid::EquidistantGrid, r::Int) + refine(g::EquidistantGrid, r::Int) Refines `grid` by a factor `r`. The factor is applied to the number of intervals which is 1 less than the size of the grid. See also: [`coarsen`](@ref) """ -function refine(grid::EquidistantGrid, r::Int) - sz = size(grid) - new_sz = (sz .- 1).*r .+ 1 - return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper) +function refine(g::EquidistantGrid, r::Int) + new_sz = (length(g) - 1)*r + 1 + return EquidistantGrid(change_length(g.points, new_sz)) end - """ coarsen(grid::EquidistantGrid, r::Int) @@ -178,14 +65,68 @@ See also: [`refine`](@ref) """ -function coarsen(grid::EquidistantGrid, r::Int) - sz = size(grid) - - if !all(n -> (n % r == 0), sz.-1) +function coarsen(g::EquidistantGrid, r::Int) + if (length(g)-1)%r != 0 throw(DomainError(r, "Size minus 1 must be divisible by the ratio.")) end - new_sz = (sz .- 1).÷r .+ 1 + new_sz = (length(g) - 1)÷r + 1 + + return EquidistantGrid(change_length(g.points, new_sz)) +end + + +""" + equidistant_grid(size::Dims, limit_lower, limit_upper) + +Construct an equidistant grid with corners at the coordinates `limit_lower` and +`limit_upper`. + +The length of the domain sides are given by the components of +`limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and +`limit_upper=(1,2)` the domain is defined as `(-1,1)x(0,2)`. The side lengths +of the grid are not allowed to be negative. + +The number of equidistantly spaced points in each coordinate direction are given +by the tuple `size`. + +Note: If `limit_lower` and `limit_upper` are integers and `size` would allow a +completely integer grid, `equidistant_grid` will still return a floating point +grid. This simlifies the implementation and avoids certain surprise +behaviours. +""" +function equidistant_grid(size::Dims, limit_lower, limit_upper) + gs = map(equidistant_grid, size, limit_lower, limit_upper) + return TensorGrid(gs...) +end - return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper) +""" + equidistant_grid(size::Int, limit_lower::T, limit_upper::T) + +Constructs a 1D equidistant grid. +""" +function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T + if any(size .<= 0) + throw(DomainError("size must be postive")) + end + + if any(limit_upper.-limit_lower .<= 0) + throw(DomainError("side length must be postive")) + end + return EquidistantGrid(range(limit_lower, limit_upper, length=size)) # TBD: Should it use LinRange instead? end + +CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary? + + +""" + change_length(::AbstractRange, n) + +Change the length of a range to `n`, keeping the same start and stop. +""" +function change_length end + +change_length(r::UnitRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n)) +change_length(r::StepRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n)) +change_length(r::StepRangeLen, n) = range(r[begin], r[end], n) +change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)
--- a/src/Grids/grid.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/Grids/grid.jl Tue May 02 20:14:39 2023 +0200 @@ -1,27 +1,59 @@ """ - Grid + Grid{T,D} + +The top level type for grids. +TODO: Should implement - Base.ndims(grid::Grid) - points(grid::Grid) + * interfaces for iteration and indexing +""" +abstract type Grid{T,D} end + +Base.ndims(::Grid{T,D}) where {T,D} = D +Base.eltype(::Type{<:Grid{T}}) where T = T +target_manifold_dim(::Grid{T}) where T = _ncomponents(T) # TBD: Name of this function?! +component_type(::Grid{T}) where T = eltype(T) """ -abstract type Grid end -function points end +# TODO +""" +function refine end + +""" +# TODO +""" +function coarsen end """ - dims(grid::Grid) +# TODO +""" +function boundary_identifiers end -A range containing the dimensions of `grid` +""" +# TODO """ -dims(grid::Grid) = 1:ndims(grid) +function boundary_grid end +# TBD Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff? + + +# TODO: Make sure that all grids implement all of the above. + """ - evalOn(grid::Grid, f::Function) +TODO: -Evaluate function `f` on `grid` +* Mention map(f,g) if you want a concrete array """ -function evalOn(grid::Grid, f::Function) - F(x) = f(x...) - return F.(points(grid)) +eval_on(g::Grid, f) = eval_on(g, f, Base.IteratorSize(g)) # TBD: Borde f vara först som i alla map, sum, och dylikt +function eval_on(g::Grid, f, ::Base.HasShape) + if hasmethod(f, (Any,)) + return LazyTensors.LazyFunctionArray((I...)->f(g[I...]), size(g)) + else + return LazyTensors.LazyFunctionArray((I...)->f(g[I...]...), size(g)) + end end + + +# TODO: Explain how and where these are intended to be used +_ncomponents(::Type{<:Number}) = 1 +_ncomponents(T::Type{<:SVector}) = length(T)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/tensor_grid.jl Tue May 02 20:14:39 2023 +0200 @@ -0,0 +1,94 @@ +""" + TensorGrid{T,D} <: Grid{T,D} + +* Only supports HasShape grids at the moment +""" +struct TensorGrid{T,D,GT<:NTuple{N,Grid} where N} <: Grid{T,D} + grids::GT + + function TensorGrid(gs...) + # T = combined_coordinate_vector_type(eltype.(gs)...) + T = mapreduce(eltype, combined_coordinate_vector_type, gs) + D = sum(ndims, gs) + + return new{T,D,typeof(gs)}(gs) + end +end + +# Indexing interface +function Base.getindex(g::TensorGrid, I...) + szs = ndims.(g.grids) + + Is = LazyTensors.split_tuple(I, szs) + ps = map((g,I)->SVector(g[I...]), g.grids, Is) + + return vcat(ps...) +end + +function Base.eachindex(g::TensorGrid) + szs = LazyTensors.concatenate_tuples(size.(g.grids)...) + return CartesianIndices(szs) +end + +# Iteration interface +Base.iterate(g::TensorGrid) = iterate(Iterators.product(g.grids...)) |> _iterate_combine_coords +Base.iterate(g::TensorGrid, state) = iterate(Iterators.product(g.grids...), state) |> _iterate_combine_coords +_iterate_combine_coords(::Nothing) = nothing +_iterate_combine_coords((next,state)) = combine_coordinates(next...), state + +Base.IteratorSize(::Type{<:TensorGrid{<:Any, D}}) where D = Base.HasShape{D}() +Base.eltype(::Type{<:TensorGrid{T}}) where T = T +Base.length(g::TensorGrid) = sum(length, g.grids) +Base.size(g::TensorGrid) = LazyTensors.concatenate_tuples(size.(g.grids)...) + + +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) + +""" +# TODO: +""" +struct TensorGridBoundary{N, BID} <: BoundaryIdentifier end +grid_id(::TensorGridBoundary{N, BID}) where {N, BID} = N +boundary_id(::TensorGridBoundary{N, BID}) where {N, BID} = BID() + +""" + boundary_identifiers(::TensorGrid) + +Returns a tuple containing the boundary identifiers for the grid. +""" +function boundary_identifiers(g::TensorGrid) + per_grid = map(eachindex(g.grids)) do i + return map(bid -> TensorGridBoundary{i, typeof(bid)}(), boundary_identifiers(g.grids[i])) + end + return LazyTensors.concatenate_tuples(per_grid...) +end + + +""" + boundary_grid(grid::TensorGrid, id::TensorGridBoundary) + +The grid for the boundary specified by `id`. +""" +function boundary_grid(g::TensorGrid, bid::TensorGridBoundary) + local_boundary_grid = boundary_grid(g.grids[grid_id(bid)], boundary_id(bid)) + new_grids = Base.setindex(g.grids, local_boundary_grid, grid_id(bid)) + return TensorGrid(new_grids...) +end + + +function combined_coordinate_vector_type(coordinate_types...) + combined_coord_length = mapreduce(_ncomponents, +, coordinate_types) + combined_coord_type = mapreduce(eltype, promote_type, coordinate_types) + + if combined_coord_length == 1 + return combined_coord_type + else + return SVector{combined_coord_length, combined_coord_type} + end +end + +function combine_coordinates(coords...) + return mapreduce(SVector, vcat, coords) + # return SVector(coords...) +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/zero_dim_grid.jl Tue May 02 20:14:39 2023 +0200 @@ -0,0 +1,26 @@ +""" + ZeroDimGrid{T} <: Grid{T,0} +# TODO +""" +struct ZeroDimGrid{T} <: Grid{T,0} + point::T +end + +# Indexing interface +Base.getindex(g::ZeroDimGrid) = g.point +Base.eachindex(g::ZeroDimGrid) = CartesianIndices(()) + +# Iteration interface +Base.iterate(g::ZeroDimGrid) = (g.point, nothing) +Base.iterate(g::ZeroDimGrid, ::Any) = nothing + +Base.IteratorSize(::Type{<:ZeroDimGrid}) = Base.HasShape{0}() +Base.length(g::ZeroDimGrid) = 1 +Base.size(g::ZeroDimGrid) = () + + +refine(g::ZeroDimGrid, ::Int) = g +coarsen(g::ZeroDimGrid, ::Int) = g + +boundary_identifiers(g::ZeroDimGrid) = () +boundary_grid(g::ZeroDimGrid, ::Any) = throw(ArgumentError("ZeroDimGrid has no boundaries"))
--- a/src/LazyTensors/lazy_array.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/LazyTensors/lazy_array.jl Tue May 02 20:14:39 2023 +0200 @@ -36,7 +36,7 @@ function Base.getindex(lfa::LazyFunctionArray{F,T,D}, I::Vararg{Int,D}) where {F,T,D} @boundscheck checkbounds(lfa, I...) - return @inbounds lfa.f(I...) + return @inbounds @inline lfa.f(I...) end
--- a/src/SbpOperators/boundaryops/boundary_operator.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/boundaryops/boundary_operator.jl Tue May 02 20:14:39 2023 +0200 @@ -13,12 +13,12 @@ end """ - BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil, region) + BoundaryOperator(grid::EquidistantGrid, closure_stencil, region) -Constructs the BoundaryOperator with stencil `closure_stencil` for a one-dimensional `grid`, restricting to +Constructs the BoundaryOperator with stencil `closure_stencil` for a `EquidistantGrid` `grid`, restricting to to the boundary specified by `region`. """ -function BoundaryOperator(grid::EquidistantGrid{1}, closure_stencil::Stencil{T,N}, region::Region) where {T,N} +function BoundaryOperator(grid::EquidistantGrid, closure_stencil::Stencil{T,N}, region::Region) where {T,N} return BoundaryOperator{T,typeof(region),N}(closure_stencil,size(grid)[1]) end
--- a/src/SbpOperators/boundaryops/boundary_restriction.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/boundaryops/boundary_restriction.jl Tue May 02 20:14:39 2023 +0200 @@ -1,25 +1,23 @@ """ - boundary_restriction(grid, closure_stencil::Stencil, boundary) + boundary_restriction(g, closure_stencil::Stencil, boundary) Creates boundary restriction operators `e` as `LazyTensor`s on `boundary` `e` is the restriction of a grid function to `boundary` using a `Stencil` `closure_stencil`. `e'` is the prolongation of a grid function on `boundary` to the whole grid using the same `closure_stencil`. -On a one-dimensional `grid`, `e` is a `BoundaryOperator`. On a multi-dimensional `grid`, `e` is the inflation of +On a one-dimensional grid, `e` is a `BoundaryOperator`. On a multi-dimensional grid, `e` is the inflation of a `BoundaryOperator`. See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref). """ -function boundary_restriction(grid, closure_stencil, boundary) - converted_stencil = convert(Stencil{eltype(grid)}, closure_stencil) - - op = BoundaryOperator(restrict(grid, dim(boundary)), converted_stencil, region(boundary)) - return LazyTensors.inflate(op, size(grid), dim(boundary)) +#TODO: Check docstring +function boundary_restriction(g::TensorGrid, stencil_set::StencilSet, boundary::TensorGridBoundary) + op = boundary_restriction(g.grids[grid_id(boundary)], stencil_set, boundary_id(boundary)) + return LazyTensors.inflate(op, size(g), grid_id(boundary)) end -""" - boundary_restriction(grid, stencil_set, boundary) - -Creates a `boundary_restriction` operator on `grid` given a `stencil_set`. -""" -boundary_restriction(grid, stencil_set::StencilSet, boundary) = boundary_restriction(grid, parse_stencil(stencil_set["e"]["closure"]), boundary) +function boundary_restriction(g::EquidistantGrid, stencil_set::StencilSet, boundary) + closure_stencil = parse_stencil(stencil_set["e"]["closure"]) + converted_stencil = convert(Stencil{eltype(g)}, closure_stencil) + return BoundaryOperator(g, converted_stencil, boundary) +end
--- a/src/SbpOperators/boundaryops/normal_derivative.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/boundaryops/normal_derivative.jl Tue May 02 20:14:39 2023 +0200 @@ -1,5 +1,5 @@ """ - normal_derivative(grid, closure_stencil::Stencil, boundary) + normal_derivative(g, closure_stencil::Stencil, boundary) Creates the normal derivative boundary operator `d` as a `LazyTensor` @@ -10,17 +10,16 @@ See also: [`BoundaryOperator`](@ref), [`LazyTensors.inflate`](@ref). """ -function normal_derivative(grid, closure_stencil, boundary) - direction = dim(boundary) - h_inv = inverse_spacing(grid)[direction] - - op = BoundaryOperator(restrict(grid, dim(boundary)), scale(closure_stencil,h_inv), region(boundary)) - return LazyTensors.inflate(op, size(grid), dim(boundary)) +#TODO: Check docstring +function normal_derivative(g::TensorGrid, stencil_set::StencilSet, boundary::TensorGridBoundary) + op = normal_derivative(g.grids[grid_id(boundary)], stencil_set, boundary_id(boundary)) + return LazyTensors.inflate(op, size(g), grid_id(boundary)) end -""" - normal_derivative(grid, stencil_set, boundary) +function normal_derivative(g::EquidistantGrid, stencil_set::StencilSet, boundary) + closure_stencil = parse_stencil(stencil_set["d1"]["closure"]) + h_inv = inverse_spacing(g) -Creates a `normal_derivative` operator on `grid` given a `stencil_set`. -""" -normal_derivative(grid, stencil_set::StencilSet, boundary) = normal_derivative(grid, parse_stencil(stencil_set["d1"]["closure"]), boundary) + scaled_stencil = scale(closure_stencil,h_inv) + return BoundaryOperator(g, scaled_stencil, boundary) +end
--- a/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/constant_interior_scaling_operator.jl Tue May 02 20:14:39 2023 +0200 @@ -18,7 +18,7 @@ end end -function ConstantInteriorScalingOperator(grid::EquidistantGrid{1}, interior_weight, closure_weights) +function ConstantInteriorScalingOperator(grid::EquidistantGrid, interior_weight, closure_weights) return ConstantInteriorScalingOperator(interior_weight, Tuple(closure_weights), size(grid)[1]) end
--- a/src/SbpOperators/volumeops/derivatives/dissipation.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/derivatives/dissipation.jl Tue May 02 20:14:39 2023 +0200 @@ -10,30 +10,31 @@ Artificial Dissipation,” Journal of Scientific Computing, vol. 21, no. 1, pp. 57–79, Aug. 2004" """ -function undivided_skewed04(g::EquidistantGrid, p, direction) +function undivided_skewed04(g::TensorGrid, p, direction) + op = undivided_skewed04(g.grids[direction], p) + return LazyTensors.inflate(op, size(g), direction) +end + +function undivided_skewed04(g::EquidistantGrid, p) T = eltype(g) interior_weights = T.(dissipation_interior_weights(p)) - D = stencil_operator_distinct_closures( + D = StencilOperatorDistinctClosures( g, dissipation_interior_stencil(interior_weights), dissipation_lower_closure_stencils(interior_weights), dissipation_upper_closure_stencils(interior_weights), - direction, ) - Dᵀ = stencil_operator_distinct_closures( + Dᵀ = StencilOperatorDistinctClosures( g, dissipation_transpose_interior_stencil(interior_weights), dissipation_transpose_lower_closure_stencils(interior_weights), dissipation_transpose_upper_closure_stencils(interior_weights), - direction, ) return D, Dᵀ end -undivided_skewed04(g::EquidistantGrid{1}, p) = undivided_skewed04(g, p, 1) - function dissipation_interior_weights(p) if p == 0 return (1,)
--- a/src/SbpOperators/volumeops/derivatives/first_derivative.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/derivatives/first_derivative.jl Tue May 02 20:14:39 2023 +0200 @@ -1,48 +1,36 @@ """ - first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) + first_derivative(g::TensorGrid, stencil_set, direction) Creates the first-derivative operator `D1` as a `LazyTensor` -`D1` approximates the first-derivative d/dξ on `grid` along the coordinate dimension specified by +`D1` approximates the first-derivative d/dξ on `g` along the coordinate dimension specified by `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` for the points in the closure regions. -On a one-dimensional `grid`, `D1` is a `VolumeOperator`. On a multi-dimensional `grid`, `D1` is the inflation of -a `VolumeOperator`. - See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref). """ -function first_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) - h_inv = inverse_spacing(grid)[direction] - - D₁ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv), scale.(closure_stencils,h_inv), odd) - return LazyTensors.inflate(D₁, size(grid), direction) +function first_derivative(g::TensorGrid, stencil_set, direction) + D₁ = first_derivative(g.grids[direction], stencil_set) + return LazyTensors.inflate(D₁, size(g), direction) end - -""" - first_derivative(grid, inner_stencil, closure_stencils) - -Creates a `first_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`. -""" -first_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = first_derivative(grid, inner_stencil, closure_stencils, 1) - - """ - first_derivative(grid, stencil_set::StencilSet, direction) + first_derivative(g::EquidistantGrid, stencil_set) -Creates a `first_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`. +Creates a `first_derivative` operator on a `EquidistantGrid` given a `StencilSet`. """ -function first_derivative(grid::EquidistantGrid, stencil_set::StencilSet, direction) +function first_derivative(g::EquidistantGrid, stencil_set::StencilSet) inner_stencil = parse_stencil(stencil_set["D1"]["inner_stencil"]) closure_stencils = parse_stencil.(stencil_set["D1"]["closure_stencils"]) - first_derivative(grid,inner_stencil,closure_stencils,direction); + return first_derivative(g, inner_stencil, closure_stencils); end - """ - first_derivative(grid, stencil_set) + first_derivative(g::EquidistantGrid, inner_stencil, closure_stencils) -Creates a `first_derivative` operator on a 1D `grid` given a `stencil_set`. +Creates a `first_derivative` operator on a `EquidistantGrid` given an `inner_stencil` and a`closure_stencils`. """ -first_derivative(grid::EquidistantGrid{1}, stencil_set::StencilSet) = first_derivative(grid, stencil_set, 1) +function first_derivative(g::EquidistantGrid, inner_stencil::Stencil, closure_stencils) + h⁻¹ = inverse_spacing(g) + return VolumeOperator(g, scale(inner_stencil,h⁻¹), scale.(closure_stencils,h⁻¹), odd) +end
--- a/src/SbpOperators/volumeops/derivatives/second_derivative.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/derivatives/second_derivative.jl Tue May 02 20:14:39 2023 +0200 @@ -1,48 +1,36 @@ """ - second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) + second_derivative(g::EquidistantGrid, inner_stencil, closure_stencils, direction) Creates the second-derivative operator `D2` as a `LazyTensor` -`D2` approximates the second-derivative d²/dξ² on `grid` along the coordinate dimension specified by +`D2` approximates the second-derivative d²/dξ² on `g` along the coordinate dimension specified by `direction`, using the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` for the points in the closure regions. -On a one-dimensional `grid`, `D2` is a `VolumeOperator`. On a multi-dimensional `grid`, `D2` is the inflation of -a `VolumeOperator`. - See also: [`VolumeOperator`](@ref), [`LazyTensors.inflate`](@ref). """ -function second_derivative(grid::EquidistantGrid, inner_stencil, closure_stencils, direction) - h_inv = inverse_spacing(grid)[direction] - - D₂ = VolumeOperator(restrict(grid, direction), scale(inner_stencil,h_inv^2), scale.(closure_stencils,h_inv^2), even) - return LazyTensors.inflate(D₂, size(grid), direction) +function second_derivative(g::TensorGrid, stencil_set, direction) + D₂ = second_derivative(g.grids[direction], stencil_set) + return LazyTensors.inflate(D₂, size(g), direction) end - -""" - second_derivative(grid, inner_stencil, closure_stencils) - -Creates a `second_derivative` operator on a 1D `grid` given `inner_stencil` and `closure_stencils`. -""" -second_derivative(grid::EquidistantGrid{1}, inner_stencil::Stencil, closure_stencils) = second_derivative(grid, inner_stencil, closure_stencils,1) - - """ - second_derivative(grid, stencil_set, direction) + second_derivative(g, stencil_set) -Creates a `second_derivative` operator on `grid` along coordinate dimension `direction` given a `stencil_set`. +Creates a `second_derivative` operator on a 1D `g` given a `stencil_set`. """ -function second_derivative(grid::EquidistantGrid, stencil_set::StencilSet, direction) +function second_derivative(g::EquidistantGrid, stencil_set::StencilSet) inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) - second_derivative(grid,inner_stencil,closure_stencils,direction); + return second_derivative(g, inner_stencil, closure_stencils) end - """ - second_derivative(grid, stencil_set) + second_derivative(g, inner_stencil, closure_stencils) -Creates a `second_derivative` operator on a 1D `grid` given a `stencil_set`. +Creates a `second_derivative` operator on a 1D `g` given `inner_stencil` and `closure_stencils`. """ -second_derivative(grid::EquidistantGrid{1}, stencil_set::StencilSet) = second_derivative(grid, stencil_set, 1) \ No newline at end of file +function second_derivative(g::EquidistantGrid, inner_stencil::Stencil, closure_stencils) + h⁻¹ = inverse_spacing(g) + return VolumeOperator(g, scale(inner_stencil,h⁻¹^2), scale.(closure_stencils,h⁻¹^2), even) +end
--- a/src/SbpOperators/volumeops/inner_products/inner_product.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inner_product.jl Tue May 02 20:14:39 2023 +0200 @@ -15,32 +15,17 @@ See also: [`ConstantInteriorScalingOperator`](@ref). """ -function inner_product(grid::EquidistantGrid, interior_weight, closure_weights) - Hs = () - - for i ∈ dims(grid) - Hs = (Hs..., inner_product(restrict(grid, i), interior_weight, closure_weights)) - end - - return foldl(⊗, Hs) +function inner_product(tg::TensorGrid, stencil_set::StencilSet) + return mapreduce(g->inner_product(g,stencil_set), ⊗, tg.grids) end -function inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights) - h = spacing(grid)[1] +function inner_product(g::EquidistantGrid, stencil_set::StencilSet) + interior_weight = parse_scalar(stencil_set["H"]["inner"]) + closure_weights = parse_tuple(stencil_set["H"]["closure"]) - H = SbpOperators.ConstantInteriorScalingOperator(grid, h*interior_weight, h.*closure_weights) - return H + h = spacing(g) + return SbpOperators.ConstantInteriorScalingOperator(g, h*interior_weight, h.*closure_weights) end -inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityTensor{eltype(grid)}() - -""" - inner_product(grid, stencil_set) +inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}() -Creates a `inner_product` operator on `grid` given a `stencil_set`. -""" -function inner_product(grid, stencil_set::StencilSet) - inner_stencil = parse_scalar(stencil_set["H"]["inner"]) - closure_stencils = parse_tuple(stencil_set["H"]["closure"]) - return inner_product(grid, inner_stencil, closure_stencils) -end
--- a/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/inner_products/inverse_inner_product.jl Tue May 02 20:14:39 2023 +0200 @@ -12,31 +12,16 @@ See also: [`ConstantInteriorScalingOperator`](@ref). """ -function inverse_inner_product(grid::EquidistantGrid, interior_weight, closure_weights) - H⁻¹s = () - - for i ∈ dims(grid) - H⁻¹s = (H⁻¹s..., inverse_inner_product(restrict(grid, i), interior_weight, closure_weights)) - end - - return foldl(⊗, H⁻¹s) +function inverse_inner_product(tg::TensorGrid, stencil_set::StencilSet) + return mapreduce(g->inverse_inner_product(g,stencil_set), ⊗, tg.grids) end -function inverse_inner_product(grid::EquidistantGrid{1}, interior_weight, closure_weights) - h⁻¹ = inverse_spacing(grid)[1] - H⁻¹ = SbpOperators.ConstantInteriorScalingOperator(grid, h⁻¹*1/interior_weight, h⁻¹./closure_weights) - return H⁻¹ +function inverse_inner_product(g::EquidistantGrid, stencil_set::StencilSet) + interior_weight = parse_scalar(stencil_set["H"]["inner"]) + closure_weights = parse_tuple(stencil_set["H"]["closure"]) + + h⁻¹ = inverse_spacing(g) + return SbpOperators.ConstantInteriorScalingOperator(g, h⁻¹*1/interior_weight, h⁻¹./closure_weights) end -inverse_inner_product(grid::EquidistantGrid{0}, interior_weight, closure_weights) = IdentityTensor{eltype(grid)}() - -""" - inverse_inner_product(grid, stencil_set) - -Creates a `inverse_inner_product` operator on `grid` given a `stencil_set`. -""" -function inverse_inner_product(grid, stencil_set::StencilSet) - inner_stencil = parse_scalar(stencil_set["H"]["inner"]) - closure_stencils = parse_tuple(stencil_set["H"]["closure"]) - return inverse_inner_product(grid, inner_stencil, closure_stencils) -end +inverse_inner_product(g::ZeroDimGrid, stencil_set::StencilSet) = IdentityTensor{component_type(g)}()
--- a/src/SbpOperators/volumeops/laplace/laplace.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/laplace/laplace.jl Tue May 02 20:14:39 2023 +0200 @@ -3,7 +3,7 @@ Implements the Laplace operator, approximating ∑d²/xᵢ² , i = 1,...,`Dim` as a `LazyTensor`. Additionally `Laplace` stores the `StencilSet` -used to construct the `LazyTensor `. +used to construct the `LazyTensor`. """ struct Laplace{T, Dim, TM<:LazyTensor{T, Dim, Dim}} <: LazyTensor{T, Dim, Dim} D::TM # Difference operator @@ -17,11 +17,9 @@ See also [`laplace`](@ref). """ -function Laplace(grid::EquidistantGrid, stencil_set::StencilSet) - inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) - closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) - Δ = laplace(grid, inner_stencil,closure_stencils) - return Laplace(Δ,stencil_set) +function Laplace(g::Grid, stencil_set::StencilSet) + Δ = laplace(g, stencil_set) + return Laplace(Δ, stencil_set) end LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D) @@ -32,24 +30,26 @@ # Base.show(io::IO, L::Laplace) = ... """ - laplace(grid::EquidistantGrid, inner_stencil, closure_stencils) - -Creates the Laplace operator operator `Δ` as a `LazyTensor` + laplace(g::Grid, stencil_set) -`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `grid`, using -the stencil `inner_stencil` in the interior and a set of stencils `closure_stencils` -for the points in the closure regions. +Creates the Laplace operator operator `Δ` as a `LazyTensor` on the given grid -On a one-dimensional `grid`, `Δ` is equivalent to `second_derivative`. On a -multi-dimensional `grid`, `Δ` is the sum of multi-dimensional `second_derivative`s -where the sum is carried out lazily. +`Δ` approximates the Laplace operator ∑d²/xᵢ² , i = 1,...,`Dim` on `g`. The +approximation depends on the type of grid and the stencil set. See also: [`second_derivative`](@ref). """ -function laplace(grid::EquidistantGrid, inner_stencil, closure_stencils) - Δ = second_derivative(grid, inner_stencil, closure_stencils, 1) - for d = 2:ndims(grid) - Δ += second_derivative(grid, inner_stencil, closure_stencils, d) +function laplace end +function laplace(g::TensorGrid, stencil_set) + # return mapreduce(+, enumerate(g.grids)) do (i, gᵢ) + # Δᵢ = laplace(gᵢ, stencil_set) + # LazyTensors.inflate(Δᵢ, size(g), i) + # end + + Δ = LazyTensors.inflate(laplace(g.grids[1], stencil_set), size(g), 1) + for d = 2:ndims(g) + Δ += LazyTensors.inflate(laplace(g.grids[d], stencil_set), size(g), d) end return Δ end +laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
--- a/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/stencil_operator_distinct_closures.jl Tue May 02 20:14:39 2023 +0200 @@ -1,22 +1,3 @@ -""" - stencil_operator_distinct_closures( - grid::EquidistantGrid, - inner_stencil, - lower_closure, - upper_closure, - direction - ) - -Creates a multi-dimensional `StencilOperatorDistinctClosures` acting on grid -functions of `grid`. - -See also: [`StencilOperatorDistinctClosures`](@ref) -""" -function stencil_operator_distinct_closures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure, direction) - op = StencilOperatorDistinctClosures(restrict(grid, direction), inner_stencil, lower_closure, upper_closure) - return LazyTensors.inflate(op, size(grid), direction) -end - """ StencilOperatorDistinctClosures{T,K,N,M,L} <: LazyTensor{T,1} @@ -37,7 +18,7 @@ size::Tuple{Int} end -function StencilOperatorDistinctClosures(grid::EquidistantGrid{1}, inner_stencil, lower_closure, upper_closure) +function StencilOperatorDistinctClosures(grid::EquidistantGrid, inner_stencil, lower_closure, upper_closure) return StencilOperatorDistinctClosures(inner_stencil, Tuple(lower_closure), Tuple(upper_closure), size(grid)) end
--- a/src/SbpOperators/volumeops/volume_operator.jl Mon May 01 11:37:09 2023 +0200 +++ b/src/SbpOperators/volumeops/volume_operator.jl Tue May 02 20:14:39 2023 +0200 @@ -10,7 +10,7 @@ parity::Parity end -function VolumeOperator(grid::EquidistantGrid{1}, inner_stencil, closure_stencils, parity) +function VolumeOperator(grid::EquidistantGrid, inner_stencil, closure_stencils, parity) return VolumeOperator(inner_stencil, Tuple(closure_stencils), size(grid), parity) end
--- a/test/Grids/equidistant_grid_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/Grids/equidistant_grid_test.jl Tue May 02 20:14:39 2023 +0200 @@ -1,39 +1,114 @@ using Sbplib.Grids using Test using Sbplib.RegionIndices +using Sbplib.LazyTensors @testset "EquidistantGrid" begin - @test EquidistantGrid(4,0.0,1.0) isa EquidistantGrid - @test EquidistantGrid(4,0.0,8.0) isa EquidistantGrid - # constuctor - @test_throws DomainError EquidistantGrid(0,0.0,1.0) - @test_throws DomainError EquidistantGrid(1,1.0,1.0) - @test_throws DomainError EquidistantGrid(1,1.0,-1.0) - @test EquidistantGrid(4,0.0,1.0) == EquidistantGrid((4,),(0.0,),(1.0,)) + @test EquidistantGrid(0:0.1:10) isa EquidistantGrid + @test EquidistantGrid(range(0,1,length=10)) isa EquidistantGrid + @test EquidistantGrid(LinRange(0,1,11)) isa EquidistantGrid + + @testset "Indexing Interface" begin + g = EquidistantGrid(0:0.1:10) + @test g[1] == 0.0 + @test g[5] == 0.4 + @test g[101] == 10.0 + + @test g[begin] == 0.0 + @test g[end] == 10.0 + + @test all(eachindex(g) .== 1:101) + end + + @testset "Iterator interface" begin + @test eltype(EquidistantGrid(0:10)) == Int + @test eltype(EquidistantGrid(0:2:10)) == Int + @test eltype(EquidistantGrid(0:0.1:10)) == Float64 + @test size(EquidistantGrid(0:10)) == (11,) + @test size(EquidistantGrid(0:0.1:10)) == (101,) + + @test collect(EquidistantGrid(0:0.1:0.5)) == [0.0, 0.1, 0.2, 0.3, 0.4, 0.5] + + @test Base.IteratorSize(EquidistantGrid{Float64, StepRange{Float64}}) == Base.HasShape{1}() + end @testset "Base" begin - @test eltype(EquidistantGrid(4,0.0,1.0)) == Float64 - @test eltype(EquidistantGrid((4,3),(0,0),(1,3))) == Int - @test size(EquidistantGrid(4,0.0,1.0)) == (4,) - @test size(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3) - @test ndims(EquidistantGrid(4,0.0,1.0)) == 1 - @test ndims(EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0))) == 2 + @test ndims(EquidistantGrid(0:10)) == 1 end @testset "spacing" begin - @test [spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(1. /3,)...] atol=5e-13 - @test [spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(0.5, 1.)...] atol=5e-13 + @test spacing(EquidistantGrid(0:10)) == 1 + @test spacing(EquidistantGrid(0:0.1:10)) == 0.1 end @testset "inverse_spacing" begin - @test [inverse_spacing(EquidistantGrid(4,0.0,1.0))...] ≈ [(3.,)...] atol=5e-13 - @test [inverse_spacing(EquidistantGrid((5,3), (0.0,-1.0), (2.0,1.0)))...] ≈ [(2, 1.)...] atol=5e-13 + @test inverse_spacing(EquidistantGrid(0:10)) == 1 + @test inverse_spacing(EquidistantGrid(0:0.1:10)) == 10 + end + + @testset "boundary_identifiers" begin + g = EquidistantGrid(0:0.1:10) + @test boundary_identifiers(g) == (Lower(), Upper()) + @inferred boundary_identifiers(g) + end + + @testset "boundary_grid" begin + g = EquidistantGrid(0:0.1:1) + @test boundary_grid(g, Lower()) == ZeroDimGrid(0.0) + @test boundary_grid(g, Upper()) == ZeroDimGrid(1.0) + end + + @testset "refine" begin + g = EquidistantGrid(0:0.1:1) + @test refine(g, 1) == g + @test refine(g, 2) == EquidistantGrid(0:0.05:1) + @test refine(g, 3) == EquidistantGrid(0:(0.1/3):1) end - @testset "points" begin - g = EquidistantGrid((5,3), (-1.0,0.0), (0.0,7.11)) - gp = points(g); + @testset "coarsen" begin + g = EquidistantGrid(0:1:10) + @test coarsen(g, 1) == g + @test coarsen(g, 2) == EquidistantGrid(0:2:10) + + g = EquidistantGrid(0:0.1:1) + @test coarsen(g, 1) == g + @test coarsen(g, 2) == EquidistantGrid(0:0.2:1) + + g = EquidistantGrid(0:10) + @test coarsen(g, 1) == EquidistantGrid(0:1:10) + @test coarsen(g, 2) == EquidistantGrid(0:2:10) + + @test_throws DomainError(3, "Size minus 1 must be divisible by the ratio.") coarsen(g, 3) + end +end + + +@testset "equidistant_grid" begin + @test equidistant_grid(4,0.0,1.0) isa EquidistantGrid + @test equidistant_grid((4,3),(0.0,0.0),(8.0,5.0)) isa TensorGrid + + # constuctor + @test_throws DomainError equidistant_grid(0,0.0,1.0) + @test_throws DomainError equidistant_grid(1,1.0,1.0) + @test_throws DomainError equidistant_grid(1,1.0,-1.0) + + @test_throws DomainError equidistant_grid((0,0),(0.0,0.0),(1.0,1.0)) + @test_throws DomainError equidistant_grid((1,1),(1.0,1.0),(1.0,1.0)) + @test_throws DomainError equidistant_grid((1,1),(1.0,1.0),(-1.0,-1.0)) + + @testset "Base" begin + @test eltype(equidistant_grid(4,0.0,1.0)) == Float64 + @test eltype(equidistant_grid((4,3),(0,0),(1,3))) <: AbstractVector{Float64} + @test size(equidistant_grid(4,0.0,1.0)) == (4,) + @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3) + @test ndims(equidistant_grid(4,0.0,1.0)) == 1 + @test ndims(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == 2 + end + + @testset "getindex" begin + g = equidistant_grid((5,3), (-1.0,0.0), (0.0,7.11)) + gp = collect(g); p = [(-1.,0.) (-1.,7.11/2) (-1.,7.11); (-0.75,0.) (-0.75,7.11/2) (-0.75,7.11); (-0.5,0.) (-0.5,7.11/2) (-0.5,7.11); @@ -43,83 +118,18 @@ @test [gp[i]...] ≈ [p[i]...] atol=5e-13 end end +end - @testset "restrict" begin - g = EquidistantGrid((5,3), (0.0,0.0), (2.0,1.0)) - @test restrict(g, 1) == EquidistantGrid(5,0.0,2.0) - @test restrict(g, 2) == EquidistantGrid(3,0.0,1.0) - g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) - @test restrict(g, 1) == EquidistantGrid(2,0.0,2.0) - @test restrict(g, 2) == EquidistantGrid(5,0.0,1.0) - @test restrict(g, 3) == EquidistantGrid(3,0.0,3.0) - @test restrict(g, 1:2) == EquidistantGrid((2,5),(0.0,0.0),(2.0,1.0)) - @test restrict(g, 2:3) == EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) - @test restrict(g, [1,3]) == EquidistantGrid((2,3),(0.0,0.0),(2.0,3.0)) - @test restrict(g, [2,1]) == EquidistantGrid((5,2),(0.0,0.0),(1.0,2.0)) - end - - @testset "boundary_identifiers" begin - g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) - bids = (CartesianBoundary{1,Lower}(),CartesianBoundary{1,Upper}(), - CartesianBoundary{2,Lower}(),CartesianBoundary{2,Upper}(), - CartesianBoundary{3,Lower}(),CartesianBoundary{3,Upper}()) - @test boundary_identifiers(g) == bids - @inferred boundary_identifiers(g) - end +@testset "change_length" begin + @test Grids.change_length(0:20, 21) == 0:20 + @test Grids.change_length(0:20, 11) == 0:2:20 + @test Grids.change_length(0:2:20, 21) == 0:20 - @testset "boundary_grid" begin - @testset "1D" begin - g = EquidistantGrid(5,0.0,2.0) - (id_l, id_r) = boundary_identifiers(g) - @test boundary_grid(g,id_l) == EquidistantGrid{Float64}() - @test boundary_grid(g,id_r) == EquidistantGrid{Float64}() - @test_throws DomainError boundary_grid(g,CartesianBoundary{2,Lower}()) - @test_throws DomainError boundary_grid(g,CartesianBoundary{0,Lower}()) - end - @testset "2D" begin - g = EquidistantGrid((5,3),(0.0,0.0),(1.0,3.0)) - (id_w, id_e, id_s, id_n) = boundary_identifiers(g) - @test boundary_grid(g,id_w) == restrict(g,2) - @test boundary_grid(g,id_e) == restrict(g,2) - @test boundary_grid(g,id_s) == restrict(g,1) - @test boundary_grid(g,id_n) == restrict(g,1) - @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}()) - end - @testset "3D" begin - g = EquidistantGrid((2,5,3), (0.0,0.0,0.0), (2.0,1.0,3.0)) - (id_w, id_e, - id_s, id_n, - id_t, id_b) = boundary_identifiers(g) - @test boundary_grid(g,id_w) == restrict(g,[2,3]) - @test boundary_grid(g,id_e) == restrict(g,[2,3]) - @test boundary_grid(g,id_s) == restrict(g,[1,3]) - @test boundary_grid(g,id_n) == restrict(g,[1,3]) - @test boundary_grid(g,id_t) == restrict(g,[1,2]) - @test boundary_grid(g,id_b) == restrict(g,[1,2]) - @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}()) - end - end + @test Grids.change_length(range(0,1,length=10), 10) == range(0,1,length=10) + @test Grids.change_length(range(0,1,length=10), 5) == range(0,1,length=5) + @test Grids.change_length(range(0,1,length=10), 20) == range(0,1,length=20) - @testset "refine" begin - @test refine(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}() - @test refine(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}() - - g = EquidistantGrid((10,5),(0.,1.),(2.,3.)) - @test refine(g, 1) == g - @test refine(g, 2) == EquidistantGrid((19,9),(0.,1.),(2.,3.)) - @test refine(g, 3) == EquidistantGrid((28,13),(0.,1.),(2.,3.)) - end - - @testset "coarsen" begin - @test coarsen(EquidistantGrid{Float64}(), 1) == EquidistantGrid{Float64}() - @test coarsen(EquidistantGrid{Float64}(), 2) == EquidistantGrid{Float64}() - - g = EquidistantGrid((7,13),(0.,1.),(2.,3.)) - @test coarsen(g, 1) == g - @test coarsen(g, 2) == EquidistantGrid((4,7),(0.,1.),(2.,3.)) - @test coarsen(g, 3) == EquidistantGrid((3,5),(0.,1.),(2.,3.)) - - @test_throws DomainError(4, "Size minus 1 must be divisible by the ratio.") coarsen(g, 4) == EquidistantGrid((3,5),(0.,1.),(2.,3.)) - end + @test Grids.change_length(LinRange(1,2,10),10) == LinRange(1,2,10) + @test Grids.change_length(LinRange(1,2,10),15) == LinRange(1,2,15) end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/grid_test.jl Tue May 02 20:14:39 2023 +0200 @@ -0,0 +1,66 @@ +using Test +using Sbplib.Grids +using Sbplib.LazyTensors +using StaticArrays + +@testset "Grid" begin + struct DummyGrid{T,D} <: Grid{T,D} end + + @test eltype(DummyGrid{Int, 2}) == Int + @test eltype(DummyGrid{Int, 2}()) == Int + + @test ndims(DummyGrid{Int, 2}()) == 2 + + @test target_manifold_dim(DummyGrid{Int, 1}()) == 1 + @test target_manifold_dim(DummyGrid{SVector{3,Float64}, 2}()) == 3 + + @testset "component_type" begin + @test component_type(DummyGrid{Int,1}()) == Int + @test component_type(DummyGrid{Float64,1}()) == Float64 + @test component_type(DummyGrid{Rational,1}()) == Rational + + @test component_type(DummyGrid{SVector{3,Int},2}()) == Int + @test component_type(DummyGrid{SVector{2,Float64},3}()) == Float64 + @test component_type(DummyGrid{SVector{4,Rational},4}()) == Rational + end +end + +@testset "eval_on" begin + @test eval_on(ZeroDimGrid(@SVector[1.,2.]), x̄->x̄[1]+x̄[2]) isa LazyArray + @test eval_on(ZeroDimGrid(@SVector[1.,2.]), x̄->x̄[1]+x̄[2]) == fill(3.) + @test eval_on(ZeroDimGrid(@SVector[3.,2.]), x̄->x̄[1]+x̄[2]) == fill(5.) + + @test eval_on(ZeroDimGrid(1.), x̄->2x̄) isa LazyArray + @test eval_on(ZeroDimGrid(1.), x̄->2x̄) == fill(2.) + + @test eval_on(EquidistantGrid(range(0,1,length=4)), x->2x) isa LazyArray + @test eval_on(EquidistantGrid(range(0,1,length=4)), x->2x) == 2 .* range(0,1,length=4) + + + g = equidistant_grid((5,3), (0.0,0.0), (2.0,1.0)) + + @test eval_on(g, x̄ -> 0.) isa LazyArray + @test eval_on(g, x̄ -> 0.) == fill(0., (5,3)) + + @test eval_on(g, x̄ -> sin(x̄[1])*cos(x̄[2])) == map(x̄->sin(x̄[1])*cos(x̄[2]), g) + + # Vector valued function + @test eval_on(g, x̄ -> @SVector[x̄[2], x̄[1]]) isa LazyArray{SVector{2,Float64}} + @test eval_on(g, x̄ -> @SVector[x̄[2], x̄[1]]) == map(x̄ -> @SVector[x̄[2], x̄[1]], g) + + # Multi-argument functions + f(x,y) = sin(x)*cos(y) + @test eval_on(g, f) == map(x̄->f(x̄...), g) + + #TODO: inference test! +end + +@testset "_ncomponents" begin + @test Grids._ncomponents(Int) == 1 + @test Grids._ncomponents(Float64) == 1 + @test Grids._ncomponents(Rational) == 1 + + @test Grids._ncomponents(SVector{3,Int}) == 3 + @test Grids._ncomponents(SVector{2,Float64}) == 2 + @test Grids._ncomponents(SVector{4,Rational}) == 4 +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/tensor_grid_test.jl Tue May 02 20:14:39 2023 +0200 @@ -0,0 +1,126 @@ +using Test +using Sbplib.Grids +using StaticArrays +using Sbplib.RegionIndices + +@testset "TensorGrid" begin + g₁ = EquidistantGrid(range(0,1,length=11)) + g₂ = EquidistantGrid(range(2,3,length=6)) + g₃ = EquidistantGrid(1:10) + g₄ = ZeroDimGrid(@SVector[1,2]) + + @test TensorGrid(g₁, g₂) isa TensorGrid + @test TensorGrid(g₁, g₂) isa Grid{SVector{2,Float64}, 2} + @test TensorGrid(g₃, g₃) isa Grid{SVector{2,Int}, 2} + @test TensorGrid(g₁, g₂, g₃) isa Grid{SVector{3,Float64}, 3} + @test TensorGrid(g₁, g₄) isa Grid{SVector{3,Float64}, 1} + @test TensorGrid(g₁, g₄, g₂) isa Grid{SVector{4,Float64}, 2} + + @testset "Indexing Interface" begin + @test TensorGrid(g₁, g₂)[1,1] isa SVector{2,Float64} + @test TensorGrid(g₁, g₂)[1,1] == [0.0,2.0] + @test TensorGrid(g₁, g₂)[3,5] == [0.2,2.8] + @test TensorGrid(g₁, g₂)[10,6] == [0.9,3.0] + + @test TensorGrid(g₁, g₃)[1,1] isa SVector{2,Float64} + @test TensorGrid(g₁, g₃)[1,1] == [0.0,1.0] + + @test TensorGrid(g₁, g₂, g₃)[3,4,5] isa SVector{3,Float64} + @test TensorGrid(g₁, g₂, g₃)[3,4,5] == [0.2, 2.6, 5.0] + + @test TensorGrid(g₁, g₄)[3] isa SVector{3,Float64} + @test TensorGrid(g₁, g₄)[3] == [0.2, 1., 2.] + + @test TensorGrid(g₁, g₄, g₂)[3,2] isa SVector{4,Float64} + @test TensorGrid(g₁, g₄, g₂)[3,2] == [0.2, 1., 2., 2.2] + + @test eachindex(TensorGrid(g₁, g₂)) == CartesianIndices((11,6)) + @test eachindex(TensorGrid(g₁, g₃)) == CartesianIndices((11,10)) + @test eachindex(TensorGrid(g₁, g₂, g₃)) == CartesianIndices((11,6,10)) + @test eachindex(TensorGrid(g₁, g₄)) == CartesianIndices((11,)) + @test eachindex(TensorGrid(g₁, g₄, g₂)) == CartesianIndices((11,6)) + end + + @testset "Iterator interface" begin + @test eltype(TensorGrid(g₁, g₂)) == SVector{2,Float64} + @test eltype(TensorGrid(g₁, g₃)) == SVector{2,Float64} + @test eltype(TensorGrid(g₁, g₂, g₃)) == SVector{3,Float64} + @test eltype(TensorGrid(g₁, g₄)) == SVector{3,Float64} + @test eltype(TensorGrid(g₁, g₄, g₂)) == SVector{4,Float64} + + @test size(TensorGrid(g₁, g₂)) == (11,6) + @test size(TensorGrid(g₁, g₃)) == (11,10) + @test size(TensorGrid(g₁, g₂, g₃)) == (11,6,10) + @test size(TensorGrid(g₁, g₄)) == (11,) + @test size(TensorGrid(g₁, g₄, g₂)) == (11,6) + + @test Base.IteratorSize(TensorGrid(g₁, g₂)) == Base.HasShape{2}() + @test Base.IteratorSize(TensorGrid(g₁, g₃)) == Base.HasShape{2}() + @test Base.IteratorSize(TensorGrid(g₁, g₂, g₃)) == Base.HasShape{3}() + @test Base.IteratorSize(TensorGrid(g₁, g₄)) == Base.HasShape{1}() + @test Base.IteratorSize(TensorGrid(g₁, g₄, g₂)) == Base.HasShape{2}() + + @test iterate(TensorGrid(g₁, g₂))[1] isa SVector{2,Float64} + @test iterate(TensorGrid(g₁, g₃))[1] isa SVector{2,Float64} + @test iterate(TensorGrid(g₁, g₂, g₃))[1] isa SVector{3,Float64} + @test iterate(TensorGrid(g₁, g₄))[1] isa SVector{3,Float64} + @test iterate(TensorGrid(g₁, g₄, g₂))[1] isa SVector{4,Float64} + + @test collect(TensorGrid(g₁, g₂)) == [@SVector[x,y] for x ∈ range(0,1,length=11), y ∈ range(2,3,length=6)] + @test collect(TensorGrid(g₁, g₃)) == [@SVector[x,y] for x ∈ range(0,1,length=11), y ∈ 1:10] + @test collect(TensorGrid(g₁, g₂, g₃)) == [@SVector[x,y,z] for x ∈ range(0,1,length=11), y ∈ range(2,3,length=6), z ∈ 1:10] + @test collect(TensorGrid(g₁, g₄)) == [@SVector[x,1,2] for x ∈ range(0,1,length=11)] + @test collect(TensorGrid(g₁, g₄, g₂)) == [@SVector[x,1,2,y] for x ∈ range(0,1,length=11), y ∈ range(2,3,length=6)] + end + + @testset "refine" begin + g1(n) = EquidistantGrid(range(0,1,length=n)) + g2(n) = EquidistantGrid(range(2,3,length=n)) + + @test refine(TensorGrid(g1(11), g2(6)),1) == TensorGrid(g1(11), g2(6)) + @test refine(TensorGrid(g1(11), g2(6)),2) == TensorGrid(g1(21), g2(11)) + @test refine(TensorGrid(g1(11), g2(6)),3) == TensorGrid(g1(31), g2(16)) + @test refine(TensorGrid(g1(11), g₄), 1) == TensorGrid(g1(11), g₄) + @test refine(TensorGrid(g1(11), g₄), 2) == TensorGrid(g1(21), g₄) + end + + @testset "coarsen" begin + g1(n) = EquidistantGrid(range(0,1,length=n)) + g2(n) = EquidistantGrid(range(2,3,length=n)) + + @test coarsen(TensorGrid(g1(11), g2(6)),1) == TensorGrid(g1(11), g2(6)) + @test coarsen(TensorGrid(g1(21), g2(11)),2) == TensorGrid(g1(11), g2(6)) + @test coarsen(TensorGrid(g1(31), g2(16)),3) == TensorGrid(g1(11), g2(6)) + @test coarsen(TensorGrid(g1(11), g₄), 1) == TensorGrid(g1(11), g₄) + @test coarsen(TensorGrid(g1(21), g₄), 2) == TensorGrid(g1(11), g₄) + end + + @testset "boundary_identifiers" begin + @test boundary_identifiers(TensorGrid(g₁, g₂)) == map((n,id)->TensorGridBoundary{n,id}(), (1,1,2,2), (Lower,Upper,Lower,Upper)) + @test boundary_identifiers(TensorGrid(g₁, g₄)) == (TensorGridBoundary{1,Lower}(),TensorGridBoundary{1,Upper}()) + end + + @testset "boundary_grid" begin + @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂) + @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄) + end +end + +@testset "combined_coordinate_vector_type" begin + @test Grids.combined_coordinate_vector_type(Float64) == Float64 + @test Grids.combined_coordinate_vector_type(Float64, Int) == SVector{2,Float64} + @test Grids.combined_coordinate_vector_type(Float32, Int16, Int32) == SVector{3,Float32} + + @test Grids.combined_coordinate_vector_type(SVector{2,Float64}) == SVector{2,Float64} + @test Grids.combined_coordinate_vector_type(SVector{2,Float64}, SVector{1,Float64}) == SVector{3,Float64} + @test Grids.combined_coordinate_vector_type(SVector{2,Float64}, SVector{1,Int}, SVector{3, Float32}) == SVector{6,Float64} +end + +@testset "combine_coordinates" begin + @test Grids.combine_coordinates(1,2,3) isa SVector{3, Int} + @test Grids.combine_coordinates(1,2,3) == [1,2,3] + @test Grids.combine_coordinates(1,2.,3) isa SVector{3, Float64} + @test Grids.combine_coordinates(1,2.,3) == [1,2,3] + @test Grids.combine_coordinates(1,@SVector[2.,3]) isa SVector{3, Float64} + @test Grids.combine_coordinates(1,@SVector[2.,3]) == [1,2,3] +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/zero_dim_grid_test.jl Tue May 02 20:14:39 2023 +0200 @@ -0,0 +1,44 @@ +using Test +using Sbplib.Grids +using StaticArrays + +@testset "ZeroDimGrid" begin + @test ZeroDimGrid(1) isa ZeroDimGrid{Int} + @test ZeroDimGrid([1,2,3]) isa ZeroDimGrid{Vector{Int}} + @test ZeroDimGrid(@SVector[1.0,2.0]) isa ZeroDimGrid{SVector{2,Float64}} + + @testset "Indexing Interface" begin + g = ZeroDimGrid(@SVector[1,2]) + + @test g[] == [1,2] + @test eachindex(g) == CartesianIndices(()) + end + + @testset "Iterator interface" begin + g = ZeroDimGrid(@SVector[1,2]) + + @test Base.IteratorSize(g) == Base.HasShape{0}() + @test eltype(g) == SVector{2,Int} + @test length(g) == 1 + @test size(g) == () + @test collect(g) == fill(@SVector[1,2]) + end + + @testset "refine" begin + @test refine(ZeroDimGrid(@SVector[1.0,2.0]),1) == ZeroDimGrid(@SVector[1.0,2.0]) + @test refine(ZeroDimGrid(@SVector[1.0,2.0]),2) == ZeroDimGrid(@SVector[1.0,2.0]) + end + + @testset "coarsen" begin + @test coarsen(ZeroDimGrid(@SVector[1.0,2.0]),1) == ZeroDimGrid(@SVector[1.0,2.0]) + @test coarsen(ZeroDimGrid(@SVector[1.0,2.0]),2) == ZeroDimGrid(@SVector[1.0,2.0]) + end + + @testset "boundary_identifiers" begin + @test boundary_identifiers(ZeroDimGrid(@SVector[1.0,2.0])) == () + end + + @testset "boundary_grid" begin + @test_throws ArgumentError("ZeroDimGrid has no boundaries") boundary_grid(ZeroDimGrid(@SVector[1.0,2.0]), :bid) + end +end
--- a/test/Manifest.toml Mon May 01 11:37:09 2023 +0200 +++ b/test/Manifest.toml Tue May 02 20:14:39 2023 +0200 @@ -2,7 +2,7 @@ julia_version = "1.8.5" manifest_format = "2.0" -project_hash = "23260eda65ade7d11fffed313a68520d0bc053fc" +project_hash = "f2b0634c12bbed93a17efc88d466604d5a07c465" [[deps.ArgTools]] uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f" @@ -250,6 +250,17 @@ uuid = "276daf66-3868-5448-9aa4-cd146d93841b" version = "2.1.7" +[[deps.StaticArrays]] +deps = ["LinearAlgebra", "Random", "StaticArraysCore", "Statistics"] +git-tree-sha1 = "2d7d9e1ddadc8407ffd460e24218e37ef52dd9a3" +uuid = "90137ffa-7385-5640-81b9-e52037218182" +version = "1.5.16" + +[[deps.StaticArraysCore]] +git-tree-sha1 = "6b7ba252635a5eff6a0b0664a41ee140a1c9e72a" +uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" +version = "1.4.0" + [[deps.Statistics]] deps = ["LinearAlgebra", "SparseArrays"] uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
--- a/test/Project.toml Mon May 01 11:37:09 2023 +0200 +++ b/test/Project.toml Tue May 02 20:14:39 2023 +0200 @@ -2,6 +2,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
--- a/test/SbpOperators/boundaryops/boundary_operator_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/boundaryops/boundary_operator_test.jl Tue May 02 20:14:39 2023 +0200 @@ -10,8 +10,7 @@ @testset "BoundaryOperator" begin closure_stencil = Stencil(2.,1.,3.; center = 1) - g_1D = EquidistantGrid(11, 0.0, 1.0) - g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) + g_1D = EquidistantGrid(range(0,1,length=11)) @testset "Constructors" begin @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T @@ -30,7 +29,7 @@ end @testset "Application" begin - v = evalOn(g_1D,x->1+x^2) + v = eval_on(g_1D,x->1+x^2) u = fill(3.124) @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3] @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2] @@ -38,7 +37,7 @@ @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]] - v = evalOn(g_1D, x->1. +x*im) + v = eval_on(g_1D, x->1. +x*im) @test (op_l*v)[] isa ComplexF64 u = fill(1. +im)
--- a/test/SbpOperators/boundaryops/boundary_restriction_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/boundaryops/boundary_restriction_test.jl Tue May 02 20:14:39 2023 +0200 @@ -9,27 +9,24 @@ @testset "boundary_restriction" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4) e_closure = parse_stencil(stencil_set["e"]["closure"]) - g_1D = EquidistantGrid(11, 0.0, 1.0) - g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) + g_1D = equidistant_grid(11, 0.0, 1.0) + g_2D = equidistant_grid((11,15), (0.0, 0.0), (1.0,1.0)) @testset "boundary_restriction" begin @testset "1D" begin - e_l = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Lower}()) - @test e_l == boundary_restriction(g_1D,stencil_set,CartesianBoundary{1,Lower}()) + e_l = boundary_restriction(g_1D,stencil_set,Lower()) @test e_l == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Lower()) @test e_l isa BoundaryOperator{T,Lower} where T @test e_l isa LazyTensor{T,0,1} where T - e_r = boundary_restriction(g_1D,e_closure,CartesianBoundary{1,Upper}()) - @test e_r == boundary_restriction(g_1D,stencil_set,CartesianBoundary{1,Upper}()) + e_r = boundary_restriction(g_1D,stencil_set,Upper()) @test e_r == BoundaryOperator(g_1D,Stencil{Float64}(e_closure),Upper()) @test e_r isa BoundaryOperator{T,Upper} where T @test e_r isa LazyTensor{T,0,1} where T end @testset "2D" begin - e_w = boundary_restriction(g_2D,e_closure,CartesianBoundary{1,Upper}()) - @test e_w == boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}()) + e_w = boundary_restriction(g_2D,stencil_set,CartesianBoundary{1,Upper}()) @test e_w isa InflatedTensor @test e_w isa LazyTensor{T,1,2} where T end @@ -37,8 +34,8 @@ @testset "Application" begin @testset "1D" begin - e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(e_closure), boundary_identifiers(g_1D)) - v = evalOn(g_1D,x->1+x^2) + e_l, e_r = boundary_restriction.(Ref(g_1D), Ref(stencil_set), boundary_identifiers(g_1D)) + v = eval_on(g_1D,x->1+x^2) u = fill(3.124) @test (e_l*v)[] == v[1] @@ -47,7 +44,7 @@ end @testset "2D" begin - e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(e_closure), boundary_identifiers(g_2D)) + e_w, e_e, e_s, e_n = boundary_restriction.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) v = rand(11, 15) u = fill(3.124)
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl Tue May 02 20:14:39 2023 +0200 @@ -7,24 +7,23 @@ import Sbplib.SbpOperators.BoundaryOperator @testset "normal_derivative" begin - g_1D = EquidistantGrid(11, 0.0, 1.0) - g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0)) + g_1D = equidistant_grid(11, 0.0, 1.0) + g_2D = equidistant_grid((11,12), (0.0, 0.0), (1.0,1.0)) @testset "normal_derivative" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - d_closure = parse_stencil(stencil_set["d1"]["closure"]) @testset "1D" begin - d_l = normal_derivative(g_1D, d_closure, CartesianBoundary{1,Lower}()) - @test d_l == normal_derivative(g_1D, stencil_set, CartesianBoundary{1,Lower}()) + d_l = normal_derivative(g_1D, stencil_set, Lower()) + @test d_l == normal_derivative(g_1D, stencil_set, Lower()) @test d_l isa BoundaryOperator{T,Lower} where T @test d_l isa LazyTensor{T,0,1} where T end @testset "2D" begin - d_w = normal_derivative(g_2D, d_closure, CartesianBoundary{1,Lower}()) - d_n = normal_derivative(g_2D, d_closure, CartesianBoundary{2,Upper}()) + d_w = normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) + d_n = normal_derivative(g_2D, stencil_set, CartesianBoundary{2,Upper}()) Ix = IdentityTensor{Float64}((size(g_2D)[1],)) Iy = IdentityTensor{Float64}((size(g_2D)[2],)) - d_l = normal_derivative(restrict(g_2D,1),d_closure,CartesianBoundary{1,Lower}()) - d_r = normal_derivative(restrict(g_2D,2),d_closure,CartesianBoundary{1,Upper}()) + d_l = normal_derivative(g_2D.grids[1], stencil_set, Lower()) + d_r = normal_derivative(g_2D.grids[2], stencil_set, Upper()) @test d_w == normal_derivative(g_2D, stencil_set, CartesianBoundary{1,Lower}()) @test d_w == d_l⊗Iy @test d_n == Ix⊗d_r @@ -33,14 +32,13 @@ end end @testset "Accuracy" begin - v = evalOn(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y) - v∂x = evalOn(g_2D, (x,y)-> 2*x + y) - v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x) + v = eval_on(g_2D, (x,y)-> x^2 + (y-1)^2 + x*y) + v∂x = eval_on(g_2D, (x,y)-> 2*x + y) + v∂y = eval_on(g_2D, (x,y)-> 2*(y-1) + x) # TODO: Test for higher order polynomials? @testset "2nd order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - d_closure = parse_stencil(stencil_set["d1"]["closure"]) - d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(d_closure), boundary_identifiers(g_2D)) + d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) @test d_w*v ≈ -v∂x[1,:] atol = 1e-13 @test d_e*v ≈ v∂x[end,:] atol = 1e-13 @@ -50,8 +48,7 @@ @testset "4th order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - d_closure = parse_stencil(stencil_set["d1"]["closure"]) - d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(d_closure), boundary_identifiers(g_2D)) + d_w, d_e, d_s, d_n = normal_derivative.(Ref(g_2D), Ref(stencil_set), boundary_identifiers(g_2D)) @test d_w*v ≈ -v∂x[1,:] atol = 1e-13 @test d_e*v ≈ v∂x[end,:] atol = 1e-13
--- a/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/constant_interior_scaling_operator_test.jl Tue May 02 20:14:39 2023 +0200 @@ -33,7 +33,7 @@ @test_throws DomainError ConstantInteriorScalingOperator(4,(2,3), 3) @testset "Grid constructor" begin - g = EquidistantGrid(11, 0., 2.) + g = equidistant_grid(11, 0., 2.) @test ConstantInteriorScalingOperator(g, 3., (.1,.2)) isa ConstantInteriorScalingOperator{Float64} end end
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl Tue May 02 20:14:39 2023 +0200 @@ -27,7 +27,7 @@ end @testset "undivided_skewed04" begin - g = EquidistantGrid(20, 0., 11.) + g = equidistant_grid(20, 0., 11.) D,Dᵀ = undivided_skewed04(g, 1) @test D isa LazyTensor{Float64,1,1} @@ -35,18 +35,20 @@ @testset "Accuracy conditions" begin N = 20 - g = EquidistantGrid(N, 0//1,2//1) + g = equidistant_grid(N, 0//1,2//1) h = only(spacing(g)) @testset "D_$p" for p ∈ [1,2,3,4] D,Dᵀ = undivided_skewed04(g, p) @testset "x^$k" for k ∈ 0:p - v = evalOn(g, x->monomial(x,k)) - vₚₓ = evalOn(g, x->monomial(x,k-p)) + v = eval_on(g, x->monomial(x,k)) + vₚₓ = eval_on(g, x->monomial(x,k-p)) @test D*v == h^p * vₚₓ end end + + # TODO: Add 2D tests end @testset "transpose equality" begin @@ -67,7 +69,7 @@ return Dmat end - g = EquidistantGrid(11, 0., 1.) + g = equidistant_grid(11, 0., 1.) @testset "D_$p" for p ∈ [1,2,3,4] D,Dᵀ = undivided_skewed04(g, p)
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl Tue May 02 20:14:39 2023 +0200 @@ -24,73 +24,71 @@ @testset "Constructors" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - g₁ = EquidistantGrid(11, 0., 1.) - g₂ = EquidistantGrid((11,14), (0.,1.), (1.,3.)) + g₁ = equidistant_grid(11, 0., 1.) + g₂ = equidistant_grid((11,14), (0.,1.), (1.,3.)) - @test first_derivative(g₁, stencil_set, 1) isa LazyTensor{Float64,1,1} + @test first_derivative(g₁, stencil_set) isa LazyTensor{Float64,1,1} @test first_derivative(g₂, stencil_set, 2) isa LazyTensor{Float64,2,2} - @test first_derivative(g₁, stencil_set, 1) == first_derivative(g₁, stencil_set) interior_stencil = CenteredStencil(-1,0,1) closure_stencils = [Stencil(-1,1, center=1)] - @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa LazyTensor{Float64,1,1} - @test first_derivative(g₁, interior_stencil, closure_stencils, 1) isa VolumeOperator - @test first_derivative(g₁, interior_stencil, closure_stencils, 1) == first_derivative(g₁, interior_stencil, closure_stencils) - @test first_derivative(g₂, interior_stencil, closure_stencils, 2) isa LazyTensor{Float64,2,2} + @test first_derivative(g₁, interior_stencil, closure_stencils) isa LazyTensor{Float64,1,1} end @testset "Accuracy conditions" begin N = 20 - g = EquidistantGrid(N, 0//1,2//1) + g = equidistant_grid(N, 0//1,2//1) @testset for order ∈ [2,4] stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) - D₁ = first_derivative(g, stencil_set, 1) + D₁ = first_derivative(g, stencil_set) @testset "boundary x^$k" for k ∈ 0:order÷2 - v = evalOn(g, x->monomial(x,k)) + v = eval_on(g, x->monomial(x,k)) @testset for i ∈ 1:closure_size(D₁) - x, = points(g)[i] + x, = g[i] @test (D₁*v)[i] == monomial(x,k-1) end @testset for i ∈ (N-closure_size(D₁)+1):N - x, = points(g)[i] + x, = g[i] @test (D₁*v)[i] == monomial(x,k-1) end end @testset "interior x^$k" for k ∈ 0:order - v = evalOn(g, x->monomial(x,k)) + v = eval_on(g, x->monomial(x,k)) - x, = points(g)[10] + x, = g[10] @test (D₁*v)[10] == monomial(x,k-1) end end end @testset "Accuracy on function" begin - # 1D - g = EquidistantGrid(30, 0.,1.) - v = evalOn(g, x->exp(x)) - @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)] - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) - D₁ = first_derivative(g, stencil_set, 1) + @testset "1D" begin + g = equidistant_grid(30, 0.,1.) + v = eval_on(g, x->exp(x)) + @testset for (order, tol) ∈ [(2, 6e-3),(4, 2e-4)] + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) + D₁ = first_derivative(g, stencil_set) - @test D₁*v ≈ v rtol=tol + @test D₁*v ≈ v rtol=tol + end end - # 2D - g = EquidistantGrid((30,60), (0.,0.),(1.,2.)) - v = evalOn(g, (x,y)->exp(0.8x+1.2*y)) - @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)] - stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) - Dx = first_derivative(g, stencil_set, 1) - Dy = first_derivative(g, stencil_set, 2) + @testset "2D" begin + g = equidistant_grid((30,60), (0.,0.),(1.,2.)) + v = eval_on(g, (x,y)->exp(0.8x+1.2*y)) + @testset for (order, tol) ∈ [(2, 6e-3),(4, 3e-4)] + stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order) + Dx = first_derivative(g, stencil_set, 1) + Dy = first_derivative(g, stencil_set, 2) - @test Dx*v ≈ 0.8v rtol=tol - @test Dy*v ≈ 1.2v rtol=tol + @test Dx*v ≈ 0.8v rtol=tol + @test Dy*v ≈ 1.2v rtol=tol + end end end end
--- a/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/derivatives/second_derivative_test.jl Tue May 02 20:14:39 2023 +0200 @@ -15,24 +15,18 @@ closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) Lx = 3.5 Ly = 3. - g_1D = EquidistantGrid(121, 0.0, Lx) - g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly)) + g_1D = equidistant_grid(121, 0.0, Lx) + g_2D = equidistant_grid((121,123), (0.0, 0.0), (Lx, Ly)) @testset "Constructors" begin @testset "1D" begin - Dₓₓ = second_derivative(g_1D,inner_stencil,closure_stencils,1) - @test Dₓₓ == second_derivative(g_1D,inner_stencil,closure_stencils) - @test Dₓₓ == second_derivative(g_1D,stencil_set,1) - @test Dₓₓ == second_derivative(g_1D,stencil_set) - @test Dₓₓ isa VolumeOperator + Dₓₓ = second_derivative(g_1D, stencil_set) + @test Dₓₓ == second_derivative(g_1D, inner_stencil, closure_stencils) + @test Dₓₓ isa LazyTensor{Float64,1,1} end @testset "2D" begin - Dₓₓ = second_derivative(g_2D,inner_stencil,closure_stencils,1) - D2 = second_derivative(g_1D,inner_stencil,closure_stencils,1) - I = IdentityTensor{Float64}(size(g_2D)[2]) - @test Dₓₓ == D2⊗I - @test Dₓₓ == second_derivative(g_2D,stencil_set,1) - @test Dₓₓ isa LazyTensor{T,2,2} where T + Dₓₓ = second_derivative(g_2D,stencil_set,1) + @test Dₓₓ isa LazyTensor{Float64,2,2} end end @@ -45,10 +39,10 @@ maxOrder = 4; for i = 0:maxOrder-1 f_i(x) = 1/factorial(i)*x^i - monomials = (monomials...,evalOn(g_1D,f_i)) + monomials = (monomials...,eval_on(g_1D,f_i)) end - v = evalOn(g_1D,x -> sin(x)) - vₓₓ = evalOn(g_1D,x -> -sin(x)) + v = eval_on(g_1D,x -> sin(x)) + vₓₓ = eval_on(g_1D,x -> -sin(x)) # 2nd order interior stencil, 1nd order boundary stencil, # implies that L*v should be exact for monomials up to order 2. @@ -77,15 +71,15 @@ end @testset "2D" begin - l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2)); + l2(v) = sqrt(prod(spacing.(g_2D.grids))*sum(v.^2)); binomials = () maxOrder = 4; for i = 0:maxOrder-1 f_i(x,y) = 1/factorial(i)*y^i + x^i - binomials = (binomials...,evalOn(g_2D,f_i)) + binomials = (binomials...,eval_on(g_2D,f_i)) end - v = evalOn(g_2D, (x,y) -> sin(x)+cos(y)) - v_yy = evalOn(g_2D,(x,y) -> -cos(y)) + v = eval_on(g_2D, (x,y) -> sin(x)+cos(y)) + v_yy = eval_on(g_2D,(x,y) -> -cos(y)) # 2nd order interior stencil, 1st order boundary stencil, # implies that L*v should be exact for binomials up to order 2. @@ -94,7 +88,7 @@ Dyy = second_derivative(g_2D,stencil_set,2) @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 - @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 + @test Dyy*binomials[3] ≈ eval_on(g_2D,(x,y)->1.) atol = 5e-9 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2 end @@ -107,8 +101,8 @@ # due to accumulation of round-off errors/cancellation errors? @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 - @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 - @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9 + @test Dyy*binomials[3] ≈ eval_on(g_2D,(x,y)->1.) atol = 5e-9 + @test Dyy*binomials[4] ≈ eval_on(g_2D,(x,y)->y) atol = 5e-9 @test Dyy*v ≈ v_yy rtol = 5e-4 norm = l2 end end
--- a/test/SbpOperators/volumeops/inner_products/inner_product_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inner_product_test.jl Tue May 02 20:14:39 2023 +0200 @@ -10,47 +10,39 @@ Lx = π/2. Ly = Float64(π) Lz = 1. - g_1D = EquidistantGrid(77, 0.0, Lx) - g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) - g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz)) - integral(H,v) = sum(H*v) + g_1D = equidistant_grid(77, 0.0, Lx) + g_2D = equidistant_grid((77,66), (0.0, 0.0), (Lx,Ly)) + g_3D = equidistant_grid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz)) @testset "inner_product" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) @testset "0D" begin - H = inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure) - @test H == inner_product(EquidistantGrid{Float64}(), stencil_set) - @test H == IdentityTensor{Float64}() + H = inner_product(ZeroDimGrid(0.), stencil_set) @test H isa LazyTensor{T,0,0} where T end @testset "1D" begin - H = inner_product(g_1D, quadrature_interior, quadrature_closure) - @test H == inner_product(g_1D, stencil_set) - @test H isa ConstantInteriorScalingOperator + H = inner_product(g_1D, stencil_set) @test H isa LazyTensor{T,1,1} where T end @testset "2D" begin - H = inner_product(g_2D, quadrature_interior, quadrature_closure) - H_x = inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure) - H_y = inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure) - @test H == inner_product(g_2D, stencil_set) + H = inner_product(g_2D, stencil_set) + H_x = inner_product(g_2D.grids[1], stencil_set) + H_y = inner_product(g_2D.grids[2], stencil_set) @test H == H_x⊗H_y @test H isa LazyTensor{T,2,2} where T end + + # TBD: Should there be more tests? end @testset "Sizes" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) @testset "1D" begin - H = inner_product(g_1D, quadrature_interior, quadrature_closure) + H = inner_product(g_1D, stencil_set) @test domain_size(H) == size(g_1D) @test range_size(H) == size(g_1D) end @testset "2D" begin - H = inner_product(g_2D, quadrature_interior, quadrature_closure) + H = inner_product(g_2D, stencil_set) @test domain_size(H) == size(g_2D) @test range_size(H) == size(g_2D) end @@ -61,52 +53,44 @@ v = () for i = 0:4 f_i(x) = 1/factorial(i)*x^i - v = (v...,evalOn(g_1D,f_i)) + v = (v...,eval_on(g_1D,f_i)) end - u = evalOn(g_1D,x->sin(x)) + u = eval_on(g_1D,x->sin(x)) @testset "2nd order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_1D, quadrature_interior, quadrature_closure) + H = inner_product(g_1D, stencil_set) for i = 1:2 - @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 + @test sum(H*v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 end - @test integral(H,u) ≈ 1. rtol = 1e-4 + @test sum(H*u) ≈ 1. rtol = 1e-4 end @testset "4th order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_1D, quadrature_interior, quadrature_closure) + H = inner_product(g_1D, stencil_set) for i = 1:4 - @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 + @test sum(H*v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 end - @test integral(H,u) ≈ 1. rtol = 1e-8 + @test sum(H*u) ≈ 1. rtol = 1e-8 end end @testset "2D" begin b = 2.1 v = b*ones(Float64, size(g_2D)) - u = evalOn(g_2D,(x,y)->sin(x)+cos(y)) + u = eval_on(g_2D,(x,y)->sin(x)+cos(y)) @testset "2nd order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_2D, quadrature_interior, quadrature_closure) - @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 - @test integral(H,u) ≈ π rtol = 1e-4 + H = inner_product(g_2D, stencil_set) + @test sum(H*v) ≈ b*Lx*Ly rtol = 1e-13 + @test sum(H*u) ≈ π rtol = 1e-4 end @testset "4th order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_2D, quadrature_interior, quadrature_closure) - @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 - @test integral(H,u) ≈ π rtol = 1e-8 + H = inner_product(g_2D, stencil_set) + @test sum(H*v) ≈ b*Lx*Ly rtol = 1e-13 + @test sum(H*u) ≈ π rtol = 1e-8 end end end
--- a/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/inner_products/inverse_inner_product_test.jl Tue May 02 20:14:39 2023 +0200 @@ -9,29 +9,22 @@ @testset "Diagonal-stencil inverse_inner_product" begin Lx = π/2. Ly = Float64(π) - g_1D = EquidistantGrid(77, 0.0, Lx) - g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) + g_1D = equidistant_grid(77, 0.0, Lx) + g_2D = equidistant_grid((77,66), (0.0, 0.0), (Lx,Ly)) @testset "inverse_inner_product" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) @testset "0D" begin - Hi = inverse_inner_product(EquidistantGrid{Float64}(), quadrature_interior, quadrature_closure) - @test Hi == inverse_inner_product(EquidistantGrid{Float64}(), stencil_set) - @test Hi == IdentityTensor{Float64}() + Hi = inverse_inner_product(ZeroDimGrid(1.), stencil_set) @test Hi isa LazyTensor{T,0,0} where T end @testset "1D" begin - Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure) - @test Hi == inverse_inner_product(g_1D, stencil_set) - @test Hi isa ConstantInteriorScalingOperator + Hi = inverse_inner_product(g_1D, stencil_set) @test Hi isa LazyTensor{T,1,1} where T end @testset "2D" begin - Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure) - Hi_x = inverse_inner_product(restrict(g_2D,1), quadrature_interior, quadrature_closure) - Hi_y = inverse_inner_product(restrict(g_2D,2), quadrature_interior, quadrature_closure) - @test Hi == inverse_inner_product(g_2D, stencil_set) + Hi = inverse_inner_product(g_2D, stencil_set) + Hi_x = inverse_inner_product(g_2D.grids[1], stencil_set) + Hi_y = inverse_inner_product(g_2D.grids[2], stencil_set) @test Hi == Hi_x⊗Hi_y @test Hi isa LazyTensor{T,2,2} where T end @@ -39,15 +32,13 @@ @testset "Sizes" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) @testset "1D" begin - Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure) + Hi = inverse_inner_product(g_1D, stencil_set) @test domain_size(Hi) == size(g_1D) @test range_size(Hi) == size(g_1D) end @testset "2D" begin - Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure) + Hi = inverse_inner_product(g_2D, stencil_set) @test domain_size(Hi) == size(g_2D) @test range_size(Hi) == size(g_2D) end @@ -55,45 +46,37 @@ @testset "Accuracy" begin @testset "1D" begin - v = evalOn(g_1D,x->sin(x)) - u = evalOn(g_1D,x->x^3-x^2+1) + v = eval_on(g_1D,x->sin(x)) + u = eval_on(g_1D,x->x^3-x^2+1) @testset "2nd order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_1D, quadrature_interior, quadrature_closure) - Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure) + H = inner_product(g_1D, stencil_set) + Hi = inverse_inner_product(g_1D, stencil_set) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end @testset "4th order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_1D, quadrature_interior, quadrature_closure) - Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure) + H = inner_product(g_1D, stencil_set) + Hi = inverse_inner_product(g_1D, stencil_set) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end end @testset "2D" begin - v = evalOn(g_2D,(x,y)->sin(x)+cos(y)) - u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) + v = eval_on(g_2D,(x,y)->sin(x)+cos(y)) + u = eval_on(g_2D,(x,y)->x*y + x^5 - sqrt(y)) @testset "2nd order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_2D, quadrature_interior, quadrature_closure) - Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure) + H = inner_product(g_2D, stencil_set) + Hi = inverse_inner_product(g_2D, stencil_set) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end @testset "4th order" begin stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=4) - quadrature_interior = parse_scalar(stencil_set["H"]["inner"]) - quadrature_closure = parse_tuple(stencil_set["H"]["closure"]) - H = inner_product(g_2D, quadrature_interior, quadrature_closure) - Hi = inverse_inner_product(g_2D, quadrature_interior, quadrature_closure) + H = inner_product(g_2D, stencil_set) + Hi = inverse_inner_product(g_2D, stencil_set) @test Hi*H*v ≈ v rtol = 1e-15 @test Hi*H*u ≈ u rtol = 1e-15 end
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl Tue May 02 20:14:39 2023 +0200 @@ -4,40 +4,40 @@ using Sbplib.Grids using Sbplib.LazyTensors -# Default stencils (4th order) -operator_path = sbp_operators_path()*"standard_diagonal.toml" -stencil_set = read_stencil_set(operator_path; order=4) -inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) -closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) -g_1D = EquidistantGrid(101, 0.0, 1.) -g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) +@testset "Laplace" begin + # Default stencils (4th order) + operator_path = sbp_operators_path()*"standard_diagonal.toml" + stencil_set = read_stencil_set(operator_path; order=4) + g_1D = equidistant_grid(101, 0.0, 1.) + g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) -@testset "Laplace" begin @testset "Constructors" begin @testset "1D" begin - Δ = laplace(g_1D, inner_stencil, closure_stencils) - @test Laplace(g_1D, stencil_set) == Laplace(Δ, stencil_set) - @test Laplace(g_1D, stencil_set) isa LazyTensor{T,1,1} where T + @test Laplace(g_1D, stencil_set) == Laplace(laplace(g_1D, stencil_set), stencil_set) + @test Laplace(g_1D, stencil_set) isa LazyTensor{Float64,1,1} end @testset "3D" begin - Δ = laplace(g_3D, inner_stencil, closure_stencils) - @test Laplace(g_3D, stencil_set) == Laplace(Δ,stencil_set) - @test Laplace(g_3D, stencil_set) isa LazyTensor{T,3,3} where T + @test Laplace(g_3D, stencil_set) == Laplace(laplace(g_3D, stencil_set),stencil_set) + @test Laplace(g_3D, stencil_set) isa LazyTensor{Float64,3,3} end end # Exact differentiation is measured point-wise. In other cases # the error is measured in the l2-norm. @testset "Accuracy" begin - l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2)); + l2(v) = sqrt(prod(spacing.(g_3D.grids))*sum(v.^2)); polynomials = () maxOrder = 4; for i = 0:maxOrder-1 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i) - polynomials = (polynomials...,evalOn(g_3D,f_i)) + polynomials = (polynomials...,eval_on(g_3D,f_i)) end - v = evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z)) - Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z)) + # v = eval_on(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z)) + # Δv = eval_on(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z)) + + v = eval_on(g_3D, x̄ -> sin(x̄[1]) + cos(x̄[2]) + exp(x̄[3])) + Δv = eval_on(g_3D, x̄ -> -sin(x̄[1]) - cos(x̄[2]) + exp(x̄[3])) + @inferred v[1,2,3] # 2nd order interior stencil, 1st order boundary stencil, # implies that L*v should be exact for binomials up to order 2. @@ -67,19 +67,24 @@ end @testset "laplace" begin + operator_path = sbp_operators_path()*"standard_diagonal.toml" + stencil_set = read_stencil_set(operator_path; order=4) + g_1D = equidistant_grid(101, 0.0, 1.) + g_3D = equidistant_grid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) + @testset "1D" begin - Δ = laplace(g_1D, inner_stencil, closure_stencils) - @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils, 1) - @test Δ isa LazyTensor{T,1,1} where T + Δ = laplace(g_1D, stencil_set) + @test Δ == second_derivative(g_1D, stencil_set) + @test Δ isa LazyTensor{Float64,1,1} end @testset "3D" begin - Δ = laplace(g_3D, inner_stencil, closure_stencils) - @test Δ isa LazyTensor{T,3,3} where T - Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1) - Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2) - Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3) + Δ = laplace(g_3D, stencil_set) + @test Δ isa LazyTensor{Float64,3,3} + Dxx = second_derivative(g_3D, stencil_set, 1) + Dyy = second_derivative(g_3D, stencil_set, 2) + Dzz = second_derivative(g_3D, stencil_set, 3) @test Δ == Dxx + Dyy + Dzz - @test Δ isa LazyTensor{T,3,3} where T + @test Δ isa LazyTensor{Float64,3,3} end end
--- a/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/stencil_operator_distinct_closures_test.jl Tue May 02 20:14:39 2023 +0200 @@ -6,40 +6,9 @@ import Sbplib.SbpOperators.Stencil import Sbplib.SbpOperators.StencilOperatorDistinctClosures -import Sbplib.SbpOperators.stencil_operator_distinct_closures - -@testset "stencil_operator_distinct_closures" begin - lower_closure = ( - Stencil(-1,1, center=1), - ) - - inner_stencil = Stencil(-2,2, center=1) - - upper_closure = ( - Stencil(-3,3, center=1), - Stencil(-4,4, center=2), - ) - - g₁ = EquidistantGrid(5, 0., 1.) - g₂ = EquidistantGrid((5,5), (0.,0.), (1.,1.)) - h = 1/4 - - A₁ = stencil_operator_distinct_closures(g₁, inner_stencil, lower_closure, upper_closure, 1) - A₂¹ = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 1) - A₂² = stencil_operator_distinct_closures(g₂, inner_stencil, lower_closure, upper_closure, 2) - - v₁ = evalOn(g₁, x->x) - - u = [1., 2., 2., 3., 4.]*h - @test A₁*v₁ == u - - v₂ = evalOn(g₂, (x,y)-> x + 3y) - @test A₂¹*v₂ == repeat(u, 1, 5) - @test A₂²*v₂ == repeat(3u', 5, 1) -end @testset "StencilOperatorDistinctClosures" begin - g = EquidistantGrid(11, 0., 1.) + g = equidistant_grid(11, 0., 1.) lower_closure = ( Stencil(-1,1, center=1),
--- a/test/SbpOperators/volumeops/volume_operator_test.jl Mon May 01 11:37:09 2023 +0200 +++ b/test/SbpOperators/volumeops/volume_operator_test.jl Tue May 02 20:14:39 2023 +0200 @@ -14,7 +14,7 @@ @testset "VolumeOperator" begin inner_stencil = CenteredStencil(1/4, 2/4, 1/4) closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2)) - g = EquidistantGrid(11,0.,1.) + g = equidistant_grid(11,0.,1.) @testset "Constructors" begin op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)