Programming Homework 05

Submission Deadline: 08.07.2026 - 23:55

The homework solution has to be handed in as a group solution via Moodle.
The notebooks must run without errors in our environment on the RWTH-Cluster. See here for more information on how to get started.

Participant list

Please add your names and student ID (Matrikelnummer) here:

Name, ID
Name, ID
Name, ID
Name, ID

Learning goals

Understand the basic structure of a mixed finite element method for Darcy flow
Investigate the role of permeability
Understand the meaning of the Eigenstructure of the permeability matrix

Before you start

Make sure that you update the library folder in your CMM directory. Download link: library.tar.gz
Make sure that you update your dolfinx_latest.sif container image.

You can use the following command on the RWTH cluster/Linux:

    wget https://mbd.pages.rwth-aachen.de/courses/cmm/content/exercises/notebooks/homework05.out.ipynb
    rm -rf library
    wget https://mbd.pages.rwth-aachen.de/courses/cmm/content/exercises/notebooks/library.tar.gz
    tar -xvzf library.tar.gz
    rm library.tar.gz
    rm dolfinx_latest.sif
    apptainer pull oras://registry.git.rwth-aachen.de/mbd/courses/containers/dolfinx:latest

Todos

Make sure that all tests in the notebook run without errors. The following tasks need to be completed:

update the participant list
Task 1 — implement the piecewise permeability k_ufl for the multilayer problem (P1)
Task 2 — implement the Darcy term d(u, v) in the weak form of the multilayer problem (P1)
Task 3 — implement the analytical solution u_exact_multilayer for the multilayer problem (P1)
Task 4 — compute the permeability matrix K from eigenvectors and eigenvalues (P2)
Task 5 — implement the Darcy term d_anisotropic(u, v) for the anisotropic problem (P2)

Imports

import os
import sys
import numpy as np

from mpi4py import MPI
import pyvista
import tqdm.autonotebook

import dolfinx.plot
from  dolfinx.fem import petsc
from petsc4py import PETSc
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, la

from dolfinx.fem import (
    Constant,
    Function,
    dirichletbc,
    assemble_scalar,
    extract_function_spaces,
    form,
    functionspace,
    locate_dofs_topological,
)

from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile, VTXWriter
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner

import library.plot
import library.helpers
import library.boundary 

CMM_DIR = os.getcwd()
os.environ["PYTHONPATH"] = f"{CMM_DIR}:{os.environ.get('PYTHONPATH', '')}"

# overload print function to print only on rank 0
def print_rank0(*args, **kwargs):
    if MPI.COMM_WORLD.rank == 0:
        print(*args, **kwargs)

output_dir = ""
if MPI.COMM_WORLD.rank == 0:
    output_dir = library.helpers.make_unique_dir(os.path.join(CMM_DIR, 'results_hw05'))

# mpi share output_dir to all ranks
output_dir = MPI.COMM_WORLD.bcast(output_dir, root=0)

# Check if we are running interactively (e.g. Jupyter notebook with visual frontend)
# or as a script: https://stackoverflow.com/a/64523765
interactive_mode_is_on = hasattr(sys, 'ps1')

if os.getenv('DEVCONTAINER') and interactive_mode_is_on:
    pyvista.start_xvfb()

Mixed finite elements for Darcian flow

We start by writing the Darcy problem as

\[ \begin{aligned} - \mu\,\mathbf{K}^{-1} \mathbf{u} - \boldsymbol{\nabla} p + \mathbf{f} &= \mathbf{0}, \text{ on } \Omega, \\ \boldsymbol{\nabla}\cdot\mathbf{u} &= \psi \text{ on } \Omega, \end{aligned} \]

with boundary conditions

\[ \begin{aligned} \mathbf{u} \cdot \mathbf{n} = u_n \text{ on } \Gamma_E, \\ p\,\mathbf{n} = \mathbf{g} \text{ on } \Gamma_N, \end{aligned} \]

Here, \(\mathbf{u} \in \mathbb{R}^d\) is the Darcy velocity, \(p \in \mathbb{R}\) is the pressure, \(\mathbf{f} \in \mathbb{R}^d\) is a body force vector.

In the simplest case where permeability is a scalar field with values \(k\), the permeability matrix \(\mathbf{K}\) can be written as \(\mathbf{K} = k \, \mathbf{I}\), where \(k\) is the scalar permeability and \(\mathbf{I}\) is the identity matrix.

In this case, we can also write the first term as \(- \frac{\mu}{k}\, \mathbf{u}\).

The weak formulation

Here, we consider the following dual mixed formulation of the problem (see for example Eqns. (9) and (10) of www.doi.org/10.1016/j.jhydrol.2011.12.024):

The weak formulation of the Darcy problem reads as follows: Find \((\mathbf{u}, p) \in \mathbb{V} \times \mathbb{Q}\) such that

\[ \begin{aligned} \int_{\Omega} \mu \, \mathbf{K}^{-1} \, \mathbf{u} \cdot \mathbf{v} \, \text{d}x - \int_{\Omega} p \, \nabla \cdot \mathbf{v} \, \text{d}x &= \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, \text{d}x + \int_{\Gamma^N} \mathbf{g} \cdot \mathbf{v} \, \text{d}S \quad \forall \, \mathbf{v} \in \mathbb{V}, \\ -\int_{\Omega} \nabla \cdot \mathbf{u} \, q \, \text{d}x &= -\int_{\Omega} \psi \, q \, \text{d}x \quad \forall \, q \in \mathbb{Q}. \end{aligned} \]

which can be written as

\[ \begin{aligned} d(\mathbf{u}, \mathbf{v}) + b(\mathbf{v}, p) &= L(\mathbf{v}) \quad \forall \, \mathbf{v} \in \mathbb{V}, \\ b(\mathbf{u}, q) &= S(q) \quad \forall \, q \in \mathbb{Q} \\ \end{aligned} \]

with

\[ \begin{aligned} d \,:\, \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\,, \quad & d(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \mu \, \mathbf{K}^{-1} \, \mathbf{u} \cdot \mathbf{v} \, \text{d}x \\ b \,:\, \mathbb{R}^d \times \mathbb{R} \to \mathbb{R}\,, \quad & b(\mathbf{u}, q) = -\int_{\Omega} \nabla \cdot \mathbf{u} \, q \, \text{d}x, \\ L \,:\, \mathbb{R}^d \to \mathbb{R}\,, \quad & L(\mathbf{v}) = \int_{\Omega} \mathbf{f} \cdot \mathbf{v} \, \text{d}x + \int_{\Gamma^N} \mathbf{g} \cdot \mathbf{v} \,\text{d}S, \\ S \,:\, \mathbb{R} \to \mathbb{R}\,, \quad & S(q) = -\int_{\Omega} \psi \, q \, \text{d}x. \end{aligned} \]

The discrete computational domain

We start by creating a discrete approximation of the computational domain \(\Omega\). Here we use triangular elements to discretize the domain.

mesh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [30, 30], CellType.triangle
)

The finite element approximation space

First we choose the finite-element approximation space, i.e. the space in which we search for an approximate solution of the problem.

In the following we set up the discrete approximation spaces \(\mathbb{V}_h\) and \(\mathbb{Q}_h\) for the velocity and pressure fields, respectively.

from ufl import MixedFunctionSpace

k = 0
# Function spaces for the velocity and for the pressure
V_h = fem.functionspace(mesh, ("Raviart-Thomas", k + 1))
Q_h = fem.functionspace(mesh, ("Discontinuous Lagrange", k))
VQ = MixedFunctionSpace(V_h, Q_h)

# Function space for visualising the velocity field
gdim = mesh.geometry.dim
W = fem.functionspace(mesh, ("Discontinuous Lagrange", k + 1, (gdim,)))

u_vis = Function(W, name="u_visualization")
u_vis.name = "darcy_velocity"

from ufl import (
    FacetNormal,
    Identity,
    TestFunction,
    TrialFunction,
    div,
    dot,
    ds,
    dS,
    dx,
    inner,
    lhs,
    nabla_grad,
    avg,
    rhs,
    sym,
    CellDiameter,
    outer,
)

# Mesh element surface normal:
n = FacetNormal(mesh)
# External body force source term:
f_body = Constant(mesh, PETSc.ScalarType((0, 0)))
# Dynamic viscosity:
mu = Constant(mesh, PETSc.ScalarType(1e-8))

P1 — Multilayer Darcy Flow

We consider the multilayer problem from Table 2 of www.doi.org/10.1016/j.jhydrol.2011.12.024. The unit square domain is \(\Omega = [0, 1]^2\) and the permeability field is given by 5 equally sized vertical layers with different scalar permeability values.

The values of the permeability field are given as \[ k(y) = \begin{cases} 8.0 & \text{if } y < 0.2 \\ 3.0 & \text{if } 0.2 <= y < 0.4 \\ 0.5 & \text{if } 0.4 <= y < 0.6 \\ 5.0 & \text{if } 0.6 <= y < 0.8 \\ 1.0 & \text{if } y >= 0.8 \end{cases} \quad \cdot 10^{-13} \text{m}^2 \]

The flow is driven by a pressure difference \(\Delta p = 10^5\,\text{Pa}\) between the left inlet and right outlet of the domain (\(\Gamma^{N,\text{in}}\) and \(\Gamma^{N,\text{out}}\)).

Although the domain is 2D, the solution is effectively 1D by symmetry: the pressure gradient acts purely in the \(x\)-direction and the permeability varies only with \(y\), so \(u_y = 0\) everywhere and only \(u_x\) is non-trivial.

Analytical solution hint: For pressure-driven flow in each layer, Darcy’s law gives \(u_x^{(i)} = \frac{k_i}{\mu}\,\frac{\Delta p}{L}\). With \(\mu = 10^{-8}\,\text{Pa}\cdot\text{s}\), \(\Delta p = 10^5\,\text{Pa}\), and \(L = 1\,\text{m}\), the expected velocity in layer \(i\) is \(k_i \times 10^{-13} \times 10^5 / 10^{-8} = k_i\,\text{m/s}\), i.e. the prefactors 8.0, 3.0, 0.5, 5.0, 1.0 are the expected \(u_x\) values directly.

1. Permeability field

k1 = Constant(mesh, PETSc.ScalarType(1e-13))
k2 = Constant(mesh, PETSc.ScalarType(5e-13))
k3 = Constant(mesh, PETSc.ScalarType(0.5e-13))
k4 = Constant(mesh, PETSc.ScalarType(3e-13))
k5 = Constant(mesh, PETSc.ScalarType(8e-13))

STUDENT TODO — Task 1: Piecewise permeability function

In the following cell, we implement the piecewise permeability function given above. The numerical values are given in the variables k_i with \(i \in \{1, 2, 3, 4, 5\}\). It remains to write a UFL expression that returns the correct scalar permeability for a given \(y\) coordinate.

Hint 1: We can get the \(x\) or \(y\) coordinate of a UFL spatial coordinate x_ufl in the domain by using x_ufl[0] and x_ufl[1] respectively.

Hint 2: We can use conditional statements to implement piecewise functions in UFL. The syntax is conditional(condition, true_value, false_value), where condition is a boolean expression, true_value is the value returned if the condition is true, and false_value is the value returned if the condition is false.

Hint 3: conditional statements can be nested, e.g. the true_value and false_value can again be conditionals.

For example

from ufl import conditional
from ufl import SpatialCoordinate

x_ufl = SpatialCoordinate(mesh)
expression = conditional(x_ufl[0] < 0.5, 1.0, 0.0)

returns 1.0 if the \(x\) coordinate is less than 0.5, and 0.0 otherwise.

Write an expression that returns the correct stepwise permeability defined above as a UFL expression stored in k_ufl:

from ufl import conditional, SpatialCoordinate

x_ufl = SpatialCoordinate(mesh)
k_ufl = 0.0 # YOURCODE HERE

2. Pressure boundary conditions

As mentioned earlier, the flow is driven by a pressure gradient, imposed through the pressure boundary conditions on the left and right boundaries of the domain \(\Gamma^{N,in}\) and \(\Gamma^{N,out}\), respectively. The pressure difference is given by \(\Delta p = 10^5\,Pa\).

The pressure boundary term appears as a natural boundary condition in the weak form of the problem, which can be seen by performing integration by parts of the pressure term:

\[ -\int_{\Omega} \nabla p \cdot \mathbf{v} \, \text{d}x = \int_{\Omega} p \, \nabla \cdot \mathbf{v} \, \text{d}x - \int_{\Gamma_N} \underline{p} \, \mathbf{v} \cdot \mathbf{n} \, \text{d}S \quad . \]

We can prescribe pressure boundary conditions by inserting them into the right-hand side boundary term for the underlined pressure \(\underline{p}\) in the above equation.

is_left = conditional(x_ufl[0] < 1e-6, 1.0, 0.0)
is_right = conditional(1.0 - x_ufl[0] < 1e-6, 1.0, 0.0)
is_natural_pressure_boundary =  is_left + is_right

pressure_left = Constant(mesh, PETSc.ScalarType(1e5))
pressure_right = Constant(mesh, PETSc.ScalarType(0))

pressure_value_expression_multilayer = is_left*pressure_left + is_right*pressure_right

3. Weak form

STUDENT TODO — Task 2: Darcy term in the weak form

Given our definition of the permeability function k_ufl above, implement the weak form Darcy term

\[d(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \mu \, \mathbf{K}^{-1} \, \mathbf{u} \cdot \mathbf{v} \, \text{d}x\]

in the following cell.

d = lambda u, v: 0.0 # YOURCODE HERE

The rest of the weak form is then assembled below as usual:

from ufl import extract_blocks

b = lambda v, p: -inner(p, div(v)) * dx
L = lambda f, psi, v: inner(f, v) * dx - inner(psi, q_h) * dx
pressure_bc_term_multilayer = lambda v: - inner(v, n) * pressure_value_expression_multilayer * is_natural_pressure_boundary * ds 

# Define variational problem
u_h, p_h = ufl.TrialFunctions(VQ)
v_h, q_h = ufl.TestFunctions(VQ)

f = Function(W) # Initialize zero body force source term
psi = Function(Q_h) # Initialize zero mass source term

weak_form_lhs = d(u_h,v_h) + b(v_h, p_h) + b(u_h, q_h) 
weak_form_rhs = L(f, psi, v_h) + pressure_bc_term_multilayer(v_h) 
# Add this term just to make the extract_blocks function work on the correct block structure
weak_form_rhs += inner(fem.Constant(mesh, default_real_type(0.0)), q_h) * dx

weak_form_darcy_multilayer_lhs = fem.form(extract_blocks(weak_form_lhs))
weak_form_darcy_multilayer_rhs = fem.form(extract_blocks(weak_form_rhs))

4. Boundary conditions and solver

The top and bottom boundaries are no-outflow boundaries. In the mixed dual formulation, the essential boundary condition is on the normal component of the velocity field \(\mathbf{u} \cdot \mathbf{n}\), which we set to zero on the top and bottom boundaries \(\Gamma^{E}\).

def no_outflow_boundary_marker(x):
    return np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)

The following syntax is similar to the one used before. The actual implementation has one quirk, which is due to the way FEniCSx handles Dirichlet boundary conditions for some particular finite element spaces, such as the Raviart-Thomas space used here.

It is a bit counterintuitive because the velocity \(u_D\) is a 2D zero vector but only the normal component is set to zero. This is because for the Raviart-Thomas element, the degree of freedom corresponds to the normal component of the velocity field \(\mathbf{u} \cdot \mathbf{n}\) (and not to approximate velocity values as in previous examples using Taylor-Hood elements). We forego these technicalities here but if one is interested in details you can start here.

 # dimension of the mesh facets (1D lines for a 2D mesh)
fdim = mesh.topology.dim - 1

u_D = Function(V_h, name="u_D")
facets_no_inflow = locate_entities_boundary(mesh, fdim, no_outflow_boundary_marker)
sub_dofs = dolfinx.fem.locate_dofs_topological(V_h, fdim, facets_no_inflow)
no_normal_flow_bc = dolfinx.fem.dirichletbc(u_D, sub_dofs)
velocity_bc_multilayer = [no_normal_flow_bc]

Our solver function

Similar to previous exercises, the following function takes the weak form and boundary conditions and solves the linear system. For convenience, it returns the sub-components of the solution vector belonging to the velocity and pressure at the degrees of freedom separately.

from pathlib import Path

def solve_mixed_darcy_problem(simulation_name, 
                              weak_form_lhs, 
                              weak_form_rhs, 
                              essential_boundary_conditions=None):
    """
    Assemble and solve the linear system resulting from the weak form of the mixed Darcy problem.
    """ 

    folder = Path(output_dir +"/" +simulation_name)
    folder.mkdir(exist_ok=True, parents=True)
    vtk_file_velocity_abs_path = folder / str("darcy_velocity.pvd")
    vtk_file_pressure_abs_path = folder / str("pressure.pvd")

    if essential_boundary_conditions is None:
        A = assemble_matrix_block(weak_form_lhs)
        b = assemble_vector_block(weak_form_rhs, weak_form_lhs)
    else:
        A = assemble_matrix_block(weak_form_lhs, bcs=essential_boundary_conditions)
        b = assemble_vector_block(weak_form_rhs, weak_form_lhs, bcs=essential_boundary_conditions)

    A.assemble()

    # Create and configure solver
    ksp = PETSc.KSP().create(mesh.comm)  # type: ignore
    ksp.setOperators(A)
    ksp.setType("preonly")
    ksp.getPC().setType("lu")
    ksp.getPC().setFactorSolverType("mumps")
    opts = PETSc.Options()  # type: ignore
    opts["mat_mumps_icntl_14"] = 80  # Increase MUMPS working memory
    opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix (pressure nullspace)
    opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix (pressure nullspace)
    opts["ksp_error_if_not_converged"] = 1
    ksp.setFromOptions()

    x = A.createVecRight()
    try:
        ksp.solve(b, x)
    except PETSc.Error as e:  # type: ignore
        if e.ierr == 92:
            print("The required PETSc solver/preconditioner is not available. Exiting.")
            print(e)
            exit(0)
        else:
            raise e

    # Split the solution
    u_sol = fem.Function(V_h)
    p_sol = fem.Function(Q_h)
    
    p_sol.name = "p"
    offset = V_h.dofmap.index_map.size_local * V_h.dofmap.index_map_bs
    u_sol.x.array[:offset] = x.array_r[:offset]
    u_sol.x.scatter_forward()
    p_sol.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
    p_sol.x.scatter_forward()

    # Subtract the average of the pressure if pressure is only determined up to a constant
    # p_h.x.array[:] -= domain_average(mesh, p_h)

    u_vis = fem.Function(W)
    u_vis.name = "u"
    u_vis.interpolate(u_sol)

    vtk_writer_u = dolfinx.io.VTKFile(mesh.comm, vtk_file_velocity_abs_path, "w+")
    vtk_writer_p = dolfinx.io.VTKFile(mesh.comm, vtk_file_pressure_abs_path, "w+")
    vtk_writer_u.write_function([u_vis], t=0.0)
    vtk_writer_p.write_function([p_sol], t=0.0)

    return u_sol, p_sol, u_vis

5. Solution

Solve the steady multilayer problem

u_sol, p_sol, u_vis = solve_mixed_darcy_problem("multilayer_flow",
                                                weak_form_lhs=weak_form_darcy_multilayer_lhs, 
                                                weak_form_rhs=weak_form_darcy_multilayer_rhs, 
                                                essential_boundary_conditions=velocity_bc_multilayer)

if interactive_mode_is_on:
    library.plot.contour_plot(u_vis, geometry_dim=2, create_geometry_from_space=W)
    library.plot.contour_plot(p_sol, geometry_dim=2, create_geometry_from_space=W)

6. Verification

STUDENT TODO — Task 3: Analytical solution

To verify that the numerical solution is correct, we compare it with the analytical solution of the multilayer problem.

Implement the UFL expression that returns a vector evaluating to the analytical solution at a given point \((x,y)\) in the domain \(\Omega = [0,1]^2\).

Hint 1: We can use the as_vector([analytical_velocity_x, analytical_velocity_y]) method to create a vector of UFL expressions. We need to implement the expressions analytical_velocity_x and analytical_velocity_y for the velocity components \(u_x\) and \(u_y\).

Hint 2: Recall that we can use UFL conditional expressions to implement piecewise functions of the \(x\) coordinate x_ufl[0] and \(y\) coordinate x_ufl[1]. The conditional expression has the form conditional(condition, true_value, false_value), where condition is a boolean expression that evaluates to True or False, and true_value and false_value are the values returned accordingly.

from ufl import as_vector

ux_exact = 0.0 # REPLACE WITH YOURCODE
uy_exact = 0.0 # REPLACE WITH YOURCODE
u_exact_multilayer = as_vector([ux_exact, uy_exact])

Verification

The cell below computes the \(L^2\) error between the numerical and analytical solutions and checks that it is below the expected tolerance.

from dolfinx.fem import assemble_scalar, form

u_ex = u_exact_multilayer
L2_error = form(dot(u_sol - u_ex, u_sol - u_ex) * dx)
error_L2 = np.sqrt(mesh.comm.allreduce(assemble_scalar(L2_error), op=MPI.SUM))

def test_multilayer_solution():
    
    if (mesh.comm.rank == 0):
        if np.isclose(error_L2, 0.0, atol=1e-10):
            print_rank0("Steady multilayer test passed.")
        else:
            print_rank0("Steady multilayer test failed.")
            print_rank0(f"Expected error_L2: 0.0, got {error_L2}")

test_multilayer_solution()

Steady multilayer test passed.

P2 — Anisotropic Permeability

In P1, permeability was a scalar field. We now extend to a case where the permeability is a tensor. For the 2D flow example, the permeability will be a \(2 \times 2\) matrix. In principle this can again be a function of the spatial coordinates, but for simplicity we use a constant permeability tensor.

Note that P1 is a special case of P2: scalar permeability \(k\) corresponds to the isotropic permeability matrix \(\mathbf{K} = k\,\mathbf{I}\), where \(\mathbf{I}\) is the identity matrix. In P2 we generalise to a full \(2\times2\) matrix with off-diagonal entries.

1. Eigenvectors and eigenvalues

In particular, we are trying to understand how the permeability “acts” on the velocity field.

Let a basis of eigenvectors \(\mathbf{e}_1, \mathbf{e}_2\) of the permeability tensor \(\mathbf{K}\) be given

\[ \mathbf{e}_1 = \begin{pmatrix} 0.90630779 \\ 0.42261826 \end{pmatrix}, \quad \mathbf{e}_2 = \begin{pmatrix} -0.42261826 \\ 0.90630779 \end{pmatrix} \]

with corresponding eigenvalues \(\lambda_1=3\cdot 10^{-1}\) and \(\lambda_2 = 1\cdot 10^{-4}\).

# Matrix where each column is an eigenvector
eigen_vectors = np.array([[0.90630779, -0.42261826], [0.42261826, 0.90630779]])
eigen_values = np.array([3.0e-1, 1.0e-4])

We can quickly visualize the normalized eigenvectors:

eig_vec1 = eigen_vectors[:, 0]
eig_vec2 = eigen_vectors[:, 1]
library.plot.plot_eigenvectors_2d(eig_vec1, eig_vec2, scale=0.5)

2. Permeability matrix

STUDENT TODO — Task 4: Permeability matrix

From the given eigenvectors and eigenvalues, compute the permeability matrix \(\mathbf{K}\) as a NumPy array. We look for a matrix \(\mathbf{K}\) that has the eigenvectors \(\mathbf{e}_1\) and \(\mathbf{e}_2\) as columns and the eigenvalues \(\lambda_1\) and \(\lambda_2\).

Note: The matrix should be programmatically computed and not hardcoded, such that our way of obtaining \(\mathbf{K}\) is traceable in the code.

# YOUR CODE
K = np.array([[0, 0], [0, 0]]) # REPLACE WITH YOURCODE

Perform the following check to see if our matrix is correct:

from numpy import linalg as LA
def test_permeability_matrix(your_K):
    """ 
    Check that we correctly computed the permeability matrix K 
    by going backwards and check if your matrix really has the eigenvalues and eigenvectors we expect.
    """
    your_matrix_eigen_values, your_matrix_eigen_vectors = LA.eig(your_K)

    assert np.allclose(your_matrix_eigen_values, eigen_values), "Eigenvalues do not match!"
    assert np.allclose(your_matrix_eigen_vectors, eigen_vectors), "Eigenvectors do not match!"

    print("TEST PASSED: Permeability matrix K is correct!")

print("Permeability Matrix:")
print(K)
test_permeability_matrix(K)

Permeability Matrix:
[[0.246436   0.11486836]
 [0.11486836 0.053664  ]]
TEST PASSED: Permeability matrix K is correct!

from ufl import as_matrix

k11_ufl = Constant(mesh, PETSc.ScalarType(K[0, 0]))
k12_ufl = Constant(mesh, PETSc.ScalarType(K[0, 1]))
k21_ufl = Constant(mesh, PETSc.ScalarType(K[1, 0]))
k22_ufl = Constant(mesh, PETSc.ScalarType(K[1, 1]))

K_ufl = as_matrix([[k11_ufl, k12_ufl], [k21_ufl, k22_ufl]])

def inverse_2x2_ufl(matrix):
    matrix_ufl = as_matrix(matrix)
    assert matrix_ufl.ufl_shape == (2, 2), "inverse_2x2_ufl only works on 2x2 matrices!"

    a = matrix_ufl[0, 0]
    b = matrix_ufl[0, 1]
    c = matrix_ufl[1, 0]
    d = matrix_ufl[1, 1]

    return (1 / (a * d - b * c)) * as_matrix([[d, -b], [-c, a]])

3. Weak form

STUDENT TODO — Task 5: Darcy term for matrix permeability

We now generalise the weak form Darcy term to the case of a matrix permeability. Given the matrix K_ufl implemented above, implement the weak form Darcy term

\[ d_{\text{anisotropic}}(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \mu \, \mathbf{K}^{-1} \, \mathbf{u} \cdot \mathbf{v} \, \text{d}x \]

in the following cell.

We can use the inverse_2x2_ufl function above to compute the inverse of K_ufl.

K_inv = 0.0 # YOURCODE
d_anisotropic = lambda u, v: 0.0 # YOURCODE

For this numerical example we work in dimensionless units and rescale pressure and viscosity to \(\mathcal{O}(1)\) to avoid numerical issues:

pressure_left.value = 1.0
pressure_right.value = 0.0
mu = Constant(mesh, PETSc.ScalarType(1.0))

# LHS
weak_form_left = d_anisotropic(u_h,v_h) + b(v_h, p_h) + b(u_h, q_h)
weak_form_darcy_anisotropic_lhs = fem.form(extract_blocks(weak_form_left))
weak_form_darcy_anisotropic_rhs = fem.form(extract_blocks(weak_form_rhs))

4. Solution

Let’s compute the numerical solution:


u_sol_anisotropic, p_sol_anisotropic, u_vis_anisotropic = solve_mixed_darcy_problem("example_anisotropic",
                                                                                    weak_form_lhs=weak_form_darcy_anisotropic_lhs, 
                                                                                    weak_form_rhs=weak_form_darcy_multilayer_rhs,
                                                                                    essential_boundary_conditions=velocity_bc_multilayer)

if interactive_mode_is_on:
    library.plot.contour_plot_2d_with_eigenvectors(u_vis_anisotropic,
                                               geometry_dim=2,
                                               eigen_vectors=eigen_vectors,
                                               vector_loc = [0.5, 0.5],
                                               scale_glyph=0.15,
                                               create_geometry_from_space=W
                                               )
    library.plot.contour_plot(p_sol_anisotropic, geometry_dim=2, create_geometry_from_space=W)

5. Verification

The eigenvalue ratio \(\lambda_1/\lambda_2 = 3000\) means flow is 3000 times easier along \(\mathbf{e}_1\) than along \(\mathbf{e}_2\). This is why the resulting Darcy velocity is strongly aligned with \(\mathbf{e}_1\), even though the imposed pressure gradient drives flow in the \(x\)-direction.

The following test checks if our numerical solution DOF vector matches the L2 norm of the correct solution, so you can double check if you implemented the last weak form term correctly:

Verification

The cell below checks that the \(L^2\) norms of the velocity and pressure DOF vectors match the expected values for the correct anisotropic solution.

def test_anisotropic_solution():
    l2_norm_velocity = np.linalg.norm(u_sol_anisotropic.x.array)
    l2_norm_pressure = np.linalg.norm(p_sol_anisotropic.x.array)
    
    if (mesh.comm.rank == 0):
        if (
            np.isclose(l2_norm_pressure, 30.78219412159688, atol=1e-10) 
            and 
            np.isclose(l2_norm_velocity, 0.11438640312926375, atol=1e-10)
        ):
            print_rank0("Matrix permeability case test passed.")
        else:
            print_rank0("Matrix permeability case test failed. The solution seems to be incorrect :(")

test_anisotropic_solution()

Matrix permeability case test passed.

Prepare for your presentation

During the presentation, be prepared to answer the following questions/tasks. We will ask 1-2 questions per person.

For questions that specifically ask for equations, please be ready to write them down by hand without additional notes

What was the exercise about? Give a short summary and name the tasks that you had to solve.
Show how the pressure boundary condition is a natural boundary condition for the mixed Darcy problem by performing integration by parts
Sketch the setup of the 1D multilayer problem and explain the boundary conditions. Why do we expect a purely 1D flow in this example?
Give the general weak form of the Darcy term (the term involving the permeability matrix \(\mathbf{K}\)) in the mixed dual formulation shown in this exercise. How can we simplify this for a scalar permeability? How is the case of a scalar permeability related to the general permeability matrix?
For the anisotropic permeability case, (roughly) sketch the normalized eigenvectors of the permeability matrix \(\mathbf{K}\). How are the eigenvectors related to the resulting flow field?
What is the role of the eigenvalues of the permeability matrix \(\mathbf{K}\) in this context? Explain using the Eigenvalues and Eigenvectors of the permeability matrix \(\mathbf{K}\) that we computed above.