How to Implement Non-Linear Problem in Firedrake

Question

I am new to firedrake / fenics. I believe I have a relatively simple non-linear toy problem to solve: 1D, 1-component diffusion reaction equation at steady-state. Boundary conditions are no flux at x = 1 and a constant concentration at x = 0.

When the reaction is first order m = 1, this problem solves and gives a reasonable solution. However, I cannot figure out how to implement this problem for reaction orders where m is not equal to 1 or 0.

I have tried various formulations for rxn_func and I have documented the various errors I receive when using these different formulations.

I have consulted the firedrake documentation, specifically: this unsteady non-linear problem. I tried to mimic the formulation there where they use inner(dot(u,nabla_grad(u)), v) to evaluate ((u⋅∇)u)⋅v.

I am pretty sure there is a simple solution, but I cannot decode the error messages I am receiving. Any help is appreciated.

'''
    Simple Diffusion Reaction Problem
    
    PDE / Governing Equation
    D₁ ∂²c₁/∂x² - kᵣₓₙ⋅(c₁)ᵐ = 0
    
    Boundary Conditions
    c₁(x=0) = c₁_∞ 
    ₁(x=L) = 0
'''

import firedrake as fd
import firedrake_adjoint as fda
from firedrake import inner, grad, dot, dx, dS, ds, jump, sqrt
import matplotlib.pyplot as plt
import numpy as np

import os

max_iter = 200

L      = 1.0                     # cm
Diff_1 = fd.Constant(1.0e-6)     # cm² / s
k_rxn = fd.Constant(3.0e-6)      # s⁻¹
c_1_bulk = fd.Constant(1.0e-3)   # mol / cm³

Flux_1_x_eq_L = fd.Constant(0.0)

# Create mesh
n = 100
mesh = fd.IntervalMesh(n, L) # (ncells, length)
x = fd.SpatialCoordinate(mesh)

# Define discrete function spaces and functions
V = fd.FunctionSpace(mesh, "CG", 1)
W = V

u, v = fd.TrialFunction(W), fd.TestFunction(W)
u_sol = fd.Function(W)

print(f"DOFs: {W.dim()}")

rxn_func = k_rxn * u         # works - reaction order = 1
rxn_func = k_rxn * dot(u, u) # NotImplementedError: Cannot take length of non-vector expression.
rxn_func = k_rxn * inner(u, u) # NotImplementedError: Cannot take length of non-vector expression.
rxn_func = k_rxn * pow(u, 2) # ufl.algorithms.check_arities.ArityMismatch: Applying nonlinear operator Power to expression depending on form argument v_1.

# Define and solve the PDE equations
A  = (-Diff_1 * inner(grad(u), grad(v)) - inner(rxn_func, v) ) * dx

# Add in the boundary conditions
# weakly imposed
L = -v * Flux_1_x_eq_L * ds(2)             # ₁(x=L) = Flux_1_x_eq_L

# strongly imposed
bc_c1_0 = fd.DirichletBC(W.sub(0), c_1_bulk, 1)    # c₁(x=0) = c₁_∞

# solve the system of equations
fd.solve(A == L, u_sol, bcs=[bc_c1_0])

c1 = u_sol

x_coord = mesh.coordinates.dat.data

# plot the concentration profile
fig, axes = plt.subplots(nrows=1, ncols=1)
axes.plot(x_coord, c1.dat.data[:] * 1e3)
axes.set_title('Concentration Profiles')
axes.set_xlabel('Position (cm)')
axes.set_ylabel('Concentration (mol / L)')

fig.tight_layout()
fig.savefig('OneComponent_Diffusion_Reaction.png', format='png')

How to Implement Non-Linear Problem in Firedrake

0 Answers0