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