Mercurial > repos > public > sbplib_julia

changeset 1651:707fc9761c2b feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2024 12:47:26 +0200
parents b7dcd3dd3181 (current diff) 8250cf5a3ce9 (diff)
children 65b2d2c72fbc
files Manifest.toml Project.toml src/Grids/Grids.jl src/SbpOperators/volumeops/laplace/laplace.jl test/SbpOperators/volumeops/laplace/laplace_test.jl
diffstat 23 files changed, 589 insertions(+), 128 deletions(-) [+]
line wrap: on
line diff
diff -r b7dcd3dd3181 -r 707fc9761c2b Manifest.toml
--- a/Manifest.toml	Wed Jun 26 12:36:41 2024 +0200
+++ b/Manifest.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.2"
+julia_version = "1.10.4"
 manifest_format = "2.0"
 project_hash = "32fac879810099260f177c27318d3f26de4a00cc"
 
@@ -16,14 +16,15 @@
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d"
+git-tree-sha1 = "ed2ec3c9b483842ae59cd273834e5b46206d6dda"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.9.0"
+version = "7.11.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceCUDSSExt = "CUDSS"
     ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
     ArrayInterfaceReverseDiffExt = "ReverseDiff"
@@ -34,6 +35,7 @@
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
     ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
@@ -45,9 +47,9 @@
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80"
+git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.14.0"
+version = "4.15.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -56,7 +58,7 @@
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.0+0"
+version = "1.1.1+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -75,9 +77,9 @@
 uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621"
+git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.13.0"
+version = "1.14.0"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -145,9 +147,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2"
+git-tree-sha1 = "6e00379a24597be4ae1ee6b2d882e15392040132"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.3"
+version = "1.9.5"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -158,9 +160,9 @@
     Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
diff -r b7dcd3dd3181 -r 707fc9761c2b Notes.md
--- a/Notes.md	Wed Jun 26 12:36:41 2024 +0200
+++ b/Notes.md	Wed Jun 26 12:47:26 2024 +0200
@@ -1,5 +1,41 @@
 # Notes
 
+## How to dispatch for different operators
+We have a problem in how dispatch for different operators work.
+ * We want to keep the types simple and flat (Awkward to forward `apply`)
+ * We want to dispatch SATs on the parameters of the continuous operator. (a * div for example)
+ * We want to allow keeping the same stencil_set across different calls. (maybe not so bad for the user to be responsible)
+
+Could remove the current opset idea and introduce a description of continuous operators
+ ```julia
+abstract type DifferentialOperator end
+
+struct Laplace <: DifferentialOperator end
+struct Advection <: DifferentialOperator
+    v
+end
+
+difference_operator(::Laplace, grid, stencil_set) = ... # Returns a plain LazyTensor. Replaces the current `laplace()` function.
+sat_tensors(::Laplace, grid, stencil_set, bc) = ...
+
+sat(::DifferentialOperator, grid, stencil_set, bc) = ...
+ ```
+
+
+### Update 2024-06-26
+We will run into trouble if we start assuming things about the coupling
+between the continuous and discrete setting. We could add representations of
+continuous operators but we will also need representations of discrete
+operators. Ideally it should be possible to ignore the continuous
+representations and only work with the discrete operators without losing
+functionality. The discrete representations does not have to be LazyTensors.
+The could be used as inputs to methods for `sat`, `difference_operator` and so
+on.
+
+To see need for a fully functional discrete layer we can consider the
+optimization of material parameters or something similar. In this case we do
+not necessarily want to handle continuous objects.
+
 ## Reading operators
 
 Jonatan's suggestion is to add methods to `Laplace`, `SecondDerivative` and
diff -r b7dcd3dd3181 -r 707fc9761c2b Project.toml
--- a/Project.toml	Wed Jun 26 12:36:41 2024 +0200
+++ b/Project.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -1,7 +1,7 @@
 name = "Sbplib"
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
 authors = ["Jonatan Werpers <jonatan@werpers.com>", "Vidar Stiernström <vidar.stiernstrom@it.uu.se>, and contributors"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff -r b7dcd3dd3181 -r 707fc9761c2b TODO.md
--- a/TODO.md	Wed Jun 26 12:36:41 2024 +0200
+++ b/TODO.md	Wed Jun 26 12:47:26 2024 +0200
@@ -21,6 +21,10 @@
        See (https://docs.julialang.org/en/v1/manual/types/#man-custom-pretty-printing)
  - [ ] Samla noggrannhets- och SBP-ness-tester för alla operatorer på ett ställe
  - [ ] Move export statements to top of each module
+ - [ ] Implement apply_transpose for
+      - [ ] ElementwiseTensorOperation
+      - [ ] VolumeOperator
+      - [ ] Laplace
 
 
  - [ ] Gå igenom alla typ parametrar och kolla om de är motiverade. Både i signaturer och typer, tex D i VariableSecondDerivative. Kan vi använda promote istället?
diff -r b7dcd3dd3181 -r 707fc9761c2b benchmark/Manifest.toml
--- a/benchmark/Manifest.toml	Wed Jun 26 12:36:41 2024 +0200
+++ b/benchmark/Manifest.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.2"
+julia_version = "1.10.4"
 manifest_format = "2.0"
 project_hash = "25bba7b4a00465d5a2b00b589eb10e3301c31f2a"
 
@@ -25,14 +25,15 @@
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "44691067188f6bd1b2289552a23e4b7572f4528d"
+git-tree-sha1 = "ed2ec3c9b483842ae59cd273834e5b46206d6dda"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.9.0"
+version = "7.11.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceCUDSSExt = "CUDSS"
     ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
     ArrayInterfaceReverseDiffExt = "ReverseDiff"
@@ -43,6 +44,7 @@
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
     ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
@@ -63,9 +65,9 @@
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "c955881e3c981181362ae4088b35995446298b80"
+git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.14.0"
+version = "4.15.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -74,7 +76,7 @@
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.1.0+0"
+version = "1.1.1+0"
 
 [[deps.DataAPI]]
 git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
@@ -185,9 +187,9 @@
 version = "1.2.0"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621"
+git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.13.0"
+version = "1.14.0"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -268,7 +270,7 @@
 deps = ["StaticArrays", "TOML", "TiledIteration"]
 path = ".."
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-version = "0.1.0"
+version = "0.1.1"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -300,9 +302,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "bf074c045d3d5ffd956fa0a461da38a44685d6b2"
+git-tree-sha1 = "6e00379a24597be4ae1ee6b2d882e15392040132"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.3"
+version = "1.9.5"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -313,9 +315,9 @@
     Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
diff -r b7dcd3dd3181 -r 707fc9761c2b benchmark/benchmarks.jl
--- a/benchmark/benchmarks.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/benchmark/benchmarks.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -15,9 +15,9 @@
 ll(d) = ntuple(i->0., d)
 lu(d) = ntuple(i->1., d)
 
-g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)[1])
-g2 = equidistant_grid(ll(2), lu(2), sz(2))
-g3 = equidistant_grid(ll(3), lu(3), sz(3))
+g1 = equidistant_grid(ll(1)[1], lu(1)[1], sz(1)...)
+g2 = equidistant_grid(ll(2), lu(2), sz(2)...)
+g3 = equidistant_grid(ll(3), lu(3), sz(3)...)
 
 v1 = rand(sz(1)...)
 v2 = rand(sz(2)...)
diff -r b7dcd3dd3181 -r 707fc9761c2b docs/Manifest.toml
--- a/docs/Manifest.toml	Wed Jun 26 12:36:41 2024 +0200
+++ b/docs/Manifest.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.4"
 manifest_format = "2.0"
 project_hash = "4f0756199bb5f6739a5f4697152617efc4e0705c"
 
@@ -10,15 +10,15 @@
 version = "0.0.1"
 
 [[deps.AbstractTrees]]
-git-tree-sha1 = "faa260e4cb5aba097a73fab382dd4b5819d8ec8c"
+git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
-version = "0.4.4"
+version = "0.4.5"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "0fb305e0253fd4e833d486914367a2ee2c2e78d0"
+git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.1"
+version = "4.0.4"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -29,16 +29,19 @@
 version = "1.1.1"
 
 [[deps.ArrayInterface]]
-deps = ["Adapt", "LinearAlgebra", "Requires", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "bbec08a37f8722786d87bedf84eae19c020c4efa"
+deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "ed2ec3c9b483842ae59cd273834e5b46206d6dda"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.7.0"
+version = "7.11.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceBandedMatricesExt = "BandedMatrices"
     ArrayInterfaceBlockBandedMatricesExt = "BlockBandedMatrices"
     ArrayInterfaceCUDAExt = "CUDA"
+    ArrayInterfaceCUDSSExt = "CUDSS"
+    ArrayInterfaceChainRulesExt = "ChainRules"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    ArrayInterfaceReverseDiffExt = "ReverseDiff"
     ArrayInterfaceStaticArraysCoreExt = "StaticArraysCore"
     ArrayInterfaceTrackerExt = "Tracker"
 
@@ -46,7 +49,10 @@
     BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
     BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+    CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
+    ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
@@ -56,11 +62,17 @@
 [[deps.Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[deps.CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.4"
+
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
-git-tree-sha1 = "75bd5b6fc5089df449b5d35fa501c846c9b6549b"
+git-tree-sha1 = "b1c55339b7c6c350ee89f2c1604299660525b248"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "4.12.0"
+version = "4.15.0"
 weakdeps = ["Dates", "LinearAlgebra"]
 
     [deps.Compat.extensions]
@@ -69,7 +81,7 @@
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.1+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -82,10 +94,10 @@
 version = "0.9.3"
 
 [[deps.Documenter]]
-deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "Test", "Unicode"]
-git-tree-sha1 = "2613dbec8f4748273bbe30ba71fd5cb369966bac"
+deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"]
+git-tree-sha1 = "5461b2a67beb9089980e2f8f25145186b6d34f91"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "1.2.1"
+version = "1.4.1"
 
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
@@ -94,30 +106,30 @@
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "4558ab818dcceaab612d1bb8c19cee87eda2b83c"
+git-tree-sha1 = "1c6317308b9dc757616f0b5cb379db10494443a7"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.5.0+0"
+version = "2.6.2+0"
 
 [[deps.FileWatching]]
 uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
 
 [[deps.Git]]
 deps = ["Git_jll"]
-git-tree-sha1 = "51764e6c2e84c37055e846c516e9015b4a291c7d"
+git-tree-sha1 = "04eff47b1354d702c3a85e8ab23d539bb7d5957e"
 uuid = "d7ba0133-e1db-5d97-8f8c-041e4b3a1eb2"
-version = "1.3.0"
+version = "1.3.1"
 
 [[deps.Git_jll]]
 deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "b30c473c97fcc1e1e44fab8f3e88fd1b89c9e9d1"
+git-tree-sha1 = "d18fb8a1f3609361ebda9bf029b60fd0f120c809"
 uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
-version = "2.43.0+0"
+version = "2.44.0+2"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
-git-tree-sha1 = "8b72179abc660bfab5e28472e019392b97d0985c"
+git-tree-sha1 = "b6d6bfdd7ce25b0f9b2f6b3dd56b2673a66c8770"
 uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
-version = "0.2.4"
+version = "0.2.5"
 
 [[deps.IfElse]]
 git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
@@ -212,9 +224,9 @@
 version = "1.2.0"
 
 [[deps.OffsetArrays]]
-git-tree-sha1 = "6a731f2b5c03157418a20c12195eb4b74c8f8621"
+git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e"
 uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-version = "1.13.0"
+version = "1.14.0"
 weakdeps = ["Adapt"]
 
     [deps.OffsetArrays.extensions]
@@ -223,13 +235,13 @@
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "cc6e1927ac521b659af340e0ca45828a3ffc748f"
+git-tree-sha1 = "a028ee3cb5641cccc4c24e90c36b0a4f7707bdf5"
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.0.12+0"
+version = "3.0.14+0"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -249,15 +261,15 @@
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.1"
+version = "1.4.3"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -291,7 +303,7 @@
 deps = ["StaticArrays", "TOML", "TiledIteration"]
 path = ".."
 uuid = "5a373a26-915f-4769-bcab-bf03835de17b"
-version = "0.1.0"
+version = "0.1.1"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -306,9 +318,9 @@
 
 [[deps.Static]]
 deps = ["IfElse"]
-git-tree-sha1 = "f295e0a1da4ca425659c57441bcb59abb035a4bc"
+git-tree-sha1 = "d2fdac9ff3906e27f7a618d47b676941baa6c80c"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "0.8.8"
+version = "0.8.10"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "Requires", "SparseArrays", "Static", "SuiteSparse"]
@@ -323,9 +335,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5"
+git-tree-sha1 = "6e00379a24597be4ae1ee6b2d882e15392040132"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.1"
+version = "1.9.5"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -336,9 +348,9 @@
     Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.SuiteSparse]]
 deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
@@ -369,6 +381,15 @@
 uuid = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 version = "0.5.0"
 
+[[deps.TranscodingStreams]]
+git-tree-sha1 = "d73336d81cafdc277ff45558bb7eaa2b04a8e472"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.10.10"
+weakdeps = ["Random", "Test"]
+
+    [deps.TranscodingStreams.extensions]
+    TestExt = ["Test", "Random"]
+
 [[deps.UUIDs]]
 deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
diff -r b7dcd3dd3181 -r 707fc9761c2b ext/SbplibMakieExt.jl
--- a/ext/SbplibMakieExt.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/ext/SbplibMakieExt.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -2,6 +2,7 @@
 
 using Sbplib.Grids
 using Makie
+using StaticArrays
 
 
 function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
@@ -33,39 +34,58 @@
 end
 
 
-function make_plot(g,gf)
-    v,f,c = verticies_and_faces_and_values(g,gf)
-    mesh(v,f,color=c,
+## Grids
+
+Makie.convert_arguments(::Type{<:Scatter}, g::Grid) = (reshape(map(Point,g),:),) # (map(Point,collect(g)[:]),)
+function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,2})
+    M = collect(g)
+
+    function cat_with_NaN(a,b)
+        vcat(a,[@SVector[NaN,NaN]],b)
+    end
+
+    xlines = reduce(cat_with_NaN, eachrow(M))
+    ylines = reduce(cat_with_NaN, eachcol(M))
+
+    return (cat_with_NaN(xlines,ylines),)
+end
+
+Makie.plot!(plot::Plot(Grid{<:Any,2})) = lines!(plot, plot.attributes, plot[1])
+
+
+## Grid functions
+
+### 1D
+function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
+    (collect(g), gf)
+end
+
+function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
+    (collect(g), gf)
+end
+
+Makie.plot!(plot::Plot(Grid{<:Any,1}, AbstractArray{<:Any,1})) = lines!(plot, plot.attributes, plot[1], plot[2])
+
+### 2D
+function Makie.convert_arguments(::Type{<:Surface}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    (getindex.(g,1), getindex.(g,2), gf)
+end
+
+function Makie.plot!(plot::Plot(Grid{<:Any,2},AbstractArray{<:Any, 2}))
+    r = @lift verticies_and_faces_and_values($(plot[1]), $(plot[2]))
+    v,f,c = (@lift $r[1]), (@lift $r[2]), (@lift $r[3])
+    mesh!(plot, plot.attributes, v, f;
+        color=c,
         shading = NoShading,
     )
 end
-
-# scatter(collect(g)[:])
+# TBD: Can we define `mesh` instead of the above function and then forward plot! to that?
 
-function Makie.surface(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2}; kwargs...)
-    surface(getindex.(g,1), getindex.(g,2), gf;
-        shading = NoShading,
-        kwargs...,
-    )
+function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    ps = map(g,gf) do (x,y), z
+        @SVector[x,y,z]
+    end
+    (reshape(ps,:),)
 end
 
-function Makie.mesh(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2}; kwargs...)
-    v,f,c = verticies_and_faces_and_values(g, gf)
-    mesh(v,f,color=c,
-        shading = NoShading,
-        kwargs...,
-    )
 end
-
-function Makie.plot!(plot::Plot(Grid{<:Any,2},AbstractArray{<:Any, 2}))
-    # TODO: How to handle kwargs?
-    # v,f,c = verticies_and_faces_and_values(plot[1], plot[2])
-    r = @lift verticies_and_faces_and_values($(plot[1]), $(plot[2]))
-    v,f,c = (@lift $r[1]), (@lift $r[2]), (@lift $r[3])
-    mesh!(plot, v, f, color=c,
-        shading = NoShading,
-    )
-end
-
-Makie.convert_arguments(::Type{<:Scatter}, g::Grid) = (map(Tuple,collect(g)[:]),)
-end
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Grids/Grids.jl
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Grids/equidistant_grid.jl
--- a/src/Grids/equidistant_grid.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -107,6 +107,9 @@
 behaviors.
 """
 function equidistant_grid(limit_lower, limit_upper, dims::Vararg{Int})
+    if !(length(limit_lower) == length(limit_upper) == length(dims))
+        throw(ArgumentError("All arguments must be of the same length"))
+    end
     gs = map(equidistant_grid, limit_lower, limit_upper, dims)
     return TensorGrid(gs...)
 end
@@ -129,6 +132,7 @@
 
 
 equidistant_grid(hb::HyperBox, dims::Vararg{Int}) = equidistant_grid(hb.a, hb.b, dims...)
+# TODO: One dimensional grids shouldn't have vector eltype right?, Change here or in HyperBox?
 
 function equidistant_grid(c::Chart, dims::Vararg{Int})
     lg = equidistant_grid(parameterspace(c), dims...)
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Grids/grid.jl
--- a/src/Grids/grid.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/Grids/grid.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -74,7 +74,7 @@
     end
 end
 
-Base.size(cv) = size(cv.v)
+Base.size(cv::ArrayComponentView) = size(cv.v)
 Base.getindex(cv::ArrayComponentView, i::Int) = cv.v[i][cv.component_index...]
 Base.getindex(cv::ArrayComponentView, I::Vararg{Int}) = cv.v[I...][cv.component_index...]
 IndexStyle(::Type{<:ArrayComponentView{<:Any,<:Any,AT}}) where AT = IndexStyle(AT)
@@ -107,6 +107,8 @@
 """
 function boundary_identifiers end
 
+# TBD: Boundary identifiers for charts and atlases?
+
 """
     boundary_grid(g::Grid, id::BoundaryIdentifier)
 
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Grids/manifolds.jl
--- a/src/Grids/manifolds.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/Grids/manifolds.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -19,6 +19,7 @@
 """
 abstract type ParameterSpace{D} end
 Base.ndims(::ParameterSpace{D}) where D = D
+# TBD:  Should implement domain_dim?
 
 struct HyperBox{T,D} <: ParameterSpace{D}
     a::SVector{D,T}
@@ -98,8 +99,11 @@
 which will both allow calling `jacobian(c,ξ)`.
 """
 jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ)
+# TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian?
 
 
+# TBD: Should Charts, parameterspaces have boundary names?
+
 """
     Atlas
 
@@ -160,6 +164,16 @@
 (c::LineSegment)(s) = (1-s)*c.a + s*c.b
 
 
+function linesegments(ps...)
+    return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1]
+end
+
+
+function polygon_edges(ps...)
+    n = length(ps)
+    return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(Ps)]
+end
+
 struct Circle{T,PT} <: Curve
     c::PT
     r::T
@@ -189,8 +203,5 @@
     s(ξ̄...)
 end
 
+# TODO: Implement jacobian() for the different mapping helpers
 
-function polygon_sides(Ps...)
-    n = length(Ps)
-    return [t->line(t,Ps[i],Ps[mod1(i+1,n)]) for i ∈ eachindex(Ps)]
-end
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Grids/mapped_grid.jl
--- a/src/Grids/mapped_grid.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -14,6 +14,7 @@
 
 Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d)
 Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d)
+# TODO: axes(...)?
 
 # Iteration interface
 
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Grids/tensor_grid.jl
diff -r b7dcd3dd3181 -r 707fc9761c2b src/LazyTensors/lazy_tensor_operations.jl
--- a/src/LazyTensors/lazy_tensor_operations.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/LazyTensors/lazy_tensor_operations.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -121,6 +121,7 @@
 
 Base.:*(a::T, tm::LazyTensor{T}) where T = TensorComposition(ScalingTensor{T,range_dim(tm)}(a,range_size(tm)), tm)
 Base.:*(tm::LazyTensor{T}, a::T) where T = a*tm
+Base.:-(tm::LazyTensor) = (-one(eltype(tm)))*tm
 
 """
     InflatedTensor{T,R,D} <: LazyTensor{T,R,D}
diff -r b7dcd3dd3181 -r 707fc9761c2b src/SbpOperators/SbpOperators.jl
--- a/src/SbpOperators/SbpOperators.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/SbpOperators/SbpOperators.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -11,7 +11,6 @@
 export sbp_operators_path
 
 # Operators
-export boundary_quadrature
 export boundary_restriction
 export inner_product
 export inverse_inner_product
@@ -22,20 +21,34 @@
 export second_derivative
 export second_derivative_variable
 export undivided_skewed04
-
-using Sbplib.RegionIndices
-using Sbplib.LazyTensors
-using Sbplib.Grids
+export closure_size
 
 @enum Parity begin
     odd = -1
     even = 1
 end
 
-export closure_size
 
+# Boundary conditions
+export BoundaryCondition
+export NeumannCondition
+export DirichletCondition
+export discretize_data
+export boundary_data
+export boundary
+export sat
+export sat_tensors
+
+# Using
+using Sbplib.RegionIndices
+using Sbplib.LazyTensors
+using Sbplib.Grids
+
+# Includes
 include("stencil.jl")
 include("stencil_set.jl")
+include("boundary_conditions/boundary_condition.jl")
+include("boundary_conditions/sat.jl")
 include("volumeops/volume_operator.jl")
 include("volumeops/stencil_operator_distinct_closures.jl")
 include("volumeops/constant_interior_scaling_operator.jl")
diff -r b7dcd3dd3181 -r 707fc9761c2b src/SbpOperators/boundary_conditions/boundary_condition.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundary_conditions/boundary_condition.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -0,0 +1,62 @@
+"""
+    BoundaryCondition
+
+Description of a boundary condition. Implementations describe the kind of
+boundary condition, what boundary the condition applies to, and any associated
+data. Should implement [`boundary`](@ref) and may implement
+[`boundary_data`](@ref) if applicable.
+
+For examples see [`DirichletCondition`](@ref) and [`NeumannCondition`](@ref)
+"""
+abstract type BoundaryCondition end
+
+"""
+    boundary(::BoundaryCondition)
+
+The boundary identifier of the BoundaryCondition.
+"""
+function boundary end
+
+"""
+    boundary_data(::BoundaryCondition)
+
+If implemented, the data associated with the BoundaryCondition.
+"""
+function boundary_data end
+
+"""
+    discretize_data(grid, bc::BoundaryCondition)
+
+The data of `bc` as a lazily evaluated grid function on the boundary grid
+specified by `boundary(bc)`.
+"""
+function discretize_data(grid, bc::BoundaryCondition)
+    return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc))
+end
+
+"""
+    DirichletCondition{DT,BID}
+
+A Dirichlet condition with `data::DT` on the boundary
+specified by the boundary identifier `BID`.
+"""
+struct DirichletCondition{DT,BID} <: BoundaryCondition
+    data::DT
+    boundary::BID
+end
+boundary_data(bc::DirichletCondition) = bc.data
+boundary(bc::DirichletCondition) = bc.boundary
+
+"""
+    NeumannCondition{DT,BID}
+
+A Neumann condition with `data::DT` on the boundary
+specified by the boundary identifier `BID`.
+"""
+struct NeumannCondition{DT,BID} <: BoundaryCondition
+    data::DT
+    boundary::BID
+end
+boundary_data(bc::NeumannCondition) = bc.data
+boundary(bc::NeumannCondition) = bc.boundary
+
diff -r b7dcd3dd3181 -r 707fc9761c2b src/SbpOperators/boundary_conditions/sat.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/SbpOperators/boundary_conditions/sat.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -0,0 +1,25 @@
+"""
+    sat_tensors(op, grid, bc::BoundaryCondition; kwargs...)
+
+The penalty tensor and boundary operator used to construct a
+simultaneous-approximation-term for imposing `bc` related to `op`.
+
+For `penalty_tensor, L  = sat_tensors(...)` then `SAT(u,g) =
+penalty_tensor*(L*u - g)`  where `g` is the boundary data.
+"""
+function sat_tensors end
+
+
+"""
+    sat(op, grid, bc::BoundaryCondition; kwargs...)
+
+Simultaneous-Approximation-Term for a general `bc` to `op`. Returns a function
+`SAT(u,g)` weakly imposing `bc` when added to `op*u`.
+
+Internally `sat_tensors(op, grid, bc; ...)` is called to construct the
+necessary parts for the SAT.
+"""
+function sat(op, grid, bc::BoundaryCondition; kwargs...)
+    penalty_tensor, L = sat_tensors(op, grid, bc; kwargs...)
+    return SAT(u, g) = penalty_tensor*(L*u - g)
+end
diff -r b7dcd3dd3181 -r 707fc9761c2b src/SbpOperators/operators/standard_diagonal.toml
--- a/src/SbpOperators/operators/standard_diagonal.toml	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -34,11 +34,14 @@
     {s = ["1", "-2", "1"], c = 1},
 ]
 
+D2.positivity = {theta_M = "0.3636363636", theta_R = "1.000000538455350", m_b = "2"}
+
 D2variable.inner_stencil = [["1/2", "1/2", "0"],[ "-1/2", "-1", "-1/2"],["0", "1/2", "1/2"]]
 D2variable.closure_stencils = [
         {s = [["2", "-1", "0"],["-3", "1",   "0"],["1","0","0"]], c = 1},
 ]
 
+
 [[stencil_set]]
 
 order = 4
@@ -65,6 +68,8 @@
     {s = [ "-1/49",     "0",   "59/49", "-118/49", "64/49", "-4/49"], c = 4},
 ]
 
+D2.positivity = {theta_M = "0.2505765857", theta_R = "0.577587500088313", m_b = "4"}
+
 D2variable.inner_stencil = [
     ["-1/8",   "1/6", "-1/8",   "0",    "0"  ],
     [ "1/6",   "1/2",  "1/2",  "1/6",   "0"  ],
@@ -136,7 +141,6 @@
 ]
 
 
-
 [[stencil_set]]
 
 order = 6
@@ -150,3 +154,5 @@
 
 e.closure = ["1"]
 d1.closure = ["-25/12", "4", "-3", "4/3", "-1/4"]
+
+D2.positivity = {theta_M = "0.1878687080", theta_R = "0.3697", m_b = "7"}
diff -r b7dcd3dd3181 -r 707fc9761c2b src/SbpOperators/stencil_set.jl
diff -r b7dcd3dd3181 -r 707fc9761c2b src/SbpOperators/volumeops/laplace/laplace.jl
--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -51,9 +51,9 @@
     end
     return Δ
 end
+
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
-
 function laplace(grid::MappedGrid, stencil_set)
     J = jacobian_determinant(grid)
     J⁻¹ = DiagonalTensor(map(inv, J))
@@ -74,3 +74,79 @@
         end
     end
 end
+
+
+"""
+    sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+
+The operators required to construct the SAT for imposing a Dirichlet
+condition. `H_tuning` and `R_tuning` are used to specify the strength of the
+penalty.
+
+See also: [`sat`](@ref),[`DirichletCondition`](@ref), [`positivity_decomposition`](@ref).
+"""
+function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
+    id = boundary(bc)
+    set  = Δ.stencil_set
+    H⁻¹ = inverse_inner_product(g,set)
+    Hᵧ = inner_product(boundary_grid(g, id), set)
+    e = boundary_restriction(g, set, id)
+    d = normal_derivative(g, set, id)
+    B = positivity_decomposition(Δ, g, bc; H_tuning, R_tuning)
+    penalty_tensor = H⁻¹∘(d' - B*e')∘Hᵧ
+    return penalty_tensor, e
+end
+
+"""
+    sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+
+The operators required to construct the SAT for imposing a Neumann condition.
+
+See also: [`sat`](@ref), [`NeumannCondition`](@ref).
+"""
+function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+    id = boundary(bc)
+    set  = Δ.stencil_set
+    H⁻¹ = inverse_inner_product(g,set)
+    Hᵧ = inner_product(boundary_grid(g, id), set)
+    e = boundary_restriction(g, set, id)
+    d = normal_derivative(g, set, id)
+
+    penalty_tensor = -H⁻¹∘e'∘Hᵧ
+    return penalty_tensor, d
+end
+
+"""
+    positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+
+Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive
+definite with respect to the boundary quadrature. Here `d` is the normal
+derivative and `e` is the boundary restriction operator. `B` can then be used
+to form a symmetric and energy stable penalty for a Dirichlet condition. The
+parameters `H_tuning` and `R_tuning` are used to specify the strength of the
+penalty and must be greater than 1. For details we refer to
+https://doi.org/10.1016/j.jcp.2020.109294
+"""
+function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+    @assert(H_tuning ≥ 1.)
+    @assert(R_tuning ≥ 1.)
+    Nτ_H, τ_R = positivity_limits(Δ,g,bc)
+    return H_tuning*Nτ_H + R_tuning*τ_R
+end
+
+# TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
+# change bc::BoundaryCondition to id::BoundaryIdentifier
+function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)
+    h = spacing(g)
+    θ_H = parse_scalar(Δ.stencil_set["H"]["closure"][1])
+    θ_R = parse_scalar(Δ.stencil_set["D2"]["positivity"]["theta_R"])
+
+    τ_H = 1/(h*θ_H)
+    τ_R = 1/(h*θ_R)
+    return τ_H, τ_R
+end
+
+function positivity_limits(Δ::Laplace, g::TensorGrid, bc::DirichletCondition)
+    τ_H, τ_R = positivity_limits(Δ, g.grids[grid_id(boundary(bc))], bc)
+    return τ_H*ndims(g), τ_R
+end
diff -r b7dcd3dd3181 -r 707fc9761c2b src/Sbplib.jl
diff -r b7dcd3dd3181 -r 707fc9761c2b test/Grids/equidistant_grid_test.jl
--- a/test/Grids/equidistant_grid_test.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -118,6 +118,8 @@
     @test_throws DomainError equidistant_grid((1.0,1.0),(1.0,1.0), 1, 1)
     @test_throws DomainError equidistant_grid((1.0,1.0),(-1.0,-1.0), 1, 1)
 
+    @test_throws ArgumentError equidistant_grid((0.0,),(8.0,5.0), 4, 3, 4)
+
     @testset "Base" begin
         @test eltype(equidistant_grid(0.0, 1.0, 4)) == Float64
         @test eltype(equidistant_grid((0,0),(1,3), 4, 3)) <: AbstractVector{Float64}
diff -r b7dcd3dd3181 -r 707fc9761c2b test/Grids/tensor_grid_test.jl
diff -r b7dcd3dd3181 -r 707fc9761c2b test/Manifest.toml
--- a/test/Manifest.toml	Wed Jun 26 12:36:41 2024 +0200
+++ b/test/Manifest.toml	Wed Jun 26 12:47:26 2024 +0200
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.10.0"
+julia_version = "1.10.4"
 manifest_format = "2.0"
 project_hash = "f2b0634c12bbed93a17efc88d466604d5a07c465"
 
@@ -16,14 +16,14 @@
 
 [[deps.BenchmarkTools]]
 deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "f1f03a9fa24271160ed7e73051fba3c1a759b53f"
+git-tree-sha1 = "f1dff6729bc61f4d49e140da1af55dcd1ac97b2f"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.4.0"
+version = "1.5.0"
 
 [[deps.CompilerSupportLibraries_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "1.0.5+1"
+version = "1.1.1+0"
 
 [[deps.Dates]]
 deps = ["Printf"]
@@ -117,9 +117,9 @@
 
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
-git-tree-sha1 = "7d6dd4e9212aebaeed356de34ccf262a3cd415aa"
+git-tree-sha1 = "a2d09619db4e765091ee5c6ffe8872849de0feea"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
-version = "0.3.26"
+version = "0.3.28"
 
     [deps.LogExpFunctions.extensions]
     LogExpFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -163,7 +163,7 @@
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
-version = "0.3.23+2"
+version = "0.3.23+4"
 
 [[deps.OpenLibm_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -189,15 +189,15 @@
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "03b4c25b43cb84cee5c90aa9b5ea0a78fd848d2f"
+git-tree-sha1 = "5aa36f7049a63a1528fe8f7c3f2113413ffd4e1f"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.2.0"
+version = "1.2.1"
 
 [[deps.Preferences]]
 deps = ["TOML"]
-git-tree-sha1 = "00805cd429dcb4870060ff49ef443486c262e38e"
+git-tree-sha1 = "9306f6085165d270f7e3db02af26a400d580f5c6"
 uuid = "21216c6a-2e73-6563-6e65-726566657250"
-version = "1.4.1"
+version = "1.4.3"
 
 [[deps.Printf]]
 deps = ["Unicode"]
@@ -238,9 +238,9 @@
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "e2cfc4012a19088254b3950b85c3c1d8882d864d"
+git-tree-sha1 = "2f5d4697f21388cbe1ff299430dd169ef97d7e14"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.3.1"
+version = "2.4.0"
 
     [deps.SpecialFunctions.extensions]
     SpecialFunctionsChainRulesCoreExt = "ChainRulesCore"
@@ -250,9 +250,9 @@
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "f68dd04d131d9a8a8eb836173ee8f105c360b0c5"
+git-tree-sha1 = "6e00379a24597be4ae1ee6b2d882e15392040132"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.1"
+version = "1.9.5"
 
     [deps.StaticArrays.extensions]
     StaticArraysChainRulesCoreExt = "ChainRulesCore"
@@ -263,9 +263,9 @@
     Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [[deps.StaticArraysCore]]
-git-tree-sha1 = "36b3d696ce6366023a0ea192b4cd442268995a0d"
+git-tree-sha1 = "192954ef1208c7019899fbf8049e717f92959682"
 uuid = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
-version = "1.4.2"
+version = "1.4.3"
 
 [[deps.Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -293,9 +293,9 @@
 
 [[deps.TestSetExtensions]]
 deps = ["DeepDiffs", "Distributed", "Test"]
-git-tree-sha1 = "3a2919a78b04c29a1a57b05e1618e473162b15d0"
+git-tree-sha1 = "ccebd99935be339d2ad907589708ba1f0d62bab3"
 uuid = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
-version = "2.0.0"
+version = "3.0.0"
 
 [[deps.Tullio]]
 deps = ["DiffRules", "LinearAlgebra", "Requires"]
diff -r b7dcd3dd3181 -r 707fc9761c2b test/SbpOperators/boundary_conditions/boundary_condition_test.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/boundary_conditions/boundary_condition_test.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -0,0 +1,39 @@
+using Test
+
+using Sbplib.Grids
+using Sbplib.RegionIndices
+using Sbplib.SbpOperators
+
+@testset "BoundaryCondition" begin
+    grid_1d = equidistant_grid(0.0, 1.0, 11)
+    grid_2d = equidistant_grid((0.0, 0.0), (1.0,1.0), 11, 15)
+    grid_3d = equidistant_grid((0.0, 0.0, 0.0), (1.0,1.0, 1.0), 11, 15, 13)
+    (id_l,_) = boundary_identifiers(grid_1d)
+    (_,_,_,id_n) = boundary_identifiers(grid_2d)
+    (_,_,_,_,id_b,_) = boundary_identifiers(grid_3d)
+
+    g = 3.14
+    f(x,y,z) = x^2+y^2+z^2
+    @testset "Constructors" begin
+        @test DirichletCondition(g,id_l) isa DirichletCondition{Float64,Lower}
+        @test NeumannCondition(f,id_b) isa NeumannCondition{<:Function,CartesianBoundary{3,Lower}}
+    end
+
+    @testset "boundary" begin
+        @test boundary(DirichletCondition(g,id_l)) == id_l
+        @test boundary(NeumannCondition(f,id_b)) == id_b
+    end
+
+    @testset "boundary_data" begin
+        @test boundary_data(DirichletCondition(g,id_l)) == g
+        @test boundary_data(NeumannCondition(f,id_b)) == f
+    end
+
+    @testset "discretize_data" begin
+        @test fill(g) ≈ discretize_data(grid_1d,DirichletCondition(g,id_l))
+        @test g*ones(11,1) ≈ discretize_data(grid_2d,DirichletCondition(g,id_n))
+        X = repeat(0:1/10:1, inner = (1,15))
+        Y = repeat(0:1/14:1, outer = (1,11))
+        @test map((x,y)->f(x,y,0), X,Y') ≈ discretize_data(grid_3d,NeumannCondition(f,id_b))
+    end
+end
diff -r b7dcd3dd3181 -r 707fc9761c2b test/SbpOperators/boundary_conditions/sat_test.jl
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test/SbpOperators/boundary_conditions/sat_test.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -0,0 +1,71 @@
+using Test
+
+using Sbplib.Grids
+using Sbplib.LazyTensors
+using Sbplib.SbpOperators
+
+stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order = 4)
+
+struct MockOp end
+
+function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::DirichletCondition; a = 1.)
+    e = boundary_restriction(g, stencil_set, boundary(bc))
+    L = a*e
+    sat_op = e'
+    return sat_op, L
+end
+
+function SbpOperators.sat_tensors(op::MockOp, g::Grid, bc::NeumannCondition)
+    e = boundary_restriction(g, stencil_set, boundary(bc))
+    d = normal_derivative(g, stencil_set, boundary(bc))
+    L = d
+    sat_op = e'
+    return sat_op, L
+end
+
+@testset "sat" begin
+    op = MockOp()
+    @testset "1D" begin
+        grid  = equidistant_grid(0., 1., 11)
+        l, r = boundary_identifiers(grid)
+        u = eval_on(grid, x-> 1. + 2x^2)
+        dc = DirichletCondition(1.0, l)
+        g_l = discretize_data(grid, dc)
+        SAT_l = sat(op, grid, dc)
+        @test SAT_l(u, g_l) ≈ zeros((size(grid))) atol = 1e-13
+        
+        nc = NeumannCondition(4.0, r)
+        g_r = discretize_data(grid, nc)
+        SAT_r = sat(op, grid, nc)
+        @test SAT_r(u, g_r) ≈ zeros((size(grid))) atol = 1e-13
+    end
+    @testset "2D" begin
+        grid  = equidistant_grid((0.,0.), (1.,1.), 11, 13)
+        W, E, S, N = boundary_identifiers(grid)
+        u = eval_on(grid, (x,y) -> x+y^2)
+
+        dc_W = DirichletCondition(1.0, W)
+        SAT_W = sat(op, grid, dc_W)
+        g_W = discretize_data(grid, dc_W)
+        r_W = zeros(size(grid))
+        r_W[1,:] .= map(y -> (y^2-1.), range(0., 1., length=13))
+        @test SAT_W(u, g_W) ≈ r_W atol = 1e-13
+
+        dc_E = DirichletCondition(2, E)
+        SAT_E = sat(op, grid, dc_E; a = 2.)
+        g_E = discretize_data(grid, dc_E)
+        r_E = zeros(size(grid))
+        r_E[end,:] .= map(y -> (2*(1. + y^2)-2.), range(0., 1., length=13))
+        @test SAT_E(u, g_E) ≈ r_E atol = 1e-13
+
+        nc_S = NeumannCondition(.0, S)
+        SAT_S = sat(op, grid, nc_S)
+        g_S = discretize_data(grid, nc_S)
+        @test SAT_S(u, g_S) ≈ zeros(size(grid)) atol = 1e-13
+
+        nc_N = NeumannCondition(2.0, N)
+        SAT_N = sat(op, grid, nc_N)
+        g_N = discretize_data(grid, nc_N)
+        @test SAT_N(u, g_N) ≈ zeros(size(grid)) atol = 1e-13
+    end
+end
diff -r b7dcd3dd3181 -r 707fc9761c2b test/SbpOperators/stencil_set_test.jl
diff -r b7dcd3dd3181 -r 707fc9761c2b test/SbpOperators/volumeops/laplace/laplace_test.jl
--- a/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Jun 26 12:36:41 2024 +0200
+++ b/test/SbpOperators/volumeops/laplace/laplace_test.jl	Wed Jun 26 12:47:26 2024 +0200
@@ -110,3 +110,66 @@
     end
 end
 
+@testset "sat_tensors" begin
+    # TODO: The following tests should be implemented
+    #       1. Symmetry D'H == H'D (test_broken below)
+    #       2. Test eigenvalues of and/or solution to Poisson
+    #       3. Test tuning of Dirichlet conditions
+    #
+    #       These tests are likely easiest to implement once
+    #       we have support for generating matrices from tensors.
+
+    operator_path = sbp_operators_path()*"standard_diagonal.toml"
+    orders = (2,4)
+    tols = (5e-2,5e-4)
+    sz = (201,401)
+    g = equidistant_grid((0.,0.), (1.,1.), sz...)
+    
+    # Verify implementation of sat_tensors by testing accuracy and symmetry (TODO) 
+    # of the operator D = Δ + SAT, where SAT is the tensor composition of the 
+    # operators from sat_tensor. Note that SAT*u should approximate 0 for the 
+    # conditions chosen.
+
+    @testset "Dirichlet" begin
+        for (o, tol) ∈ zip(orders,tols)
+            stencil_set = read_stencil_set(operator_path; order=o)
+            Δ = Laplace(g, stencil_set)
+            H = inner_product(g, stencil_set)
+            u = collect(eval_on(g, (x,y) -> sin(π*x)sin(2*π*y)))
+            Δu = collect(eval_on(g, (x,y) -> -5*π^2*sin(π*x)sin(2*π*y)))
+            D = Δ 
+            for id ∈ boundary_identifiers(g)
+                D = D + foldl(∘, sat_tensors(Δ, g, DirichletCondition(0., id)))
+            end
+            e = D*u .- Δu
+            # Accuracy
+            @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+            # Symmetry
+            r = randn(size(u))
+            @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+        end
+    end
+
+    @testset "Neumann" begin
+        @testset "Dirichlet" begin
+            for (o, tol) ∈ zip(orders,tols)
+                stencil_set = read_stencil_set(operator_path; order=o)
+                Δ = Laplace(g, stencil_set)
+                H = inner_product(g, stencil_set)
+                u = collect(eval_on(g, (x,y) -> cos(π*x)cos(2*π*y)))
+                Δu = collect(eval_on(g, (x,y) -> -5*π^2*cos(π*x)cos(2*π*y)))
+                D = Δ 
+                for id ∈ boundary_identifiers(g)
+                    D = D + foldl(∘, sat_tensors(Δ, g, NeumannCondition(0., id)))
+                end
+                e = D*u .- Δu
+                # Accuracy
+                @test sqrt(sum(H*e.^2)) ≈ 0 atol = tol
+                # Symmetry
+                r = randn(size(u))
+                @test_broken (D'∘H - H∘D)*r .≈ 0 atol = 1e-13 # TODO: Need to implement apply_transpose for D.
+            end
+        end
+    end
+end
+