As mentioned in the comment above, the problem is related with (and essentially equivalent to, in my view) this one, but let me elaborate with an answer anyway.
The explanation for the usage of NeumannValue
can be found in Details section of document of NeumannValue
, and the tutorial Finite Element Method Usage Tips (Particularly Formal Partial Differential Equations and NeumannValue
and Formal Partial Differential Equations section) but maybeprobably a bit hard to follow. In short, the NeumannValue
is defined based on a so-called formal PDE, which can be checked with the new-in-12.2 NDSolve`FEM`GetInactivePDE
. (If you're not yet in v12.2, try the function in this post. )
Let's check the formal PDE corresponding to the 2D wave equation:
<< NDSolve`FEM`;
{eqLHS, ic} =
With[{u = u[t, x, y]}, {D[u, {t, 2}] -
Laplacian[u, {x, y}], {u == E^(-5 (x^2 + y^2)), D[u, t] == 0} /. t -> 0}];
GetInactivePDE@
First@NDSolve`ProcessEquations[{eqLHS == 0, ic},
u, {t, 0, trange}, {x, y} ∈ region]

As mentioned in Details section of document of NeumannValue
:
$∇·(-c ∇u-α u+γ)+…=f+\text{NeumannValue}[g-q u,\text{pred}]$ is used to specify the flux over the part of the boundary $∂Ω$ where pred is true, such that $\overset{\rightharpoonup }{n}·(c ∇u+α u-γ)=g-q u$ holds.
Clearly, in our case $-c ∇u-α u+γ=$-{{1, 0}, {0, 1}}.Inactive[Grad][u[t, x, y], {x, y}]
.
Next step is to determine the $g-q u$. For $x=2$, the needed 1st order ABC is $\frac{\partial u}{\partial t}+\frac{\partial u}{\partial x}=0$, $\overset{\rightharpoonup }{n}=(1,0)$, so
{1, 0} . (-(-{{1, 0}, {0, 1}} . Inactive[Grad][u[t, x, y], {x, y}])) == gMinuspu
should be equivalent to
With[{u = u[t, x, y]}, D[u, t] + D[u, x] == 0]
at $x=2$. Then it's clear gminuspu
i.e. $g-pu$ should be -D[u[t, x, y], t]
at $x=2$.
Similarly we'll find at all the boundary of the rectangle region, the needed $g-pu$ terms are all -D[u[t, x, y], t]
. So the needed NeumannValue
is
sol = NDSolveValue[{eqLHS == NeumannValue[-D[u[t, x, y], t], True], ic},
u, {t, 0, trange}, {x, y} ∈ region];
plot = Table[
DensityPlot[sol[t, x, y], {x, -range, range}, {y, -range, range},
PlotRange -> {-.5, 1}, ImageSize -> 300, PlotPoints -> 50], {t, 0, 3.8, .075}];
ListAnimate[plot]

From the animation we can see the reflection is still somewhat obvious, this is expected because we're using 1st order ABC. To reduce the reflection, you may consider 2nd order ABC as shown here. Setting up PML as shown in the tutorial Acoustics in the Time Domain is also a possible choice.