How to include nonlocal coupling in PDEs

variational formulation
1

Hey Fenics community,

I am currently trying to solve a system of a pair partial differential equations. Two of my trial functions within my system are the three-dimensional vector-valued functions \boldsymbol{x}(u) and \boldsymbol{f}(u) defined on the unit interval u\in[0,1]. We call \boldsymbol{x}(u) the centreline coordinates, \boldsymbol{f}(u) the force line density and u the material coordinate.

The two functions \boldsymbol{x}(u) and \boldsymbol{f}(u) are fully coupled through a complicated integral equation

\begin{align} \tag{1} \boldsymbol{u}(u)&=\frac{1}{8\pi\mu}\left(\boldsymbol{u}_{\text{local}}(u)+\boldsymbol{u}_{\text{nonlocal}}(u)\right), \\ \end{align}
\begin{align} \boldsymbol{u}_{\text{local}}(u)&=\ln\left(\frac{L^{2}-4s(u)^{2}}{(R_{0}\varphi(u))^{2}}\right)\left[\boldsymbol{I}+\hat{\boldsymbol{t}}(u)\otimes\hat{\boldsymbol{t}}(u)\right]\cdot\boldsymbol{f}(u)\\ \tag{2} &+\left[\boldsymbol{I}-3\hat{\boldsymbol{t}}(u)\otimes\hat{\boldsymbol{t}}(u)\right]\cdot\boldsymbol{f}(u) \end{align}
\begin{align} \boldsymbol{u}_{\text{nonlocal}}(u)&=\int_{0}^{1}\left[\left(\frac{\boldsymbol{I}+\hat{\boldsymbol{r}}(u,u')\otimes\hat{\boldsymbol{r}}(u,u')}{r(u,u')}\right)\cdot\boldsymbol{f}(u')\right.\\ \tag{3} &-\left.\frac{\boldsymbol{I}+\hat{\boldsymbol{t}}(u)\otimes\hat{\boldsymbol{t}}(u)}{\left|s(u')-s(u)\right|}\cdot\boldsymbol{f}(u)\right]\left|\boldsymbol{x}_{u}(u')\right|du', \end{align}

The functions \hat{\boldsymbol{r}}(u,u'), r(u,u'), \hat{\boldsymbol{t}}(u) and s(u) and \left|\boldsymbol{x}_{u}(u')\right| can all be expressed in terms of the centreline coordinates \boldsymbol{x}(u). The vector-valued function \boldsymbol{u}(u) on the l.h.s. is the centreline velocity which can be expressed in terms of the centreline coordinates by a finite difference approximation. Alternatively, it would also be possible to express the system of partial differential equations in terms of centreline velocity \boldsymbol{u}(u) instead of the centreline coordinates \boldsymbol{x}(u) from the get go.

To solve the system of partial equations with fenics, I choose a uniform grid on the unit interval with N grid points and I approximate the force line density and the centreline vector components as continuous piecewise linear functions.

Because equations (2) and (3) are linear in the force line density \boldsymbol{f}(u), their solutions must be a linear function of the vertex values \boldsymbol{f}_{i}. Thus, there exists a matrix \boldsymbol{M} which maps the composite force line density vector \boldsymbol{F} onto the composite centreline velocity vector \boldsymbol{U} such that

\boldsymbol{U} =\boldsymbol{M}\boldsymbol{F},\quad\boldsymbol{M}\in\mathbb{R}^{3N\times3N}

with \boldsymbol{U} and \boldsymbol{F} being defined as

\boldsymbol{U} =(\boldsymbol{u}_{1},\boldsymbol{u}_{2},\ldots,\boldsymbol{u}_{N})\in\mathbb{R}^{3N}, \boldsymbol{F} =(\boldsymbol{f}_{1},\boldsymbol{f}_{2},\ldots\boldsymbol{u}_{N})\in\mathbb{R}^{3N}.

The vectors \boldsymbol{u}_{i} and \boldsymbol{f}_{i} refer to the centreline velocity and the force line density at the i-th vertix.

My first question is, would be in principle possible to pass the equation (1-3) symbolically to fenics using UFL. A potential problem is that the integrand includes divisions by zero if u=u' because r(u,u)=0 and \left|s(u)-s(u)\right|=0. Note that the integral converges nonetheless because terms with opposite signs converge to the same value and cancel in the limes u\rightarrow u'.

As a workaround, I derived the matrix elements of \boldsymbol{M} by hand which allows me to handle the division by zero problem by setting problematic matrix elements to zero. The inverse of \boldsymbol{M} exists, i.e. I can compute

\begin{align} \tag{4} \boldsymbol{F}=\boldsymbol{M}^{-1}\boldsymbol{U} \end{align}

which is the equation I want to include in my larger system of partial differential equations. I would approximate \boldsymbol{M}^{-1} to be a function of the the centreline coordinates \boldsymbol{x}(u) from the previous time step to linearize my equations.

What would be the best way to go about this? If possible, I would like to include equation (4) symbolically for \boldsymbol{M}^{-1} given numerically. If this is not possible, then I could insert the matrix \boldsymbol{M}^{-1} as a block into the bigger matrix representing my linearized discretized system of partial differential equations. To do that, I probably would need a dummy equation to include \boldsymbol{f(u)} into my set of equations. Here is some code to show how I define the function space for \boldsymbol{f}(u) and \boldsymbol{u}(u).

Thanks for taking the time!

from fenics import *
import numpy as np


N = 129 

mesh = UnitIntervalMesh(N - 1)

P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
P3 = MixedElement(3 * [P1])

# inv_M would be computed using the solution for x from the 
# previous time step
inv_M = np.random.rand((3*N, 3*N))

u = TrialFunction(P1) # 3xN
f = TrialFunction(P1) # 3xN
2

Were you able to solve this?