Mercurial > repos > public > sbplib_julia

changeset 927:d360fc2d9620 feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 21 Feb 2022 23:36:41 +0100
parents 47425442bbc5 (current diff) de1625deb27e (diff)
children d83f685f1031
files test/SbpOperators/boundaryops/normal_derivative_test.jl
diffstat 4 files changed, 105 insertions(+), 37 deletions(-) [+]
line wrap: on
line diff
--- a/src/Grids/EquidistantGrid.jl	Mon Feb 21 23:33:29 2022 +0100
+++ b/src/Grids/EquidistantGrid.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -1,27 +1,23 @@
-"""
-    EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T})
-	EquidistantGrid{T}()
-
-`EquidistantGrid` is a grid with equidistant grid spacing per coordinat direction.
+export EquidistantGrid
+export spacing
+export inverse_spacing
+export restrict
+export boundary_identifiers
+export boundary_grid
+export refine
+export coarsen
 
-`EquidistantGrid(size, limit_lower, limit_upper)` construct the grid with the
-domain defined by the two points P1, and P2 given by `limit_lower` and
-`limit_upper`. The length of the domain sides are given by the components of
-(P2-P1). E.g for a 2D grid with P1=(-1,0) and P2=(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 `size`.
+"""
+    EquidistantGrid{Dim,T<:Real} <: AbstractGrid
 
-`EquidistantGrid{T}()` constructs a 0-dimensional grid.
-
+`Dim`-dimensional equidistant grid with coordinates of type `T`.
 """
 struct EquidistantGrid{Dim,T<:Real} <: AbstractGrid
     size::NTuple{Dim, Int}
     limit_lower::NTuple{Dim, T}
     limit_upper::NTuple{Dim, T}
 
-    # General constructor
-    function EquidistantGrid(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where Dim where 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
@@ -30,12 +26,31 @@
         end
         return new{Dim,T}(size, limit_lower, limit_upper)
     end
+end
 
-	# Specialized constructor for 0-dimensional grid
-	EquidistantGrid{T}() where T = new{0,T}((),(),())
+"""
+    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.
+
+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
-export EquidistantGrid
+
+"""
+    EquidistantGrid{T}()
 
+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)
@@ -62,18 +77,16 @@
 """
     spacing(grid::EquidistantGrid)
 
-The spacing between the grid points of the grid.
+The spacing between grid points.
 """
 spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1)
-export spacing
 
 """
     inverse_spacing(grid::EquidistantGrid)
 
-The reciprocal of the spacing between the grid points of the grid.
+The reciprocal of the spacing between grid points.
 """
 inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid)
-export inverse_spacing
 
 """
     points(grid::EquidistantGrid)
@@ -93,7 +106,7 @@
 """
     restrict(::EquidistantGrid, dim)
 
-Pick out given dimensions from the grid and return a grid for them
+Pick out given dimensions from the grid and return a grid for them.
 """
 function restrict(grid::EquidistantGrid, dim)
     size = grid.size[dim]
@@ -102,7 +115,6 @@
 
     return EquidistantGrid(size, limit_lower, limit_upper)
 end
-export restrict
 
 """
     boundary_identifiers(::EquidistantGrid)
@@ -114,7 +126,6 @@
 	 ...)
 """
 boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),dimension(g)))...)...,)
-export boundary_identifiers
 
 
 """
@@ -133,5 +144,40 @@
 	end
     return restrict(grid,orth_dims)
 end
-export boundary_grid
 boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
+
+
+"""
+    refine(grid::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{dimension(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
+end
+
+"""
+    coarsen(grid::EquidistantGrid, r::Int)
+
+Coarsens `grid` by a factor `r`. The factor is applied to the number of
+intervals which is 1 less than the size of the grid. If the number of
+intervals are not divisible by `r` an error is raised.
+
+See also: [`refine`](@ref)
+"""
+function coarsen(grid::EquidistantGrid, r::Int)
+    sz = size(grid)
+
+    if !all(n -> (n % r == 0), sz.-1)
+        throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
+    end
+
+    new_sz = (sz .- 1).÷r .+ 1
+
+    return EquidistantGrid{dimension(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
+end
--- a/src/SbpOperators/operators/standard_diagonal.toml	Mon Feb 21 23:33:29 2022 +0100
+++ b/src/SbpOperators/operators/standard_diagonal.toml	Mon Feb 21 23:36:41 2022 +0100
@@ -31,7 +31,7 @@
 ]
 
 e.closure = ["1"]
-d1.closure = {s = ["-3/2", "2", "-1/2"], c = 1}
+d1.closure = {s = ["3/2", "-2", "1/2"], c = 1}
 
 [[stencil_set]]
 
@@ -57,4 +57,4 @@
 ]
 
 e.closure = ["1"]
-d1.closure = {s = ["-11/6", "3", "-3/2", "1/3"], c = 1}
+d1.closure = {s = ["11/6", "-3", "3/2", "-1/3"], c = 1}
--- a/test/Grids/EquidistantGrid_test.jl	Mon Feb 21 23:33:29 2022 +0100
+++ b/test/Grids/EquidistantGrid_test.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -98,4 +98,26 @@
                 @test_throws DomainError boundary_grid(g,CartesianBoundary{4,Lower}())
             end
     end
+
+    @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
 end
--- a/test/SbpOperators/boundaryops/normal_derivative_test.jl	Mon Feb 21 23:33:29 2022 +0100
+++ b/test/SbpOperators/boundaryops/normal_derivative_test.jl	Mon Feb 21 23:36:41 2022 +0100
@@ -41,10 +41,10 @@
             (d_w, d_e, d_s, d_n) = 
                 map(id -> normal_derivative(g_2D, d_closure, id), boundary_identifiers(g_2D))
 
-            @test d_w*v ≈ v∂x[1,:] atol = 1e-13
-            @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
-            @test d_s*v ≈ v∂y[:,1] atol = 1e-13
-            @test d_n*v ≈ -v∂y[:,end] atol = 1e-13
+            @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
+            @test d_e*v ≈ v∂x[end,:] atol = 1e-13
+            @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
+            @test d_n*v ≈ v∂y[:,end] atol = 1e-13
         end
 
         @testset "4th order" begin
@@ -53,10 +53,10 @@
             (d_w, d_e, d_s, d_n) = 
                 map(id -> normal_derivative(g_2D, d_closure, id), boundary_identifiers(g_2D))
 
-            @test d_w*v ≈ v∂x[1,:] atol = 1e-13
-            @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
-            @test d_s*v ≈ v∂y[:,1] atol = 1e-13
-            @test d_n*v ≈ -v∂y[:,end] atol = 1e-13
+            @test d_w*v ≈ -v∂x[1,:] atol = 1e-13
+            @test d_e*v ≈ v∂x[end,:] atol = 1e-13
+            @test d_s*v ≈ -v∂y[:,1] atol = 1e-13
+            @test d_n*v ≈ v∂y[:,end] atol = 1e-13
         end
     end
 end