How to implement inhomogeneus Neumann boundary condition

1

Hello everyone,

I am using FEniCS to solve the following 1D system equation:

\begin{array}{l} \frac{{\partial \xi }}{{\partial t}} = - {{\hat L}_\sigma }\hat g\left( \xi \right) + \hat \kappa \frac{{{\partial ^2}\xi }}{{\partial t}} - {{\hat L}_\eta }\dot h\left( \xi \right)\left( {{e^{\left( {1 - \alpha } \right)\hat \phi }} - \hat c{e^{ - \alpha \hat \phi }}} \right)\\ \frac{\partial }{{\partial x}}\left[ {\hat K\left( \xi \right)\frac{{\partial \hat \phi }}{{\partial x}} + {{\hat K}_D}\left( \xi \right)\frac{1}{{\hat c}}\frac{{\partial \hat c}}{{\partial x}}} \right] = - \frac{{\partial \xi }}{{\partial t}}\\ \hat \varepsilon \left( \xi \right)\frac{{\partial \hat c}}{{\partial t}} - \frac{\partial }{{\partial x}}\left( {\hat D\left( \xi \right)\frac{{\partial \hat c}}{{\partial x}}} \right) = - \frac{{\partial \xi }}{{\partial t}} \end{array}

With the following boundaries conditions:

\begin{array}{l} \xi = 1,\,\,\,\hat \phi = - 0.05,\,\,\,\frac{{\partial \hat c}}{{\partial x}} = 0\,\,\,\,\,\,at\,\,\hat x = 0\\ \xi = 0,\,\,\,0.3\frac{{\partial \hat \phi }}{{\partial x}} + 0.4\frac{{\partial \hat c}}{{\partial x}} = 1,\,\,\,\hat c = 1\,\,\,at\,\,\hat x = 20 \end{array}

I´m trying to implement the inhomogeneus Neumann boundary condition as the manual explains and I´m not sure what I´m doing wrong.

The code regarding to the boundaries conditions is the next :

# %%Define boundary condition
lox = 20

x = sym.symbols('x[0]') 
xi, c, phi = u.split()

# g = dphi/dx =  1/0.3 - (0.4/0.3)*dc/dx

g = Constant(1/0.3) - Constant(0.4/0.3)*sym.diff(c, x).subs(x,lox)

# Define boundary conditions
boundary_conditions = {1: {'Neumann': g}}   # x = 1
                       
# Define boundary subdomains
tol = 1e-14

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], lox, tol)

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim())
#boundary_markers.set_all(9999)
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 1)

#Redefine boundary integration measure
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Collect Neumann integrals
integrals_N = []
for i in boundary_conditions:
    if 'Neumann' in boundary_conditions[i]:
        if boundary_conditions[i]['Neumann'] != 0:
            g = boundary_conditions[i]['Neumann']
            integrals_N.append(g*v_2*ds(i))

def boundary0(x, on_boundary):
    return on_boundary and near(x[0], 0)
def boundaryL(x, on_boundary):
    return on_boundary and near(x[0], lox)

bc_xi1 = DirichletBC(V.sub(0), Constant(chi_ini), boundary0)

bc_xi2 = DirichletBC(V.sub(0), Constant(chi_fin), boundaryL)

bc_c2 = DirichletBC(V.sub(1), Constant(c_fin), boundaryL)

bc_phi1 = DirichletBC(V.sub(2), Constant(phi_ini), boundary0)

bcs = [bc_xi1, bc_xi2,  bc_c2, bc_phi1] # Condiciones dirichlet

################################################################

The message error is the following:

File “”, line 1, in
runfile(’/home/josemanuel/Documents/Fenics/Newmann1D.py’, wdir=’/home/josemanuel/Documents/Fenics’)

File “/usr/lib/python3/dist-packages/spyder_kernels/customize/spydercustomize.py”, line 827, in runfile
execfile(filename, namespace)

File “/usr/lib/python3/dist-packages/spyder_kernels/customize/spydercustomize.py”, line 110, in execfile
exec(compile(f.read(), filename, ‘exec’), namespace)

File “/home/josemanuel/Documents/Fenics/Newmann1D.py”, line 154, in
g = Constant(1/0.3) - Constant(0.4/0.3)*sym.diff(c, x).subs(x,lox)

File “/usr/lib/python3/dist-packages/sympy/core/function.py”, line 2448, in diff
return Derivative(f, *symbols, **kwargs)

File “/usr/lib/python3/dist-packages/sympy/core/function.py”, line 1243, in new
expr = sympify(expr)

File “/usr/lib/python3/dist-packages/sympy/core/sympify.py”, line 328, in sympify
coerced = coerce(a)

File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py”, line 468, in float
raise RuntimeError(“Cannot convert spatially varying function to float.”)

RuntimeError: Cannot convert spatially varying function to float.

Anyone who can help me?
Thank you in advance.

2

It’s a little hard to debug without a fully-runnable example, but it looks like you’re trying to use SymPy to take derivatives. That should not be necessary, since FEniCS’s Unified Form Language (UFL) has robust symbolic differentiation capabilities.

For the case of partial derivatives with respect to spatial coordinates, you can use the dx() method of a UFL object, e.g.,

#g = Constant(1/0.3) - Constant(0.4/0.3)*sym.diff(c, x).subs(x,lox)
g = Constant(1/0.3) - Constant(0.4/0.3)*c.dx(0)

where the argument to dx specifies the index of coordinate with respect to which the partial derivative is taken. (N.b. that the dx() method is not related to the integration measure object dx.) Also, because the finite element assembly procedure will only evaluate this term on the specified boundary, there is no need to manually substitute the value of the boundary’s x coordinate.

For other symbolic derivatives, there is a diff function in UFL similar to the SymPy diff; see the discussion here for an example.

1 Like