4
$\begingroup$

I have the following code to solve the wave equation in 2D:

range=2;
trange=5;
region=Region[Rectangle[range{-1,-1},range{1,1}]];
wave2D = NDSolveValue[
      {D[u[t,x,y] ,{t,2}]- 
       Laplacian[u[t, x, y], {x, y}] == 0, 
       u[0, x, y] == E^(-5 (x^2 + y^2)), 
       Derivative[1, 0, 0][u][0, x, y] == 0
       (* absorbing condition *)
   },
   u, {t, 0, trange}, {x, y} \[Element] region] // Quiet;

I want to add an absorbing boundary condition: $\partial_t u=\partial_x u$ on the boundary and I think I have to use NeumannValue somehow but I don't understand the documentation of NeumannValue. Can someone show me how to do that and explain the idea behind NeumannValue in simple terms?

Edit: with the help of @Ulrich Neumann (perfect name for this question) I was able to make the following plot

plot = Table[
   Plot3D[wave2D[t, x, y], {x, -range, range}, {y, -range, range}, 
    PlotRange -> {-.5, 1}, ImageSize -> 300], {t, 0, 3.8, .075}];
ListAnimate[plot]

enter image description here

Improve this question
$\endgroup$
6

2 Answers 2

Reset to default
4
$\begingroup$

With NeumannValue the absorbing boundary conditions might be formulated for x==2and x==-2. Roughly NeumannValue describes the flux perpendicular to the boundary:

    range = 2;
trange = 5;
region = Region[Rectangle[range {-1, -1}, range {1, 1}]];
wave2D = NDSolveValue[{D[u[t, x, y], {t, 2}] - 
     Laplacian[u[t, x, y], {x, y}] == 
    NeumannValue[-Derivative[1, 0, 0][u][t, x, y], x == 2]+NeumannValue[ -Derivative[1, 0, 0][u][t, x, y], x == -2], 
   u[0, x, y] == E^(-5 (x^2 + y^2)), 
   Derivative[1, 0, 0][u][0, x, y] == 0
   (*absorbing condition*)}, 
  u, {t, 0, trange}, {x, y} \[Element] region] 

Manipulate[Plot3D[wave2D[t, x, y], Element[{x, y}, region],PlotRange -> {-1, 1}], {{t, 2 }, 0, 5, Appearance -> "Labeled"}]

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Very nice! I want it along all boundaries but that works by setting the predicate to true. What I still think is weird is that you don't need the term laplacian u == 0 anymore, does NeumannValue evaluate to zero for the region that is not on the boundary? $\endgroup$
    – AccidentalTaylorExpansion
    Apr 30, 2022 at 14:29
  • $\begingroup$ @AccidentalTaylorExpansion Laplacian[u[t, x, y], {x, y}]is still part of the pde! $\endgroup$
    – Ulrich Neumann
    Apr 30, 2022 at 14:31
  • $\begingroup$ I meant to say $\partial_{tt}u-\nabla^2 u=0$ does not occur on its own anymore, only in conjunction with NeumannValue, which would imply for me that NeumannValue evaluates to zero inside the bulk and to something else (in this case $\partial_t u$ on the boundary). It's still a bit confusing how it works but I guess it works now $\endgroup$
    – AccidentalTaylorExpansion
    Apr 30, 2022 at 15:06
2
$\begingroup$

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]

enter image description here

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]

enter image description here

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ (+1) do you have suggestions for how to improve the documentation to get rid of the 'probaby hard to follow'? $\endgroup$
    – user21
    May 3, 2022 at 16:53
  • $\begingroup$ @user21 One thing I can think of is, perhaps the document should emphasize a bit more about the importance of formal PDE/Inactive PDE. Adding example showing how to determine the NeumannValue with the help of GetInactivePDE (for certain strange PDE) may improve the readability of document, too. As we've seen, checking the corresponding formal PDE is almost unavoidable for setting NeumannValue correctly. $\endgroup$
    – xzczd
    May 3, 2022 at 17:04
  • $\begingroup$ Thanks, I added some more information to the documentation. I feel that at the end of the day one has to think through the derivation of the FEM like it is hinted add in the NeumannValue ref page or the FEM tutorials to really understand the meaning of NeumannValue and how to use that. I am not sure there is really a short cut. But again, I tried to add more information to the documentation and made GetInactivePDE more prominent in conjunction with the sections on understanding PDE parsing and NeumannValue. $\endgroup$
    – user21
    May 4, 2022 at 5:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.