# weak formulation of coupled pdes for fenics

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.

• Regarding the actual implementation, this is a FEniCS question -- please ask on the FEniCS mailing list. Regarding the weak formulation: what have you done so far? Getting the weak formulation isn't hard, in principle: you multiply by a test function and then integrate by parts where necessary. – Wolfgang Bangerth Jul 7 '13 at 18:50
• I would add to Wolfgang's comment that there are two places here where I think recasting the terms in divergence form would simplify and improve your formulation. – Bill Barth Jul 7 '13 at 20:51
• @BillBarth,what do you mean by recasting the terms in divergence? and where exactly? – Andrei Jul 7 '13 at 21:34
• @WolfgangBangerth, the difficulty comes in the last part of your statement, "where necessary" - how do I know where it's necessary and where not? – Andrei Jul 7 '13 at 21:42
• The typical way to finding weak forms is to multiply by a test function and integrate by parts using the divergence theorem to convert terms that looks like $\int (\nabla \cdot (k \nabla u)) v dx$ into $\int k \nabla u \cdot \nabla v dx$ plus a boundary term. You generally know where to apply integration by parts by looking for the divergence of something. It appears to me that the pair of terms in parens in the first equation and the whole term in square brackets in the second come from the divergence of something. Recasting this way should lead to knowledge of where to integrate by parts. – Bill Barth Jul 7 '13 at 21:57

## 1 Answer

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.