|
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