Expressions depending on spatially varying coefficients and their derivatives in advection-diffusion-PDE

Hello everyone,

I have a question regarding the implementation of spacially varying parameters for a 2d advection-diffusion-reaction-equation (ADRE). My main problem is that the PDE I want to solve does not arise from common applications, but instead as an adjoint problem to an ADRE, so I was not able to find an example of how to implement it in FEniCS.

The PDE defined an on a quadratic region \Omega=[0,X]\times[0,Z] looks like this:

\partial_t u = \alpha (\partial_{xx}u + \partial_{zz}u) + (\partial_x\alpha - \beta_x) \partial_{x} u + (\partial_z\alpha - \beta_z) \partial_{z} u + \gamma u

with Robin BCs (also depending on \alpha(x,z) and \beta(x,z)) on the left and right side
-\partial_x u(0,z,t) = \frac{\beta_x(0,z,t)}{\alpha(0,z,t)}u(0,z,t)
-\partial_x u(X,z,t) = \frac{\beta_x(X,z,t)}{\alpha(X,z,t)} u(X,z,t)

and Dirichlet BCs at the top and bottom
u(x,Z,t) = 0
u(x,0,t) = 0

with the original diffusion coefficient \alpha(x,z), flow rate \beta(x,z) (with its x-component \beta_x(x,z) and z-component \beta_z(x,z)) and reaction coefficient \gamma(x,z), all sufficiently smooth etc. and implemented as dolfin Expressions.

My main problem lies in the FEniCS implementation of the terms (\partial_x\alpha - \beta_x) and (\partial_z\alpha - \beta_z) in the PDE and \frac{\beta_x(X,z,t)}{\alpha(X,z,t)}, \frac{\beta_x(0,z,t)}{\alpha(0,z,t)} in the BCs. Most of the things I tried resulted in shape errors or similar, I am simply not aware how to properly combine Expressions. Here is a minimal example of the rest of my code, with some failed attempts at implementing it commented out:

from fenics import *
import numpy as np

T = 5.0            # final time
X = 0.5
Z = 1.0
num_steps = 400    # number of time steps
dt = T / num_steps # time step size
# Placeholder for solution
values = []

# Subdomains
class BoundaryZ1(SubDomain):
     # def inside(self, x, on_boundary):
         # return x[1] >= Z - DOLFIN_EPS and on_boundary
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1], Z, tol)

class BoundaryX1(SubDomain):
     # def inside(self, x, on_boundary):
         # return x[0] >= X - DOLFIN_EPS and on_boundary
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0], X, tol)  
         
class BoundaryX0(SubDomain):
     # def inside(self, x, on_boundary):
         # return x[0] <= DOLFIN_EPS and on_boundary
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0], 0, tol)
        
class BoundaryZ0(SubDomain):
     # def inside(self, x, on_boundary):
         # return x[1] <= DOLFIN_EPS and on_boundary
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1], 0, tol)
        # return on_boundary

# Create mesh and define function space
nx = nz = 50
mesh = RectangleMesh(Point(0,0),Point(X,Z),nx,nz)

V = FunctionSpace(mesh, 'CG',1)
Vv = VectorFunctionSpace(mesh, "CG", 1)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# Diffusion Coefficient and Velocity field as expressions

alpha = Expression(("(0.01*(x[0]+1))+(0.01*(x[1]+1))"), domain=mesh, degree = 2)
beta = Function(Vv)
beta = Expression(("(x[1])","(-x[0]-0.5)"), domain=mesh, degree = 2)
gamma = Constant('0') # Reaction coefficient set to zero for now

# Attempts at implementation
# betax = project(beta[0],V)
# betaz = project(beta[1],V)
# beta_s = Expression(("al[0]-bx","al[1]-bz"), al=alpha.dx, bx=betax, bz=betaz, domain=mesh, degree = 2)
# beta_s = Expression(("-bx","-bz"), al=alpha.dx, bx=beta, bz=beta, domain=mesh, degree = 2)

# Subdomains
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1);
boundary_markers.set_all(9999)
# Inlet
bz1 = BoundaryZ1()
bz1.mark(boundary_markers, 0)
# Left Wall
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 1)
# Right Wall
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 2)
# Outlet
bz0 = BoundaryZ0()
bz0.mark(boundary_markers, 3)
# Define boundary conditions

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Define initial value
u_Ini = Function(V)

u_Ini = Constant(1)
u_n = interpolate(u_Ini, V)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)

t = 0.0

g = Constant(0) # Neumann Value
# r = Expression(("a/b"), a=alpha, b=betax, domain=mesh, degree = 2) # ?  # Robin Value ("Diffusion Coefficient")
s = Constant(0) # Robin Value ("Ambient temp.")
u_D = Constant(0)

boundary_conditions = {0: {'Dirichlet': u_D},
                       1: {'Robin': (r, s)},
                       2: {'Robin': (r, s)},
                       3: {'Dirichlet': u_D}
                       }

bcs = []
for j in boundary_conditions:
    if 'Dirichlet' in boundary_conditions[j]:
        bc = DirichletBC(V, boundary_conditions[j]['Dirichlet'],
                         boundary_markers, j)
        bcs.append(bc)
integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        g = boundary_conditions[i]['Neumann']
        integrals_N.append(g*v*ds(i))
integrals_R = []
for i in boundary_conditions:
    if 'Robin' in boundary_conditions[i]:
        r, s = boundary_conditions[i]['Robin']
        integrals_R.append(r*(u - s)*v*ds(i))



# Time-Stepping formulation from the FEniCS Book, Chapter 31. 
# theta = 1 Backwards Euler
# theta = 0.5 Crank-Nicholson
# theta = Constant(1.0)
theta = Constant(0.5)

a0= alpha_s*inner(grad(u_n),grad(v))*dx + inner(beta_s,grad(u_n))*v*dx + gamma*u*v*dx + f*v*dx 
a1= alpha_s*inner(grad(u),grad(v))*dx + inner(beta_s,grad(u))*v*dx + gamma*u*v*dx + f*v*dx

A = (1/dt)*inner(u, v)*dx - (1/dt)*inner(u_n,v)*dx + theta*a1 + (1-theta)*a0 + sum(integrals_N) +  sum(integrals_R)

F = A
a = lhs(F)
L = rhs(F)
    
u = Function(V)

u.assign(u_n)
u.rename("u", "u")

problem = LinearVariationalProblem(a, L, u, bcs)
solver  = LinearVariationalSolver(problem)
t = 0.0
for n in range(num_steps):
    solver.solve()
    u_n.assign(u)
    t += dt

Any help would be massively appreciated!

Best regards
Folke

Please note that your code is far more than a simple example, and the following would have been sufficient to describe your problem:

from fenics import *
import numpy as np

X = 0.5
Z = 1.0
num_steps = 400    # number of time steps

# Create mesh and define function space
nx = nz = 50
mesh = RectangleMesh(Point(0,0),Point(X,Z),nx,nz)
alpha = Expression(("(0.01*(x[0]+1))+(0.01*(x[1]+1))"), domain=mesh, degree = 2)
beta = Expression(("(x[1])","(-x[0]-0.5)"), domain=mesh, degree = 2)
gamma = Constant('0') # Reaction coefficient set to zero for now

which can be solved with the following lines;

alphax = alpha.dx(0)
betaz = beta.dx(1)

Hello, and thank you for answering!

Ill admit my example is a bit involved, however I thought it might be useful to post some context to how I implement my variational formulation and boundary conditions.

I do not understand what you refer to when you say your code “solves” it: In my issue as described above I do not use the derivative of \beta, for example for the term \partial_x\alpha - \beta_x I need the difference between the derivative in x of \alpha and the x-component of the vector-valued \beta. That is what I cannot get to work in FEniCS, and I’m asking how to implement the operations.

If I try, for example:

alpha = Expression(("(0.01*(x[0]+1))+(0.01*(x[1]+1))"), domain=mesh, degree = 2)
beta = Function(Vv)
beta = Expression(("(x[1])","(-x[0]-0.5)"), domain=mesh, degree = 2)
gamma = Constant('0') # Reaction coefficient set to zero for now

betax = project(beta[0],V)
betaz = project(beta[1],V)
alphadx = alpha.dx(0)
alphadz = alpha.dx(1)
alpha_s = alpha
beta_s = Expression(("alx-bx","alz-bz"), alx=alphadx, alz=alphadz, bx=betax, bz=betaz, domain=mesh, degree = 2)
gamma_s = gamma

g = Constant(0) # Neumann Value
r = Expression(("b/a"), a=alpha, b=betax, domain=mesh, degree = 2) # Robin Value ("Diffusion Coefficient")
s = Constant(0) # Robin Value ("Ambient temp.")

The I get a sleugh of errors, in essence saying that dijitso is unable to compile my Expression for \beta_s = (\partial_x\alpha - \beta_x).

So, how would I implement it in a way that works?

As there was so much other code in your last code, I did not get what your actual issue was.
Your actual issue can be solved as follows (no need to wrap an expression inside an expression):

from fenics import *

X = 0.5
Z = 1.0
num_steps = 400    # number of time steps

# Create mesh and define function space
nx = nz = 50
mesh = RectangleMesh(Point(0,0),Point(X,Z),nx,nz)

alpha = Expression(("(0.01*(x[0]+1))+(0.01*(x[1]+1))"), domain=mesh, degree = 2)
beta = Expression(("(x[1])","(-x[0]-0.5)"), domain=mesh, degree = 2)
alphax = alpha.dx(0)
alphaz = alpha.dx(1)
beta_s = as_vector((alphax-beta[0], alphaz-beta[1]))
assemble(inner(beta_s, beta_s)*dx)
1 Like

Thank you! Your suggestion was exactly what I was looking for.

My initial problem is solved. However, I have a follow up question that popped up when I looked at the obtained solution relating to my implementation of the Robin-BCs:

My initial implementation of the Robin-BCs such as -\partial_xu(0,z,t) = \frac{\beta_x(0,z,t)}{\alpha(0,z,t)}u(0,z,t) was based on the requirement \alpha(0,z,t) \partial_x u(0,z,t) + \beta_x(0,z,t)u(0,z,t) = 0 for x=0 resp. x=X, which I tried to write in the form -\kappa\frac{\partial u}{\partial n} = r(u-s) as in Volume I of the FEniCS Tutorial book with r=\frac{\beta_x(0,z,t)}{\alpha(0,z,t)} and s=0.

The solution I obtained now shows that at x = 0 and x = X the state variable u is constant at zero while at the same time the slope is set to r=\frac{\beta_x(0,z,t)}{\alpha(0,z,t)}, which does not fulfill the requierement \alpha(0,z,t) \partial_x u(0,z,t) + \beta_x(0,z,t)u(0,z,t) = 0 at all. Has anybody experienced this when implementing Robin-BCs? Is this behaviour intended, or did I simply make a mistake implementing it?

Note that \beta/\alpha is not well defined if \alpha=0, so I would be carefull to scale your PDE by this.
I would also like to note that after doing integration by parts, the alpha will be part of the boundary integral.
Also most PDEs with spatially varying coefficients is on the form \nabla \cdot (\alpha(x) \nabla u), and not \alpha \nabla \cdot( \nabla u). Make sure that this is the PDE you would like to solve. (Edit: I guess you have expanded your differentiation operation before doing integration by parts. I would strongly suggest to not do this)