Boundary conditions for nonlinear advection-diffusion

Question

Hello, I'm working on a nonlinear advection-diffusion equation in 1D but a little unsure about the boundary conditions. The equation depends on a variable q and goes as follows:

$$ \frac{d}{dt}f(y,t) = \nabla [ J(y)f(y,t) + D \nabla f^{2-q}(y, t)] $$

For the time integration I chose an implicit Euler scheme

# velocity field
I = Expression(('(2.0-q)*sinh(x[0])',),q=q,degree=2)
def operator(u, v):
    return factor*(2.0-q)*pow(u,1.0-q)*dot(grad(u),grad(v))*dx + u*dot(I, grad(v))*dx 
    #? + f*v*dx - g*v*ds

F = inner(u-u0, v)*dx + dt*operator(u, v)

The initial condition is composed of two Gauß peaks, symetrically around zero. Since I don't want any of my 'mass' to escape the system I first imposed Dirichlet BC's. But then I discovered this link scicomp.stackexchange.com that argues I need Robin BC's which now makes sense. This would mean
that the term in square brackets from above has to vanish at the boundaries. Do I just have to set g to zero?
Also I just get a solution for values of q in $[1, 1.7]$ even though I should also be able to choose negative values or larger ones than 1.7. But doing that the Newton solver returns NaN and fails. Do my functions have to fulfill continiuty conditions in order for the solver to work? I already tried adjusting my mesh and time sizes, but I get most 3 or 4 time step integrations before it fails. Can this be connected to the boundary conditions?

Comments

From your equation, integrate the diffusive term by parts; this results in a diffusive flux term, which if you want zero diffusive flux, set it to zero. Next, if you have to impose the advective flux at the boundary. If the time-invariant velocity $J$ is zero on the boundary, there is no advective flux.

What do you mean by $f^{2-q}$? is that taking $f$ to the power of $2-q$? If so you should bring the entire function to that power, including parentheses. Also, it is nice to have a working bit of example code...

Here's some:

from fenics import *
import matplotlib.pyplot as plt

mesh = UnitIntervalMesh(1000)
P1   = FunctionSpace(mesh, 'CG', 1)

u    = TrialFunction(P1)
v    = TestFunction(P1)
u0   = Function(P1)

J    = interpolate(Expression("-pow(x[0] - 0.5, 2) + 0.25", degree=1), P1)
D    = Constant(0.1)
g    = Constant(0.0)
dt   = Constant(0.001)
dudt = (u - u0) / dt

R = + J * u * v * dx \
    - D * u.dx(0) * v.dx(0) * dx \
    + D * g * v * ds \
    - dudt * v * dx

a = lhs(R)
L = rhs(R)

u = Function(P1)
u0.vector()[500] = 100.0
u.vector()[500]  = 100.0

fig = plt.figure()
ax  = fig.add_subplot(111)
p,  = ax.plot(u.vector().array())

plt.ion()
plt.show()

t = 0.0
while t < 10:
  u0.assign(u)
  solve(a == L, u)
  t += dt(0)

  p.set_ydata(u.vector().array())
  ax.set_ylim([u.vector().min(), u.vector().max()])
  plt.draw()
  plt.pause(0.00000001)

plt.ioff()

The explicit scheme introduces errors, I think this is why the solution appears to gain "mass"...

edit The pure Neuman problem here is not well posed, the solution is only unique up to a constant. Additional constraint on the solution is required, but I am not sure what... Imposition of Dirichlet BC weakly is the only thing I can think of.

---

Thank you for your answer! I've got indeed a working version but just for limited q, which I find
strange. Yes the term is $[f(y,t)]^{2-q}$. I also tried differetiating it, giving me $(2-q) f^{1-q} \nabla f(y, t)$ , because that might be faster. If you try values for q greater than 1.8 or smaller than 1.0 the Newton solver doesn't converge, and I've got no idea why. The norm should be conserved because I use a implicit euler scheme. Maybe you can answer me also another question. One can either partially integrate the advection term and endfd up multiplying it by $\nabla v$ or you can do the differetiation ($sinh(y)$ is also affected by the $\nabla$) thus adding an additional source term.Is one mehod better than the other?

from fenics import *
import os
import numpy as np

parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

T = 2.0
Nt = 1000
dt = T/Nt
L = 8.0
Nx = 500

mesh = IntervalMesh(Nx, -L, L)
V = FunctionSpace(mesh, 'CG', 1)

u_D = Constant(0.0)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

gamma = 0.281
sigma = gamma/np.sqrt(8*np.log(2.0))
y0 = 2.915;

ic = Expression('1/sqrt(8*pi*pow(sigma, 2.0)) * exp(-pow(x[0]-y0, 2)) / (2*pow(sigma, 2.0)) \
+ 1/sqrt(8*pi*pow(sigma, 2.0)) * exp(-pow(x[0]+y0, 2)) / (2*pow(sigma, 2.0))',
                 y0=y0,sigma=sigma,pi=np.pi,degree=2)

v = TestFunction(V)
u0 = Function(V)
u = Function(V)
f = Constant(0.0)
g = Constant(0.0)

name = "test"
q = 0.95
I = Expression(('(2.0-q)*sinh(x[0])',),q=q,degree=2)
factor = 160.0/1150.0

def operator(u, v):
    return factor*(2.0-q)*pow(u,1.0-q)*dot(grad(u),grad(v))*dx 
                  + u*dot(I, grad(v))*dx 
                  + f*v*dx - g*v*ds
    #alternatively factor*dot(grad(pow(u,2.0-q)),grad(v))*dx 


F = inner(u-u0, v)*dx + dt*operator(u, v)

# normalize to one. integrate over the whole domain
u0.interpolate(ic)
S = assemble(u0*dx(mesh))
u0.vector()[:] /= S
u.assign(u0)

solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

solve(F == 0, u, bc)

vtkfile = File('clean/clean.pvd')
t = 0.0
u.rename('u', 'density')


vtkfile << (u, t)

while t < T:
    solve(F == 0, u, bc)
    vtkfile << (u, t)

    u0.assign(u)
    t += dt

Answer 1

The range of q that allows the solution to converge might have to do with the derivative (which is calculated via the chain rule for the Newton iterations) of the term $f^{2-q}$ being undefined when $f$ is zero.

I wrote about how to deal with this here.