Skip to content
This repository has been archived by the owner on Jul 19, 2023. It is now read-only.

(WIP) Safer handling of DerivativeOperator.coefficients #244

Open
wants to merge 12 commits into
base: master
Choose a base branch
from
Open

(WIP) Safer handling of DerivativeOperator.coefficients #244

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
6a228b5
removing coefficient_functions.jl and the corresponding test
dangirsh May 5, 2020
7530b13
whitespace
dangirsh May 5, 2020
965bea4
remove references to coefficient_functions.jl and related test
dangirsh May 5, 2020
5c24d99
initialize DerivativeOperator.coefficients to ones instead of undef when
dangirsh May 6, 2020
fcdf775
add function to initialize coefficients based on coeff_func
dangirsh May 6, 2020
f80886e
consolidate coefficient handling to coeffcients.jl
dangirsh May 6, 2020
514a99b
forgot two convolution function in last commit
dangirsh May 6, 2020
c0b85cd
use get_coefficient in concretization routines
dangirsh May 6, 2020
db02dba
small fixes. tests passing
dangirsh May 8, 2020
2d173b5
Merge branch 'master' into fixup-coeff-func
dangirsh Sep 1, 2020
10850b5
Merge branch 'master' into fixup-coeff-func
dangirsh Dec 14, 2020
edd7abd
small fix in Base.copyto! in concretization.jl (credit simenhu)
dangirsh Dec 14, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/DiffEqOperators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,13 +27,13 @@ include("derivative_operators/multi_dim_bc_operators.jl")

### Derivative Operators
include("derivative_operators/fornberg.jl")
include("derivative_operators/coefficients.jl")
include("derivative_operators/derivative_operator.jl")
include("derivative_operators/abstract_operator_functions.jl")
include("derivative_operators/convolutions.jl")
include("derivative_operators/concretization.jl")
include("derivative_operators/ghost_derivative_operator.jl")
include("derivative_operators/derivative_operator_functions.jl")
include("derivative_operators/coefficient_functions.jl")

### Composite Operators
include("composite_operators.jl")
Expand Down
67 changes: 0 additions & 67 deletions src/derivative_operators/coefficient_functions.jl
View file
Open in desktop

This file was deleted.

32 changes: 32 additions & 0 deletions src/derivative_operators/coefficients.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
"""
```
init_coefficients(coeff_func, len::Int)
```

Return the initial value of an operator's `coefficients` field based on the type
of `coeff_func`.
"""
init_coefficients(coeff_func::Nothing, len::Int) = nothing

init_coefficients(coeff_func::Number, len::Int) = coeff_func * ones(typeof(coeff_func), len)

function init_coefficients(coeff_func::AbstractVector{T}, len::Int) where T <: Number
coeff_func
end

init_coefficients(coeff_func::Function, len::Int) = ones(Float64, len)



"""
```
get_coefficient(coefficients, index)
```
"""
get_coefficient(coefficients::AbstractVector, index::Int) = coeff[i]

# FIXME: I don't think this case is used anymore
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This case was covered by a subset of convolution functions, but I don't know what the use-case is.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's for a constant coefficient, like D*A so a constant diffusion. If the coefficient is a number, then the whole setup is actually a bitstype so it'll compile to the GPU and distribute much better, and this is a common case so it's worth supporting. I think that one function is really all that's necessary to support it.

get_coefficient(coefficients::Number, index::Int) = coefficients

# FIXME: Why use "true" here for the value 1?
get_coefficient(coefficients::Nothing, index::Int) = true
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

true is a hard 1 in Julia. true will convert to any type, so it's the recommended value if you want it to promote in further operations.

38 changes: 19 additions & 19 deletions src/derivative_operators/concretization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) wh

coeff = A.coefficients
get_coeff = if coeff isa AbstractVector
i -> coeff[i]
i = get_coefficient(coeff, i)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you reference the variable i before it's defined here? Get an error julia> Array(CenteredDifference{1}(2, 3, 0.01, 100, c_squared)) ERROR: UndefVarError: i not defined Stacktrace: [1] copyto!(::Array{Float64,2}, ::DerivativeOperator{Float64,1,false,Float64,StaticArrays.SArray{Tuple{5},Float64,1,5},StaticArrays.SArray{Tuple{1},StaticArrays.SArray{Tuple{5},Float64,1,5},1,1},Array{Float64,1},Array{Float64,1}}, ::Int64) at /home/simen/Documents/school/masteroppgave/wave_simulation/dev/DiffEqOperators/src/derivative_operators/concretization.jl:16 [2] Array(::DerivativeOperator{Float64,1,false,Float64,StaticArrays.SArray{Tuple{5},Float64,1,5},StaticArrays.SArray{Tuple{1},StaticArrays.SArray{Tuple{5},Float64,1,5},1,1},Array{Float64,1},Array{Float64,1}}, ::Int64) at /home/simen/Documents/school/masteroppgave/wave_simulation/dev/DiffEqOperators/src/derivative_operators/concretization.jl:45 (repeats 2 times) [3] top-level scope at REPL[2]:1 when using a locally merged version of this pull request.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Changed te = to -> in line 16 as the following lines, which works for me.

elseif coeff isa Number
i -> coeff
else
Expand Down Expand Up @@ -314,7 +314,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh
stencils = A.stencil_coefs

for i in 1:bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
cur_stencil = stencils
L[i,i+1:i+stl] = cur_coeff*cur_stencil
Expand All @@ -325,7 +325,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh
end

for i in bpc+1:len-bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
cur_stencil = stencils
cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
if cur_coeff >= 0
Expand All @@ -336,7 +336,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) wh
end

for i in len-bpc+1:len
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

very clean. Beautiful.

if cur_coeff < 0
cur_stencil = stencils
cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
Expand All @@ -358,7 +358,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len)
coeff = A.coefficients

for i in 1:bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
else
Expand All @@ -367,7 +367,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len)
end

for i in bpc+1:len-bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc]
else
Expand All @@ -376,7 +376,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len)
end

for i in len-bpc+1:len
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff < 0
L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc]
else
Expand All @@ -403,7 +403,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=
stencils = A.stencil_coefs

for i in 1:bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
cur_stencil = stencils
L[i,i+1:i+stl] = cur_coeff*cur_stencil
Expand All @@ -414,7 +414,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=
end

for i in bpc+1:len-bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
cur_stencil = stencils
cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
if cur_coeff >= 0
Expand All @@ -425,7 +425,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=
end

for i in len-bpc+1:len
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff < 0
cur_stencil = stencils
cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
Expand All @@ -447,7 +447,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In
coeff = A.coefficients

for i in 1:bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
else
Expand All @@ -456,7 +456,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In
end

for i in bpc+1:len-bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc]
else
Expand All @@ -465,7 +465,7 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::In
end

for i in len-bpc+1:len
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff < 0
L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc]
else
Expand All @@ -492,7 +492,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A
stencils = A.stencil_coefs

for i in 1:bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
cur_stencil = stencils
L[i,i+1:i+stl] = cur_coeff*cur_stencil
Expand All @@ -503,7 +503,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A
end

for i in bpc+1:len-bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
cur_stencil = stencils
cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
if cur_coeff >= 0
Expand All @@ -514,7 +514,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A
end

for i in len-bpc+1:len
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff < 0
cur_stencil = stencils
cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
Expand All @@ -537,7 +537,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
L = BandedMatrix{T}(Zeros(len, len+2), (stl-2, stl))

for i in 1:bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
else
Expand All @@ -546,7 +546,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
end

for i in bpc+1:len-bpc
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff >= 0
L[i,i+1:i+stl] = cur_coeff * A.stencil_coefs[1,i-bpc]
else
Expand All @@ -555,7 +555,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int
end

for i in len-bpc+1:len
cur_coeff = coeff[i]
cur_coeff = get_coefficient(coeff, i)
if cur_coeff < 0
L[i,i-stl+2:i+1] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc]
else
Expand Down
Loading