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

variational formulation
1

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, 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:

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:

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, 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.

2

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$.