Solve a PDE on a domain $\Omega$ with given boundary conditions

Question

I'm starting to study the behavior of some PDEs and I would like to run simulations in mathematica to help me visualize solutions.

For example, a prime example that I would like to study is $$ \left\{ \begin{array}{cc} -\Delta\phi + u\cdot\nabla\phi=1 & \text{in }\Omega\\ \phi=0& \text{on }\partial\Omega \end{array}\right. $$ for a fixed (smooth enough) $u$ which is divergence free and flows tangentially to the boundary.

For one dimensional problems, inputting boundary conditions is easy since there are typically only two bounary points at most. But I would like to study this PDE (in 2-d) with a variety of different $\Omega$'s. I can't find anything in the mathematica documentation about how to input arbitrary boundary conditions.

To make things simpler, I will restrict myself to reasonably nice $\Omega$'s which have nice boundaries. Suppose that $\partial\Omega$ is given implicitly by $F(x,y)=0$, e.g. $x^2+y^2-1=0$. What would be the easiest/best/most efficient way to input the boundary condition $\phi(x,y)=0$ when $F(x,y)=0$?

Edit: Feel free to use any reasonable discretization technique, or to assume that you have a function like isInRegion[x,y] if you need it.

Answer 1

Define your region $\Omega$ and its boundary $\partial\Omega$ by $$\left\{ \begin{array}{cc} f(x,y)>0 & \text{in }\Omega\\ f(x,y)=0& \text{on }\partial\Omega \end{array}\right.$$ And define another function $g(x,y)$ so that the ordered pair $(f(x,y),g(x,y))$ is invertible, i.e. you can write $x=x(f,g)$ and $y=y(f,g)$.

In the $x,y$ coordinate system, the level curves of $f$ form something like concentric circles, though irregularly shaped. So $f$ is something like a radial coordinate. If you choose $g$ to be the streamlines of the gradient of $f$, then the curves of $g$ are orthogonal to the curves of $f$, like the radial rays of constant angle in polar coordinates.

In the $f,g$ coordinate systems, the boundary conditions are rectangular, $0<f<F$ and $0<g<2\pi$ with periodic boundary conditions for the $g$ dimension: $\phi(f,0)=\phi(f,2\pi)$ and your given boundary condition for the $f$ dimension: $\phi(0,g)=0$ along with the physical condition $\phi(F,g)<\infty$.

You can look up expressions for the laplacian and directional derivative in general curvilinear coordinates, or derive it yourself. Your PDE will now look really complicated, but the boundaries are rectangular, so NDSolve should be able to handle them.

EDIT: The only really hard part is coming up with the coordinate transformations. You may find acceptable expressions using Assumptions with FullSimplify on the results of Solve. Here is an example

(* our new coordinates *)
c = {\[ScriptF], ℊ}; 

The boundary defines one family of curves. $F(x,y)=f$

(* Begin example *)
\[GothicF][x_, y_] := 1 - (x^2 + y^2)

The other family of coordinate curves is arbitrary, but needs to be compatible. $G(x,y)=g$

\[GothicG][x_, y_] := ArcTan[x, y]

Try to let the computer do the hard work of inverting the two equations above. $x=X(f,g)$ and $y=Y(f,g)$

{x[\[ScriptF]_, ℊ_], y[\[ScriptF]_, ℊ_]} = {x, y} /. 
  Assuming[(0 < \[ScriptF] < 1) && 
           (0 <= ℊ < 2*Pi), 
       FullSimplify@First@Solve[{
        \[GothicF][x, y] == \[ScriptF],
        \[GothicG][x, y] == ℊ}, {x, y}]]
(* End example *)

Once you've decided on your coordinates and the transformation, the rest is tedious and algorithmic. First we find the coefficients of the metric tensor and its inverse.

$$ds^2=g_{\text{ij}}dx^idx^j$$

dsSquared = Expand[(#1 . #1 & )[Dt[{x @@ c, y @@ c}]]]; 
g = FullSimplify@Table[
  (1/2 + 1/2 KroneckerDelta[i, j])Coefficient[dsSquared, Dt[c[[i]]]*Dt[c[[j]]]],
  {i, Length[c]}, {j, Length[c]}];
gI = Inverse[g];

Then we are ready to calculate the Laplace-Beltrami operator

$$\Delta \phi =\frac{1}{\sqrt{|g|}}\frac{\partial }{\partial x^i}\left(g^{\text{jk}}\sqrt{|g|}\frac{\partial \phi }{\partial x^k}\right)$$

Laplacian[ϕ_, c_, g_] := 
  FullSimplify@Expand@Sum[
    (1/Sqrt[Abs[Det[g]]])*D[Sqrt[Abs[Det[g]]]*gI[[j, k]]*D[ϕ @@ c, c[[k]]], c[[j]]],
    {j, Length[c]}, {k, Length[c]}]

To handle the directional derivative, we first need to find the components of the vector in our new coordinate system. $$\overset{\rightharpoonup }{u}=u^x\overset{\rightharpoonup }{e}_x+u^x\overset{\rightharpoonup }{e}_y=u^f\overset{\rightharpoonup }{e}_f+u^g\overset{\rightharpoonup }{e}_g$$ The basis vectors $\overset{\rightharpoonup }{e}_f,\overset{\rightharpoonup }{e}_g$

eSub[1] = {D[x @@ c, c[[1]]], D[y @@ c, c[[1]]]}; 
eSub[2] = {D[x @@ c, c[[2]]], D[y @@ c, c[[2]]]}; 

The basis covectors $\overset{\rightharpoonup }{e}^f,\overset{\rightharpoonup }{e}^g$

eSup[1] = (eSub[1]*g[[2, 2]] - eSub[2]*g[[2, 1]])/Det[g]; 
eSup[2] = (g[[1, 1]]*eSub[2] - g[[1, 2]]*eSub[1])/Det[g]; 

The components of the vector are then $$u^f=\overset{\rightharpoonup }{e}^f\cdot \overset{\rightharpoonup }{u}=\overset{\rightharpoonup }{e}^f\cdot \overset{\rightharpoonup }{e}_xu^x+\overset{\rightharpoonup }{e}^f\cdot \overset{\rightharpoonup }{e}_yu^y$$

{uf, ug} = {{eSup[1] . {1, 0}, eSup[1] . {0, 1}},
            {eSup[2] . {1, 0}, eSup[2] . {0, 1}}} . {ux, uy}; 

And the directional derivative is $$\overset{\rightharpoonup }{u}\cdot \nabla \phi =u^f\frac{\partial }{\partial f}\phi +u^g\frac{\partial }{\partial g}\phi$$

DirectionalDerivative[ϕ_, c_, u_] := Sum[
  u[[i]]*D[ϕ @@ c, c[[i]]],
  {i, Length[c]}]

Finally, we construct your differential equation

eqn = -Laplacian[ϕ, c, g] + DirectionalDerivative[ϕ, c, {uf, ug}] == 1