Deep dive: evaluate state using FEniCS¶

In this notebook we solve the model via the finite element method using FEniCSx.

import pathlib
import matplotlib.pyplot as plt
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, default_scalar_type, fem, mesh
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import plotly.graph_objects as go
import plotly.io as pio
from plotly.subplots import make_subplots

The following line allow the plot at the end

pio.renderers.default = "notebook_connected"

Choice of geometry¶

In this notebook we fix some of the parameters arbitrarily.

The final time $T$ is 2 (representing seconds).

T = 2.0

The space geometry (or domain) $\Omega$ is the 1D unit interval, i.e. $(0,1)$.

x0 = 0.0
xL = 1.0

The diffusion coefficients $\kappa_1$ and $\kappa_2$ are constant at 1.

def kappa_1(x):
    return np.ones_like(x[0])

def kappa_2(x):
    return np.ones_like(x[0])

The initial value $y_\circ$ is constant at 5.

def y0(x):
    return 5 * np.ones_like(x[0])

The control vector is $\boldsymbol{\mu} = [3, 3, 3, 3]$.

mu = np.array([3.0, 3.0, 3.0, 3.0])

As the input function $u$ we take: $$u(t) := \begin{cases} -3 & t \le 4/3 \\ +3 & 4/3 < t \le 2 \end{cases}$$

def u_func(t):
    return -3 * (t <= 4 / 3) + 3 * (t > 4 / 3) * (t <= 2.0)

Finite Element discretization¶

To solve this problem numerically we use:

linear Lagrangian elements in space on 201 node points over the interval $\Omega$.

nx = 201
degree = 1
K = 201

implicit Euler method in time in 201 time steps (counting initial and final time)

dt = T / (K - 1)

Solution with FEniCS step by step¶

We first create the discretized space and the function space for the state $(y,q)$:

Omega = mesh.create_interval(
    comm=MPI.COMM_WORLD,
    nx=nx-1,  # this function wants the number of intervals rather than the number of nodes
    points=[x0, xL],
)

P1 = element("Lagrange", Omega.basix_cell(), degree)
V = fem.functionspace(Omega, mixed_element([P1, P1]))

As $V$ is the solution space for the coupled state $(y,q)$, we also define the individual subspaces:

V_sub_0 = V.sub(0)
Vy, Vy_to_V_sub_0 = V_sub_0.collapse()

V_sub_1 = V.sub(1)
Vq, Vq_to_V_sub_1 = V_sub_1.collapse()

The solution at the current iteration $(y^n, q^n)$ is called v and the solution at the previous iteration $(y^{n-1}, q^{n-1})$ is called v_old. Here we initialize the two states in the joint function space.

v = fem.Function(V)  # current solution
v_old = fem.Function(V)  # solution from previous converged step

# Split mixed functions
y, q = ufl.split(v)
y_old, q_old = ufl.split(v_old)

The nonlinearity is defined as follow:

y_ = ufl.variable(y)
q_ = ufl.variable(q)
f = ufl.sqrt(y_) * ufl.sinh(q_)

The initial condition should only be applied to the “y-part” of the state v. Here is how this is done.

v.x.array[:] = 0.0
v.sub(0).interpolate(y0)
v.x.scatter_forward()

As in the weak formulation the Dirichlet boundary conditions only determine the space where the solution is living, here they are collected in a list of conditions bcs:

def on_left_boundary(x):
    return np.isclose(x[0], x0)


boundary_dofs = fem.locate_dofs_geometrical(Vq, on_left_boundary)

bc0 = fem.dirichletbc(default_scalar_type(0), Vq_to_V_sub_1[boundary_dofs], V_sub_1)
bcs = [bc0]

We interpolate the diffusion equations in the corresponding space:

k1 = fem.Function(Vy)
k1.interpolate(kappa_1)

k2 = fem.Function(Vq)
k2.interpolate(kappa_2)

We define the test functions and the nonlinear PDE operator F (where the two weak functions will be summed):

phi, psi = ufl.TestFunctions(V)
b = fem.Constant(Omega, ScalarType(0.0))
Fy = (
    ufl.inner(y, phi) * ufl.dx
    - ufl.inner(y_old, phi) * ufl.dx
    + dt * mu[0] * k1 * ufl.inner(ufl.grad(y), ufl.grad(phi)) * ufl.dx
    - dt * mu[1] * ufl.inner(f, phi) * ufl.dx
)
Fq = (
    mu[2] * k2 * ufl.inner(ufl.grad(q), ufl.grad(psi)) * ufl.dx
    + mu[3] * ufl.inner(f, psi) * ufl.dx
    + b * psi * ufl.ds
)
F = Fy + Fq

Here we define the solver options and we initialize a problem object. This object is going to be solved in each iteration, while updating all the functions.

# Create nonlinear problem
petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(default_real_type).eps) * 1e-2,
    "snes_atol": 0,
    "snes_rtol": 0,
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
}
problem = NonlinearProblem(
    F,
    v,
    bcs=bcs,
    petsc_options_prefix="fenics_FE_",
    petsc_options=petsc_options,
)

We collect data into the matrices Y and Q:

\[\begin{equation*} \mathrm{Y} = \left[ \begin{array}{c|c|c} y^1 & \dots & y^K \end{array} \right] \qquad \mathrm{Q} = \left[ \begin{array}{c|c|c} q^1 & \dots & q^K \end{array} \right] \end{equation*}\]

These matrices will contain the columns of the coordinate arrays (in the FE basis).

Y = np.zeros((nx, K-1))
Q = np.zeros((nx, K-1))

The solver is advanced in time from to until a terminal time is reached

t_array = np.linspace(1, T, K)  # array of time steps
for k, t in enumerate(t_array[1:]):  # we skip time=0
    v_old.x.array[:] = v.x.array
    b.value = u_func(t)
    problem.solve()
    Y[:, k] = v.x.array[Vy_to_V_sub_0]
    Q[:, k] = v.x.array[Vq_to_V_sub_1]

Next, let us plot the states in an interactive plot.

fig = go.Figure()

# Array of x values
x = Omega.geometry.x[:, 0]

# Create subplots
fig = make_subplots(
    rows=1,
    cols=2,
    column_widths=[0.5, 0.5],
)

# Add all $y^k$ and $q^k$ functions to trace
for k, t in enumerate(t_array[1:]):
    fig.add_traces(
        [
            go.Scatter(
                x=x,
                y=Y[:,k],
                visible=False,
                name=f"y(t={t:0.3f})",
            ),
            go.Scatter(
                x=x,
                y=Q[:,k],
                visible=False,
                name=f"q(t={t:0.3f})",
            ),
        ],
        rows=[1,1],
        cols=[1,2],
    )

# Make t=0 visible
fig.data[0].visible = True
fig.data[1].visible = True

# Create and add slider
steps = []
for i in range(K-1):
    step = dict(
        method="update",
        args=[{"visible": [False] * len(fig.data)},
              {"title": "Slider switched to step: " + str(i)}],  # layout attribute
    )
    step["args"][0]["visible"][2 * i] = True  # Toggle i'th trace to "visible"
    step["args"][0]["visible"][2 * i + 1] = True  # Toggle i'th trace to "visible"
    steps.append(step)

sliders = [dict(
    active=0,
    currentvalue={"prefix": "Time: "},
    pad={"t": 50},
    steps=steps
)]  

# Add slider to figure
fig.update_layout(
    sliders=sliders,
    yaxis={"range": [0, 7]},
    legend={"x": 0.5, "y":1.3, "xanchor": "center", "orientation":"h"},
)

# Show final interactive plot
fig.show()

We see that the state $y$ is “inflated” by the positive input $u$. When $u$ turns negative, $q$ changes to a more “decrescent” shape. This makes the state $y$ “deflate”.