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

Commit

Permalink
Inplace Operation
Browse files Browse the repository at this point in the history
  • Loading branch information
mjsheikh committed Feb 6, 2021
1 parent e60c3bd commit 61366b5
Show file tree
Hide file tree
Showing 3 changed files with 18 additions and 20 deletions.
6 changes: 3 additions & 3 deletions examples/Fast_Diffusion.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
2 changes: 1 addition & 1 deletion src/DiffEqOperators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
30 changes: 14 additions & 16 deletions src/derivative_operators/derivative_operator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 61366b5

Please sign in to comment.