I’m looking at the advection-diffusion reaction demo, which solves the following equations \frac{\partial u_1}{\partial t} + w \cdot \nabla u_1 - \nabla \cdot (\epsilon \nabla u_1) = f_1 - Ku_1 u_2 \frac{\partial u_2}{\partial t} + w \cdot \nabla u_2 - \nabla \cdot (\epsilon \nabla u_2) = f_2 - Ku_1 u_2 \frac{\partial u_3}{\partial t} + w \cdot \nabla u_3 - \nabla \cdot (\epsilon \nabla u_3) = f_3 + Ku_1 u_2 - Ku_3
The demo propogates the solutions by the Euler method, using the following variational form at each time step \int_{\Omega} \Delta t^{-1} (u^{n+1}_1 - u^n_1) v_1 + w \cdot \nabla u^{n+1}_1 v_1 + \epsilon \nabla u^{n+1}_1 \cdot \nabla v_1 ) dx - \int_\Omega f_1 dx - \int_{\Omega} - K u^{n+1}_1 u^{n+1}_2 v_1 dx + \text{other terms} = 0.
The demo solves this equations with homogeneous Neumann boundary conditions on the entire boundary (specifically, \partial u_i / \partial n = 0 for i = 1,2,3 and n is along the surface), so that all the boundary integrals vanish.
My question is how is this well posed? The boundary condition only specifies the derivative of the functions at the boundary, so the solutions seem only defined up to an overall constant? Fenics is able to solve the problem just fine, however, and the demo documentation interprets the solution as physical.
This demo explains how to properly treat problems with only Neumann boundary conditions, by adding additional constraints of the form \int_\Omega u dx = constant.
It would be great if somehow Fenics “luckily” finds the “physically correct” solution by how it solves the problem for reaction-diffusion problems, as adding additional constraints only makes the FEM problem more complicated. Indeed, in my own (limited) reaction-diffusion simulations I have not yet noticed an issue. However, can I be sure that there won’t be subtle errors associated with the ill-poisedness of the PDE?
Thanks in advance!
You’re forgetting the time dimension. The initial condition u(t=0) = u_0 could be interpreted as the “Dirichlet” condition you’d probably looking for. The problem is well posed.
I’m not sure I understand your answer. Isn’t the u(t = 0) = u_0 essentially just a “source” term for the solution to u(t = \Delta t) equation? And then each solution thereafter the solution to the previous step is a source term for the next step. Where does a Dirichlet boundary condition come in?
You’re conflating the discretisation with the original problem. The original problem is posed on some domain which incorporates a space domain \Omega and time domain \mathcal{T}. So in order to find u : \Omega \times \mathcal{T} \rightarrow \mathbb{R} such that
u_t - \nabla^2 u + u = f
we need some kind of boundary data. So at t=0 we enforce u(x, t=0) = u_0(x) and on (x, t) \in \partial\Omega \times \{t \neq 0\} we impose \nabla u \cdot n = 0.
I understand what you are saying now, and I agree, but it doesn’t quite answer my question. I agree that the overall problem is well-posed, but how can I be sure that the discretization is? What’s to prevent the solver from finding an incorrect solution at the first time step? A solution that is wrong by an overall constant? Or any time step?
It seems to me that the solution is only guaranteed to be correct if the following constraints are added \int_{\Omega} dx u^{n + 1}_1 = \int_\Omega (u^{n}_1 - Ku^n_1u^n_2 ) dx \int_{\Omega} dx u^{n + 1}_2 = \int_\Omega (u^{n}_2 - Ku^n_1u^n_2 ) dx \int_{\Omega} dx u^{n + 1}_3 = \int_\Omega (u^{n}_3 + Ku^n_1u^n_2 - Ku^n_3) dx.
Does the solution found by Fenics satisfy these constraints even though they don’t appear explicitly in the formulation?
Oops! I resolved this issue. Thanks @nate for helping me out with this question!
As this demo explains, the Poisson equation indeed has no unique solution with pure Neumann boundary conditions: \nabla^2 \phi = f \text{ on } \Omega, \frac{\partial \phi}{\partial n} = g \text{ on } \partial \Omega.
This is because the equations only relate derivatives of \phi to each other. However, a discretized diffusion equation is of the form \Delta t \nabla^2\phi^{N+1} = \phi^{N+1} - \phi^{N} \text{ on } \Omega \frac{\partial \phi^{N+1}}{\partial n} = 0 \text{ on } \partial \Omega.
This is equivalent to modifying the Poisson equation above to include a new term \nabla^2 \phi = f + h\phi.
with h \neq 0. Now, this equation is well-posed even with pure Neumann boundary conditions, because it relates the function \phi with it’s derivatives, and does not relate only the derivatives of \phi as in the case of the Poisson equation.
Thus, not only is the entire time-dependent equation well-poised as pointed out by @nate, but also every PDE individually for each time-step in the time discretization.