weak formulation of coupled pdes for fenics

Question

I am trying to implement the following system of time-dependent, coupled nonlinear pdes in FEniCS:

$$\partial_{t}\rho+2\left(\nabla\rho\nabla\phi+\rho\nabla^{2}\phi\right)+2\sigma\rho\left(\rho-1\right)=0$$

$$\partial_{t}\phi+\left(\nabla\phi\right)^{2}-\frac{1}{2}\left[\frac{\nabla^{2}\rho}{\rho}-\frac{1}{2}\frac{\left(\nabla\rho\right)^{2}}{\rho^{2}}\right]+\rho-1=0$$

where the unknown scalar functions $\rho(\mathbf{r},t)$ and $\phi(\mathbf{r},t)$ are real-valued (2D). Boundary conditions are homogeneous Neumann for both functions.

As I understand from the documentation of Fenics, the first step is to write the weak formulation of the problem, but the examples in the tutorial are all very basic and I would really appreciate if someone showed me how to implement the system above.

Answer 1

To find the weaks forms we use the identity $$ \nabla\cdot(\rho\,\nabla\phi) = \rho\,\nabla^2\phi + \nabla\rho \cdot \nabla\phi $$ and its vector counterpart $$ \nabla\cdot(\varphi\,\mathbf{v}) = \varphi\,(\nabla\cdot\mathbf{v}) + (\nabla\varphi) \cdot \mathbf{v} \,. $$

---

Then the first equation can be written as $$ \partial_{t}\rho+2\,\nabla\cdot(\rho\,\nabla\phi)+2\sigma\rho\left(\rho-1\right)=0 \,. $$ If we consider only the spatial terms and use a weighting function $w$, the weak form is $$ \int_\Omega (\partial_{t}\rho)\,w\,d\Omega + 2\int_\Omega \left[\nabla\cdot(\rho\,\nabla\phi)+\sigma\rho\left(\rho-1\right)\right]w\,d\Omega=0 \,. $$ Using the vector identity above, we get $$ \int_\Omega (\partial_{t}\rho)\,w\,d\Omega + 2\int_\Omega \left[\nabla\cdot(w\,\rho\,\nabla\phi)-\nabla w\cdot(\rho\,\nabla\phi) + w\sigma\rho\left(\rho-1\right)\right]\,d\Omega=0 \,. $$ Apply the divergence theorem to get $$ \int_\Omega (\partial_{t}\rho)\,w\,d\Omega + 2\int_\Gamma (w\,\rho\,\nabla\phi)\cdot\mathbf{n}\,d\Gamma - 2\int_\Omega \left[\nabla w\cdot(\rho\,\nabla\phi) - w\sigma\rho\left(\rho-1\right)\right]\,d\Omega=0 \,. $$

---

Using the same identity as before, we also have $$ \nabla\phi\cdot\nabla\phi = \nabla\cdot(\phi\,\nabla\phi) - \phi\,\nabla\cdot(\nabla\phi) $$ and $$ \nabla\cdot\left(\frac{1}{\rho}\,\nabla\rho\right) = \frac{\nabla^2\rho}{\rho} -\frac{\nabla\rho \cdot \nabla\rho}{\rho^2}\,. $$ Assuming that the extra factor of $1/2$ in the second equation is a typo, we can write that equation as $$ \partial_{t}\phi+\nabla\cdot(\phi\,\nabla\phi) - \phi\,\nabla\cdot(\nabla\phi)-\frac{1}{2}\nabla\cdot\left(\frac{1}{\rho}\,\nabla\rho\right)+\rho-1=0 $$ It is here that the appropriate choice of terms to consider at the boundaries becomes important. The choice is typically determined by the ease of application of the boundary conditions and how the two sets of equations are coupled in the numerical method.