SciML / MethodOfLines.jl

Automatic Finite Difference PDE solving with Julia SciML
https://docs.sciml.ai/MethodOfLines/stable/
MIT License
151 stars 27 forks source link

Absorbing Boundary Conditions #235

Open gladisor opened 1 year ago

gladisor commented 1 year ago

Dear MOL Developers,

I am simulating acoustic waves in 1 and 2 dimensions. My current goal is to try to simulate an open domain where waves which leave the domain do not return. My understanding is that this can be done with an absorbing boundary condition or a Perfectly Matched Layer.

In the 1 dimension case I have been successful in implementing this boundary condition.

`

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, IfElse using ModelingToolkit.Symbolics: CallWithMetadata using OrdinaryDiffEq: ODEIntegrator, ODESolution import GLMakie

@parameters speed @parameters t @variables u(..)

@parameters x [bounds = (x_min, x_max)] Dxx = Differential(x)^2 Dtt = Differential(t)^2 C(x, t) = (x, t) -> speed eq = Dtt(u(x, t)) ~ C(x, t) ^ 2 * Dxx(u(x, t)) x_min, x_max = getbounds(x)

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

bcs = [ u(x, 0.0) ~ exp(-x^2), Dt(u(x_min, t)) - speed Dx(u(x_min, t)) ~ 0., Dt(u(x_max, t)) + speed Dx(u(x_max, t)) ~ 0., Dt(u(x, 0.0)) ~ 0. ]

@named sys = PDESystem(eq, bcs, [x in (x_min, x_max), t in (0.0, 10.0)], [x, t], [u(x, t)]) disc = MOLFiniteDifference([x => 30], t)

prob = discretize(sys, disc) grid = get_discrete(sys, disc) iter = init(prob, Tsit5(), advance_to_tstop = true, saveat = 0.05)

step!(iter) `

Animating the resulting solution shows the correct behavior of the wave leaving the domain.

wave_1d

However when I attempt to extend this behavior to the two dimensional case by including a second spacial variable (y) there is an error message.

` @parameters speed @parameters t @variables u(..)

@parameters x [bounds = (x_min, x_max)] @parameters y [bounds = (y_min, y_max)]

Dxx = Differential(x)^2 Dyy = Differential(y)^2 Dtt = Differential(t)^2 C(x, y, t) = (x, y, t) -> speed eq = Dtt(u(x, y, t)) ~ C(x, y, t) ^ 2 * (Dxx(u(x, y, t)) + Dyy(u(x, y,t)) x_min, x_max = getbounds(x) y_min, y_max = getbounds(y)

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

bcs = [ u(x, y, 0.0) ~ exp(-(x^2 + y^2)), Dt(u(x_min, y, t)) - wave.speed Dx(u(x_min, y, t)) ~ 0., Dt(u(x_max, y, t)) + wave.speed Dx(u(x_max, y, t)) ~ 0., Dt(u(x, y_min, t)) - wave.speed Dy(u(x, y_min, t)) ~ 0., Dt(u(x, y_max, t)) + wave.speed Dy(u(x, y_max, t)) ~ 0., Dt(u(x, y, 0.0)) ~ 0. ]

@named sys = PDESystem(eq, bcs, [x in (x_min, x_max), y in (y_min, y_max), t in (0.0, 10.0)], [x, y, t], [u(x, y, t)]) disc = MOLFiniteDifference([x => 30, y => 30], t)

prob = discretize(sys, disc) grid = get_discrete(sys, disc) iter = init(prob, Tsit5(), advance_to_tstop = true, saveat = 0.05)

step!(iter) `

The error message

ERROR: LoadError: ArgumentError: Equations (1680), states (1664), and initial conditions (1664) are of different lengths. To allow a different number of equations than states use kwarg check_length=false.

Please let me know any advice you have to resolve this problem.

Thanks

gladisor commented 1 year ago

I just lowered the discretization threshold to 20 from 30 and it works! This issue must be related to #187