# Implementing Robin conditions in the variational formulation of a diffusion-type equation

**URL:** <https://fenicsproject.discourse.group/t/implementing-robin-conditions-in-the-variational-formulation-of-a-diffusion-type-equation/3703>\
**Category:** variational formulation\
**Created:** [July 5, 2020, 7:36pm UTC](https://fenicsproject.discourse.group/t/implementing-robin-conditions-in-the-variational-formulation-of-a-diffusion-type-equation/3703 "2020-07-05T19:36:31Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![legibility](https://avatars.discourse-cdn.com/v4/letter/l/7ab992/32.png) [@legibility](https://fenicsproject.discourse.group/u/legibility)\
**Post date:** [July 5, 2020, 7:36pm UTC](https://fenicsproject.discourse.group/t/implementing-robin-conditions-in-the-variational-formulation-of-a-diffusion-type-equation/3703/1 "2020-07-05T19:36:31Z")

</div>

Hello all.

I am trying to solve a diffusion-type equation in 1D with Robin boundary conditions on the LHS and RHS. The code itself is probably too long to post here and ever receive an answer, but I have verified that the issue is in my implementation of the boundary conditions, so hopefully what follows will suffice.

To validate this, I have a known analytical solution that also fulfils the boundary conditions. However, as the time step increases, the numerical solution begins to differ from the numerical solution at the inner boundary, with the problem propagating to make the simulation useless.

However, from inspection of the numerical solution, **it also seems to satisfy the Robin boundary condition**! I am certain that the analytical solution is correct - so what could I be doing wrong?

I have read as much as I can about this issue, and tried many different options, but without success, so I am now hoping that someone could help me here.

## Mathematical formulation

The diffusion equation is as follows.

\frac{\partial u}{\partial t} =A \frac{\partial}{\partial x} \left[B \left( \frac{\partial u}{\partial x}\right)\right]\text{ on } [0, L]\\ \frac{\partial u}{\partial x} = Cu \text{ at } x=0, L

Discretising in time, multiplying by `TestFunction` v, and integrating, I have as follows:

u -\Delta t A \frac{\partial}{\partial x} \left[B \left( \frac{\partial u}{\partial x}\right)\right] = u^n\\ \int uv\mathrm{d} x -\Delta t \int A \frac{\partial}{\partial x} \left[B \left( \frac{\partial u}{\partial x}\right)\right]v\mathrm{d} x = \int u^n v \mathrm{d} x\\ \int uv\mathrm{d} x - \Delta t \left[Av \left[B \left( \frac{\partial u}{\partial x}\right)\right]\right]\_0^L + \Delta t \int \frac{\partial \left(A v\right)}{\partial x} \left[B \left( \frac{\partial u}{\partial x}\right)\right]v\mathrm{d} x = \int u^n v \mathrm{d} x

Applying the boundary condition:

\int uv\mathrm{d} x - \Delta t \left[Av \left[BC u\right]\right]\_0^L + \Delta t \int \frac{\partial \left(A v\right)}{\partial x} \left[B \left( \frac{\partial u}{\partial x}\right)\right]v\mathrm{d} x = \int u^n v \mathrm{d} x

## UFL implementation

Translating this for FEniCS is done as usual.

From [Simple Poisson equation in 1D with von Neumann condition on the left](https://fenicsproject.discourse.group/t/simple-poisson-equation-in-1d-with-von-neumann-condition-on-the-left/2179/3), I know that I need to explicitly impose the sign on the two surface elements separately. I believe that this could be a source of my mistake, but I have tried all four possible options (i.e. + or - in front of each of the two lines).

**Each of these options resulted in the exact same problem mentioned above, down to the apparent satisfaction of the Robin boundary condition**. This is what makes me suspicious of my implementation, as surely changing those ought to have some difference?

I mark the locations x=0, L by subclassing `Subdomain` and marking with `mark_0=0` and `mark_L=1`. I have checked that the marking function is implemented correctly.

Then, with `A, B, C` as suitable `Expressions`, and `dbdx_A` another expression containing the manually differentiated \frac{\mathrm{d} A}{x}, I can write down F:

```auto
F = u*v*dx
F += - Dt*A*B*C*u*v*ds(mark_L) # This sign has been switched with +, with no difference
F += + Dt*A*B*C*u*v*ds(mark_0) #This sign has been switched with -, with no difference
F += + Dt* (A*Dx(v,0) + dbdr_A*v)*B*Dx(u,0) * dx
F += - (u_n)*v*dx

```

_Aside - I intend for A and B to eventually depend on u - hence the non-linear implementation, but as yet they do not: (_`B = Expression('x[0]', degree = 5)`_). When they do, say for B = u^{1.5} x, my hope is that I can define, for example,_ `A = (u**1.5) * Expression('x[0]', degree=5)`.

I Can then solve as follows:

```auto
u.assign(u_n) # Avoids non-convergence issues for non-linear solvers on the first step, good practice for the time-dependent problem.
for n in range(100): # Take 100 steps
    t += Dt
    solve(F==0, u, bcs=[], solver_parameters={"newton_solver":{"linear_solver": "lu",
                                                            "convergence_criterion": "incremental", # or residual
                                                            "maximum_iterations": 30,
                                                            "relative_tolerance": 1e-7,
                                                            "absolute_tolerance": 1e-7}})

    u_n.assign(u)

```

## Validation failing

To check the suitability of the solution, I have an analytical solution u\_A (x, t). However, when I set the initial condition of u as u\_A (x, t=0) the error increases in a discretisation-independent way (i.e. for varying mesh density and timestep).

As I mentioned above, if I enforce Dirichlet boundary conditions, setting u(0, t) = u\_A (0, t) and u(L, t) = u\_A (L, t) for all t via `bcs`, it works with perfect agreement and appropriate discretisation-dependent error variation. When I do so I simply remove the lines of `F` containing `ds(-)`, as (from my understanding) they are then no longer evaluated.

# Discussion

I hope that this is a problem to do with my understanding of the variational formulation, and not something more subtle to do with degrees of expressions or the `solver_parameters`. I have run with extremely fine x, t discretisation parameters and values of `Expression` `degree` up to 10 to no avail.

Since it’s a time-dependent problem, I don’t think I need any further conditions to be imposed, unlike the example from [https://fenicsproject.org/docs/dolfin/latest/python/demos/neumann-poisson/demo\_neumann-poisson.py.html](https://fenicsproject.org/docs/dolfin/latest/python/demos/neumann-poisson/demo_neumann-poisson.py.html), but please do correct me.

Many thanks!

_Aside - it appears that it is not possible to have, for example,_ `Dx(A*v, 0)`_, though I have as yet failed to find an explanation of why - but I would love to find out._

---

<div class="post-metadata">

**Author:** ![huadong\_gao](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/huadong_gao/32/191_2.png) [@huadong\_gao](https://fenicsproject.discourse.group/u/huadong_gao)\
**Post date:** [July 6, 2020, 12:26pm UTC](https://fenicsproject.discourse.group/t/implementing-robin-conditions-in-the-variational-formulation-of-a-diffusion-type-equation/3703/2 "2020-07-06T12:26:28Z")

</div>

Hi legibility

` It seems that the boundary condition given by $$\frac{\partial u}{\partial x} = Cu \text{ at } x=0, L$$ is not common. The Robin boundary condition should be like this $$C u+\frac{\partial u}{\partial \nv} =0, \text{ at } x=0, L$$, where $C>0$.`