sciml / methodoflines.jl Goto Github PK

Automatic Finite Difference PDE solving with Julia SciML

Home Page: https://docs.sciml.ai/MethodOfLines/stable/

License: MIT License

Julia 100.00%

pde finite-difference-method sciml scientific-ml partial-differential-equations julia high-performance-computing symbolic-computation parallel-computing

methodoflines.jl's People

Contributors

Stargazers

Watchers

methodoflines.jl's Issues

Second order time derivative - The Wave Equations

@xtalax Could you help me , implement this feature ? Any pointers on how to get started would be helpful

Structural simplify Nonlinear problems

Due to recent upgrades to MTK, it is now possible to do this for performance and retrieve the solution. Credit to @YingboMa.

Finalize upwinding in MOLFiniteDifference

This requires a rule application so that the substitution can look for an arbitrary a *D(u) and rewrite with a . Almost working code is there as a comment:

        forward_weights(i,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][i[j]],space[j][i[j]+1]])
        reverse_weights(i,j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][i[j]-1],space[j][i[j]]])
        upwinding_rules = [@rule(*(~~a,(Differential(nottime[j]))(u),~~b) => IfElse.ifelse(*(~~a..., ~~b...,)>0,
                                *(~~a..., ~~b..., dot(reverse_weights(i,j),depvars[k][central_neighbor_idxs(i,j)[1:2]])),
                                *(~~a..., ~~b..., dot(forward_weights(i,j),depvars[k][central_neighbor_idxs(i,j)[2:3]]))))
                                for j in 1:length(nottime), k in 1:length(pdesys.depvars)]

but this requires an additional feature in the SymbolicUtils rules for interpolation into rules.

@shashi @YingboMa

Differential w.r.t. multiple variables Any[] not allowed during MOLFiniteDifference discretize

I'm trying to discretize the following Nerst-Plank system in a single spatial dimension; however, I've received the following error during the check_equations() call:

ERROR: ArgumentError: Differential w.r.t. multiple variables Any[τ, Z] are not allowed.

Full code:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters Z τ
@variables c(..) ϕ(..) S(..)
Dz = Differential(Z)
Dzz = Differential(Z)^2
Dτ = Differential(τ)

Z_min = Z_min = τ_min = 0.0
Z_max = 1.0
τ_max = 50.0

# Parameters
iₐ = 10.
i₀ = 20.
D₁ = 1e-11
F = 96487.
R = 8.314
T = 298.15
Lₓ = 50e-6
c₀ = 1000.
Mᵥ = 1.2998e-5
Cᵣ = cᵥ = 1.

# Dimensionless Parameters
δ = iₐ*Lₓ*Cᵣ/(2*F*D₁*c₀)
tD = abs(Cᵣ)*Lₓ^2/(3600*D₁)
Da = i₀*Lₓ/(2*F*D₁*c₀)
Pe = Mᵥ*i₀*3600/(Lₓ*Cᵣ*F)
ϕ₀ = R*T/F

# PDE
eq = [
       Dz(c(Z,τ))*Dz(ϕ(Z,τ))+c(Z,τ)*(Dzz(ϕ(Z,τ)))/(1.0-S(τ))^2 ~ 0.0, 
       Dτ(c(Z,τ)) ~ (-Pe*(c(Z,τ)/c₀)^(1/2)*(-ϕ(Z,τ)/ϕ₀)*(1.0-Z)*(Dz(c(Z,τ))))/(tD*(1.0-S(τ)))+(Dzz(c(Z,τ))/(tD*(1.0-S(τ))^2))
]
    
    
domains = [Z ∈ Interval(Z_min, Z_max),
           τ ∈ Interval(τ_min, τ_max)]

# BCs
bcs =  [   
        c(Z,0) ~ 1.0, #Concentration @ τ=0
        c(Z_max,τ)*(Dz(ϕ(1,τ))) ~ -δ*(1-S(τ)), #ϕ @ Z = 1
        ϕ(Z_max,τ) ~ 0.0, #ϕ @ Z = 1
        Dz(c(Z_max,τ)) ~ -δ*(1-S(τ)), #Concentration @ Z = 1
        Dz(ϕ(Z_max,τ)) ~ 0.0, #ϕ @ Z = 0
        Dz(ϕ(Z_min,τ))/(1-S(τ)) ~ -Da*(c(Z_min,τ)^(1/2))-ϕ(Z_min,τ), #ϕ BC @ Z = 0
        Dz(c(Z_min,τ))/(1-S(τ)) ~ -Da*c(Z_min,τ)^(1/2)*(-ϕ(Z_min,τ))+cᵥ*Pe*c(Z_min,τ)^(3/2)*(-ϕ(Z_min,τ)) #Concentration @ Z = 0
       ] 

@named pdesys = PDESystem(eq,bcs,domains,[Z,τ],[c(Z,τ), ϕ(Z,τ)])

N = 10

dz = (Z_max-Z_min)/N

order = 2

discretization = MOLFiniteDifference([Z => dz], τ, approx_order=order, grid_align=center_align)

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

Full error:

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: Differential w.r.t. multiple variables Any[τ, Z] are not allowed.
Stacktrace:
 [1] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:171
 [2] ODESystem(deqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}}, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Any}, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, connections::Nothing, preface::Nothing, events::Vector{ModelingToolkit.SymbolicContinuousCallback}, tearing_state::Nothing, substitutions::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:120
 [3] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Num}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Num, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:174
 [4] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [6] top-level scope
   @ ~/Documents/Git/test/Scratch1.jl:143

S(τ) is a moving boundary with velocity defined as:

Dτ(S(τ)) ~ (Pe*ϕ(Z,τ)*c(Z,τ)^(1/2))/(τ*D₁)

that I am planning to combine and solve with the nerst-plank discretisation. If anyone has a better approach to this, I'm definitely open to suggestions.

Possible issue with test case 00a in MOL 1d linear convection

Hello @xtalax . I found a possible small error in test case 00a in the 1d linear convection tests.

Here is the link to the file and location.

https://github.com/SciML/MethodOfLines.jl/blob/master/test/pde_systems/MOL_1D_Linear_Convection.jl#L66

The problem is that the test case is defined for "Dt(u(t,x)) ~ Dx(u(t,x))" on line 66.
BUT, line 74 of the file defines the equations as

 eq  = Dt(u(t,x)) ~ -Dx(u(t,x))

So there seems to be an erroneus negative sign in front of the Dx term. I think this is probably just a copy and paste oversight from the first test case in the file 00, which does have a negative term. I just wanted to make sure that this was actually an error, or perhaps I just missed some code below where you removed that negative sign.

@testset "Test 00a: Dt(u(t,x)) ~ Dx(u(t,x))" begin

Moving over from a DiscreteSpace representation of the domain to DiscretizedVariable equations

Moving towards implementing a DiscretizedVariable with a smart recursive getindex and custom dict based indexing like Dict([x=>i, z=>k]) would allow for sampling whole expressions at once, and much greater flexibility, especially for mixed argument equations. Could contain also automatic interpolations and different underlying grid types per variable, helping with keeping grid types flexible and enabling staggered grids.

This would implement both the Array and Sym interfaces, allowing discrete equations that can be indexed into to give the equation at each interior point. This avoids the current need to wrap arguments to derivatives in a separate variable, and makes the rules much more comprehensible as well as allowing future extensibility i.e. #32.

Derivatives become the demarcation between different types of sampling and become a custom subtype of DiscretizedVariable, with special subtypes for Nonlinear laplacian/spherical and other types of derivatives with special handling.

There would be a step that recognises elements of an equation and replaces these with rules, creating a DiscretizedEquation and then the resulting equation is simply indexed into to generate the interior/BCs.

Is this a good idea? How much work is it to implement the Sym interface?

Plot recipes

Hello, thanks for creating this package. I would like to help with implementing plot recipes and I was thinking that having a data structure that stores / computes the grid + time span together with the solution data could help. I have worked on a package that used this idea (PICDataStructures.jl) and I was thinking to create a more general solution by re-implementing some of its parts that could benefit the SciML ecosystem more broadly.

More technically, I'm thinking that there is a need for a new type for a PDE solution which could be used for plot(pde_sol), as with the current types we plot the solution of the PDE discretization and transforming it back to something useful for users requires further processing (i.e. #82). This new type would need to be built with information about the discretization, so that we can map the discretization type to how the grid and timespan are represented in the solution, so that we know how to plot them.

Another aspect that, I believe, would be of interest would be extensibility: How easy would it be to extend the above mentioned solution to new discretization types besides MethodOfLines or how could one re-use this to plot solutions obtained from outside SciML, i.e. read a finite difference solution from a file and plot it (this was the motivation for my PICDataStructures.jl package.

One immediate question would be what plotting recipe system(s) will be used: RecipesBase or Makie/MakieCore.

Let me know what you think about this.

Allow arbitrary approximation order

At the moment approximation order is hardcoded to 2, this needs to be generalized to higher orders for the following:

Cartesian Finite Difference.
Spherical Finite Difference.
Nonlinear Cartesian Laplacian.
Nonlinear Spherical Laplacian.

Rearrange equations in to a set form pre discretization for consistency

Following from the conversation here.

Having constant coefficient advection terms on the same side as the time derivative causes the wrong scheme to be used from the perspective of if that term appeared on the other side of the equation.

Casting equations to a set form, nominally eq.lhs - eq.rhs ~ 0 pre discretization and discretizing from this perspective would ensure consistency.

PDE discretization error

Hi all, I’m currently working on implementing a relatively simple reaction-diffusion-convection model of a 2D hydrogen flame. The PDESystem is created successfully, but the discretization step using methodoflines.jl fails with an error that I do not understand:
(I'm using method of lines v0.2.0 and modeling toolkit v8.8.5)

ERROR: ArgumentError: collection must be non-empty
Stacktrace:
 [1] first
   @ ./abstractarray.jl:387 [inlined]
 [2] MethodOfLines.BoundaryHandler(bcs::Vector{Equation}, s::MethodOfLines.DiscreteSpace{2, 1, MethodOfLines.CenterAlignedGrid}, depvar_ops::Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}, tspan::Tuple{Float64, Float64}, derivweights::MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/bcs/boundary_types.jl:142
 [3] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/discretization/MOL_discretization.jl:41
 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/discretization/MOL_discretization.jl:120
 [5] top-level scope
   @ ./timing.jl:210 [inlined]
 [6] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame.jl:0

If you know what’s going on here I’d love to know. In case it’s useful, I’ve added the model definition below (the error occurs at the very last line)

## Define variables and derivative operators

# the two spatial dimensions and time
@parameters x y t
# the mass fractions MF (respectively H₂, O₂, and H₂O), temperature,
# and the source terms of the mass fractions and temperature
@variables MF[1:3](..) T(..) Smf[1:3](..) ST(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

# laplace operator
∇²(u) = Dxx(u) + Dyy(u)

# gradient operator
grad(u) = [Dx(u), Dy(u)]


## Define domains

xmin = ymin = 0.0 # cm
xmax = 1.8 # cm
ymax = 0.9 # cm
tmin = 0.0
tmax = 0.06 # s

domains = [
    x ∈ Interval(xmin, xmax),
    y ∈ Interval(ymin, ymax),
    t ∈ Interval(tmin, tmax)
]


## Define parameters

# diffusivity coefficient (for temperature and mass fractions)
κ =  2.0 #* cm^2/s
# constant (divergence-free) velocity field
U = [50, 0] #.* cm/s
# density of the mixture
ρ = 1.39e-3 #* g/cm^3
# molecular weights (respectively H₂, O₂, and H₂O)
W = [2.016, 31.9, 18] #.* g/mol
# stoichiometric coefficients
ν = [2, 1, 2]
# heat of the reaction
Q = 9800 #* K
# universal gas constant
R = 8.314472 * 100 #* 1e-2 J/mol/K -> because J = Nm = 100 N cm
# y coordinates of the inlet area on the left side of the domain
inlety = (0.3, 0.6) # mm

# TODO try out different values (systematically)
# pre-exponential factor in source term
A = 5.5e11 # dimensionless
# activation energy in source term
E = 5.5e13 * 100 # 1e-2 J/mol -> J = Nm = 100 N cm

## Define the model

# diffusion and convection terms for the mass fractions
# TODO check if you can vectorize these equations by joining them
eqs = [
    Dt( MF[i](x,y,t) ) ~ κ * ∇²( MF[i](x,y,t) ) - U' * grad(MF[i](x,y,t))
    for i in 1:3
]
# diffusion and convection terms for the temperature
push!(eqs, Dt(T(x,y,t)) ~ κ * ∇²( T(x,y,t) ) - U' * grad(T(x,y,t)))

# TODO continue here please


function atInlet(x,y,t)
    return (inlety[1] < y) * (y < inlety[2])
end

bcs = [
    #
    # initial conditions
    #

    T(x, y, 0) ~ 300.0, # K; around 26 C
    # domain is empty at the start
    MF[1](x,y,0) ~ 0.0,
    MF[2](x,y,0) ~ 0.0, 
    MF[3](x,y,0) ~ 0.0,

    #
    # boundary conditions
    #

    # left side
    T(xmin,y,t) ~ atInlet(xmin,y,t) * 950 + (1-atInlet(xmin,y,t)) *  300, # K
    MF[1](xmin,y,t) ~ atInlet(xmin,y,t) * 0.0282, # mass fraction
    MF[2](xmin,y,t) ~ atInlet(xmin,y,t) * 0.2259,
    MF[3](xmin,y,t) ~ 0.0,

    # right side
    Dt( T(xmax,y,t) ) ~ 0.0, # K/s
    Dt( MF[1](xmax,y,t) ) ~ 0.0, # mass fraction/s
    Dt( MF[2](xmax,y,t) ) ~ 0.0,
    Dt( MF[3](xmax,y,t) ) ~ 0.0,

    # top
    Dt( T(x,ymax,t) ) ~ 0.0, # K/s
    Dt( MF[1](x,ymax,t) ) ~ 0.0, # mass fraction/s
    Dt( MF[2](x,ymax,t) ) ~ 0.0,
    Dt( MF[3](x,ymax,t) ) ~ 0.0,

    # bottom
    Dt( T(x,ymin,t) ) ~ 0.0, # K/s
    Dt( MF[1](x,ymin,t) ) ~ 0.0, # mass fraction/s
    Dt( MF[2](x,ymin,t) ) ~ 0.0,
    Dt( MF[3](x,ymin,t) ) ~ 0.0,
]

@named pdesys = PDESystem(eqs, bcs, domains, [x,y,t],[T(x, y, t)])


## Discretize the system

N = 32 
dx = 1/N
dy = 1/N
order = 2

discretization = MOLFiniteDifference(
    [x=>dx, y=>dy], t, approx_order=order
)

# this creates an ODEProblem or a NonlinearProblem, depending on the system
problem = discretize(pdesys, discretization)

DAE/ODE Problem specified in MTK with domains and BCs

@xtalax As you suggested I submit the following DAE problem that I would like to reduce to an ODE system producing a dataset that can be handled by other SciML packages like DataDrivenDiffEq and others.

using ModelingToolkit, DomainSets
using LinearAlgebra
using DataDrivenDiffEq
using OrdinaryDiffEq, Sundials
using Plots; gr()

ModelingToolkit.@parameters begin
    ϕ₀ =0.0401                      # Phillips'curve parameter intercept 0.04/(1. - 0.04^2)
    ϕ₁ =6.41e-05                    # Phillips'curve parameter slope 0.04^3/(1. - 0.04^2)
    η  =0                           # -1 < η < inf 0=CD, inf=Leontief, 1=linear 
    k₀ =-0.0065                     # investment function, constant
    k₁ =-5.0                        # investment function, affine parameter
    k₂ =20                          # investment function exponent
    α  =0.025
    δ  =0.05
    β  =0.02
    r  =0.03
    cor=3.0
    fp =0.333
    b  =0.135
end

@variables begin
    t
    ω(t)  = 0.75
    λ(t)  = 0.90
    d₁(t) = 0.5
    𝛑ₑ(t)
    Y(t)  = 100.0
    𝛎(t)
    𝛟(t) 
    κ(t) 
    g(t)    
end

Dₜ = Differential(t)

domains = [t ∈ Interval(0.0,80.0),
           ω ∈ Interval(0.0,1.0),
           λ ∈ Interval(0.0,1.0),
           d₁ ∈ Interval(0.0,3.0)]


# Variable functions

𝛎 = (t, ω)     -> (1/fp)*((1 -  ω)/b)^η
𝛟 = (t, λ)     -> -ϕ₀ + ϕ₁/((1 - λ)^2)
κ = (t, 𝛑ₑ)    -> k₀ + exp(k₁ + k₂*𝛑ₑ)
g = (t, ω, 𝛑ₑ) -> κ(t, 𝛑ₑ)/𝛎(t, ω) - δ


sfcsys = [    
    Dₜ(ω)  ~ ω*(𝛟(t, λ) - α)
    Dₜ(λ)  ~ λ*(g(t, ω, 𝛑ₑ) - α - β)
    Dₜ(d₁) ~ d₁*(r - κ(t, 𝛑ₑ)/𝛎(t, ω)) + ω - 1 + κ(t, 𝛑ₑ)
    Dₜ(Y)  ~ Y*g(t, ω, 𝛑ₑ)
    0.    ~ 𝛑ₑ - 1 + ω + r*d₁
]

ModelingToolkit.@named primalKeen = ModelingToolkit.ODESystem(sfcsys)
tspan = (0.0, 80.0)

rdKeen = structural_simplify(primalKeen; simplify=true)
full_equations(rdKeen)

dt=0.25
odaeprob = ODAEProblem(rdKeen,[],tspan,jac=true,sparse=true)
odae_sol = solve(odaeprob,Tsit5(),abstol=1/10^14,reltol=1/10^14, saveat=dt)    # 11 ms

Variables in separate but overlapping domains

An important piece of the P2D battery model is to be able to simulate variables with separate but overlapping domains, corresponding to the different components of the battery. More specifically, one variable's domain is a subset of another's.
[EDIT: slightly more complex MWE requiring offsets]
For example, with domains 0<x<3 and 0<y<1 and 2<z<3 where y and z are subsets of x:

using ModelingToolkit
import IfElse

@parameters t x y z
@variables u(..) v(..) w(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Dy = Differential(y)
Dyy = Differential(y)^2
Dz = Differential(z)
Dzz = Differential(z)^2

sy = u(t,y) + v(t,y)
sz = u(t,z) + w(t,z)
sx = IfElse.ifelse(x<1,sy,IfElse.ifelse(x>2,sz,0))

eqs  = [Dt(u(t,x)) ~ Dxx(u(t,x)) + sx,
        Dt(v(t,y)) ~ Dyy(v(t,y)) + sy,
        Dt(w(t,z)) ~ Dzz(w(t,z)) + sz,]
bcs = [u(0,x) ~ 1,
       v(0,y) ~ 1,
       w(0,z) ~ 1,
       u(t,0) ~ 0.0, u(t,3) ~ 0.0,
       v(t,0) ~ 0.0, v(t,1) ~ 0.0,
       w(t,2) ~ 0.0, w(t,3) ~ 0.0]

domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,3.0),
           y ∈ Interval(0.0,1.0),
           z ∈ Interval(2.0,3.0)]

@named pdesys = PDESystem(eqs,bcs,domains,[t,x,y,z],[u(t,x),v(t,y),w(t,z)])

I believe there may have been a previous issue or PR somewhere in SciML for this, but can't find it

BoundaryConditionError: Could not read boundary condition 'Differential(x)(Differential(x)(u(t, 0))) ~ 0.0'

Trying to simulate a wave equation for a beam problem( similar to wave equation)..
encounter a error when calling discretize. Any idea what i am missing ?

BoundaryConditionError: Could not read boundary condition 'Differential(x)(Differential(x)(u(t, 0))) ~ 0.0'
throw_bc_err(::Equation) at MOL_discretization.jl:10
get_bcs(::Array{Equation,1}, ::IntervalDomain{Float64}, ::IntervalDomain{Float64}) at MOL_discretization.jl:40
discretize(::PDESystem, ::MOLFiniteDifference{Float64}) at MOL_discretization.jl:293
top-level scope at tut1.jl:41

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)             # Required in ICs and equation
Dtt = Differential(t)^2           # required in equation
Dxxxx = Differential(x)^4       # required in equation
Dxx = Differential(x)^2         # required in BCs
# some parameters
EI  = 291.6667;
m = 1.3850;
c = 0.01;
p = 1.0;
L = 5.0;

# 1D PDE and boundaru conditions
eq  = m*Dtt(u(t,x)) + c*Dt(u(t,x)) + EI*Dxxxx(u(t,x)) ~0

ic_bc = [u(0,x) ~ (p*x*(x^3 + L^3 -2*L*x^2)/(24*EI)), #for all 0 < u < L
        Dt(u(0,x)) ~ 0.,        # for all 0 < u < L
        u(t,0) ~ 0.,            # for all t > 0,, displacement zero at u=0
        u(t,5) ~ 0.,            # for all t > 0,, displacement zero at u=L
        Dxx(u(t,0)) ~ 0.,       # for all t > 0,, curvature zero at u=0
        Dxx(u(t,5)) ~ 0.]       # for all t > 0,, curvature zero at u=L

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

dt = 0.1   # dt related to saving the data.. not actual dt
# PDE sustem
pdesys = PDESystem(eq,ic_bc,domains,[t,x],[u])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference(dx,order)

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

Solving FPE with boundary condition at t=T rather than t=0

First of all, thank you for creating such a nice set of tools for solving FPEs and ODEs!

I ran into an issue when trying to solve an FPE (arising from an HJB equation) that has a boundary condition at t=T rather than at t=0. The FPE has the known solution, given in the code by trueJ, but I wanted to see how well Julia's MOL approximate this solution (as preparation for solving more complex problems). Unfortunately, I don't seem to get the following code to work:

using ModelingToolkit, MethodOfLines
import ModelingToolkit: Interval, infimum, supremum

# define the problem with indep. variables x & t and dependent variable J
@parameters x t
@variables J(..)
Dx = Differential(x)
Dt = Differential(t)
# parameters
T=5.0
α=1.0
# FPE boundaries
X=5.0
trueJ(x, t) = 0.5*x^2 / (T - t + 1/α)
# problem equations
eqs = [ Dt(J(x,t)) ~ 0.5*Dx(J(x,t))^2 ]
bcs = [ J(x,T) ~ trueJ(x, T),
        J(X,t) ~ trueJ(X, t),
        J(-X,t) ~ trueJ(-X, t)]
domains = [ t ∈ Interval(0.0, T),
            x ∈ Interval(-X, X)]
@named pde_system = PDESystem(eqs, bcs, domains, [x, t], [J(x, t)])

# discretized system
N = 4
dx = 2*X/N
order = 2
discretization = MOLFiniteDifference([x=>dx], t, approx_order=order, grid_align=center_align)
@time prob = discretize(pde_system, discretization)

Running the above yields

The system of equations is:
Equation[Differential(t)(J[2](t)) ~ 0.5((0.3J[3](t) - 0.3J[2](t))^2), Differential(t)(J[3](t)) ~ 0.5((0.3J[4](t) - 0.3J[3](t))^2), J[4](t) ~ 12.5 / (6.0 - t), J[1](t) ~ 12.5 / (6.0 - t)][Base.OneTo(4)]

Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

Term{Real, Base.ImmutableDict{DataType, Any}}[J[1](t), J[2](t), J[3](t), J[4](t)][Base.OneTo(4)]

ArgumentError: Term{Real, Base.ImmutableDict{DataType, Any}}[J[2](t), J[3](t)] are missing from the variable map.

Stacktrace:
  [1] throw_missingvars(vars::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/variables.jl:89
  [2] varmap_to_vars(varmap::Vector{Pair}, varlist::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/variables.jl:42
  [3] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:has_difference,), Tuple{Bool}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:566
  [4] (ODEProblem{true})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Float64, Float64}, parammap::SciMLBase.NullParameters; callback::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:657
  [5] (ODEProblem{true})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Float64, Float64}, parammap::SciMLBase.NullParameters) (repeats 2 times)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:656
  [6] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:635
  [7] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:635
  [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:152
  [9] top-level scope
    @ ./timing.jl:220 [inlined]
 [10] top-level scope
    @ ./In[86]:0
 [11] eval
    @ ./boot.jl:373 [inlined]
 [12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196

As you can see, I have boundary conditions on some chosen x range, as well as on t=T. That's the natural thing to do for the problem I am interested in. If I instead change the boundary condition to t=0 by using

bcs = [ J(x,0.0) ~ trueJ(x, 0.0),
        J(X,t) ~ trueJ(X, t),
        J(-X,t) ~ trueJ(-X, t)]

the code seems to work without problems. Are boundary conditions on t=T simply not supported? Or is something else going on here?

Hermite finite differences

https://www.colorado.edu/amath/sites/default/files/attached-files/2020_f_hermite-fd_ima_j_num_anal.pdf

Typos in the Brusselator PDE tutorial

Missing comma and wrong parameter name while in discretization:

"discretization = MOLFiniteDifference([x=>dx, y=>dy], t, approx_order=order grid_type=center_align)"

Which should be :

"discretization = MOLFiniteDifference([x=>dx, y=>dy], t, approx_order=order, grid_align=center_align)"

Generalize Neumann handling to higher dimensions

In SciML/DiffEqOperators.jl#349 there is a portion:

        derivars = [[dot(left_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(2),1)[1:2][2]],depvars[j][central_neighbor_idxs(CartesianIndex(2),1)[1:2][1]]]),
                    dot(right_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][1]],depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][2]]])]
                    for j in 1:length(pdesys.depvars)]

This is hardcoded to only be the inward facing derivative in a one dimensional sense because CartesianIndex(2) is going to fail against a higher dimensional point. Some starter code is:

edgeindices = [indices[e...] for e in edges]
depvarmaps = reduce(vcat,[substitute.((pdesys.depvars[i],),edgevals) .=> edgevars[i] for i in 1:length(pdesys.depvars)])
left_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][1],space[j][2]])
right_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, 0.0, [space[j][end-1],space[j][end]])
central_neighbor_idxs(i,j) = [i+CartesianIndex([ifelse(l==j,-1,0) for l in 1:length(nottime)]...),i,i+CartesianIndex([ifelse(l==j,1,0) for l in 1:length(nottime)]...)]
derivars = [[[[dot(left_weights(j),[depvars[j][central_neighbor_idxs(edgeindices[i][k],1)[1:2][2]],depvars[j][central_neighbor_idxs(edgeindices[i][k],1)[1:2][1]]]),
             dot(right_weights(j),[depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][1]],depvars[j][central_neighbor_idxs(CartesianIndex(length(space[1])-1),1)[end-1:end][2]]])]
             for j in 1:length(pdesys.depvars)] for k in 1:length(edgeindices[i])] for i in 1:length(edgeindices)]

Periodic wraparound

Periodic boundaries are currently implemented with a relation u[1] ~ u[end], and u[end] is on the interior. The boundary upwinding leads to instability on the upper index with the brusselator, so we need to move to wrapping around at the boundary so that the centered stencil is always used.

I'm doing this by creating a Val{Bool} ecapsulating whether a dimension is periodic in #45 and indexing through a function that dispatches on Val{true} and Val{false}, hopefully leading to the compiler removing this check.

Register

@JuliaRegistrator register()

2D variable connected to 1D variable via boundary

An important (eponymous) part of the P2D battery model is that a 2D variable interacts with a 1D variable on one of the boundaries of its domain. For example:

@parameters t x r 
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
Dr = Differential(r)
Drr = Differential(r)^2

s = u(t,x) + v(t,x,1)

eqs  = [Dt(u(t,x)) ~ Dxx(u(t,x)) + s,
        Dt(v(t,x,r)) ~ Drr(v(t,x,r))]
bcs = [u(0,x) ~ 1,
       v(0,x,r) ~ 1,
       Dx(u(t,0)) ~ 0.0, Dx(u(t,1)) ~ 0.0,
       Dr(v(t,x,0)) ~ 0.0, Dr(v(t,x,1)) ~ s]

domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0),
           r ∈ Interval(0.0,1.0)]

@named pdesys = PDESystem(eqs,bcs,domains,[t,x,r],[u(t,x),v(t,x,r)])

In this example, u and v are connected via s, a function of t and x only, which is a source term for u and a boundary condition for v.
Note that there is no spatial derivative for v in the x direction, so v does not take boundary conditions at x=0 and x=1.

PDE definition cannot handle simple multiplication by a constant?

this gives an error, with D₀::Float64.

	diff_eq = ∂t(c(x, t)) ~ D₀ * ∂x(
		 ∂x(c(x, t)) / (1.0 + exp(α * (c(x, t) - χ)))
		)

my hack to get it to work:

	diff_eq = ∂t(c(x, t)) ~ ∂x(
		 ∂x(c(x, t)) / (1.0/D₀ + exp(α * (c(x, t) - χ))/D₀)
		)

the error is:

MethodError: no method matching collect_ivs_from_nested_differential!(::Set{Any}, ::SymbolicUtils.Div{Real, SymbolicUtils.Term{Real, Nothing}, SymbolicUtils.Add{Real, Float64, Dict{Any, Number}, Nothing}, Nothing})

Closest candidates are:

collect_ivs_from_nested_differential!(::Any, !Matched::SymbolicUtils.Term)

Can `MethodOfLines` solve stochastic PDEs and PDEs with jumps?

Say, I was just wondering if MethodOfLines can solve stochastic PDEs and PDEs with jumps? I looked at the list of known limitations and it did say that these types of problems are excluded. Also, since MOL uses the DiffEq solvers--which support jumps and stochastic ODES, I thought that some of this might carry over to MethodOfLines.

I looked through the test cases and did not seen any examples of a stochastic PDE. I imagine I could provide the noise model in the definitition of the SDEProblem.

Before I start all this, I figured I would just ask if these types of problems are supported yet. Thanks.

discretize only generates N-1 variables for two-sided Dirichlet BC of Laplace equation

The following code (hopes to) implement a discretization of the Laplace equation with Dirichlet BCs on the left and right end of the Interval into 5 grid points, but discretize constructs a problem that only represents four points. Is that by design, due to Dirichlet BC handling when grid_align = center_align?

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x t
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

∇²(u) = Dxx(u)

x_min = 0.0
x_max = 1.0

t_min = 0.0
t_max = 11.5

α = 10.

eq = [Dt(u(x,t)) ~  α*∇²(u(x,t))]

domains = [x ∈ Interval(x_min, x_max),
           t ∈ Interval(t_min, t_max)]


bcs = [u(x,0) ~ 20,
       Dx(u(0,t)) ~ 0,
       Dx(u(x_max,t)) ~ 0] 

@named pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])

N = 5
dx = 1/N
order = 2
discretization = MOLFiniteDifference([x=>dx], t, approx_order=order, grid_align=center_align)

# despite N = 5, prob only has 4 states.
prob = discretize(pdesys,discretization)

Complex Number support

These are currently unsupported, see #80 for example code

`ExtraVariablesSystemException` when using a Robin BC

hi! solving the diffusion equation with (i) concentration-dependent diffusion coefficient and (ii) Robin BC.

I get an error:

ExtraVariablesSystemException: The system is unbalanced. There are 102 highest order derivative variables and 101 equations.

More variables than equations, here are the potential extra variable(s):

here is our PDE:

	@parameters x t

	@variables c(..)

	∂t  = Differential(t)
	∂x  = Differential(x)
	∂²x = Differential(x) ^ 2

	diff_eq = ∂t(c(x, t)) ~ ∂x((D₀ + α * c(x, t)) * ∂x(c(x, t)))

	bcs = [c(x, 0) 	 ~ 0.0,  	       # initial condition
                  R * (D₀ + α * c(0, t)) * ∂x(c(0, t)) ~ c(0, t) - cₑ, # Robin
                  ∂x(c(ℓ, t)) ~ 0.0]   # no flux

so my student found that, if we re-arrange the Robin BC to instead look like:

	   R * ∂x(c(0, t)) + (cₑ - c(0, t)) / (D₀ + α * c(0, t)) ~ 0.0,

it works! maybe something to add to the docs?

(for novices like me, would be great to have c[1](t) and c[end](t) inferred by default instead of cutting off the solution at the boundaries.) thanks!
EDIT: nevermind, it doesn't start at c[3](t). I had an additional BC in there by accident. sorry. I can figure out how to infer c[1](t) from c[2](t). I'll edit this post when I figure it out. :)

Long compile time Brusselator example

I was trying out the Brusselator example from the docs.
http://methodoflines.sciml.ai/dev/tutorials/brusselator/

It executed just fine but it took extremely long.
On first execution the discretization step took 4min, and the solve 20min.
On second execution this turned into 40sec for discretization and 50sec for the solve.

Then i tried N = 64 instead of N = 32.
This gave a discretization time of 3min. But the solve however didnt finish within 3 hours. At which point i stopped the program.

When it is taking this long the MOL.jl package is not usable for me.

Julia version: "1.6.3"
Windows 10 PRO

Full package versions:
[c3fe647b] AbstractAlgebra v0.25.1
[6e4b80f9] BenchmarkTools v1.3.1
[764a87c0] BoundaryValueDiffEq v2.7.1
[052768ef] CUDA v3.8.5
[861a8166] Combinatorics v1.0.2
[864edb3b] DataStructures v0.18.11
[39dd38d3] Dierckx v0.5.2
[9fdde737] DiffEqOperators v4.41.0
[b552c78f] DiffRules v1.10.0
[0c46a032] DifferentialEquations v7.1.0
[31c24e10] Distributions v0.25.51
[5b8099bc] DomainSets v0.5.9
[6a86dc24] FiniteDiff v2.11.0
[26cc04aa] FiniteDifferences v0.12.24
[615f187c] IfElse v0.1.1
[c8e1da08] IterTools v1.4.0
[2c470bb0] Kronecker v0.5.1
[b964fa9f] LaTeXStrings v1.3.0
[23fbe1c1] Latexify v0.15.14
[5078a376] LazyArrays v0.22.6
[c03570c3] Memoize v0.4.4
[94925ecb] MethodOfLines v0.1.0
[961ee093] ModelingToolkit v8.5.2
[8913a72c] NonlinearSolve v0.3.16
[1dea7af3] OrdinaryDiffEq v6.7.1
[91a5bcdd] Plots v1.27.1
[92933f4c] ProgressMeter v1.7.1
[1277b4bf] ShiftedArrays v1.0.0
[90137ffa] StaticArrays v1.4.4
[2913bbd2] StatsBase v0.33.16
[0999239e] SugarBLAS v0.1.0
[d1185830] SymbolicUtils v0.19.7
[0c5d862f] Symbolics v4.3.0
[64499a7a] WriteVTK v1.14.2

MOLFiniteDifference return of the full PDE variables

Right now it returns each of the pieces of the solution, but it would be nice to get the observed variable handling working so that sol[u] returns the whole (shaped?) object for u. This requires upstream array symbolics and observed handling in ModelingToolkit.jl

Support discretization when lhs/rhs equal to zero

When the equation is setup like so:

Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0

an InvalidSystemException is thrown:

ERROR: InvalidSystemException: The derivative variable must be isolated to the left-hand side of the equation like `Differential(t)(u₂(t)) ~ ...`.
 Got 0 ~ 137.50000000000009u₃(t) - (96.87499999999991u₂(t)) - (33.3333333333334u₄(t)) - (7.2916666666667105u₅(t)) - Differential(t)(u₂(t)) - (u₂(t)*(3.750000000000008u₃(t) + 0.20833333333333368u₅(t) + 0.2499999999999997(sech(5.0 + t)^2)*(3tanh(5.0 + t) - 1)*(7.5 + 7.5tanh(5.0 + t)) - (2.7083333333333344u₂(t)) - (1.250000000000007u₄(t)))) - (19.583333333333307(sech(5.0 + t)^2)*(3tanh(5.0 + t) - 1)*(7.5 + 7.5tanh(5.0 + t))).

The current work around is re-writing the equation as:

Dt(u(x,t)) ~ -u(x,t)*Dx(u(x,t)) - α*Dx2(u(x,t)) - β*Dx3(u(x,t)) - γ*Dx4(u(x,t))

Here is the full script to reproduce:

    using ModelingToolkit,DiffEqOperators,LinearAlgebra,OrdinaryDiffEq

    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0

    du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2

    bcs = [u(x,0) ~ u_analytic(x,0),
           Dx(u(-10,t)) ~ du(-10,t),
           Dx(u(10,t)) ~ du(10,t)]

    # Space and time domains
    domains = [x ∈ IntervalDomain(-10.0,10.0),
               t ∈ IntervalDomain(0.0,0.5)]
    # Discretization
    dx = 0.4; dt = 0.2

    discretization = MOLFiniteDifference([x=>dx],t;centered_order=4)
    pdesys = PDESystem(eq,bcs,domains,[x,t],[u(x,t)])
    prob = discretize(pdesys,discretization) # <-- fails here.

Error using discretize for Neumann and Robin boundary conditions, but not Dirichlet

Running examples from the heat equation tutorial and my own problems, I get

prob = discretize(pdesys,discretization)

ERROR: AssertionError: Number of offside points should not exceed the primary wind points

using Neumann and Robin BCs, but it succeeds using Dirichlet.

(also a small note: the PDE system specification line in these docs should update from
pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])
to
@named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

Combining ODE and PDE

I would like to create an advection-reaction simulation that combines LaGrangian and Eulerian frames of reference, which means combining ODEs and PDEs together in one system:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets
using Plots

@parameters x y t
@variables so4(..)
@variables xx(t) yy(t) so2_puff(t)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)

x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 7.0

N = 2

dx = (x_max-x_min)/N
dy = (y_max-y_min)/N

# Interaction between Eulerian and Lagrangian frames of reference
puff_conc(x,y,t) = IfElse(x-dx/2 < xx < x+dx/2 && y-dy/2 < yy < y+dy/2, so2_puff(t), 0.0)
@register puff_conc(x, y, t)

# Circular winds.
θ(x,y) = atan(y.-0.5, x.-0.5)
u(x,y) =  -sin(θ(x,y))
v(x,y) = cos(θ(x,y))

k=0.01 # k is reaction rate
eq = [
    # Lagrangian puff model.
    Dt(xx) ~ u(xx,yy),
    Dt(yy) ~ v(xx,yy),
    Dt(so2_puff) ~ -k*so2_puff,

    # Eulerian model.
    Dt(so4(x,y,t)) ~ -u(x,y)*Dx(so4(x,y,t)) - v(x,y)*Dy(so4(x,y,t)) + k*puff_conc(x,y,t),
]

domains = [x ∈ Interval(x_min, x_max),
              y ∈ Interval(y_min, y_max),
              t ∈ Interval(t_min, t_max)]

# Periodic BCs
bcs = [so4(x,y,t_min) ~ 0.0,
       so4(x_min,y,t) ~ so4(x_max,y,t),
       so4(x,y_min,t) ~ so4(x,y_max,t),
]

@named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[so4(x,y,t), xx, yy, so2_puff])

discretization = MOLFiniteDifference([x=>dx, y=>dy], t, approx_order=2, grid_align=center_align)
@time prob = discretize(pdesys,discretization)

I guess there might be a couple of problems with this, but the first one that shows up is this:

The system of equations is:
Equation[Differential(t)(xx(t)) ~ -sin(atan(yy(t) - 0.5, xx(t) - 0.5)), Differential(t)(yy(t)) ~ cos(atan(yy(t) - 0.5, xx(t) - 0.5)), Differential(t)(so2_puff(t)) ~ -0.01so2_puff(t), Differential(t)(so4[2, 2](t)) ~ 0.01puff_conc(0.5, 0.5, t) + 2.0so4[2, 3](t) - 2.0so4[2, 2](t), Differential(t)(so4[3, 2](t)) ~ 0.01puff_conc(1.0, 0.5, t) + 2.0so4[3, 3](t) - 2.0so4[3, 2](t), Differential(t)(so4[2, 3](t)) ~ 0.01puff_conc(0.5, 1.0, t) + 1.2246467991473532e-16so4[2, 2](t) + 2.0so4[3, 3](t) - 2.0so4[2, 3](t), Differential(t)(so4[3, 3](t)) ~ 0.01puff_conc(1.0, 1.0, t) + 1.414213562373095so4[2, 3](t) + 1.4142135623730951so4[3, 2](t) - 2.82842712474619so4[3, 3](t), so4[2, 1](t) ~ so4[2, 3](t), so4[3, 1](t) ~ so4[3, 3](t), so4[1, 2](t) ~ so4[3, 2](t), so4[1, 3](t) ~ so4[3, 3](t), so4[1, 1](t) ~ 0][Base.OneTo(12)]

Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

Term{Real, Base.ImmutableDict{DataType, Any}}[so2_puff(t), so4[1, 1](t), so4[2, 1](t), so4[3, 1](t), so4[1, 2](t), so4[2, 2](t), so4[3, 2](t), so4[1, 3](t), so4[2, 3](t), so4[3, 3](t), xx(t), yy(t)][Base.OneTo(12)]
ERROR: ArgumentError: Term{Real, Base.ImmutableDict{DataType, Any}}[xx(t), yy(t), so2_puff(t)] are missing from the variable map.
Stacktrace:
...

Thanks in advance for any information you might have regarding whether this type of simulation should be expected to work, or any suggestions for going about fixing it. Thanks!

discretize fails on Robin Boundary condition applied to Laplace equation in spherical domain

Solving a simple heat conduction problem on a sphere (1D in radial direction, assumed rotational symmetry) with symmetry boundary condition on r=0 and convective heat flux boundary condition on r=R (boundary condition of Robin type) fails with error Boundary condition Differential(r)(T(0.01, t)) ~ 100 - T(0.01, t) is not on a boundary of the domain, or is not a valid boundary condition. MethodOfLines documentation claims that Robin BCs should work.

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters r t
@variables T(..)
Dt = Differential(t)
Dr = Differential(r)
Drr = Differential(r)^2

t_min = 0.0
t_max = 11.5

R = 10e-3 # Sphere radius
T_ambient = 100 # ambient temperature


T0(r,t) = 20 # initial temperature field

# Laplace equation in spherical coordinates, describing heat conduction in sphere
eq = [Dt(T(r,t)) ~ 1/r^2 * Dr(r^2 * Dr(T(r,t))) ]


domains = [r ∈ Interval(0, R),
           t ∈ Interval(t_min, t_max)]

# BCs
bcs = [T(r,0) ~ T0(r,0), # initial values
       Dr(T(0,t)) ~ 0,   # symmetry at r = 0
       Dr(T(R,t)) ~ T_ambient - T(R,t)]  # Robin type BC describing convective heat transport on sphere surface
       


@named pdesys = PDESystem(eq,bcs,domains,[r,t],[T(r,t)])

N = 32
dr = 1/N
order = 2

discretization = MOLFiniteDifference([r=>dr], t, approx_order=order, upwind_order = 1, grid_align = center_align)


# fails on third Boundary Condition (Robin one) with error:
#  Boundary condition Differential(r)(T(0.01, t)) ~ 100 - T(0.01, t) is not on a boundary of the domain, or is not a valid boundary condition
prob = discretize(pdesys,discretization)

Plotting example for Nonlinear problems

I tried to simulate a 1D wave on a string with fixed ends and plotted the resulting animation and I'm posing the code here so that it could be refined into a tutorial.
The code is:

using ModelingToolkit
using MethodOfLines
using DomainSets
using NonlinearSolve
using Plots
# using SciMLNLSolve

@parameters t x
@variables E(..) u(..)

∂t = Differential(t)
∂tt = ∂t^2
∂xx = Differential(x)^2

c = 1

wave_eq1D = ∂tt(u(x, t)) ~ c^2 * ∂xx(u(x, t))

f(x) = sin(x)
g(x) = 0

bcs = [
    u(0, t) ~ 0,
    u(3, t) ~ 0,
    u(x, 0) ~ f(x),
    ∂t(u(x, 0)) ~ g(x)
]

t_min = 0.0
t_max = 2.0
x_min = 0.0
x_max = 3.0

# Space and time domains
domains = [
    t ∈ Interval(t_min, t_max),
    x ∈ Interval(x_min, x_max)
]

@named pde_system = PDESystem(wave_eq1D, bcs, domains, [t, x], [u(x, t)])

# Method of lines discretization
dx = 0.03
dt = 0.04

Nx = floor(Int64, (x_max - x_min) / dx) + 1

order = 2
discretization = MOLFiniteDifference([x => dx, t => dt], approx_order=2)

# Convert the PDE problem into an nonlinear problem
prob = discretize(pde_system, discretization)

sol = solve(prob, NewtonRaphson(), tol=1e-5)

anim = @animate for i in eachindex(t_min:dt:t_max)
    scatter(x_min:dx:x_max, sol[(i-1)*Nx+1:i*Nx], y_lims=(-1,1))
end

gif(anim, "wave.gif", fps=15)

The code needs a bit of cleanup and I'll also have to check if the chosen dt and dx make sense given the stability criteria. Moreover, the problem is a bit big for an example / tutorial since the string length is quite big with respect to the natural wavelength of the problem. I'll try to find more appropriate parameters over the following days.

Thanks to @xtalax for helping me with getting the above to work.

heat equation with variable material properties

In accordance with this forum entry about solving the heat equation with variable material the problem as bug entry with full code. Two approaches are included

when trying the original equation

eqs = [Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x)*Dx(T(t,x))))]

I get this:

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: MethodError: no method matching collect_ivs_from_nested_operator!(::Set{Any}, ::SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, ::Type{Differential})
Closest candidates are:
  collect_ivs_from_nested_operator!(::Any, ::Term, ::Any) at ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:178
Stacktrace:
 [1] collect_ivs_from_nested_operator!(ivs::Set{Any}, x::Term{Real, Nothing}, target_op::Type)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:182
 [2] collect_ivs(eqs::Vector{Equation}, op::Type)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:155
 [3] collect_ivs
   @ ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:149 [inlined]
 [4] check_equations(eqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/utils.jl:169
 [5] ODESystem(deqs::Vector{Equation}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}}, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Any}, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{ODESystem}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, connections::Nothing, preface::Nothing, events::Vector{ModelingToolkit.SymbolicContinuousCallback}, tearing_state::Nothing, substitutions::Nothing; checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:120
 [6] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Num}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{ODESystem}, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Num, Float64}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, checks::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/f218S/src/systems/diffeqs/odesystem.jl:174
 [7] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [9] top-level scope

With auxilary variables

eqs  = [Dt(e(t,x)) ~ Dx(q(t,x)),
        e(t,x) ~ ρ(T(t,x))*h(T(t,x)),
        q(t,x) ~ λ(T(t,x))*Dx(T(t,x))]

the error is:

ERROR: AssertionError: Variable q(t, x) does not appear in any equation, therefore cannot be solved for
Stacktrace:
 [1] build_variable_mapping(m::Matrix{Int64}, vars::Vector{Any}, pdes::Vector{Equation})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/interiormap.jl:69
 [2] MethodOfLines.InteriorMap(pdes::Vector{Equation}, boundarymap::Dict{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}, Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Vector{Any}}}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/interiormap.jl:20
 [3] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:52
 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [5] top-level scope
   @ ~/Julia/PCMfluiddynamics/heateq.jl:165

Full code:

using IfElse
using ModelingToolkit, MethodOfLines
using ModelingToolkit: Interval
using Symbolics

"""clamp function"""
function clamp(x,edge1,edge2) 
    return (x - edge1) / ( edge2 - edge1)
end

"""derivative of clamp function"""
function dclamp(x,edge1,edge2)
    return 1.0 / ( edge2 - edge1)
end

"""7-th order smoothstep function"""
function smoothstep(x;edge1=0.,edge2=1.)
    xx = clamp(x,edge1,edge2)
    return IfElse.ifelse(xx < 0., 0.,
        IfElse.ifelse(xx > 1., 1.0, -20.0*xx^7.0+70.0*xx^6.0-84.0*xx^5.0+35.0*xx^4.0))
end

"""derivative of 7-th order smoothstep function"""
function dsmoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
dx = dclamp(x,edge1,edge2)
return IfElse.ifelse(xx < 0., 0.,
    IfElse.ifelse(xx > 1., 0., (-140.0*xx^6.0+420.0*xx^5.0-420.0*xx^4.0+140.0*xx^3.0)*dx))
end

"""integral of 7-th order smoothstep function"""
function ismoothstep(x;edge1=0.,edge2=1.)
xx = clamp(x,edge1,edge2)
de = edge2 - edge1
iFct(c,j,de) = de*j^(c+1)/(c+1)
return IfElse.ifelse(xx < 0., 0.,
        IfElse.ifelse(xx > 1., (x-edge2)-20.0*de/8.0+70.0*de/7.0-84.0*de/6.0+35.0*de/5.0,
            -20.0*iFct(7.0,xx,de)+70.0*iFct(6.0,xx,de)-84.0*iFct(5.0,xx,de)+35.0*iFct(4.0,xx,de)))
end

"""Heaviside (Step) Funktion"""
function H(t)
       0.5 * (sign(t) + 1)
end

T₀ = 40+273.15
T₁ = 42+273.15
ξ(T) = smoothstep(T;edge1=T₀,edge2=T₁)
dξ(T) = dsmoothstep(T;edge1=T₀,edge2=T₁)
iξ(T) = ismoothstep(T;edge1=T₀,edge2=T₁)

"""density"""
function ρ(T;ξ=ξ)
	ρₛ = 0.86
	ρₗ(T) =  -0.000636906389320033 * T + 0.975879174796585
	return ξ(T) * ρₗ(T) + (1-ξ(T)) * ρₛ
end

"""derivitative of density"""
function dρ(T;ξ=ξ)
	ρₛ = 0.86
	ρₗ(T) =  -0.000636906389320033 * T + 0.975879174796585
	dρₗ(T) =  -0.000636906389320033
    return dξ(T) * ρₗ(T) + dρₗ(T) * ξ(T) - dξ(T) * ρₛ
end

"""thermal conductivity"""
function λ(T;ξ=ξ)
	λₗ(T) = 0.00000252159374600214 * T^2 - 0.00178215201558874 * T + 0.461994765195229
	λₛ(T) = -0.00000469854522539654 * T^2 + 0.00207972452683292 * T + 0.00000172718360302859
	return ξ(T) * λₗ(T) + (1-ξ(T)) * λₛ(T)
end

"""derivative of thermal conductivity"""
function dλ(T;ξ=ξ)
	λₗ(T) = 0.00000252159374600214 * T^2 - 0.00178215201558874 * T + 0.461994765195229
	dλₗ(T) = 2*0.00000252159374600214 * T - 0.00178215201558874
	λₛ(T) = -0.00000469854522539654 * T^2 + 0.00207972452683292 * T + 0.00000172718360302859
	dλₛ(T) = -2*0.00000469854522539654 * T + 0.00207972452683292
    return dξ(T)*λₗ(T) + dλₗ(T)*ξ(T) - dξ(T)*λₛ(T) + (1-ξ(T))*dλₛ(T)
end

"""specific enthalpie cp = const."""
function h(T;ξ=ξ,href=0.,Tref=35+273.15)
	cₗ = 2.0e3
	cₛ = 2.0e3
	Δh = 200e3
	return href + iξ(T)*cₗ+((T-Tref)-iξ(T))*cₛ+ξ(T)*Δh
end

"""derivative specific enthalpie cp = const."""
function dh(T;ξ=ξ)
	cₗ = 2.0e3
	cₛ = 2.0e3
	Δh = 200e3
    return ξ(T)*cₗ + cₛ - ξ(T)*cₛ + dξ(T)*Δh
end

# variables, parameters, and derivatives
@parameters t x

# 
@variables T(..)
# when using of auxilary variables
@variables T(..) q(..) e(..)

Dt = Differential(t)
Dx = Differential(x)
DT = Differential(T)
Dxx = Differential(x)^2

@register ρ(T)
@register λ(T)
@register h(T)
@register dρ(T)
@register dλ(T)
@register dh(T)

Symbolics.derivative(::typeof(ρ), args::NTuple{1,Any}, ::Val{1}) = dρ(args[1])
Symbolics.derivative(::typeof(λ), args::NTuple{1,Any}, ::Val{1}) = dλ(args[1])
Symbolics.derivative(::typeof(h), args::NTuple{1,Any}, ::Val{1}) = dh(args[1])

# domain edges
t_min, t_max = (0.,200.0)
x_min, x_max = (0.,5.e-3)
# discretization parameters
dx=0.1e-3
order=2

# original equation
eqs = [Dt(ρ(T(t,x))*h(T(t,x))) ~ Dx(λ(T(t,x)*Dx(T(t,x))))]
# usage of auxiliary variables
eqs  = [Dt(e(t,x)) ~ Dx(q(t,x)),
        e(t,x) ~ ρ(T(t,x))*h(T(t,x)),
        q(t,x) ~ λ(T(t,x))*Dx(T(t,x))]

Tᵢ = 273.15+30.0
Tⱼ = 273.15+50.0
# when using auxilary variables
eᵢ = ρ(Tᵢ)*h(Tᵢ)
qᵢ = 0.

# initial and boundary conditions
bcs = [T(t_min,x) ~ Tᵢ,
       T(t,x_min) ~ Tᵢ+(Tⱼ-Tᵢ)*H(t-10.0),
       Dx(T(t,x_max)) ~ 0.]
# when using auxilary variables
bcs = [T(t_min,x) ~ Tᵢ,
       T(t,x_min) ~ Tᵢ+(Tⱼ-Tᵢ)*H(t-10.0),
       Dx(T(t,x_max)) ~ 0.,
       e(t_min,x) ~ eᵢ,
       q(t_min,x) ~ 0.]

# space and time domains
domains = [t ∈ Interval(t_min, t_max),
           x ∈ Interval(x_min, x_max)]

# PDE system
@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [T(t,x)])
# when using auxilary variables
@named pdesys = PDESystem(eqs, bcs, domains, [t,x], [T(t,x),e(t,x),q(t,x)])

# method of lines discretization
discretization = MOLFiniteDifference([x=>dx], t; approx_order=order)
prob = ModelingToolkit.discretize(pdesys, discretization)
sol = solve(prob,Tsit5(),abstol=1e-8,saveat=1)

MOLFiniteDifference periodic boundaries

https://discourse.julialang.org/t/periodic-boundary-conditions-for-pde-system-using-modellingtoolkit-and-diffeqbase/51827/3

Intermediate matrix form in `MOLFiniteDifference`

The code in MOLFiniteDifference is getting pretty messy as we add more functionality for nonlinear laplacian, upwinding, etc.

One way to simplify the code would be to have an intermediate matrix form - which then gets scalarized by the compiler - instead of directly going to the scalar form by looping over all indices. Ideally this would use DerivativeOperator to stop this package diverging into two almost entirely distinct packages (the two examples in the readme). Can DerivativeOperator be used like this? If so,

can it already do nonlinear laplacian?
how can this then be scalarized?
Should boundary conditions remain as they are or use the *bc convolution?

I realize this is a little bit of a step backwards because that is closer to how it was done before the refactoring of MOLFiniteDifference, so worth discussing the pros and cons of each approach here a little bit before changing anything

Vector calculus operators for the symbolic form

(Opening a separate issue from SciML/DiffEqOperators.jl#153 SciML/DiffEqOperators.jl#146 which are about the lazy operators)

Add vector calculus operators (Laplacian, Divergence, Gradient) to supplement the Differential operator.

This would enable writing equations in generic form, with the geometry then specified later by the domains (also need to consider boundary conditions in this case). Some sample code for Poisson:

@parameters x y
@variables u(..)
Δ = Laplacian(x,y)

# 2D PDE
eq  = Δ(u(x,y)) ~ -sin(pi*x)*sin(pi*y)

One possible way to implement this is to have a @rule which replaces the vector calculus operators with the appropriate Differential operators (based on the type of the domains) and then reuse the existing MOLFiniteDifference code.
For non-Cartesian geometries, this requires at least SciML/DiffEqOperators.jl#354 , and probably a bunch of other stuff (e.g. how do you implement finite differences at r=0).

Combining ODE and PDE (part 2)

Hello,

Continuing from #83, and starting with this test:

MethodOfLines.jl/test/pde_systems/MOL_1D_Linear_Diffusion.jl

Lines 742 to 789 in 17f4cbf

 @testset "Test 13: one linear diffusion with mixed BCs, one ODE" begin 

 # Method of Manufactured Solutions 

 u_exact = (x,t) -> exp.(-t) * sin.(x) 

 v_exact = (t) -> exp.(-t) 

 @parameters t x 

 @variables u(..) v(..) 

 Dt = Differential(t) 

 Dx = Differential(x) 

 Dxx = Dx^2 

 # 1D PDE and boundary conditions 

 eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)), 

 Dt(v(t)) ~ -v(t)] 

 bcs = [u(0,x) ~ sin(x), 

 v(0) ~ 1, 

 u(t,0) ~ 0, 

 Dx(u(t,1)) ~ exp(-t) * cos(1)] 

 # Space and time domains 

 domains = [t ∈ Interval(0.0,1.0), 

 x ∈ Interval(0.0,1.0)] 

 # PDE system 

 @named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t)]) 

 # Method of lines discretization 

 l = 100 

 dx = range(0.0,1.0,length=l) 

 dx_ = dx[2]-dx[1] 

 order = 2 

 discretization = MOLFiniteDifference([x=>dx_],t) 

 # Convert the PDE problem into an ODE problem 

 prob = discretize(pdesys,discretization) 

 # Solve ODE problem 

 sol = solve(prob,Tsit5(),saveat=0.1) 

 x_sol = dx[2:end-1] 

 t_sol = sol.t 

 # Test against exact solution 

 for i in 1:length(sol) 

 @test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:length(x_sol)], atol=0.01)) 

 @test v_exact(t_sol[i]) ≈ sol.u[i][end] atol=0.01 

 end 

 end

I change the equation Dt(u(t,x)) ~ Dxx(u(t,x)) to be Dt(u(t,x)) ~ Dxx(u(t,x)) + v(t):

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters t x
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Dx^2

# 1D PDE and boundary conditions
eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)) + v(t), # This is the only line that is significantly changed from the test.
        Dt(v(t)) ~ -v(t)]

bcs = [u(0,x) ~ sin(x),
        v(0) ~ 1,
        u(t,0) ~ 0,
        Dx(u(t,1)) ~ exp(-t) * cos(1)]
domains = [t ∈ Interval(0.0,1.0),
            x ∈ Interval(0.0,1.0)]
@named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t)])
discretization = MOLFiniteDifference([x=>0.01],t)
prob = discretize(pdesys,discretization)

I was expecting that this would add a value to that equation that varied in t but was constant across all x. Instead, I get the following error: ERROR: ArgumentError: reducing over an empty collection is not allowed.

Is this a bug, or is it the expected behavior? Thanks for your help!

Second order in time

Implement recognition of du0 for second order in time PDEs

how to get the spatial discretization? + solution at spatial boundary not included right?

regarding the heat eqn examples:

can we get the spatial discretization from MOLFiniteDifference somehow instead of xs = 0:Δx:ℓ each time? to warrant this I do:

	Nₓ = length(prob.u0) + 2
	@assert (length(0:Δx:ℓ) == Nₓ)

[think this should be in the docs or explained on this page] am I right that the solution to the PDE on the points on the spatial boundary are not included in the solution (regardless of the boundary condition)? right now I manually infer these so I can do an integration over the spatial domain.

THANK YOU for this package.

Feature request: support for integral terms

Help! All the PDE models I'm interested in integrate over one of the dimensions, and NeuralPDE.jl is too slow for many of them. As MOL is the standard way to discretize these models in the literature, would it be possible to add in support for integrals?

Solving nonlinear PDE

I am trying to solve this nonlinear PDE du/dt + u du/dx = ..... but gives error that u is not defined. I am posting the code and the error here.

using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve, DifferentialEquations

# Define the system of equations
@parameters t x 
@variables u(..) 

Dt = Differential(t)
Dx = Differential(x)


tmin = 0.;
tmax = 6000.;

xmin = 0.;
xmax = 200.;


dx = 10.0;
dt = 3.0;
order = 2; # Order of the finite element approximation, must be even



eqs = [#Dt(u(t,x)) + u(t,x) * Dx(u(t,x)) ~ 0.0
        Dt(u(t,x)) + 0.5 * Dx(u(t,x)^2) ~ 0.0]    

## Note that we are defining hu(t,x) because we cannot do 
domain = [x ∈ Interval(xmin,xmax), 
          t ∈ Interval(tmin,tmax)]

# Define the BC function 

bcs = [u(tmin, x) ~ 0.,
       u(t, xmin) ~ 0.,
       Dt(u(t, xmax))+ 3 *  Dx(u(t,xmax)) ~ 0.0]


@named pdesys = PDESystem(eqs, bcs, domain, [t, x], [u(t,x)])

# Create discretization


discretization = MOLFiniteDifference([x => dx], t, approx_order = order, grid_align = center_align) # edge_align gives better results with Neumann BCs, though uses extra grid points

prob = ModelingToolkit.discretize(pdesys, discretization)

sol = NonlinearSolve.solve(prob, Tsit5())

Finite differences for Laplacian in a spherical domain

https://web.mit.edu/braatzgroup/analysis_of_finite_difference_discretization_schemes_for_diffusion_in_spheres_with_variable_diffusivity.pdf

Includes nonlinear diffusion

inplement finite difference WENO schemes

Hi,

I'm wondering if it is possible to add finite volumes for conservative PDEs in the future?

Thank you for your great work btw

Utilities tests

@testset "count differentials 1D" begin
    @parameters t x
    @variables u(..)
    Dt = Differential(t)

    Dx = Differential(x)
    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))

    Dxx = Differential(x)^2
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))

    Dxxxx = Differential(x)^4
    eq  = Dt(u(t,x)) ~ -Dxxxx(u(t,x))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 4
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
end

@testset "count differentials 2D" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))
    @test first(DiffEqOperators.differential_order(eq.rhs, x.val)) == 2
    @test first(DiffEqOperators.differential_order(eq.rhs, y.val)) == 2
    @test isempty(DiffEqOperators.differential_order(eq.rhs, t.val))
    @test first(DiffEqOperators.differential_order(eq.lhs, t.val)) == 1
    @test isempty(DiffEqOperators.differential_order(eq.lhs, x.val))
    @test isempty(DiffEqOperators.differential_order(eq.lhs, y.val))
end

@testset "count with mixed terms" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dx = Differential(x)
    Dy = Differential(y)
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + Dx(Dy(u(t,x,y)))
    @test DiffEqOperators.differential_order(eq.rhs, x.val) == Set([2, 1])
    @test DiffEqOperators.differential_order(eq.rhs, y.val) == Set([2, 1])
end

@testset "Kuramoto–Sivashinsky equation" begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0
    @test DiffEqOperators.differential_order(eq.lhs, x.val) == Set([4, 3, 2, 1])
end

It looks like those weren't copied over.

Defining an array of variables causes error with ModelingToolkit.discretize()

The code below generates an error message at ModelingToolkit.discretize()

begin
	# Dependencies
	using DiffEqOperators, MethodOfLines, OrdinaryDiffEq, ModelingToolkit, DomainSets, NonlinearSolve

	N = 6 # number of dependent variables
	
	# Variables, parameters, and derivatives
	@parameters x
	@variables u[1:N](..)
	Dx = Differential(x)
	Dxx = Differential(x)^2
	
	# Domain edges
	x_min= 0.
	x_max = 1.

	# Discretization parameters
	dx = 0.1
	order = 2
	
	# Equations
	eqs  = Vector{ModelingToolkit.Equation}(undef, N)
	for i = 1:N
		eqs[i] = Dxx(u[i](x)) ~ u[i](x)
	end

	# Initial and boundary conditions
	bcs = Vector{ModelingToolkit.Equation}(undef, 2*N)
	for i = 1:N
		bcs[i] = Dx(u[i](x_min)) ~ 0.
	end	
	
	for i = 1:N
		bcs[i+N] = u[i](x_max) ~ rand()
	end
	
	# Space and time domains
	domains = [x ∈ Interval(x_min, x_max)]
	
	# PDE system
	@named pdesys = PDESystem(eqs, bcs, domains, [x], [u[i](x) for i = 1:N])
	
	# Method of lines discretization
	discretization = MOLFiniteDifference([x=>dx], nothing, approx_order=order)
	prob = ModelingToolkit.discretize(pdesys,discretization)
	
	# # Solution of the ODE system
	sol = NonlinearSolve.solve(prob,NewtonRaphson())
end

This is the full error message

MethodError: no method matching nameof(::SymbolicUtils.Term{SymbolicUtils.FnType{Tuple, Real}, Nothing})

Closest candidates are:

nameof(!Matched::SymbolicUtils.Sym) at C:\Users\User\.julia\packages\SymbolicUtils\v2ZkM\src\types.jl:144

nameof(!Matched::ModelingToolkit.AbstractSystem) at C:\Users\User\.julia\packages\ModelingToolkit\Uaky4\src\systems\abstractsystem.jl:140

nameof(!Matched::DataType) at C:\Users\User\AppData\Local\Programs\Julia-1.7.0\share\julia\base\reflection.jl:223

...

(::MethodOfLines.var"#12#27"{Nothing})(::SymbolicUtils.Term{Real, Base.ImmutableDict{DataType, Any}})@discretize_vars.jl:43
[email protected]:47[inlined]
_collect(::Vector{Any}, ::Base.Generator{Vector{Any}, MethodOfLines.var"#12#27"{Nothing}}, ::Base.EltypeUnknown, ::Base.HasShape{1})@array.jl:744
collect_similar(::Vector{Any}, ::Base.Generator{Vector{Any}, MethodOfLines.var"#12#27"{Nothing}})@array.jl:653
map(::Function, ::Vector{Any})@abstractarray.jl:2849
MethodOfLines.DiscreteSpace(::Vector{Symbolics.VarDomainPairing}, ::Vector{Any}, ::Vector{SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}}, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})@discretize_vars.jl:41
symbolic_discretize(::ModelingToolkit.PDESystem, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})@MOL_discretization.jl:37
discretize(::ModelingToolkit.PDESystem, ::MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})@MOL_discretization.jl:120
top-level scope@Local: 45[inlined]

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Circular domains

I want to learn to solve simple PDE problems by Julia, using the finite difference method; I want to change the rectangular domain to a circular domain in the example I found, but it doesn't work.

PDE examples are found here: https://diffeqoperators.sciml.ai/stable/symbolic_tutorials/mol_heat/
Here are my various tests. Is there any information and books to get a beginner like me started? I did not find any useful information in the documentation.
Thank you everyone!!!

# Initial
domains = [x ∈ Interval(0.0, 1.0),
    y ∈ Interval(0.0, 1.0)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0.0..1.0), Symbolics.VarDomainPairing(y, 0.0..1.0)]


domains = [(x, y) ∈ UnitDisk()]
@parameters x y
domains = [x ∈ UnitCircle(), y ∈ UnitCircle()]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, UnitCircle()), Symbolics.VarDomainPairing(y, UnitCircle())]

domains = [(x, y) ∈ (0 .. 1) × (0 .. 1)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing((x, y), (0 .. 1) × (0 .. 1))]

domains = [x ∈ (0 .. 1) × (0 .. 1), y ∈ (0 .. 1) × (0 .. 1)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, (0 .. 1) × (0 .. 1)), Symbolics.VarDomainPairing(y, (0 .. 1) × (0 .. 1))]

# it works!
domains = [x ∈ (0 .. 1), y ∈ (0 .. 1)]
@show domains
# domains = Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0..1), Symbolics.VarDomainPairing(y, 0..1)]

using ModelingToolkit, DiffEqOperators, LinearAlgebra, DomainSets
using ModelingToolkit: Differential
using DifferentialEquations

@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
dx = 0.1
dy = 0.1

bcs = [u(0, y) ~ 0.0,
    u(1, y) ~ y,
    u(x, 0) ~ 0.0,
    u(x, 1) ~ x]

# Space and time domains
domains = [x ∈ Interval(0.0, 1.0),
    y ∈ Interval(0.0, 1.0)]

@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])

# Note that we pass in `nothing` for the time variable `t` here since we
# are creating a stationary problem without a dependence on time, only space.
discretization = MOLFiniteDifference([x => dx, y => dy], nothing, centered_order = 2)

prob = discretize(pdesys, discretization)
sol = solve(prob)

using Plots
xs, ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]
u_sol = reshape(sol.u, (length(xs), length(ys)))
plot(xs, ys, u_sol, linetype = :contourf, title = "solution")

One-dimensional NLS using MethodOfLines ‘type Array has no field lhs`

I’m trying to solve the one-dimensional nonlinear Schrodinger equation using MethodOfLines.jl. I looked at the code on the docs: http://methodoflines.sciml.ai/dev/tutorials/brusselator/ 3, and adapted it as:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x z
@variables A(..)
Dx   = Differential(x)
Dz   = Differential(z)
Dzz  = Differential(z)^2

xmin = 0.
xmax = 1e-1
zmax = 10.
zmin = -zmax

c0 = 1.
A0(x,z) = c0*sech(c0*z/sqrt(2))*exp(im*c0^2*x/2)

domains = [x ∈ Interval(xmin,xmax), z ∈ Interval(zmin,zmax)]

eq = [im*Dx(A(x,z))+Dzz(A(x,z)) ~ -abs2(A(x,z))*A(x,z)]
bcs = [A(xmin,z) ~ A0(xmin,z), 
        A(x,zmin) ~ 0, 
        A(x,zmax) ~ 0]

@named pdesys = PDESystem(eq,bcs,domains,[x,z],[A(x,z)])

N       = 100
dz      = 1/N
order   = 2

discretization = MOLFiniteDifference([z=>dz], x)

@time prob = discretize(pdesys,discretization)

However this won’t run, giving error type Array has no field lhs. I think the only main difference between my code and the tutorial is that my PDE is for a single function.

I'm running Julia v1.7.2, with package versions:

DiffEqOperators v4.41.0
MethodOfLines v0.2.0
~~ModelingToolkit v6.7.1~~ ModelingToolkit v8.5.5
DomainSets v0.5.9

The full error is:

ERROR: type Array has no field lhs
Stacktrace:
  [1] getproperty
    @ .\Base.jl:42 [inlined]
  [2] (::MethodOfLines.var"#65#67"{Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}})(x::Vector{Equation})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\MOL_utils.jl:64
  [3] mapreduce_first(f::MethodOfLines.var"#65#67"{Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}}, op::Function, x::Vector{Equation})
    @ Base .\reduce.jl:394
  [4] _mapreduce(f::MethodOfLines.var"#65#67"{Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}}, op::typeof(union), #unused#::IndexLinear, A::Vector{Vector{Equation}})
    @ Base .\reduce.jl:405
  [5] _mapreduce_dim(f::Function, op::Function, #unused#::Base._InitialValue, A::Vector{Vector{Equation}}, #unused#::Colon)
    @ Base .\reducedim.jl:330
  [6] #mapreduce#725
    @ .\reducedim.jl:322 [inlined]
  [7] mapreduce(f::Function, op::Function, A::Vector{Vector{Equation}})
    @ Base .\reducedim.jl:322
  [8] get_all_depvars(pdesys::PDESystem, depvar_ops::Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\MOL_utils.jl:64
  [9] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\discretization\MOL_discretization.jl:18
 [10] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines C:\Users\ruvil\.julia\packages\MethodOfLines\9MLPx\src\discretization\MOL_discretization.jl:146
 [11] top-level scope
    @ .\timing.jl:220

SumnationbyParts finite differences

We probably want sumnation by parts: See DEO #506

This could be implemented relatively quickly by leveraging this package

PDE discretization fails when a parameter is used

Hi all, I've working on a 2D hydrogen flame model that I've described in #59. I've tried to change the diffusivity parameter kappa by using @parameters x y t kappa and defining the PDESystem using

@named pdesys = PDESystem(
    eqs, # partial differential equations
    bcs, # initial/boundary conditions
    domains, # domain of the independent variables (i.e. spatial/time)
    [x,y,t], # independent variables
    [YF(x,y,t), YO(x,y,t), YP(x,y,t), T(x,y,t)], # dependent variables
    [κ], # parameters
    )

The error that I get is this:

Discretization failed at structural_simplify, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: Dict(kv): kv needs to be an iterator of tuples or pairs
Stacktrace:
 [1] Dict(kv::Vector{Any})
   @ Base ./dict.jl:132
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [4] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame2.jl:142

caused by: BoundsError: attempt to access Num at index [2]
Stacktrace:
 [1] indexed_iterate(I::Num, i::Int64, state::Nothing)
   @ Base ./tuple.jl:95
 [2] grow_to!(dest::Dict{Num, Float64}, itr::Vector{Any}, st::Int64)
   @ Base ./dict.jl:153
 [3] grow_to!(dest::Dict{Any, Any}, itr::Vector{Any})
   @ Base ./dict.jl:145
 [4] dict_with_eltype
   @ ./abstractdict.jl:539 [inlined]
 [5] Dict(kv::Vector{Any})
   @ Base ./dict.jl:129
 [6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:124
 [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [8] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame2.jl:142

	@testset "Test 13: one linear diffusion with mixed BCs, one ODE" begin
	# Method of Manufactured Solutions
	u_exact = (x,t) -> exp.(-t) * sin.(x)
	v_exact = (t) -> exp.(-t)

	@parameters t x
	@variables u(..) v(..)
	Dt = Differential(t)
	Dx = Differential(x)
	Dxx = Dx^2

	# 1D PDE and boundary conditions
	eqs = [Dt(u(t,x)) ~ Dxx(u(t,x)),
	Dt(v(t)) ~ -v(t)]
	bcs = [u(0,x) ~ sin(x),
	v(0) ~ 1,
	u(t,0) ~ 0,
	Dx(u(t,1)) ~ exp(-t) * cos(1)]

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

	# PDE system
	@named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t)])

	# Method of lines discretization
	l = 100
	dx = range(0.0,1.0,length=l)
	dx_ = dx[2]-dx[1]
	order = 2
	discretization = MOLFiniteDifference([x=>dx_],t)

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

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

	x_sol = dx[2:end-1]
	t_sol = sol.t

	# Test against exact solution
	for i in 1:length(sol)
	@test all(isapprox.(u_exact(x_sol, t_sol[i]), sol.u[i][1:length(x_sol)], atol=0.01))
	@test v_exact(t_sol[i]) ≈ sol.u[i][end] atol=0.01
	end
	end

sciml / methodoflines.jl Goto Github PK

methodoflines.jl's People

Contributors

Stargazers

Watchers

Forkers

methodoflines.jl's Issues

Recommend Projects

Recommend Topics

Recommend Org