multiphenicsx

https://multiphenics.github.io/
Admin
multiphenicsx is a Python library providing tools for FEniCSx that enables easy prototyping of...
Tokens:78,753
Snippets:693
Trust Score:9.7
Update:5 months ago
Show doc for...
Tokens:
Raw
### Install Nanobind and Scikit-build-core via Apt

Source: https://multiphenics.github.io/installing.html

Install build dependencies python3-nanobind and python3-scikit-build-core using apt. This is required if dolfinx was installed via apt.

```bash
apt install python3-nanobind python3-scikit-build-core
```

--------------------------------

### Install Multiphenicsx with Tutorial Dependencies

Source: https://multiphenics.github.io/installing.html

Install the multiphenicsx library from PyPI. Use the '[tutorials]' extra to include dependencies needed for running the tutorials.

```bash
python3 -m pip install --no-build-isolation 'multiphenicsx[tutorials]'
```

--------------------------------

### Install Nanobind and Scikit-build-core via Pip

Source: https://multiphenics.github.io/installing.html

Install build dependencies nanobind and scikit-build-core using pip. This is required if dolfinx was installed via Docker or built from source.

```bash
python3 -m pip install nanobind scikit-build-core[pyproject]
```

--------------------------------

### Initialize Functions for Domain Decomposition

Source: https://multiphenics.github.io/tutorials/05_singular_poisson/tutorial_create_nullspace.html

Initializes DOLFINx Functions for different subdomains in a domain decomposition setup. These functions will hold the solution components.

```python
(u1_dd, u2_dd, l_dd) = (dolfinx.fem.Function(V1), dolfinx.fem.Function(V2), dolfinx.fem.Function(M))
```

--------------------------------

### Clone Multiphenicsx Repository and Checkout Version

Source: https://multiphenics.github.io/installing.html

Clone the multiphenicsx repository from GitHub and checkout the tag corresponding to the installed version. This is a prerequisite for running the tutorials.

```bash
git clone https://github.com/multiphenics/multiphenicsx.git
cd multiphenicsx
MULTIPHENICSX_VERSION=$(python3 -c "import multiphenicsx; print(multiphenicsx.__version__)")
git checkout ${MULTIPHENICSX_VERSION}
```

--------------------------------

### Conditional Skip for Optimal Control Problem Setup

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_8_navier_stokes_neumann_control.html

This cell is marked to be skipped due to a previous test failure. It contains the setup for the main optimal control problem, including the nonlinear problem definition.

```python
"Skip cell due to a previously xfailed cell.

problem = multiphenicsx.fem.petsc.NonlinearProblem(
    F, (v, p, u, z, b), bcs=bc, J=dF,
    petsc_options_prefix=\"tutorial_8_navier_stokes_neumann_control_\", petsc_options=petsc_options,
    kind=\"mpi\", restriction=restriction
)
problem.solve()
del problem
"
```

--------------------------------

### Install Nanobind and Scikit-build-core via Conda

Source: https://multiphenics.github.io/installing.html

Install build dependencies nanobind and scikit-build-core using conda. This is required if dolfinx was installed via conda.

```bash
conda install -c conda-forge nanobind scikit-build-core
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_1_eigenvalue_laplacian.html

Imports all required libraries for the tutorial, including DOLFINx, Gmsh, MPI, NumPy, SciPy, SLEPc, UFL, and VisKex.

```python
import typing
```

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.la
import dolfinx.la.petsc
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import scipy.special
import slepc4py.SLEPc
import ufl
import viskex
```

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Define Optimality Conditions and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2a_advection_diffusion_reaction.html

Sets up the bilinear and linear forms for the state and adjoint equations, along with boundary conditions. This is crucial for defining the optimization problem.

```python
state_operator = (epsilon * ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx
                  + ufl.inner(ufl.dot(beta, ufl.grad(y)), q) * ufl. dx + sigma * ufl.inner(y, q) * ufl.dx)
adjoint_operator = (epsilon * ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx
                    - ufl.inner(ufl.dot(beta, ufl.grad(p)), z) * ufl.dx + sigma * ufl.inner(p, z) * ufl.dx)
a = [[ufl.inner(y, z) * ufl.dx, None, adjoint_operator],
     [None, alpha * ufl.inner(u, v) * ufl.dx, - ufl.inner(p, v) * ufl.dx],
     [state_operator, - ufl.inner(u, q) * ufl.dx, None]]
f = [ufl.inner(y_d, z) * ufl.dx,
     None,
     ufl.inner(ff, q) * ufl.dx]
a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * ufl.dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * ufl.dx
bdofs_Y_1234 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_1234)
bdofs_Q_1234 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_1234)
bc = [dolfinx.fem.dirichletbc(zero, bdofs_Y_1234, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_1234, Q)]
```

--------------------------------

### Assemble Discrete System

Source: https://multiphenics.github.io/tutorials/05_singular_poisson/tutorial_create_nullspace.html

Assemble the stiffness matrix 'A' and the load vector 'F' from the defined weak forms using PETSc backend.

```python
# Assemble the discrete system
A = dolfinx.fem.petsc.assemble_matrix(a_cpp)
A.assemble()
F = dolfinx.fem.petsc.assemble_vector(f_cpp)
F.ghostUpdate(
    addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)  # type: ignore[arg-type]
```

--------------------------------

### Define Problem Data

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_3b_advection_diffusion_reaction_neumann_control.html

Set up problem parameters, spatial coordinates, advection velocity field (beta), and constants for the PDE.

```python
alpha = 0.07
y_d = 2.5
epsilon = 1. / 12.
x = ufl.SpatialCoordinate(mesh)
beta = ufl.as_vector((x[1] * (1 - x[1]), 0))
zero =  petsc4py.PETSc.ScalarType(0)  # type: ignore[operator]
one =  petsc4py.PETSc.ScalarType(1)  # type: ignore[operator]
sigma = dolfinx.fem.Constant(mesh, zero)
ff = dolfinx.fem.Constant(mesh, zero)
```

--------------------------------

### Define and Solve Linear Problem with PETSc Options

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_3_lagrange_multipliers_interface.html

Configures and solves a linear problem using PETSc options for performance tuning. Ensure PETSc is installed and configured for use.

```python
petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "mat_mumps_icntl_7": 4,
    "ksp_error_if_not_converged": True,
}
problem = multiphenicsx.fem.petsc.LinearProblem(
    a, f, bcs=bcs,
    petsc_options_prefix="tutorial_1_lagrange_multipliers_interface_", petsc_options=petsc_options,
    kind="mpi", restriction=restriction
)
u1, u2, l = problem.solve()
del problem
```

--------------------------------

### Define problem forms and source terms

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_1_lagrange_multipliers_linear.html

Sets up the bilinear form 'a' and linear form 'f' for the weak formulation, including the boundary condition function 'g'.

```python
# Define problem block forms
g = dolfinx.fem.Function(V)
g.interpolate(lambda x: np.sin(3 * x[0] + 1) * np.sin(3 * x[1] + 1))
a = [[ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, ufl.inner(l, v) * ufl.ds],
     [ufl.inner(u, m) * ufl.ds, None]]
f = [ufl.inner(1, v) * ufl.dx, ufl.inner(g, m) * ufl.ds]
```

--------------------------------

### Define Mesh Tags on Specific Facets

Source: https://multiphenics.github.io/tutorials/02_understanding_restrictions/tutorial_understanding_restrictions.html

Creates a mesh tag for a subset of facets based on a geometric condition. This example tags facets where the y-coordinate plus the square root of 3 times the x-coordinate is close to zero.

```python
facet_entities_subset = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim - 1, lambda x: np.fabs(x[1] + np.sqrt(3) * x[0]) < 0.01)
facet_values_subset = np.full(facet_entities_subset.shape, 1, dtype=np.int32)
facet_restriction_subset = dolfinx.mesh.meshtags(
    mesh, mesh.topology.dim - 1, facet_entities_subset, facet_values_subset)
facet_restriction_subset.name = "facet_restriction_subset"

```

--------------------------------

### Verify DOF Count and Locations for CG_2 Restriction

Source: https://multiphenics.github.io/tutorials/02_understanding_restrictions/tutorial_understanding_restrictions.html

Asserts the total number of DOFs for the CG_2 restriction and verifies specific DOF locations using polar coordinates. This confirms the correct setup of a DofMapRestriction for the second CG space.

```python
CG_2 = CG[1]
dofmap_restriction_CG_2 = dofmap_restriction_CG[facet_restriction_subset][1]
assert count_dofs(dofmap_restriction_CG_2, mesh.comm) == 5
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_2, dofmap_restriction_CG_2)) == 1
assert len(locate_dofs_by_polar_coordinates(
    0.5, 2 * np.pi / 3, CG_2, dofmap_restriction_CG_2)) == 1
```

--------------------------------

### Define Problem Data and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_8_navier_stokes_neumann_control.html

Sets up the problem parameters, including viscosity, source terms, and desired velocity field. It also defines symbolic variables for defining the desired velocity and interpolates it into the velocity space.

```python
alpha = 1.e-5
x, y = sympy.symbols("x[0], x[1]")
psi_d = 10 * (1 - sympy.cos(0.8 * np.pi * x)) * (1 - sympy.cos(0.8 * np.pi * y)) * (1 - x)**2 * (1 - y)**2
v_d_x = sympy.lambdify([x, y], psi_d.diff(y, 1))
v_d_y = sympy.lambdify([x, y], - psi_d.diff(x, 1))
v_d = dolfinx.fem.Function(Y_velocity)
v_d.interpolate(lambda x: np.stack((v_d_x(x[0], x[1]), v_d_y(x[0], x[1])), axis=0))
nu = 0.1
zero_scalar =  petsc4py.PETSc.ScalarType(0)  # type: ignore[operator]
zero_vector = np.zeros((2, ), dtype=petsc4py.PETSc.ScalarType)
ff = dolfinx.fem.Constant(mesh, zero_vector)

```

--------------------------------

### Initialize Solution Functions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2a_advection_diffusion_reaction.html

Declares DOLFINx Functions for the state (y), control (u), and adjoint (p) variables. These will store the computed solutions.

```python
(y, u, p) = (dolfinx.fem.Function(Y), dolfinx.fem.Function(U), dolfinx.fem.Function(Q))
```

--------------------------------

### Verify DOF Count and Locations for CG_1 Restriction

Source: https://multiphenics.github.io/tutorials/02_understanding_restrictions/tutorial_understanding_restrictions.html

Asserts the total number of DOFs for the CG_1 restriction and verifies specific DOF locations using polar coordinates. This is used to confirm the correct setup of a DofMapRestriction for the first CG space.

```python
CG_1 = CG[0]
dofmap_restriction_CG_1 = dofmap_restriction_CG[facet_restriction_subset][0]
assert count_dofs(dofmap_restriction_CG_1, mesh.comm) == 3
assert len(locate_dofs_by_polar_coordinates(
    0, 0, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 2 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
assert len(locate_dofs_by_polar_coordinates(
    1, 5 * np.pi / 3, CG_1, dofmap_restriction_CG_1)) == 1
```

--------------------------------

### Initialize Solution Functions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_5_stokes.html

Initializes DOLFINx Functions for the state variables (velocity 'v' and pressure 'p') and the control variable 'u', as well as adjoint variables (velocity 'z' and pressure 'b').

```python
(v, p) = (dolfinx.fem.Function(Y_velocity), dolfinx.fem.Function(Y_pressure))
u = dolfinx.fem.Function(U)
(z, b) = (dolfinx.fem.Function(Q_velocity), dolfinx.fem.Function(Q_pressure))
```

--------------------------------

### Define Optimality Conditions and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_4b_poisson_neumann_control_boundary_observation.html

Sets up the bilinear and linear forms for the optimality conditions and defines Dirichlet boundary conditions for the state and adjoint variables. Ensure that `zero` is a `dolfinx.fem.Constant` initialized to 0.

```python
a = [[ufl.inner(y, z) * ds(2), None, ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx],
     [None, alpha * ufl.inner(u, v) * ds(2), - ufl.inner(p, v) * ds(2)],
     [ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx, - ufl.inner(u, q) * ds(2), None]]
f = [ufl.inner(y_d, z) * ds(2),
     None,
     ufl.inner(ff, q) * ufl.dx]
a[0][0] += dolfinx.fem.Constant(mesh, zero) * ufl.inner(y, z) * ufl.dx
a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * ufl.dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * ufl.dx
bdofs_Y_4 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_4)
bdofs_Q_4 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_4)
bc = [dolfinx.fem.dirichletbc(zero, bdofs_Y_4, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_4, Q)]
```

--------------------------------

### Solve Optimal Control Problem

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_3b_advection_diffusion_reaction_neumann_control.html

Set up and solve the full optimal control problem using multiphenicsx's LinearProblem, which handles the coupled state-adjoint-control system.

```python
problem = multiphenicsx.fem.petsc.LinearProblem(
    a, f, bcs=bc, u=(y, u, p),
    petsc_options_prefix="tutorial_3b_advection_diffusion_reaction_", petsc_options=petsc_options,
    kind="mpi", restriction=restriction
)
problem.solve()
del problem
```

--------------------------------

### Define Optimality Conditions and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_4a_poisson_dirichlet_control.html

Sets up the bilinear and linear forms for the optimality conditions and defines Dirichlet boundary conditions for the state and adjoint variables.

```python
a = [[ufl.inner(y, z) * ufl.dx, None, ufl.inner(l, z) * ds(2), ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx],
     [None, alpha * ufl.inner(u, v) * ds(2), - ufl.inner(l, v) * ds(2), None],
     [ufl.inner(y, m) * ds(2), - ufl.inner(u, m) * ds(2), None, None],
     [ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx, None, None, None]]
f = [ufl.inner(y_d, z) * ufl.dx,
     None,
     None,
     ufl.inner(ff, q) * ufl.dx]
a[3][3] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * ufl.dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * ufl.dx
f[2] = ufl.inner(dolfinx.fem.Constant(mesh, zero), m) * ufl.dx
bdofs_Y_4 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_4)
bdofs_Q_24 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_24)
bc = [dolfinx.fem.dirichletbc(zero, bdofs_Y_4, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_24, Q)]
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_3a_advection_diffusion_reaction_neumann_control.html

Imports essential libraries for FEM, PETSc, MPI, UFL, and visualization.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import ufl
import viskex
```

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/05_singular_poisson/tutorial_create_nullspace.html

Imports essential libraries for FEM, PETSc, Gmsh, MPI, NumPy, and visualization.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.la
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import ufl
import viskex
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_3_lagrange_multipliers_interface.html

Imports essential libraries for FEM, PETSc, Gmsh, MPI, NumPy, UFL, and viskex.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
```

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Define Solution Variables

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_4a_poisson_dirichlet_control.html

Initializes the state, control, Lagrange multiplier, and adjoint variables as DOLFINx Functions.

```python
(y, u, l, p) = (dolfinx.fem.Function(Y), dolfinx.fem.Function(U), dolfinx.fem.Function(L), dolfinx.fem.Function(Q))
```

--------------------------------

### Define problem parameters and desired state

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_6_navier_stokes.html

Sets up physical parameters such as viscosity (nu) and penalization (alpha), and defines the desired velocity field (v_d) using symbolic expressions and UFL. It also initializes forcing terms.

```python
alpha = 1.e-5
epsilon = 1.e-5
x, y = sympy.symbols("x[0], x[1]")
psi_d = 10 * (1 - sympy.cos(0.8 * np.pi * x)) * (1 - sympy.cos(0.8 * np.pi * y)) * (1 - x)**2 * (1 - y)**2
v_d_x = sympy.lambdify([x, y], psi_d.diff(y, 1))
v_d_y = sympy.lambdify([x, y], - psi_d.diff(x, 1))
v_d = dolfinx.fem.Function(Y_velocity)
v_d.interpolate(lambda x: np.stack((v_d_x(x[0], x[1]), v_d_y(x[0], x[1])), axis=0))
nu = 0.01
zero_scalar =  petsc4py.PETSc.ScalarType(0)  # type: ignore[operator]
zero_vector = np.zeros((2, ), dtype=petsc4py.PETSc.ScalarType)
ff = dolfinx.fem.Constant(mesh, zero_vector)
```

--------------------------------

### Assemble matrices A and B

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_3_infsup_lagrange_multipliers.html

Assembles the stiffness matrix 'A' and the mass matrix 'B' for the eigenvalue problem using the defined forms and restrictions.

```python
# Assemble lhs and rhs matrices
A = multiphenicsx.fem.petsc.assemble_matrix(
    dolfinx.fem.form(a), bcs=[], kind="mpi", restriction=(restriction, restriction))
A.assemble()
B = multiphenicsx.fem.petsc.assemble_matrix(
    dolfinx.fem.form(b), bcs=[], kind="mpi", restriction=(restriction, restriction))
B.assemble()
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_5_stokes.html

Imports essential libraries for FEM, PETSc, MPI, NumPy, SymPy, UFL, and VisKex.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import sympy
import ufl
import viskex
```

--------------------------------

### Solve Optimal Control Problem

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2b_advection_diffusion_reaction.html

Set up and solve the full optimal control problem, including state, adjoint, and control equations, using PETSc. All boundary conditions are applied.

```python
problem = dolfinx.fem.petsc.LinearProblem(
    a, f, bcs=[*bc_state, *bc_adjoint], u=(y, u, p),
    petsc_options_prefix="tutorial_2b_advection_diffusion_reaction_", petsc_options=petsc_options,
    kind="mpi"
)
problem.solve()
del problem
```

--------------------------------

### Configure MUMPS for Singular Systems

Source: https://multiphenics.github.io/tutorials/05_singular_poisson/tutorial_create_nullspace.html

Configure a direct solver using MUMPS to solve singular linear systems. Set options `icntl=24` to detect null pivots and `icntl=25` to avoid recomputing the null space.

```python
solution = dolfinx.fem.Function(V)
ksp = petsc4py.PETSc.KSP()
ksp.create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
ksp.getPC().setFactorSetUpSolverType()
ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)  # detect null pivots
ksp.getPC().getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)  # do not compute null space again
ksp.setFromOptions()
ksp.solve(F, solution.x.petsc_vec)
solution.x.petsc_vec.ghostUpdate(
    addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)
ksp.destroy()
```

--------------------------------

### Import DOLFINx and other libraries

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2a_advection_diffusion_reaction.html

Imports necessary modules from DOLFINx, MPI, NumPy, PETSc, UFL, and Viskex for finite element computations and visualization.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import ufl
import viskex
```

--------------------------------

### Import MultiphenicsX FEM Module

Source: https://multiphenics.github.io/tutorials/02_understanding_restrictions/tutorial_understanding_restrictions.html

Imports the finite element method module from the MultiphenicsX library.

```python
import multiphenicsx.fem
```

--------------------------------

### Solve the Optimal Control Problem

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_1a_poisson.html

Sets up and solves the full optimal control problem, including the state, control, and adjoint variables, using PETSc.

```python
problem = dolfinx.fem.petsc.LinearProblem(
    a, f, bcs=bc, u=(y, u, p),
    petsc_options_prefix="tutorial_1a_poisson_", petsc_options=petsc_options,
    kind="mpi"
)
problem.solve()
del problem
```

--------------------------------

### Set Problem Parameters

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2a_advection_diffusion_reaction.html

Initializes the parameters for the advection-diffusion-reaction problem, including penalization parameter (alpha), desired state (y_d), diffusion coefficient (epsilon), advection field (beta), reaction coefficient (sigma), forcing term (ff), and zero.

```python
alpha = 1.e-5
    y_d = 1.
    epsilon = 1.e-1
    beta = ufl.as_vector((-1., -2.))
    sigma = 1.
    ff = 1.
    zero =  petsc4py.PETSc.ScalarType(0)  # type: ignore[operator]
```

--------------------------------

### Run Monolithic Stokes Inf-Sup Calculation

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_2_infsup_stokes.html

Sets up mixed function spaces, variational forms, and solves the generalized eigenvalue problem for the Stokes flow inf-sup constant using a monolithic approach. Requires DOLFINx, PETSc, and SLEPc.

```python
def run_monolithic() -> tuple[  # type: ignore[valid-type]
    petsc4py.PETSc.RealType,
    dolfinx.fem.Function[petsc4py.PETSc.ScalarType],
    dolfinx.fem.Function[petsc4py.PETSc.ScalarType],
    dolfinx.fem.Function[petsc4py.PETSc.ScalarType]
]:
    """Run standard formulation using a mixed function space."""
    # Function spaces
    W_element = basix.ufl.mixed_element([V_element, Q_element])
    W = dolfinx.fem.functionspace(mesh, W_element)

    # Test and trial functions: monolithic
    vq = ufl.TestFunction(W)
    (v, q) = ufl.split(vq)
    up = ufl.TrialFunction(W)
    (u, p) = ufl.split(up)

    # Variational forms
    lhs = (ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - ufl.inner(p, ufl.div(v)) * ufl.dx
           - ufl.inner(ufl.div(u), q) * ufl.dx)
    rhs = - ufl.inner(p, q) * ufl.dx

    # Define restriction for DOFs associated to homogenous Dirichlet boundary conditions
    dofs_W = np.arange(0, W.dofmap.index_map.size_local + W.dofmap.index_map.num_ghosts)
    W0 = W.sub(0)
    V, _ = W0.collapse()
    bdofs_V = dolfinx.fem.locate_dofs_topological((W0, V), mesh.topology.dim - 1, boundary_facets)[0]
    restriction = multiphenicsx.fem.DofMapRestriction(W.dofmap, np.setdiff1d(dofs_W, bdofs_V))

    # Assemble lhs and rhs matrices
    A = multiphenicsx.fem.petsc.assemble_matrix(
        dolfinx.fem.form(lhs), restriction=(restriction, restriction))
    A.assemble()
    B = multiphenicsx.fem.petsc.assemble_matrix(
        dolfinx.fem.form(rhs), restriction=(restriction, restriction))
    B.assemble()

    # Solve
    eps = slepc4py.SLEPc.EPS().create(mesh.comm)
    eps.setOperators(A, B)
    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)
    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.TARGET_REAL)
    eps.setTarget(1.e-5)
    eps.getST().setType(slepc4py.SLEPc.ST.Type.SINVERT)
    eps.getST().getKSP().setType("preonly")
    eps.getST().getKSP().getPC().setType("lu")
    eps.getST().getKSP().getPC().setFactorSolverType("mumps")
    eps.solve()
    assert eps.getConverged() >= 1

    # Extract leading eigenvalue and eigenvector
    vr = multiphenicsx.fem.petsc.create_vector(W, restriction=restriction)
    vi = multiphenicsx.fem.petsc.create_vector(W, restriction=restriction)
    eigv = eps.getEigenpair(0, vr, vi)
    r, i = eigv.real, eigv.imag
    assert abs(i) < 1.e-10
    assert r > 0., "r = " + str(r) + " is not positive"
    print("Inf-sup constant (monolithic): ", np.sqrt(r))

    # Transform eigenvector into eigenfunction
    r_fun = dolfinx.fem.Function(W)
    vr.ghostUpdate(
        addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)  # type: ignore[arg-type]
    with r_fun.x.petsc_vec.localForm() as r_fun_local, \
            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, W.dofmap, restriction) as vr_wrapper:
        r_fun_local[:] = vr_wrapper
    u_fun = r_fun.sub(0).collapse()
    (u_fun_1, u_fun_2) = (u_fun.sub(0).collapse(), u_fun.sub(1).collapse())
    p_fun = r_fun.sub(1).collapse()
    normalize(u_fun_1, u_fun_2, p_fun)

    eps.destroy()
    return (r, u_fun_1, u_fun_2, p_fun)
```

--------------------------------

### Plot solution u (strong BCs)

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_1_lagrange_multipliers_linear.html

Visualizes the computed solution 'u_ex' with strong Dirichlet boundary conditions applied, using viskex.

```python
viskex.dolfinx.plot_scalar_field(u_ex, "u")
```

--------------------------------

### Initialize Gmsh model

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_3_infsup_lagrange_multipliers.html

Initializes the Gmsh meshing tool and creates a new model named 'mesh'.

```python
gmsh.initialize()
gmsh.model.add("mesh")
```

--------------------------------

### Define Trial and Test Functions for UFL

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_3_lagrange_multipliers_interface.html

Sets up trial and test functions using the UFL (Unified Form Language) for defining variational forms in finite element analysis.

```python
# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
```

--------------------------------

### Import multiphenicsx libraries

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_3b_advection_diffusion_reaction_neumann_control.html

Imports specific FEM and PETSc modules from the multiphenicsx library.

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Configure and solve the eigenvalue problem

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_3_infsup_lagrange_multipliers.html

Sets up and solves the generalized eigenvalue problem (GEP) using SLEPc. It configures the solver for target real eigenvalues, uses an inverse subspace iteration method with MUMPS direct solver.

```python
# Solve
eps = slepc4py.SLEPc.EPS().create(mesh.comm)
eps.setOperators(A, B)
eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GNHEP)
eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)
eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.TARGET_REAL)
eps.setTarget(1.e-5)
eps.getST().setType(slepc4py.SLEPc.ST.Type.SINVERT)
eps.getST().getKSP().setType("preonly")
eps.getST().getKSP().getPC().setType("lu")
eps.getST().getKSP().getPC().setFactorSolverType("mumps")
eps.solve()
assert eps.getConverged() >= 1
```

--------------------------------

### Define Variational Forms and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_1b_poisson.html

Sets up the bilinear and linear forms for the optimality conditions and the cost functional. Also defines Dirichlet boundary conditions for the state and adjoint variables.

```python
a = [[ufl.inner(y, z) * dx, None, ufl.inner(ufl.grad(p), ufl.grad(z)) * dx],
     [None, alpha * ufl.inner(u, v) * dx, - ufl.inner(p, v) * dx],
     [ufl.inner(ufl.grad(y), ufl.grad(q)) * dx, - ufl.inner(u, q) * dx, None]]
f = [ufl.inner(y_d_1, z) * dx(1) + ufl.inner(y_d_2, z) * dx(2),
     None,
     ufl.inner(ff, q) * dx]
a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * dx
bdofs_Y_1 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_1)
bdofs_Q_1 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_1)
bc = [dolfinx.fem.dirichletbc(one, bdofs_Y_1, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_1, Q)]
```

--------------------------------

### Solve Optimal Control Problem

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_3a_advection_diffusion_reaction_neumann_control.html

Sets up and solves the full optimal control problem using Multiphenicsx's LinearProblem. This involves the coupled state, control, and adjoint equations with specified PETSc options and MPI restrictions.

```python
problem = multiphenicsx.fem.petsc.LinearProblem(
    a, f, bcs=bc, u=(y, u, p),
    petsc_options_prefix="tutorial_3a_advection_diffusion_reaction_", petsc_options=petsc_options,
    kind="mpi", restriction=restriction
)
problem.solve()
del problem
```

--------------------------------

### Import Core Libraries

Source: https://multiphenics.github.io/tutorials/02_understanding_restrictions/tutorial_understanding_restrictions.html

Imports essential libraries for finite element analysis, meshing, plotting, and numerical operations.

```python
import basix
import basix.ufl
import dolfinx.fem
import dolfinx.io
import dolfinx.mesh
import dolfinx.typing
import gmsh
import matplotlib as mpl
import matplotlib.collections
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import mpi4py.MPI
import mpl_toolkits.axes_grid1
import numpy as np
import numpy.typing as npt
```

--------------------------------

### Define problem data

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_1a_poisson.html

Sets up the problem parameters including the penalization parameter alpha, spatial coordinates, the desired state y_d, and the forcing term ff. The desired state is defined as a quadratic function of spatial coordinates.

```python
alpha = 1.e-5
    x = ufl.SpatialCoordinate(mesh)
    y_d = 10 * x[0] * (1 - x[0]) * x[1] * (1 - x[1])
    zero =  petsc4py.PETSc.ScalarType(0)  # type: ignore[operator]
    ff = dolfinx.fem.Constant(mesh, zero)
```

--------------------------------

### Import necessary libraries for DOLFINx

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2b_advection_diffusion_reaction.html

Imports essential libraries for finite element methods, mesh handling, I/O, linear algebra, and visualization within the DOLFINx framework.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import gmsh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import ufl
import viskex

```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_1_lagrange_multipliers_linear.html

Imports essential libraries for FEM, meshing, and parallel computing.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import ufl
import viskex
```

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Extract State Equations and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_4a_poisson_dirichlet_control.html

Extracts the state equation (a_state, f_state) and corresponding boundary conditions (bc_state) from the optimality conditions for solving the uncontrolled state.

```python
# Extract state forms from the optimality conditions
a_state = ufl.replace(a[3][0], {q: z})
f_state = ufl.replace(f[3], {q: z})
bdofs_Y_24 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_24)
bc_state = [dolfinx.fem.dirichletbc(zero, bdofs_Y_24, Y)]
```

--------------------------------

### Solve the Optimal Control Problem

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_2a_advection_diffusion_reaction.html

Sets up and solves the full optimal control problem, including state, control, and adjoint equations, using a PETSc LinearProblem.

```python
problem = dolfinx.fem.petsc.LinearProblem(
    a, f, bcs=bc, u=(y, u, p),
    petsc_options_prefix="tutorial_2a_advection_diffusion_reaction_", petsc_options=petsc_options,
    kind="mpi"
)
problem.solve()
del problem
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_3_infsup_lagrange_multipliers.html

Imports essential libraries for finite element analysis, meshing, parallel computing, and eigenvalue solvers.

```python
import dolfinx.fem
import dolfinx.io
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import slepc4py.SLEPc
import ufl
import viskex
```

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Assemble Block Matrix and Vector

Source: https://multiphenics.github.io/tutorials/05_singular_poisson/tutorial_create_nullspace.html

Assembles the block linear system matrices and vectors using DOLFINx forms and MultiPhenicsX PETSc wrappers. Ensure correct restrictions are applied for distributed computations.

```python
zero = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(0))  # type: ignore[operator]
a_dd = [[ufl.inner(ufl.grad(u1), ufl.grad(v1)) * dx(1), None, ufl.inner(l("-"), v1("-")) * dS],
        [None, ufl.inner(ufl.grad(u2), ufl.grad(v2)) * dx(2), - ufl.inner(l("+"), v2("+")) * dS],
        [ufl.inner(u1("-"), m("-")) * dS, - ufl.inner(u2("+"), m("+")) * dS, None]]
f_dd = [ufl.inner(f1, v1) * dx(1) + ufl.inner(g, v1) * ds(1),
        ufl.inner(f2, v2) * dx(2) + ufl.inner(g, v2) * ds(2),
        ufl.inner(zero, m("-")) * dS]
a_dd_cpp = dolfinx.fem.form(a_dd)
f_dd_cpp = dolfinx.fem.form(f_dd)
```

```python
# Assemble the block linear system
A_dd = multiphenicsx.fem.petsc.assemble_matrix(a_dd_cpp, kind="mpi", restriction=(restriction, restriction))
A_dd.assemble()
F_dd = multiphenicsx.fem.petsc.assemble_vector(f_dd_cpp, kind="mpi", restriction=restriction)
F_dd.ghostUpdate(
    addv=petsc4py.PETSc.InsertMode.ADD, mode=petsc4py.PETSc.ScatterMode.REVERSE)  # type: ignore[arg-type]
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_2_infsup_stokes.html

Imports essential libraries for numerical computations, FEM, and parallel processing.

```python
import typing
```

```python
import basix.ufl
import dolfinx.fem
import dolfinx.mesh
import mpi4py.MPI
import numpy as np
import numpy.typing as npt
import petsc4py.PETSc
import slepc4py.SLEPc
import ufl
import viskex
```

```python
import multiphenicsx.fem
import multiphenicsx.fem.petsc
```

--------------------------------

### Solve Optimal Control Problem

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_1b_poisson.html

Sets up and solves the full optimal control problem, which involves the state, control, and adjoint equations simultaneously. Uses PETSc with MUMPS for direct solving.

```python
problem = dolfinx.fem.petsc.LinearProblem(
    a, f, bcs=bc, u=(y, u, p),
    petsc_options_prefix="tutorial_1b_poisson_", petsc_options=petsc_options,
    kind="mpi"
)
problem.solve()
del problem
```

--------------------------------

### Define Variational Forms and Boundary Conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_3a_advection_diffusion_reaction_neumann_control.html

Defines the state operator, adjoint operator, bilinear form 'a', linear form 'f', and boundary conditions for the optimal control problem. This includes setting up Dirichlet boundary conditions for specific boundary facets.

```python
state_operator = (epsilon * ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx
                  + ufl.inner(ufl.dot(beta, ufl.grad(y)), q) * ufl.dx + sigma * ufl.inner(y, q) * ufl.dx)
adjoint_operator = (epsilon * ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx
                    - ufl.inner(ufl.dot(beta, ufl.grad(p)), z) * ufl.dx + sigma * ufl.inner(p, z) * ufl.dx)
a = [[ufl.inner(y, z) * ufl.dx, None, adjoint_operator],
     [None, alpha * ufl.inner(u, v) * ds(2), - ufl.inner(p, v) * ds(2)],
     [state_operator, - ufl.inner(u, q) * ds(2), None]]
f = [ufl.inner(y_d, z) * ufl.dx,
     None,
     ufl.inner(ff, q) * ufl.dx]
a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * ufl.dx
f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * ufl.dx
bdofs_Y_4 = dolfinx.fem.locate_dofs_topological(Y, mesh.topology.dim - 1, boundaries_4)
bdofs_Q_4 = dolfinx.fem.locate_dofs_topological(Q, mesh.topology.dim - 1, boundaries_4)
bc = [dolfinx.fem.dirichletbc(zero, bdofs_Y_4, Y),
      dolfinx.fem.dirichletbc(zero, bdofs_Q_4, Q)]
```

--------------------------------

### Define trial and test functions

Source: https://multiphenics.github.io/tutorials/04_eigenvalue_problems/tutorial_3_infsup_lagrange_multipliers.html

Defines the trial and test functions (u, l) and (v, m) for the weak form of the problem.

```python
# Define trial and test functions
(u, l) = (ufl.TrialFunction(V), ufl.TrialFunction(V))
(v, m) = (ufl.TestFunction(V), ufl.TestFunction(V))
```

--------------------------------

### Import necessary libraries

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_1b_poisson.html

Imports essential libraries for finite element analysis, meshing, and parallel computing.

```python
import dolfinx.fem
import dolfinx.fem.petsc
import dolfinx.io
import dolfinx.mesh
import gmsh
import mpi4py.MPI
import numpy as np
import petsc4py.PETSc
import ufl
import viskex
```

--------------------------------

### Define Function Spaces for Domain Decomposition

Source: https://multiphenics.github.io/tutorials/01_lagrange_multipliers/tutorial_3_lagrange_multipliers_interface.html

Sets up function spaces for different subdomains (V1, V2) and for the Lagrange multiplier (M) on the interface. This is the initial step for problems involving domain decomposition.

```python
# Define function spaces
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
V1 = V.clone()
V2 = V.clone()
M = V.clone()

```

--------------------------------

### Initialize Solution Functions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_8_navier_stokes_neumann_control.html

Initializes the solution functions for velocity and pressure, as well as the control variable. These functions will store the computed solution during the iterative process.

```python
(v, p) = (dolfinx.fem.Function(Y_velocity), dolfinx.fem.Function(Y_pressure))
u = dolfinx.fem.Function(U)
(z, b) = (dolfinx.fem.Function(Q_velocity), dolfinx.fem.Function(Q_pressure))

```

--------------------------------

### Formulate optimality conditions

Source: https://multiphenics.github.io/tutorials/06_optimal_control/tutorial_1a_poisson.html

Constructs the bilinear form 'a' and the linear form 'f' representing the optimality conditions for the optimal control problem. This includes the state equation, adjoint equation, and control update.

```python
a = [[ufl.inner(y, z) * ufl.dx, None, ufl.inner(ufl.grad(p), ufl.grad(z)) * ufl.dx],
         [None, alpha * ufl.inner(u, v) * ufl.dx, - ufl.inner(p, v) * ufl.dx],
         [ufl.inner(ufl.grad(y), ufl.grad(q)) * ufl.dx, - ufl.inner(u, q) * ufl.dx, None]]
    f = [ufl.inner(y_d, z) * ufl.dx,
         None,
         ufl.inner(ff, q) * ufl.dx]
    a[2][2] = dolfinx.fem.Constant(mesh, zero) * ufl.inner(p, q) * ufl.dx
    f[1] = ufl.inner(dolfinx.fem.Constant(mesh, zero), v) * ufl.dx
```
Tokens:
Raw