Hi, I am working on the Navier-Stokes equations coupled with a concentration-diffusion equation. I have been checking the convergence rates / error calculations, and the velocity seems to be doing very well in certain cases (e-14 error), but for the concentration the results are poor. Refining the mesh leads to a steady decrease in the error for the velocity case, but not for the concentration, and I am trying to understand the reason. Essentially I use
where u is the solution function for the velocity, c is the unknown concentration and d is the concentration test function. On the right-hand side I use some Expression f for the forcing term. ds(0) is the boundary integral on the entire boundary.
I was wondering if it is somehow the boundary integrals that I should not use like this? Because I tried, and this seems to be not equivalent to:
inner(u, grad(c)) * d * dx - (c.dx(0)).dx(0) * d * dx - (c.dx(1)).dx(1) * d * dx,
which I don’t understand, as simply by integration by parts these should be the same, right?
Hi,
Which function space are you using to compute the concentration? If you use a first order function space div(grad(u))v dx != inner(n,grad(u))ds + inner(grad(u),grad(v))dx
To expand on @dokken’s point, see my answer here:
Also, note that including the boundary term from integration-by-parts (assuming there isn’t a DirichletBC, in which case d would be zero) is essentially leaving the boundary condition unspecified. Without choosing a specific boundary condition, there are many possible solutions satisfying the interior equation, so the problem for c will be ill-posed.
Thank you so much for your answers. The problem occurs similarly in first order and second order space choices for C. I am trying to understand what you mean by “including the boundary term” in the sense of, what else can I do? I mean, mathematically it is simply incorrect not to include the boundary integral isn’t it, when you switch from (1) \Delta u * v * dx to (2) \nabla u * \nabla v * dx? Or in fenics the way I am supposed to do this is that I can simply use (2) instead of (1) and I have to specify the boundary conditions everywhere?
But what I really don’t understand is, if the boundary value is not zero let’s say on \Gamma, and you use the form (2) without the boundary integral, then your scheme will miss the nonzero boundary integral… so I don’t get it how can it be correct?
What is essentially a good combination of specifying BCs, including boundary integrals in the F form, and using approach (2) instead of approach (1)?
For most diffusion-type problems, you would divide the full boundary into two disjoint parts, the Dirichlet boundary, where you set u=g for some given g (which is usually implemented in FEniCS using a DirichletBC) and the Neumann boundary, where you enforce a Neumann boundary condition of the form \nabla u\cdot\mathbf{n} = h, for some given h (or, more generally, replacing \nabla u with the diffusive flux, if there is some diffusion tensor not equal to identity). The Neumann condition is enforced by a boundary term in which \nabla u\cdot\mathbf{n} is replaced with h (so it no longer depends on u, and becomes part of the linear form or “right-hand side”). There are many problems of physical interest where h=0 (e.g., an insulated boundary in heat conduction), and the boundary term is therefore omitted; this case is sometimes referred to as the “do-nothing boundary condition”.
Because you have both advection and diffusion, it is not always stable to enforce a Neumann boundary condition on just the diffusive flux, in particular on regions of the boundary where the flow velocity is entering the domain (i.e., the “inflow boundary”). The clearest explanation I have seen of this is found in Section 2 of the following technical report, which I’d recommend reading: https://www.oden.utexas.edu/media/reports/2004/0431.pdf
Also, if the advection velocity is large relative to the diffusion coefficient, then additional stabilization is needed to get reasonable solutions on coarse meshes; there is a FEniCS undocumented demo on one of the most popular stabilization schemes, called “SUPG”, which is also discussed in the linked report.
I appreciate your answer, thanks a lot. What is still unclear to me is: in general, for an unknown quantity s, for example if we know that on the boundary section \Gamma_1 we have s(0, y) = \sin(y), so a nonzero boundary value, and in fenics we set DirichletBC(S, Expression('sin(x[1])'), Gamma_1), then in this particular case, do we still need to include also the boundary integral in the form F on the \Gamma_1 part of the boundary? So
is it enough to simply just set the DirichletBC, or
do I also need to add something like inner(div(s), p)*ds(1) in the form F, where ds(1) is the boundary integral on \Gamma_1?
Maybe this setting sounds artificial but I feel like it would help me understand what happens. Thank you again!
Nora
On the Dirichlet boundary (if the boundary condition is strongly-enforced, i.e., built directly into the trial solution space, as is done by DirichletBC in FEniCS), then the test function goes to zero at the boundary, even if the solution is set to nonzero boundary data. So, you could include the boundary term from integration-by-parts, but it would be zero, and not affect the result.
So do I understand well that in FEniCS the test functions go to zero on and only on the boundary section for which we have DirichletBC defined? I understand they do go to zero on parts where we have DirichletBCs, but I guess they do not go to zero on parts where we have a Neumann BC, otherwise (with \nabla u \cdot n = h) the h * v * ds(1) boundary integral term in the form F would vanish.
So essentially, boundary integral terms in the F form should be used exactly for those boundary sections where we do not have a DirichletBC? And among these, where we do not have a Dirichlet condition, we wither have a NeumannBC, or, we risk the problem being ill-posed, right?