space-dependent coefficients

[ignace-computing]

Hello. Could someone please comment on using space-dependent coefficients for convection, diffusion and convection-diffusion equations? Sorry if this is not the right place to ask such things. If not yet available, I might want to help with this.

[ChrisRackauckas]

If you multiply by an array of coefficients it should work.

[michaelsachs]

Hi, I thought about using space-dependent coefficients to restrict some processes to certain regions in space but could not get past the discretization step:

When multiplying by an array of coefficients I am not sure how to handle broadcasting within the symbolic equation and end up with ERROR: LoadError: InvalidSystemException: The system is unbalanced. There are 101 highest order derivative variables and 9803 equations.

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x
@variables n(..) 
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# space-dependent coefficient 
coeff = [ones(39); zeros(60)]

# 1D PDE and boundary conditions
eqs  = @. Dt(n(t,x)) ~ Dxx(n(t,x)) - coeff*n(t,x)

bcs = [n(0,x) ~ 2x,
        Dx(n(t,0)) ~ 0,
        Dx(n(t,10)) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,100.0),
           x ∈ IntervalDomain(0.0,10.0)]

# PDE system
pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x)],[coeff])

# Method of lines discretization
Δx = 0.1
discretization = MOLFiniteDifference([x=>Δx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

Alternatively, when adding coeff(x) as a space-dependent variable the resulting error is ERROR: LoadError: ArgumentError: Dict(kv): kv needs to be an iterator of tuples or pairs. I guess this suggests that the actual values for coeff are needed when building the ODESystem in MOL_discretization.jl, but I am not sure where and how to pass them in?

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x
@variables n(..) coeff(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eqs  = Dt(n(t,x)) ~ Dxx(n(t,x)) - coeff(x)*n(t,x)

bcs = [n(0,x) ~ 2x,
        Dx(n(t,0)) ~ 0,
        Dx(n(t,10)) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,100.0),
           x ∈ IntervalDomain(0.0,10.0)]

# PDE system
pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x)],[coeff(x)])

# Method of lines discretization
Δx = 0.1
discretization = MOLFiniteDifference([x=>Δx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

[ChrisRackauckas]

oh with MOLFiniteDifference just make the coefficient be a function: eqs = @. Dt(n(t,x)) ~ Dxx(n(t,x)) - f(x)*n(t,x)

[michaelsachs]

Thanks! That makes sense and indeed works.

[ignace-computing]

Hello. Please give me a clue how one can deal with coefficients inside the spatial derivative operator. I have tried, among others, the following:

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, Test, Plots 

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
g = x -> (x-5.)^2
f = x -> 1.
eq = Dt(u(t,x)) ~ Dxx(f(x)*u(t,x)) - g(x)*u(t,x)
bcs = [u(0,x) ~ cos(x),
       Dx(u(t,0)) ~ 0,
       Dx(u(t,Float64(pi))) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,1.0),
           x ∈ IntervalDomain(0.0,Float64(pi))]

# PDE system
pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

# Method of lines discretization
dx = range(0.0,Float64(π),length=300)
order = 2
discretization = MOLFiniteDifference([x=>dx],t)
discretization_edge = MOLFiniteDifference([x=>dx],t;grid_align=center_align)

# Convert the PDE problem into an ODE problem
p1 = plot()
sol = []
for disc in [discretization, discretization_edge]
    prob = discretize(pdesys,disc)

    # Solve ODE problem
    sol = solve(prob,Tsit5(),saveat=0.1)

    if disc.grid_align == center_align
        x_sol = dx[2:end-1]
    else
        x_sol = (dx[1:end-1]+dx[2:end])/2
    end
    t_sol = sol.t

    p1 = plot(x_sol, sol.u[end])
end

p1

which does not throw any errors.

However, when one changes the definition of f as follows:

f = x -> 1. + x

it throws an error: ERROR: LoadError: InvalidSystemException: The system is unbalanced. There are 598 highest order derivative variables and 300 equations.

[ChrisRackauckas]

That should be supported by https://github.com/SciML/DiffEqOperators.jl/pull/371 ?