From 1d635a97b1df9d9472915b3c37c3d5ae610f9625 Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Fri, 24 Dec 2021 02:50:20 +0530 Subject: [PATCH 1/6] hermite-based finite difference fornberg(2020) --- src/discretization/fornberg.jl | 45 +++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 12 deletions(-) diff --git a/src/discretization/fornberg.jl b/src/discretization/fornberg.jl index 0b24ae706..7074e4df9 100644 --- a/src/discretization/fornberg.jl +++ b/src/discretization/fornberg.jl @@ -1,20 +1,23 @@ -############################################################# -# Fornberg algorithm - -# This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0) -# to obtain Finite Difference weights over arbitrary points to arbitrary order. - -function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real - #= +#= Fornberg algorithm +This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0) +and hermite-based finite difference Fornberg(2020) algorithm (https://doi.org/10.1093/imanum/draa006) +to obtain Finite Difference weights over arbitrary points to arbitrary order. +Inputs: order: The derivative order for which we need the coefficients x0 : The point in the array 'x' for which we need the coefficients x : A dummy array with relative coordinates, e.g., central differences need coordinates centred at 0 while those at boundaries need coordinates starting from 0 to the end point + dfdx : optional argument to consider weights of the first-derivative of function or not + if dfdx == false (default kwarg), implies Fornberg(1988) + dfdx == true, implies Fornberg(2020) + Outputs: + if dfdx == false (default kwarg), _C : weights to approximate derivative of required order using function values only. + else, _D,_E : weights to approximate derivative of required order using function and its first- + derivative values respectively. +=# - The approximation order of the stencil is automatically determined from - the number of requested stencil points. - =# +function calculate_weights(order::Int, x0::T, x::AbstractVector; dfdx::Bool = false) where T<:Real N = length(x) @assert order < N "Not enough points for the requested order." M = order @@ -56,5 +59,23 @@ function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real =# _C = C[:,end] _C[div(N,2)+1] -= sum(_C) - return _C + if dfdx == false + return _C + else + A = x .- x'; + s = sum(1 ./ (A + I(N)), dims = 1) .- 1; + cp = factorial.(0:M); + cc = C./cp' + c̃ = zeros(N, M+2); + for k in 1:M+1 + c̃[:,k+1] = sum(cc[:,1:k].*cc[:,k:-1:1], dims = 2); + end + E = c̃[:,1:M+1] - (x .- x0).*c̃[:,2:M+2]; + D = c̃[:,2:M+2] + 2*E.*s'; + D = D.*cp'; + E = E.*cp'; + + _D = D[:,end]; _E = E[:,end] + return _D, _E + end end From 43b99f2858eacf379d42aa159920ddd5f62ef369 Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Fri, 24 Dec 2021 02:55:26 +0530 Subject: [PATCH 2/6] Create MOLfornberg_weights.jl tests for fornberg updated fornberg algorithm --- test/pde_systems/MOLfornberg_weights.jl | 30 +++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 test/pde_systems/MOLfornberg_weights.jl diff --git a/test/pde_systems/MOLfornberg_weights.jl b/test/pde_systems/MOLfornberg_weights.jl new file mode 100644 index 000000000..4ad93be96 --- /dev/null +++ b/test/pde_systems/MOLfornberg_weights.jl @@ -0,0 +1,30 @@ +using Test, LinearAlgebra +using MethodOfLines + +@testset "finite-difference weights from fornberg(1988) & fornberg(2020)" begin + order = 2; z = 0.0; x = [-1, 0, 1.0]; + @test MethodOfLines.calculate_weights(order, z, x) == [1,-2,1] # central difference of second-derivative with unit-step + + order = 1; z = 0.0; x = [-1., 1.0]; + @test MethodOfLines.calculate_weights(order, z, x) == [-0.5,0.5] # central difference of first-derivative with unit step + + order = 1; z = 0.0; x = [0, 1]; + @test MethodOfLines.calculate_weights(order, z, x) == [-1, 1] # forward difference + + order = 1; z = 1.0; x = [0, 1]; + @test MethodOfLines.calculate_weights(order, z, x) == [-1, 1] # backward difference + + # forward-diff of third derivative with order of accuracy == 3 + order = 3; z = 0.0; x = [0,1,2,3,4,5] + @test MethodOfLines.calculate_weights(order, z, x) == [-17/4, 71/4 ,−59/2, 49/2, −41/4, 7/4] + + order = 3; z = 0.0; x = collect(-3:3) + d, e = MethodOfLines.calculate_weights(order, z, x;dfdx = true) + @test d ≈ [-167/18000, -963/2000, -171/16,0,171/16,963/2000,167/18000] + @test e ≈ [-1/600,-27/200,-27/8,-49/3,-27/8,-27/200,-1/600] + + order = 3; z = 0.0; x = collect(-4:4) + d, e = MethodOfLines.calculate_weights(order, z, x;dfdx = true) + @test d ≈ [-2493/5488000, -12944/385875, -87/125 ,-1392/125,0,1392/125,87/125,12944/385875,2493/5488000] + @test e ≈ [-3/39200,-32/3675,-6/25,-96/25,-205/12, -96/25, -6/25,-32/3675,-3/39200] +end From cfa33c419be36ac095ea666b3d81bb03433a7723 Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Fri, 24 Dec 2021 02:57:55 +0530 Subject: [PATCH 3/6] Update runtests.jl runtests for fornberg --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 20176bd9f..0aa24b867 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,6 +10,7 @@ const is_TRAVIS = haskey(ENV,"TRAVIS") if GROUP == "All" || GROUP == "MOLFiniteDifference" @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end + @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLfornberg_weights.jl") end #@time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("pde_systems/MOL_1D_Linear_Convection.jl") end @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin include("pde_systems/MOL_1D_Linear_Diffusion.jl") end @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin include("pde_systems/MOL_1D_NonLinear_Diffusion.jl") end From 9f9f446a81f9645ed374cf310b68d7e761fc5948 Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Wed, 2 Feb 2022 19:45:22 +0530 Subject: [PATCH 4/6] runtests for fornberg file --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index ed5d7fa10..ee37a07ca 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,7 @@ const is_TRAVIS = haskey(ENV,"TRAVIS") if GROUP == "All" || GROUP == "Integration Tests" @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLtest.jl") end + @time @safetestset "MOLFiniteDifference Interface" begin include("pde_systems/MOLfornberg_weights.jl") end #@time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("pde_systems/MOL_1D_Linear_Convection.jl") end @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin include("pde_systems/MOL_1D_Linear_Diffusion.jl") end @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin include("pde_systems/MOL_1D_NonLinear_Diffusion.jl") end From 094a4f965f0830493547be6c2158f777ca66ce9e Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Wed, 2 Feb 2022 19:48:43 +0530 Subject: [PATCH 5/6] hermite-based weights of Finite difference(2020) --- src/discretization/fornberg.jl | 45 +++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 12 deletions(-) diff --git a/src/discretization/fornberg.jl b/src/discretization/fornberg.jl index 1406d8bbf..c90339842 100644 --- a/src/discretization/fornberg.jl +++ b/src/discretization/fornberg.jl @@ -1,20 +1,23 @@ -############################################################# -# Fornberg algorithm - -# This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0) -# to obtain Finite Difference weights over arbitrary points to arbitrary order. - -function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real - #= +#= Fornberg algorithm +This implements the Fornberg (1988) algorithm (https://doi.org/10.1090/S0025-5718-1988-0935077-0) +and hermite-based finite difference Fornberg(2020) algorithm (https://doi.org/10.1093/imanum/draa006) +to obtain Finite Difference weights over arbitrary points to arbitrary order. +Inputs: order: The derivative order for which we need the coefficients x0 : The point in the array 'x' for which we need the coefficients x : A dummy array with relative coordinates, e.g., central differences need coordinates centred at 0 while those at boundaries need coordinates starting from 0 to the end point + dfdx : optional argument to consider weights of the first-derivative of function or not + if dfdx == false (default kwarg), implies Fornberg(1988) + dfdx == true, implies Fornberg(2020) + Outputs: + if dfdx == false (default kwarg), _C : weights to approximate derivative of required order using function values only. + else, _D,_E : weights to approximate derivative of required order using function and its first- + derivative values respectively. +=# - The approximation order of the stencil is automatically determined from - the number of requested stencil points. - =# +function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real N = length(x) @assert order < N "Not enough points for the requested order." M = order @@ -58,5 +61,23 @@ function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real if order != 0 _C[div(N,2)+1] -= sum(_C) end - return _C + if dfdx == false + return _C + else + A = x .- x'; + s = sum(1 ./ (A + I(N)), dims = 1) .- 1; + cp = factorial.(0:M); + cc = C./cp' + c̃ = zeros(N, M+2); + for k in 1:M+1 + c̃[:,k+1] = sum(cc[:,1:k].*cc[:,k:-1:1], dims = 2); + end + E = c̃[:,1:M+1] - (x .- x0).*c̃[:,2:M+2]; + D = c̃[:,2:M+2] + 2*E.*s'; + D = D.*cp'; + E = E.*cp'; + + _D = D[:,end]; _E = E[:,end] + return _D, _E + end end From d71125ef9befa961f9d1315c5b8bd0974b47c249 Mon Sep 17 00:00:00 2001 From: Nikhil Yewale Date: Wed, 2 Feb 2022 20:03:58 +0530 Subject: [PATCH 6/6] adding kwarg dfdx::Bool forgot earlier --- src/discretization/fornberg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/discretization/fornberg.jl b/src/discretization/fornberg.jl index c90339842..3aac12116 100644 --- a/src/discretization/fornberg.jl +++ b/src/discretization/fornberg.jl @@ -17,7 +17,7 @@ Inputs: derivative values respectively. =# -function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real +function calculate_weights(order::Int, x0::T, x::AbstractVector; dfdx::Bool = false) where T<:Real N = length(x) @assert order < N "Not enough points for the requested order." M = order