Coefficients sign of InitializePDECoefficients for non-linear PDE

Question

I was trying to solve a nonlinear PDE

$-\Delta u(x,y) = u(x,y)^2$

in the domain $[0,1] \times[0,1]$ with Dirichlet boundary condition $u(x,y) =1$ on the boundary.

The linearized PDE is

$\nabla\cdot(-\nabla du-\nabla u^k) = (u^k+du)^2 \approx (u^k)^2+2u^k du.$

Rearrange the above equation into standard form according to the details of InitializePDECoefficients in the Mathematica 11.0 documentation,

$\nabla\cdot(-\nabla du-\nabla u^k) + (-2u^k) du - (u^k)^2 =0.$

According to the details of InitializePDECoefficients,

the corresponding Coefficients are:

LoadCoefficients $f_1 \to (u^k)^2$

ReactionCoefficients $a_{11} \to -2u^k$

LoadDerivativeCoefficients $\gamma_1 \to \nabla u^k \, (i.e.,\{\frac{\partial u^k}{\partial x},\frac{\partial u^k}{\partial y}\})$.

The numerical solution blows up if I use the above coefficients. It will converge if I change the sign of the LoadDerivativeCoefficients to:

LoadDerivativeCoefficients $\gamma_1 \to -\nabla u^k \, (i.e.,\{-\frac{\partial u^k}{\partial x},-\frac{\partial u^k}{\partial y}\})$.

I'm confused why it works by changing the sign of the LoadDerivativeCoefficients. I also do not understand how the boundary condition is applied in the code. Can anyone explain it to me?

-----------------------------------------------------

According to @user21, there is a typo in the ref page of InitializePDECoefficients. The sign of $\alpha$ and $\gamma$ should be switched.

My understanding of the second question according to @user21's ANSWER: To implement the linearized equation using FEM in Mathematica, we break it into two parts, that is, the linear part and the non-linear ($u^k$-dependent) part:

$\nabla\cdot(-\nabla du)$ = 0 -- (1)

and

$\nabla\cdot(-\nabla u^k) + (-2u^k) du - (u^k)^2 =0.$ -- (2)

We create the load vector and stiffness matrix of the (1) outside the Newton loop to avoid repeated calculations of this part. Inside the loop, we update the load vector and stiffness matrix of the (2) for previous solution $u^k$. The total load vector and stiffness matrix are obtained by adding up the two corresponding parts. Since the increment $du$ should satisfy the Dirichlet BC that $du=0$ on the boundary, we then need to modify the load vector and stiffness matrix accordingly. It is important that the initial guess $u^0$ should satisfy the Dirichlet BC that $u=0$ on the boundary.

Answer 1

That is an unfortunate typo in the documentation in V11.2. The correct form is given here where you see that the "LoadDerivativeCoefficients" has an opposite sign than what is given in the InitializePDECoefficients ref page. So you changing the sign is the right thing to do. Note that $\alpha$ also has the wrong sign, but that is irrelevant in your example. I fixed the documentation a few weeks back when I realized the typo.

Concerning the boundary condition, note that we solve for the increment dU (not U) The boundary conditions, however, are made for applying to U. So we proceed by setting up the initial guess in such a way that it has the boundary values set, and at the same time modify the stiffness matrix to respect that in every iteration. So the increment of the boundary condition is 0.