<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Spherical coordinates in shenfun
The Helmholtz equation is given as
$$
-\nabla^2 u + \alpha u = f.
$$
In this notebook we will solve this equation on a unitsphere, using spherical coordinates. To verify the implementation we use a spherical harmonics function as manufactured solution. 
We start the implementation by importing necessary functionality from shenfun and sympy
Step1:  Define spherical coordinates $(r, \theta, \phi)$
$$
\begin{align}
x &= r \sin \theta \cos \phi \
y &= r \sin \theta \sin \phi \
z &= r \cos \theta
\end{align}
$$
using sympy. The radius r will be constant r=1. We create the three-dimensional position vector rv as a function of the two new coordinates $(\theta, \phi)$.
Step2:  We define bases with the domains $\theta \in [0, \pi]$ and $\phi \in [0, 2\pi]$. Also define a tensorproductspace, test- and trialfunction. Note that the new coordinates and the position vector are fed to the TensorProductSpace and not the individual spaces
Step3:  Use one spherical harmonic function as a manufactured solution
Step4:  Compute the right hand side on the quadrature mesh and take the scalar product
Step5:  Note that we can use the shenfun operators div and grad on a trialfunction u, and then switch the trialfunction for a sympy function ue. The operators will then make use of sympy's derivative method on the function ue. Here (-div(grad(u))+alpha*u) corresponds to the equation we are trying to solve
Step6:  Evaluated with u=ue and you get the exact right hand side f.
Tensor product matrices that make up the Helmholtz equation are then assembled as
Step7:  And the linear system of equations can be solved using the generic SolverGeneric1ND, that can be used for any problem that only has non-periodic boundary conditions in one dimension.
Step8:  Transform back to real space and compute the error.
Step9:  Postprocessing
Since we used quite few quadrature points in solving this problem, we refine the solution for a nicer plot. Note that refine simply pads Functions with zeros, which gives exactly the same accuracy, but more quadrature points in real space. u_hat has NxM quadrature points, here we refine using 3 times as many points along both dimensions
Step10:  The periodic solution does not contain the periodic points twice, i.e., the computational mesh contains $0$, but not $2\pi$. It looks better if we wrap the periodic dimension all around to $2\pi$, and this is achieved with
Step11:  Biharmonic equation
A biharmonic equation is given as
$$
\nabla^4 u + \alpha u = f.
$$
This equation is extremely messy in spherical coordinates. I cannot even find it posted anywhere. Nevertheless, we can solve it trivially with shenfun, and we can also see what it looks like
Step12:  Remember that this equation uses constant radius r=1. We now solve the equation using the same manufactured solution as for the Helmholtz equation.
Step13:  Want to see what the regular 3-dimensional biharmonic equation looks like in spherical coordinates? This is extremely tedious to derive by hand, but in shenfun you can get there with the following few lines of code
Step14:  I don't know if this is actually correct, because I haven't derived it by hand and I haven't seen it printed anywhere, but at least I know the Cartesian equation is correct
<ASSISTANT_TASK:>
Python Code:
from shenfun import *
from shenfun.la import SolverGeneric1ND
import sympy as sp
Explanation: Spherical coordinates in shenfun
The Helmholtz equation is given as
$$
-\nabla^2 u + \alpha u = f.
$$
In this notebook we will solve this equation on a unitsphere, using spherical coordinates. To verify the implementation we use a spherical harmonics function as manufactured solution. 
We start the implementation by importing necessary functionality from shenfun and sympy:
End of explanation
r = 1
theta, phi = psi = sp.symbols('x,y', real=True, positive=True)
rv = (r*sp.sin(theta)*sp.cos(phi), r*sp.sin(theta)*sp.sin(phi), r*sp.cos(theta))
Explanation: Define spherical coordinates $(r, \theta, \phi)$
$$
\begin{align}
x &= r \sin \theta \cos \phi \
y &= r \sin \theta \sin \phi \
z &= r \cos \theta
\end{align}
$$
using sympy. The radius r will be constant r=1. We create the three-dimensional position vector rv as a function of the two new coordinates $(\theta, \phi)$.
End of explanation
N, M = 64, 64
L0 = FunctionSpace(N, 'C', domain=(0, np.pi))
F1 = FunctionSpace(M, 'F', dtype='d')
T = TensorProductSpace(comm, (L0, F1), coordinates=(psi, rv, sp.Q.positive(sp.sin(theta))))
v = TestFunction(T)
u = TrialFunction(T)
Explanation: We define bases with the domains $\theta \in [0, \pi]$ and $\phi \in [0, 2\pi]$. Also define a tensorproductspace, test- and trialfunction. Note that the new coordinates and the position vector are fed to the TensorProductSpace and not the individual spaces:
End of explanation
#sph = sp.functions.special.spherical_harmonics.Ynm
#ue = sph(6, 3, theta, phi)
ue = sp.cos(8*(sp.sin(theta)*sp.cos(phi) + sp.sin(theta)*sp.sin(phi) + sp.cos(theta)))
Explanation: Use one spherical harmonic function as a manufactured solution
End of explanation
alpha = 1000
g = (-div(grad(u))+alpha*u).tosympy(basis=ue, psi=psi)
gj = Array(T, buffer=g*T.coors.sg)
g_hat = Function(T)
g_hat = inner(v, gj, output_array=g_hat)
Explanation: Compute the right hand side on the quadrature mesh and take the scalar product
End of explanation
from IPython.display import Math
Math((-div(grad(u))+alpha*u).tolatex(funcname='u', symbol_names={theta: '\\theta', phi: '\\phi'}))
#Math((grad(u)).tolatex(funcname='u', symbol_names={theta: '\\theta', phi: '\\phi'}))
Explanation: Note that we can use the shenfun operators div and grad on a trialfunction u, and then switch the trialfunction for a sympy function ue. The operators will then make use of sympy's derivative method on the function ue. Here (-div(grad(u))+alpha*u) corresponds to the equation we are trying to solve:
End of explanation
mats = inner(v, (-div(grad(u))+alpha*u)*T.coors.sg)
mats[3].mats
Explanation: Evaluated with u=ue and you get the exact right hand side f.
Tensor product matrices that make up the Helmholtz equation are then assembled as
End of explanation
u_hat = Function(T)
Sol1 = SolverGeneric1ND(mats)
u_hat = Sol1(g_hat, u_hat)
Explanation: And the linear system of equations can be solved using the generic SolverGeneric1ND, that can be used for any problem that only has non-periodic boundary conditions in one dimension.
End of explanation
uj = u_hat.backward()
uq = Array(T, buffer=ue)
print('Error =', np.sqrt(inner(1, (uj-uq)**2)), np.linalg.norm(uj-uq))
np.linalg.norm(u_hat - u_hat.backward().forward())
import matplotlib.pyplot as plt
%matplotlib inline
plt.spy(Sol1.solvers1D[1].mat, markersize=0.2)
#raise RuntimeError
Explanation: Transform back to real space and compute the error.
End of explanation
u_hat2 = u_hat.refine([N*3, M*3])
Explanation: Postprocessing
Since we used quite few quadrature points in solving this problem, we refine the solution for a nicer plot. Note that refine simply pads Functions with zeros, which gives exactly the same accuracy, but more quadrature points in real space. u_hat has NxM quadrature points, here we refine using 3 times as many points along both dimensions
End of explanation
surf3D(u_hat2, wrapaxes=[1])
Explanation: The periodic solution does not contain the periodic points twice, i.e., the computational mesh contains $0$, but not $2\pi$. It looks better if we wrap the periodic dimension all around to $2\pi$, and this is achieved with
End of explanation
Math((div(grad(div(grad(u))))+alpha*u).tolatex(funcname='u', symbol_names={theta: '\\theta', phi: '\\phi'}))
Explanation: Biharmonic equation
A biharmonic equation is given as
$$
\nabla^4 u + \alpha u = f.
$$
This equation is extremely messy in spherical coordinates. I cannot even find it posted anywhere. Nevertheless, we can solve it trivially with shenfun, and we can also see what it looks like
End of explanation
g = (div(grad(div(grad(u))))+alpha*u).tosympy(basis=ue, psi=psi)
gj = Array(T, buffer=g)
# Take scalar product
g_hat = Function(T)
g_hat = inner(v, gj, output_array=g_hat)
mats = inner(v, div(grad(div(grad(u)))) + alpha*u)
# Solve
u_hat = Function(T)
Sol1 = SolverGeneric1ND(mats)
u_hat = Sol1(g_hat, u_hat)
# Transform back to real space.
uj = u_hat.backward()
uq = Array(T, buffer=ue)
print('Error =', np.sqrt(dx((uj-uq)**2)))
Explanation: Remember that this equation uses constant radius r=1. We now solve the equation using the same manufactured solution as for the Helmholtz equation.
End of explanation
r, theta, phi = psi = sp.symbols('x,y,z', real=True, positive=True)
rv = (r*sp.sin(theta)*sp.cos(phi), r*sp.sin(theta)*sp.sin(phi), r*sp.cos(theta))
L0 = FunctionSpace(20, 'L', domain=(0, 1))
F1 = FunctionSpace(20, 'L', domain=(0, np.pi))
F2 = FunctionSpace(20, 'F', dtype='D')
T = TensorProductSpace(comm, (L0, F1, F2), coordinates=(psi, rv, sp.Q.positive(sp.sin(theta))))
p = TrialFunction(T)
Math((div(grad(div(grad(p))))).tolatex(funcname='u', symbol_names={r: 'r', theta: '\\theta', phi: '\\phi'}))
q = TestFunction(T)
A = inner(div(grad(q)), div(grad(p)))
Explanation: Want to see what the regular 3-dimensional biharmonic equation looks like in spherical coordinates? This is extremely tedious to derive by hand, but in shenfun you can get there with the following few lines of code
End of explanation
L0 = FunctionSpace(8, 'C', domain=(0, np.pi))
F1 = FunctionSpace(8, 'F', dtype='D')
F2 = FunctionSpace(8, 'F', dtype='D')
T = TensorProductSpace(comm, (L0, F1, F2))
p = TrialFunction(T)
Math((div(grad(div(grad(p))))).tolatex(funcname='u'))
Explanation: I don't know if this is actually correct, because I haven't derived it by hand and I haven't seen it printed anywhere, but at least I know the Cartesian equation is correct:
End of explanation
<END_TASK>