DiffEqOperators.jl有限差分法求解PDE的域定义问题

我想尝试用Julia求解偏微分方程,有限差分法;

例子是在这里找到的: https://diffeqoperators.sciml.ai/stable/symbolic_tutorials/mol_heat/
我想把原文中的矩形域换成其他的形式,比如圆形域,因为我是初学者,搞来搞去也不知道怎么弄,从官方的文档里面也没有找到有用的信息?
域的定义和DomainSets.jl应该怎样设置?有Julia处理偏微分方程相关的教程或资料推荐吗?
我想读这个包的源码,查看Julia包源码的正确姿势应该是什么?谢谢诸位!

# 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")

域设置参照 DomainSets.jl 的文档,看起来圆形域有自己的形式:
Circles and Spheres - JuliaApproximation/DomainSets.jl

可以去英文社区搜搜。搜到了给大家分享一下。

  • 确定具体的方法
  • 直接 @edit
  • 单步调试跟踪
julia> /(rand(2,2), rand(2,2))
2×2 Matrix{Float64}:
 1.47324  -1.00924
 1.59493  -1.35406

julia> methods(/, (Matrix{Float64},Matrix{Float64}))
# 1 method for generic function "/":
[1] /(A::AbstractVecOrMat, B::AbstractVecOrMat) in LinearAlgebra at C:\Users\woclass\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\LinearAlgebra\src\generic.jl:1148


julia> @edit /(rand(2,2), rand(2,2))
# 这里会打开编辑器


julia> using Debugger

julia> @enter /(rand(2,2), rand(2,2))
In /(A, B) at C:\Users\woclass\AppData\Local\Programs\Julia-1.7.1\share\julia\stdlib\v1.7\LinearAlgebra\src\generic.jl:1148
 1148  function (/)(A::AbstractVecOrMat, B::AbstractVecOrMat)
>1149      size(A,2) != size(B,2) && throw(DimensionMismatch("Both inputs should have the same number of columns"))
 1150      return copy(adjoint(adjoint(B) \ adjoint(A)))
 1151  end

About to run: (size)([0.39390451265431703 0.5561702006589889; 0.09287387055408114 0.4244147302394803], 2)
1|debug> q

谢谢,我在ModelingToolkit.jl的issue下也问了这个问题,回答说:This isn’t implemented for this library for now.