|
6 | 6 | # extension: .py
|
7 | 7 | # format_name: light
|
8 | 8 | # format_version: '1.5'
|
9 |
| -# jupytext_version: 1.14.7 |
| 9 | +# jupytext_version: 1.15.2 |
10 | 10 | # kernelspec:
|
11 | 11 | # display_name: Python 3 (DOLFINx complex)
|
12 | 12 | # language: python
|
|
19 | 19 | #
|
20 | 20 | # Many PDEs, such as the [Helmholtz equation](https://docs.fenicsproject.org/dolfinx/v0.4.1/python/demos/demo_helmholtz.html) require complex-valued fields.
|
21 | 21 | #
|
22 |
| -# For simplicity, let us consider a Poisson equation of the form: |
| 22 | +# For simplicity, let us consider a Poisson equation of the form: |
23 | 23 | #
|
24 | 24 | # $$-\Delta u = f \text{ in } \Omega,$$
|
25 |
| -# $$ f = -1 - 2j \text{ in } \Omega,$$ |
| 25 | +# $$ f = -1 - 2j \text{ in } \Omega,$$ |
26 | 26 | # $$ u = u_{exact} \text{ on } \partial\Omega,$$
|
27 | 27 | # $$u_{exact}(x, y) = \frac{1}{2}x^2 + 1j\cdot y^2,$$
|
28 | 28 | #
|
29 | 29 | # As in [Solving the Poisson equation](./fundamentals) we want to express our partial differential equation as a weak formulation.
|
30 | 30 | #
|
31 |
| -# We start by defining our discrete function space $V_h$, such that $u_h\in V_h$ and $u_h = \sum_{i=1}^N c_i \phi_i(x, y)$ where $\phi_i$ are **real valued** global basis functions of our space $V_h$, $c_i \in \mathcal{C}$ are the **complex valued** degrees of freedom. |
| 31 | +# We start by defining our discrete function space $V_h$, such that $u_h\in V_h$ and $u_h = \sum_{i=1}^N c_i \phi_i(x, y)$ where $\phi_i$ are **real valued** global basis functions of our space $V_h$, and $c_i \in \mathcal{C}$ are the **complex valued** degrees of freedom. |
32 | 32 | #
|
33 |
| -# Next, we choose a test function $v\in \hat V_h$ where $\hat V_h\subset V_h$ such that $v\vert_{\partial\Omega}=0$, as done in the first tutorial. |
34 |
| -# We now need to define our inner product space. We choose the $L^2$ inner product spaces, which is a *[sesquilinear](https://en.wikipedia.org/wiki/Sesquilinear_form) 2-form*, Meaning that $\langle u, v\rangle$ is a map from $V_h\times V_h\mapsto K$, and $\langle u, v \rangle = \int_\Omega u \cdot \bar v ~\mathrm{d} x$. As it is sesquilinear, we have the following properties: |
| 33 | +# Next, we choose a test function $v\in \hat V_h$ where $\hat V_h\subset V_h$ such that $v\vert_{\partial\Omega}=0$, as done in the first tutorial. |
| 34 | +# We now need to define our inner product space. We choose the $L^2$ inner product spaces, which is a _[sesquilinear](https://en.wikipedia.org/wiki/Sesquilinear_form) 2-form_, meaning that $\langle u, v\rangle$ is a map from $V_h\times V_h\mapsto K$, and $\langle u, v \rangle = \int_\Omega u \cdot \bar v ~\mathrm{d} x$. As it is sesquilinear, we have the following properties: |
35 | 35 | #
|
36 | 36 | # $$\langle u , v \rangle = \overline{\langle v, u \rangle},$$
|
37 | 37 | # $$\langle u , u \rangle \geq 0.$$
|
|
41 | 41 | # $$\int_\Omega \nabla u_h \cdot \nabla \overline{v}~ \mathrm{dx} = \int_{\Omega} f \cdot \overline{v} ~\mathrm{d} s \qquad \forall v \in \hat{V}_h.$$
|
42 | 42 | #
|
43 | 43 | # ## Installation of FEniCSx with complex number support
|
44 |
| -# FEniCSx supports both real and complex numbers, meaning that we can create a function spaces with real valued or complex valued coefficients. |
| 44 | +# |
| 45 | +# FEniCSx supports both real and complex numbers, so we can create a function space with real valued or complex valued coefficients. |
| 46 | +# |
45 | 47 |
|
46 |
| -import dolfinx |
47 | 48 | from mpi4py import MPI
|
| 49 | +import dolfinx |
48 | 50 | import numpy as np
|
49 | 51 | mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
|
50 | 52 | V = dolfinx.fem.FunctionSpace(mesh, ("Lagrange", 1))
|
|
55 | 57 | print(u_r.x.array.dtype)
|
56 | 58 | print(u_c.x.array.dtype)
|
57 | 59 |
|
58 |
| -# However, as we would like to solve linear algebra problems on the form $Ax=b$, we need to be able to use matrices and vectors that support real and complex numbers. As [PETSc](https://petsc.org/release/) is one of the most popular interfaces to linear algebra packages, we need to be able to work with their matrix and vector structures. |
| 60 | +# However, as we would like to solve linear algebra problems of the form $Ax=b$, we need to be able to use matrices and vectors that support real and complex numbers. As [PETSc](https://petsc.org/release/) is one of the most popular interfaces to linear algebra packages, we need to be able to work with their matrix and vector structures. |
59 | 61 | #
|
60 | 62 | # Unfortunately, PETSc only supports one floating type in their matrices, thus we need to install two versions of PETSc, one that supports `float64` and one that supports `complex128`. In the [docker images](https://hub.docker.com/r/dolfinx/dolfinx) for DOLFINx, both versions are installed, and one can switch between them by calling `source dolfinx-real-mode` or `source dolfinx-complex-mode`. For the `dolfinx/lab` images, one can change the Python kernel to be either the real or complex mode, by going to `Kernel->Change Kernel...` and choose `Python3 (ipykernel)` (for real mode) or `Python3 (DOLFINx complex)` (for complex mode).
|
61 | 63 | #
|
62 | 64 | # We check that we are using the correct installation of PETSc by inspecting the scalar type.
|
| 65 | +# |
63 | 66 |
|
64 | 67 | from petsc4py import PETSc
|
65 | 68 | from dolfinx.fem.petsc import assemble_vector
|
66 | 69 | print(PETSc.ScalarType)
|
67 | 70 | assert np.dtype(PETSc.ScalarType).kind == 'c'
|
68 | 71 |
|
69 | 72 | # ## Variational problem
|
| 73 | +# |
70 | 74 | # We are now ready to define our variational problem
|
| 75 | +# |
71 | 76 |
|
72 | 77 | import ufl
|
73 | 78 | u = ufl.TrialFunction(V)
|
|
78 | 83 |
|
79 | 84 | # Note that we have used the `PETSc.ScalarType` to wrap the constant source on the right hand side. This is because we want the integration kernels to assemble into the correct floating type.
|
80 | 85 | #
|
81 |
| -# Secondly, note that we are using `ufl.inner` to describe multiplication of $f$ and $v$, even if they are scalar values. This is because `ufl.inner` takes the conjugate of the second argument, as decribed by the $L^2$ inner product. One could alternatively write this out manually |
| 86 | +# Secondly, note that we are using `ufl.inner` to describe multiplication of $f$ and $v$, even if they are scalar values. This is because `ufl.inner` takes the conjugate of the second argument, as decribed by the $L^2$ inner product. One could alternatively write this out explicitly |
82 | 87 | #
|
83 | 88 | # ### Inner-products and derivatives
|
| 89 | +# |
84 | 90 |
|
85 | 91 | L2 = f * ufl.conj(v) * ufl.dx
|
86 | 92 | print(L)
|
87 | 93 | print(L2)
|
88 | 94 |
|
89 |
| -# Similarly, if we want to use the function $ufl.derivative$ to take derivatives of functionals, we need to take some special care. As `derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to be able to use it to assemble vectors. |
| 95 | +# Similarly, if we want to use the function `ufl.derivative` to take derivatives of functionals, we need to take some special care. As `derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to in order to assemble vectors. |
| 96 | +# |
90 | 97 |
|
91 | 98 | J = u_c**2 * ufl.dx
|
92 | 99 | F = ufl.derivative(J, u_c, ufl.conj(v))
|
93 | 100 | residual = assemble_vector(dolfinx.fem.form(F))
|
94 | 101 | print(residual.array)
|
95 | 102 |
|
96 | 103 | # We define our Dirichlet condition and setup and solve the variational problem.
|
| 104 | +# |
97 | 105 | # ## Solve variational problem
|
| 106 | +# |
98 | 107 |
|
99 | 108 | mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
|
100 | 109 | boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
|
|
103 | 112 | problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc])
|
104 | 113 | uh = problem.solve()
|
105 | 114 |
|
106 |
| -# We compute the $L^2$ error and the max error |
| 115 | +# We compute the $L^2$ error and the max error. |
| 116 | +# |
107 | 117 | # ## Error computation
|
| 118 | +# |
108 | 119 |
|
109 | 120 | x = ufl.SpatialCoordinate(mesh)
|
110 | 121 | u_ex = 0.5 * x[0]**2 + 1j*x[1]**2
|
|
115 | 126 | print(global_error, max_error)
|
116 | 127 |
|
117 | 128 | # ## Plotting
|
118 |
| -# Finally, we plot the real and imaginary solution |
| 129 | +# |
| 130 | +# Finally, we plot the real and imaginary solutions. |
| 131 | +# |
119 | 132 |
|
120 | 133 | import pyvista
|
121 | 134 | pyvista.start_xvfb()
|
|
0 commit comments