From 303885307554062af85f4af732cf2cab651b275b Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Tue, 26 Jan 2021 00:57:27 +0530 Subject: [PATCH 1/3] adding a nonlinear diffusion operator to derivative_operator.jl --- 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..c934c6039 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}}, + 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 95264651055284717290d5ba34567a12652977fa Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Tue, 26 Jan 2021 12:56:37 +0530 Subject: [PATCH 2/3] syntax corrections from the last commit --- src/derivative_operators/derivative_operator.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index c934c6039..7da89483c 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -33,20 +33,21 @@ function NonLinearDiffusion(second_differential_order::Int, first_differential_o 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 + #first differential is the inner one, second is the outer differential @assert approx_order>1 "approximation_order must be greater than 1." - discretization = zeros(nknots); + 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); + 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); + 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); + 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)]; + discretization += (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)] return discretization end From fcf616d543b6edaeeba042e5af639eb6c1224310 Mon Sep 17 00:00:00 2001 From: mjsheikh Date: Sun, 31 Jan 2021 13:08:24 +0530 Subject: [PATCH 3/3] Type casting fix for autodiff --- src/derivative_operators/derivative_operator.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl index 7da89483c..fbfad008c 100644 --- a/src/derivative_operators/derivative_operator.jl +++ b/src/derivative_operators/derivative_operator.jl @@ -29,7 +29,7 @@ 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}}, + 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