Concrete Operators in DiffEqOperators.jl for time-indepedent problems

We are looking into the use of DiffEqOperators for time-independent problems. We are considering the Poisson equation on the unit interval and unit square.

First a 1D implementation was made to figure out how to use it.

DiffEqOperators provides several functions to construct finite difference approximations of a differential equation. The ‘matrix’ Ah is built up in two steps:

  1. Construct the finite difference approximation using CenteredDifference or UpwindDifference
  2. Construct the boundary conditions using one of several available functions for different types (e.g. DirichletBC for the Dirichlet boundary condition).
    The complete Ah ‘matrix’ is then constructed by multiplying the result of these two steps. For example,
    A = CenteredDifference(ord, ord_approx, h, N)
    bc = Dirichlet0BC(Float64)
    Ah = A * bc
    The word ‘matrix’ is between quotation marks because DiffEqOperators does not really return a matrix but an abstract operator. This operator can then be concretized when required. DiffEqOperators also provides an implementation of \ for solving the differential equation constructed using the operators.

u = (A * bc) \ f;
u = vcat(0, u, 0)

Next the same was tried in 2D

Dxx = CenteredDifference{1}(ord, ord_approx, h, N); Dyy = CenteredDifference{2}(ord, ord_approx, h, N); bc = Dirichlet0BC(Float64)
Qx,Qy = MultiDimBC(bc, size(F))
Axx = Dxx * Qx; Ayy = Dyy * Qy; A =Axx+Ayy;
u = A \ f;

This, however, resulted in the error message

     MethodError: Cannot `convert` an object of type
       GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, Stati
     cArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64}, MultiDimDirectionalBC{Float64,
     RobinBC{Float64, Vector{Float64}}, 1, 2, 1}} to an object of type

Please advice. Thx!

1 Like

Your post has just revealed a long standing issue with concretization in multiple dimensions that we will fix asap so thank you for that. In the meantime, can you share the system you are trying to solve and I’ll find a way to solve it with the current infrastructure?

Thx! We are trying to solve Laplacian(u) = f on the domain 0 < x,y < 1. Nothing fancy. We have a home-brewed implementation that we can share. Please let us know in case we can contribute. We will be happy to test drive your solution. Thx again.

Hi, in the meantime I suggest using our other package MethodOfLines, which uses Symbolics.jl and DiffEqOperators.jl under the hood and existing ODE/Nonlinear solvers. Here’s an example, please ask away if you have any questions. Note that this does not currently scale well to very large problems, but there are improvements in the works that will hopefully make this more reasonable.

using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve

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

Dxx = Differential(x)^2
Dyy = Differential(y)^2

xmin = 0.
xmax = 1.
ymin = 0.
ymax = 1.

f(x,y) = x^2 + y^2 #If you have array data, define a piecewise function or use an interpolation

eq = Dxx(u(x,y)) + Dyy(u(x,y)) ~ f(x,y)

domain = [x ∈ Interval(xmin,xmax), 
          y ∈ Interval(ymin,ymax)]

# Arbitrary bcs can be defined in a natural way
# Dirichlet 0 would be u(xmin,y) ~ 0.
# Neumann would be Differential(x)(u(xmax,y)) ~ α(x, y) where α is a function, can be replaced with a constant value.
# Periodic would be u(xmin,y) ~ u(xmax,y)
bcs = [u(xmin, y) ~ u(xmax, y),
       u(x, ymin) ~ u(x, ymax)]

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

# Create discretization

dx = 0.1
dy = 0.1

order = 2 # Order of the finite element approximation, must be even

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

prob = discretize(pdesys, discretization)

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

# Center aligned grid
disc_x = xmin:dx:xmax
disc_y = ymin:dy:ymax

# Edge aligned grid
# disc_x_e = xmin-dx/2:dx:xmax+dx/2
# disc_y_e = ymin-dy/2:dy:ymax+dy/2

sol_u = reshape(sol.u, (length(disc_x), length(disc_y)))

# Plot the solution
using Plots
heatmap(disc_x, disc_y, sol_u)

Thx. I currently experience a difficulty installing MethodofLines. I will try to resolve it.

What’s the issue? If it comes up for you chances are it will come up for others.

Not sure. I fail to understand the output of Pkg (see below).

julia> import Pkg; Pkg.add(“MethodOfLines”)
Resolving package versions…
ERROR: Unsatisfiable requirements detected for package Symbolics [0c5d862f]:
Symbolics [0c5d862f] log:
├─possible versions are: 0.1.0-4.3.0 or uninstalled
├─restricted to versions * by an explicit requirement, leaving only versions 0.1.0-4.3.0
├─restricted by compatibility requirements with MethodOfLines [94925ecb] to versions: 4.0.0-4.3.0
│ └─MethodOfLines [94925ecb] log:
│ ├─possible versions are: 0.1.0 or uninstalled
│ └─restricted to versions * by an explicit requirement, leaving only versions 0.1.0
└─restricted by compatibility requirements with StaticArrays [90137ffa] to versions: 0.1.0-0.1.32 or uninstalled — no versions left
└─StaticArrays [90137ffa] log:
├─possible versions are: 0.8.0-1.4.1 or uninstalled
├─restricted to versions * by an explicit requirement, leaving only versions 0.8.0-1.4.1
└─restricted by compatibility requirements with ReactionMechanismSimulator [c2d78dd2] to versions: 0.8.0-0.12.5
└─ReactionMechanismSimulator [c2d78dd2] log:
├─possible versions are: 0.1.0-0.4.0 or uninstalled
└─restricted to versions * by an explicit requirement, leaving only versions 0.1.0-0.4.0

It is telling you that StaticArrays is upper bounded to version 0.12 by ReactionMechanismSimulator, but such an ancient version of StaticArrays means nly very old verions of Symbolics can be installed. MethodOfLines requires Symbolics version 4.0 or up, while the last version of Symbolics to work with StaticArrays version 0.12 is 0.1.32.

In general ReactionMechanismSimulator appears to depend on a lot of old package versions, which will make it tricky to use this package together with other packages. I do see some activity on their repo though, so maybe they’re in the process of updating their dependencies.

@ziolai I noticed that you were interested in the Brusselator equation, you should be aware that I’m currently working on allowing this equation in MethodOfLines.

What is your application, if I can ask?

@nilshg : I removed ReactionMechanismSimulator and now have MethodsOfLines up and running.

@xtalax : my objective it to introduce solving reaction-diffusion equations in a master course that I am teaching.

Thx both!

1 Like