# The Helmholtz equation on the Torus

Helmholtz equation is given as

$$
\nabla^2 u + \alpha u = f.
$$

In this notebook we will solve this equation (without boundary conditions) on the surface of a torus, using curvilinear coordinates. The surface of the torus is parametrized by    

$$
\begin{align}
x(\theta, \phi) &= (R + r \cos \theta) \cos \phi \\
y(\theta, \phi) &= (R + r \cos \theta) \sin \phi \\
z(\theta, \phi) &= r \sin \theta
\end{align}
$$

where $\theta, \phi$ are angles which make a full circle, so that their values start and end at the same point, $R$ is the distance from the center of the tube to the center of the torus,
$r$ is the radius of the tube. Note that $\theta$ is the angle in the small circle (around its center), whereas $\phi$ is the angle of the large circle, around origo.

We start the implementation by importing necessary functionality from shenfun and sympy and then defining the coordinates of the surface of the torus. Note that `rv` represents the position vector spanning the surface.


```python
from shenfun import *
from shenfun.la import SolverGeneric1ND
import sympy as sp
from IPython.display import Math

N = 96
R = 3
r = 1
theta, phi = psi = sp.symbols('x,y', real=True, positive=True)
rv = ((R + r*sp.cos(theta))*sp.cos(phi), (R + r*sp.cos(theta))*sp.sin(phi), r*sp.sin(theta)) 
```

Now create necessary bases and function spaces. Due to the geometry of the problem, the solution will be periodic in both $\theta$ and $\phi$ directions, and we choose Fourier basis functions. The basis for the $\phi$-direction can be either real to complex or complex to complex, depending on the type of the solution. Here we assume a real solution


```python
B1 = FunctionSpace(N, 'F', dtype='D', domain=(0, 2*np.pi))
B2 = FunctionSpace(N, 'F', dtype='d', domain=(0, 2*np.pi))
T = TensorProductSpace(comm, (B1, B2), coordinates=(psi, rv, sp.Q.positive(r*sp.cos(theta)+R)))
V = VectorSpace(T)
u = TrialFunction(T)
v = TestFunction(T)
T.coors.sg
```

Note the assumption `sp.Q.positive(r*sp.cos(theta)+R))`, which is there to help sympy when computing basis vectors and scaling factors. It is not completely necessary, but try to omit it and look at what happens to the expanded Helmholtz equation below.


```python
alpha = 0
du = div(grad(u))+alpha*u
g = sp.Symbol('g', real=True, positive=True) # The Jacobian**2 (T.coors.sg**2)
#replace = [(3*sp.cos(theta)+5, g/2), (sp.cos(theta)**2+6*sp.cos(theta)+10, sp.cos(theta)**2+g), (3*sp.cos(theta)+10, g/2+5)] # to simplify the look
replace = [(sp.cos(theta)+3, g)]
Math(du.tolatex(symbol_names={r: 'r', theta: '\\theta', phi: '\\phi'}, replace=replace))
```

We now create a manufactured solution that satisfies periodicity and compute the right hand side $f$.


```python
ue = sp.sin(sp.cos(theta*4))*sp.cos(4*phi)
f = (div(grad(u))+alpha*u).tosympy(basis=ue, psi=psi)
fj = Array(T, buffer=f*T.sg)
f_hat = Function(T)
f_hat = inner(v, fj, output_array=f_hat)
```

Assemble coefficient matrix and solve problem. Note that the tensorproduct matrices along axis 0 can be non-diagonal due to the measure $\cos \theta + 3$. The matrices along the second axis will all be diagonal, so we can choose to use `SolverGeneric1ND`


```python
mats = inner(v*T.sg, (div(grad(u))+alpha*u))
#mats = inner(grad(v*T.sg), -grad(u)) # + inner(v, alpha*u)
u_hat = Function(T)
sol = SolverGeneric1ND(mats)
u_hat = sol(f_hat, u_hat, constraints=((0, 0),))
uj = u_hat.backward()
uq = Array(T, buffer=ue)
print('Error =', np.sqrt(inner(1, (uj-uq)**2)))
```


```python
mats = inner(v*T.sg, u)
```

Finally, just plot the solution using plotly.


```python
surf3D(u_hat, wrapaxes=[0, 1])
```

Inspect sparsity pattern at a given Fourier wavenumber 


```python
import matplotlib.pyplot as plt
%matplotlib notebook
plt.figure(figsize=(6,4))
plt.spy(sol.MM[6].mat, markersize=0.1)
```

Inspect the vector Laplacian


```python
uv = TrialFunction(V)
vv = TestFunction(V)
du = div(grad(uv))
Math(du.tolatex(symbol_names={r: 'r', theta: '\\theta', phi: '\\phi'}, replace=replace))
```


```python
Math(T.coors.latex_basis_vectors(symbol_names={r: 'r', theta: '\\theta', phi: '\\phi'}, replace=replace))
```


```python
%matplotlib notebook
M = inner(vv*T.sg, du)
B = BlockMatrix(M)
plt.figure(figsize=(6,4))
plt.spy(B.diags(it=(4,)), markersize=0.2)
```

8e20cb5198e31ec401a97d2ba1c3d873fa7d9bc9
ipynb
Jupyter Notebook
binder/Torus.ipynb
spectralDNS/shenfun
956633aa0f1638db5ebdc497ff68a438aa22b932
BSD-2-Clause
2017-06-17T13:30:27.000Z
2022-03-20T02:33:47.000Z
binder/Torus.ipynb
spectralDNS/shenfun
956633aa0f1638db5ebdc497ff68a438aa22b932
BSD-2-Clause
2017-05-16T06:53:04.000Z
2022-02-04T10:40:44.000Z
binder/Torus.ipynb
spectralDNS/shenfun
956633aa0f1638db5ebdc497ff68a438aa22b932
BSD-2-Clause
2018-01-31T14:37:01.000Z
2022-03-31T15:07:27.000Z
Qwen/Qwen-72B
1. YES
2. YES