Optimal Control of Poisson Equation

[MartinK-99]

Hey there, I'm working on an optimal control problem which minimizes some cost functional $J(y,u)$ under the constraint $\Delta y = u$ in the domain $Ω = (0,1) \times (0,1)$ with $y = 0$ on the boundary. Here $u$ denotes the control and $y$ the state variable. When discretizing the constraint equation using FEM, one typically defines piecewise linear basis functions $\lbrace \phi_i,\ldots,\phi_l \rbrace$ for the triangulation. The control variables are assumed to be constant on each triangle for which one defines trial functions $e_i$, $i=1,\ldots,m$ which are equal to 1 within the triangle, 0 elsewhere. Therefore,
$y(x) = \sum_{i=1}^l y_i \phi_i(x), \quad u(x) = \sum_{i=1}^m u_i e_i(x)$
and notice, that the dimension of $y$ and $u$ are not the same! Inserting this in the weak formulation for the poisson equation leads to
$\int_\Omega \sum_{i=1}^l y_i \nabla \phi_i \cdot \nabla \phi_j dx = \int_\Omega \sum_{i=1}^m u_i e_i \phi_j dx$
such that in total the linear equation system is
$K y = Bu$
with stiffness matrix $K$ and another matrix $B \in \mathbb{R}^{l\times m}$ of elements $b_{ij} = \int_\Omega \phi_i e_j dx$ [1].

My question is, how do I define the trial functions $e_i$ using Ferrite and how do I assemble the matrix $B$? My attempts so far have failed at defining the functions and from there I'm not sure how to even use the DofHandler to create a sparsity pattern for a non square matrix.

I'm grateful for your help and suggestions.

[1] Tröltzsch, Fredi. Optimal Control of Partial Differential Equations : Theory, Methods, and Applications. American Mathematical Society, 2010.

[fredrikekre]

how do I define the trial functions e i using Ferrite

If I understand correctly you should be able to use 0th order discontinuous Lagrange elements (DiscontinuousLagrange{RefTriangle, 0}() for the triangle).

[MartinK-99]

I did not know that DiscontinuousLagrange interpolation exists. Thank you!
My reference is the example "Heat Equation" from your documentation. How do I continue with the DofHandler to then create the sparsity pattern for my matrix $B$? Or do I need to use the MixedDofHandler because I now mix Lagrange and DiscontinuousLagrange interpolation?

[fredrikekre]

The mixed in MixedDofHandler is if you have different element types, for example triangles and quadrilaterals, so you don't need that.

You can add two fields to the DofHandler and you would get a block-like matrix like:

$$A_{yy} \ A_{yu}$$ $$A_{uy} \ A_{uu}$$

where $K = A_{yy}$ and $M = A_{yu}$.

There are a lot of things related to this that have changed from the last release to the master branch, and this is quite a lot nicer on the master branch so if you can use that I think that would make sense for your case. I include code for both the latest release (0.3.14) and master below.

Ferrite 0.3.14

using Ferrite

grid = generate_grid(Triangle, (2, 2))

# Interpolations
ipy = Lagrange{2, RefTetrahedron, 1}()
ipu = DiscontinuousLagrange{2, RefTetrahedron, 0}()
ipg = Lagrange{2, RefTetrahedron, 1}() # Geometric interpolation

# CellValues
qr = QuadratureRule{2, RefTetrahedron}(2)
cvy = CellScalarValues(qr, ipy, ipg)
cvu = CellScalarValues(qr, ipu, ipg)

# DofHandler
dh = DofHandler(grid)
add!(dh, :y, ipy)
add!(dh, :u, ipu)
close!(dh)

# Renumber so we have a global block matrix
renumber!(dh, DofOrder.FieldWise())

# Hack to compute global ndofs for each field
ndofs_u = getncells(grid)
ndofs_y = ndofs(dh) - ndofs_u
global_dofs_y = 1:ndofs_y
global_dofs_u = (ndofs_y + 1):ndofs(dh)

# Global matrix
A = create_sparsity_pattern(dh)

# Assembly
function assemble(A, dh, cvy, cvu)
    # This dummy rhs should not be necessary
    dummy_rhs = zeros(ndofs(dh))
    assembler = start_assemble(A, rhs)
    Ae = zeros(ndofs_per_cell(dh), ndofs_per_cell(dh))
    for cell in CellIterator(dh)
        # Reset local matrix
        fill!(Ae, 0)
        # Update cell values for current cell
        reinit!(cvy, cell)
        reinit!(cvu, cell)
        # Local dof ranges
        yr = dof_range(dh, :y)
        ur = dof_range(dh, :u)
        # Cell assembly
        for qp in 1:getnquadpoints(cvu)
            # Quadrature weight
            dΩ = getdetJdV(cvu, qp)
            # Assemble parts with y-test function
            for (i, I) in pairs(yr)
                # Assemble yy block
                ∇ϕi = shape_gradient(cvy, qp, i)
                for (j, J) in pairs(yr)
                    ∇ϕj = shape_gradient(cvy, qp, j)
                    Ae[I, J] += (∇ϕi ⋅ ∇ϕj) * dΩ
                end
                # Assemble yu block
                ϕi = shape_value(cvy, qp, i)
                for (j, J) in pairs(ur)
                    ej = shape_value(cvu, qp, j)
                    Ae[I, J] += (ϕi * ej) * dΩ
                end
            end

            # Assemble parts with u-test function
            # Same as above but just exchanging indices etc, but if
            # you only care about the upper two blocks of the matrix
            # you don't need this.
            # [...]

            # Assemble element contribution to global matrix
            dofs = celldofs(cell)
            assemble!(assembler, dofs, Ae)
        end
    end
    return A
end

# Do the assembly
A = assemble(A, dh, cvy, cvu)

# Extract K and M from A
K = A[global_dofs_y, global_dofs_y]
M = A[global_dofs_y, global_dofs_u]

Ferrite master

using Ferrite, BlockArrays

grid = generate_grid(Triangle, (2, 2))

# Interpolations
ipy = Lagrange{RefTriangle, 1}()
ipu = DiscontinuousLagrange{RefTriangle, 0}()
ipg = Lagrange{RefTriangle, 1}() # Geometric interpolation

# CellValues
qr = QuadratureRule{RefTriangle}(2)
cvy = CellValues(qr, ipy, ipg)
cvu = CellValues(qr, ipu, ipg)

# DofHandler
dh = DofHandler(grid)
add!(dh, :y, ipy)
add!(dh, :u, ipu)
close!(dh)

# Renumber so we have a global block matrix
renumber!(dh, DofOrder.FieldWise())

# Hack to compute global ndofs for each field
ndofs_u = getncells(grid)
ndofs_y = ndofs(dh) - ndofs_u
global_dofs_y = 1:ndofs_y
global_dofs_u = (ndofs_y + 1):ndofs(dh)

# Global block matrix
sp = BlockSparsityPattern([ndofs_y, ndofs_u])
add_sparsity_entries!(sp, dh)
A = allocate_matrix(BlockMatrix, sp)

# [Assembly function equivalent from the [email protected] example]

# Do the assembly
A = assemble(A, dh, cvy, cvu)

# Extract K and M from A
K = @view A[Block(1), Block(1)]
M = @view A[Block(1), Block(2)]

In particular, on the master branch you can assemble as a block matrix directly which is pretty nice.

[MartinK-99]

I'll look into this. Thank you very much for you support!