From a686de73e8e5e06610e729d6ea625f2b639690bc Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Sun, 31 Jan 2021 19:33:07 +0530 Subject: [PATCH 1/6] Adding NonLinearDiffusion operator to derivate_operator.jl returning discretized function --- src/DiffEqOperators.jl | 2 +- .../derivative_operator.jl | 25 +++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 0af00efb2..9af8bd2a5 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -59,7 +59,7 @@ end export MatrixFreeOperator export AnalyticalJacVecOperator, JacVecOperator, getops export AbstractDerivativeOperator, DerivativeOperator, - CenteredDifference, UpwindDifference + CenteredDifference, UpwindDifference, NonLinearDiffusion export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC, MultiDimDirectionalBC, ComposedMultiDimBC export compose, decompose, perpsize diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index e87535959..2115da40d 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,6 +26,31 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end +struct NonLinearDiffusion end + +function NonLinearDiffusion(second_differential_order::Int, first_differential_order::Int, approx_order::Int, + p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real}, + nknots::Int) where {T<:Real,N} + #p is given by bc1*k , k being the diffusion coefficient + #q is given by bc2*u , u being the desired function + @assert approx_order>1 "approximation_order must be greater than 1." + discretization = zeros(nknots); + + if first_differential_order > 0 + discretization += (CenteredDifference(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p); + else + discretization += q[2:(nknots + 1)].*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p); + end + + for l = 1:(second_differential_order - 1) + discretization += binomial(second_differential_order,l)*(CenteredDifference(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order - l,approx_order,dx,nknots)*p); + end + + discretization += (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]; + return discretization + +end + struct CenteredDifference{N} end function CenteredDifference{N}(derivative_order::Int, From 56c251c8ed4d1ae757d92917acaf94354426e836 Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Mon, 1 Feb 2021 01:05:42 +0530 Subject: [PATCH 2/6] Test added --- examples/Fast_Diffusion.jl | 46 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 examples/Fast_Diffusion.jl diff --git a/examples/Fast_Diffusion.jl b/examples/Fast_Diffusion.jl new file mode 100644 index 000000000..2b56a7651 --- /dev/null +++ b/examples/Fast_Diffusion.jl @@ -0,0 +1,46 @@ +# # Fast-Diffusion Problem +# +# This example demonstrates the use of 'NonLinearDiffusion' operator to solve time-dependent non-linear diffusion PDEs with coefficient having dependence on the unknown function. +# Here we consider a fast diffusion problem with Dirichlet BCs on unit interval: +# ∂ₜu = ∂ₓ(k*∂ₓu) +# k = 1/u² +# u(x=0,t) = exp(-t) +# u(x=1,t) = 1/(1.0 + exp(2t)) +# u(x, t=0) = u₀(x) + +# The analytical solution for this is given by : + +u_analytic(x, t) = 1 / sqrt(x^2 + exp(2*t)) + +# +# Reproducing it numerically +# + +using DiffEqOperators, OrdinaryDiffEq + +nknots = 100 +h = 1.0/(nknots+1) +knots = range(h, step=h, length=nknots) +n = 1 # Outer differential order +m = 1 # Inner differential order +approx_ord = 2 + +u0 = u_analytic.(knots,0.0) + +t0 = 0.0 +t1 = 1.0 + +function f(u,p,t) + bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5)) + l = bc*u + k = l.^(-2) # Diffusion Coefficient + NonLinearDiffusion(n,m,approx_ord,k,l,h,nknots) +end + +prob = ODEProblem(f, u0, (t0, t1)) +alg = KenCarp4() +sol = solve(prob,alg) + +using Test +@test u_analytic.(knots, t1) ≈ sol(t1) rtol=1e-3 + From 5a5397d85029fe8665c9360ab4c7c5e72d11fad0 Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Sun, 7 Feb 2021 00:57:51 +0530 Subject: [PATCH 3/6] Inplace Operation --- examples/Fast_Diffusion.jl | 6 ++-- src/DiffEqOperators.jl | 2 +- .../derivative_operator.jl | 30 +++++++++---------- 3 files changed, 18 insertions(+), 20 deletions(-) diff --git a/examples/Fast_Diffusion.jl b/examples/Fast_Diffusion.jl index 2b56a7651..1c67a69d2 100644 --- a/examples/Fast_Diffusion.jl +++ b/examples/Fast_Diffusion.jl @@ -26,15 +26,15 @@ m = 1 # Inner differential order approx_ord = 2 u0 = u_analytic.(knots,0.0) - +du = similar(u0) t0 = 0.0 t1 = 1.0 -function f(u,p,t) +function f(du,u,p,t) bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5)) l = bc*u k = l.^(-2) # Diffusion Coefficient - NonLinearDiffusion(n,m,approx_ord,k,l,h,nknots) + NonLinearDiffusion!(du,n,m,approx_ord,k,l,h,nknots) end prob = ODEProblem(f, u0, (t0, t1)) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 9af8bd2a5..6268f7e3b 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -59,7 +59,7 @@ end export MatrixFreeOperator export AnalyticalJacVecOperator, JacVecOperator, getops export AbstractDerivativeOperator, DerivativeOperator, - CenteredDifference, UpwindDifference, NonLinearDiffusion + CenteredDifference, UpwindDifference, NonLinearDiffusion! export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC, MultiDimDirectionalBC, ComposedMultiDimBC export compose, decompose, perpsize diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 2115da40d..5bc11f0e7 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,29 +26,27 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end -struct NonLinearDiffusion end +struct NonLinearDiffusion! end + +function NonLinearDiffusion!(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int, + p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real}, + nknots::Int) where {T<:Real, N} + #q is given by bc*u , u being the unknown function + #p is given as some function of q , p being the diffusion coefficient -function NonLinearDiffusion(second_differential_order::Int, first_differential_order::Int, approx_order::Int, - p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real}, - nknots::Int) where {T<:Real,N} - #p is given by bc1*k , k being the diffusion coefficient - #q is given by bc2*u , u being the desired function @assert approx_order>1 "approximation_order must be greater than 1." - discretization = zeros(nknots); - if first_differential_order > 0 - discretization += (CenteredDifference(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p); + du .= (CenteredDifference(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p) else - discretization += q[2:(nknots + 1)].*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p); + du .= q[2:(nknots + 1)].*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p) end - + for l = 1:(second_differential_order - 1) - discretization += binomial(second_differential_order,l)*(CenteredDifference(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order - l,approx_order,dx,nknots)*p); + du .= du .+ binomial(second_differential_order,l)*(CenteredDifference(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order - l,approx_order,dx,nknots)*p) end - - discretization += (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]; - return discretization - + + du .= du .+ (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)] + end struct CenteredDifference{N} end From b5a4510a8e21445ea682de48f6970d9503537349 Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Sun, 7 Feb 2021 08:56:52 +0530 Subject: [PATCH 4/6] Test set added --- examples/Fast_Diffusion.jl | 46 ----------------- .../derivative_operator.jl | 2 - test/Fast_Diffusion.jl | 50 +++++++++++++++++++ 3 files changed, 50 insertions(+), 48 deletions(-) delete mode 100644 examples/Fast_Diffusion.jl create mode 100644 test/Fast_Diffusion.jl diff --git a/examples/Fast_Diffusion.jl b/examples/Fast_Diffusion.jl deleted file mode 100644 index 1c67a69d2..000000000 --- a/examples/Fast_Diffusion.jl +++ /dev/null @@ -1,46 +0,0 @@ -# # Fast-Diffusion Problem -# -# This example demonstrates the use of 'NonLinearDiffusion' operator to solve time-dependent non-linear diffusion PDEs with coefficient having dependence on the unknown function. -# Here we consider a fast diffusion problem with Dirichlet BCs on unit interval: -# ∂ₜu = ∂ₓ(k*∂ₓu) -# k = 1/u² -# u(x=0,t) = exp(-t) -# u(x=1,t) = 1/(1.0 + exp(2t)) -# u(x, t=0) = u₀(x) - -# The analytical solution for this is given by : - -u_analytic(x, t) = 1 / sqrt(x^2 + exp(2*t)) - -# -# Reproducing it numerically -# - -using DiffEqOperators, OrdinaryDiffEq - -nknots = 100 -h = 1.0/(nknots+1) -knots = range(h, step=h, length=nknots) -n = 1 # Outer differential order -m = 1 # Inner differential order -approx_ord = 2 - -u0 = u_analytic.(knots,0.0) -du = similar(u0) -t0 = 0.0 -t1 = 1.0 - -function f(du,u,p,t) - bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5)) - l = bc*u - k = l.^(-2) # Diffusion Coefficient - NonLinearDiffusion!(du,n,m,approx_ord,k,l,h,nknots) -end - -prob = ODEProblem(f, u0, (t0, t1)) -alg = KenCarp4() -sol = solve(prob,alg) - -using Test -@test u_analytic.(knots, t1) ≈ sol(t1) rtol=1e-3 - diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 5bc11f0e7..871768665 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,8 +26,6 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end -struct NonLinearDiffusion! end - function NonLinearDiffusion!(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int, p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real}, nknots::Int) where {T<:Real, N} diff --git a/test/Fast_Diffusion.jl b/test/Fast_Diffusion.jl new file mode 100644 index 000000000..75ff8045c --- /dev/null +++ b/test/Fast_Diffusion.jl @@ -0,0 +1,50 @@ +# # Fast-Diffusion Problem +# +# This example demonstrates the use of 'NonLinearDiffusion' operator to solve time-dependent non-linear diffusion PDEs with coefficient having dependence on the unknown function. +# Here we consider a fast diffusion problem with Dirichlet BCs on unit interval: +# ∂ₜu = ∂ₓ(k*∂ₓu) +# k = 1/u² +# u(x=0,t) = exp(-t) +# u(x=1,t) = 1/(1.0 + exp(2t)) +# u(x, t=0) = u₀(x) +using Test +using DiffEqOperators, OrdinaryDiffEq + +@testset "1D Nonlinear fast diffusion equation" begin + # The analytical solution for this is given by : + + u_analytic(x, t) = 1 / sqrt(x^2 + exp(2*t)) + + # + # Reproducing it numerically + # + + + nknots = 100 + h = 1.0/(nknots+1) + knots = range(h, step=h, length=nknots) + n = 1 # Outer differential order + m = 1 # Inner differential order + approx_ord = 2 + + u0 = u_analytic.(knots,0.0) + du = similar(u0) + t0 = 0.0 + t1 = 1.0 + + function f(du,u,p,t) + bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5)) + l = bc*u + k = l.^(-2) # Diffusion Coefficient + NonLinearDiffusion!(du,n,m,approx_ord,k,l,h,nknots) + end + + prob = ODEProblem(f, u0, (t0, t1)) + alg = KenCarp4() + sol = solve(prob,alg) + + for t in t0:0.1:t1 + @test u_analytic.(knots, t) ≈ sol(t) rtol=1e-3 + end + +end \ No newline at end of file From ba2870ef995a1397d4eaebcc5cb7b13cb91ea455 Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Sun, 7 Feb 2021 09:28:06 +0530 Subject: [PATCH 5/6] runtest added --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index c5c8b74ea..45c2fc123 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,7 @@ if GROUP == "All" || GROUP == "Interface" @time @safetestset "Validate Boundary Padded Array Concretization" begin include("boundary_padded_array.jl") end @time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("multi_dim_bc_test.jl") end @time @safetestset "2nd order check" begin include("2nd_order_check.jl") end + @time @safetestset "Non-linear Diffusion" begin include("Fast_Diffusion.jl") end #@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails #@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end @time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end From 70944f2706476b70ffa74096f55d8bd91329cc1b Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Mon, 8 Feb 2021 00:26:26 +0530 Subject: [PATCH 6/6] out of place function added --- src/DiffEqOperators.jl | 2 +- src/derivative_operators/derivative_operator.jl | 11 ++++++++++- test/Fast_Diffusion.jl | 2 +- 3 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl index 6268f7e3b..a498b3457 100644 --- a/src/DiffEqOperators.jl +++ b/src/DiffEqOperators.jl @@ -59,7 +59,7 @@ end export MatrixFreeOperator export AnalyticalJacVecOperator, JacVecOperator, getops export AbstractDerivativeOperator, DerivativeOperator, - CenteredDifference, UpwindDifference, NonLinearDiffusion! + CenteredDifference, UpwindDifference, nonlinear_diffusion, nonlinear_diffusion! export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC, MultiDimDirectionalBC, ComposedMultiDimBC export compose, decompose, perpsize diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 871768665..1539d71f5 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -26,7 +26,7 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv coeff_func :: F end -function NonLinearDiffusion!(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int, +function nonlinear_diffusion!(du::AbstractVector{T}, second_differential_order::Int, first_differential_order::Int, approx_order::Int, p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real}, nknots::Int) where {T<:Real, N} #q is given by bc*u , u being the unknown function @@ -47,6 +47,15 @@ function NonLinearDiffusion!(du::AbstractVector{T}, second_differential_order::I end +# An out of place workaround for the mutating version +function nonlinear_diffusion(second_differential_order::Int, first_differential_order::Int, approx_order::Int, + p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T , AbstractVector{T} , Real}, + nknots::Int) where {T<:Real, N} + + du = similar(q,length(q) - 2) + return nonlinear_diffusion!(du,second_differential_order,first_differential_order,approx_order,p,q,dx,nknots) +end + struct CenteredDifference{N} end function CenteredDifference{N}(derivative_order::Int, diff --git a/test/Fast_Diffusion.jl b/test/Fast_Diffusion.jl index 75ff8045c..ca76c9c89 100644 --- a/test/Fast_Diffusion.jl +++ b/test/Fast_Diffusion.jl @@ -36,7 +36,7 @@ using DiffEqOperators, OrdinaryDiffEq bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5)) l = bc*u k = l.^(-2) # Diffusion Coefficient - NonLinearDiffusion!(du,n,m,approx_ord,k,l,h,nknots) + nonlinear_diffusion!(du,n,m,approx_ord,k,l,h,nknots) end prob = ODEProblem(f, u0, (t0, t1))