# Non-orthogonal curvilinear coordinates

In this notebook we solve Poisson's equation on a 2D wavy domain, using non-orthogonal basis vectors.


```python
import sympy as sp
import matplotlib.pyplot as plt
from IPython.display import Math
from shenfun import *
from shenfun.la import Solver2D
```

Define the non-orthogonal curvilinear coordinates. Here `psi=(u, v)` are the new coordinates and `rv` is the position vector

$$
\vec{r} = u \mathbf{i} + v\left(1- \frac{\sin 2u}{10} \right) \mathbf{j}
$$

where $\mathbf{i}$ and $\mathbf{j}$ are the Cartesian unit vectors in $x$- and $y$-directions, respectively.


```python
u = sp.Symbol('x', real=True, positive=True)
v = sp.Symbol('y', real=True)
psi = (u, v)
rv = (u, v*(1-sp.sin(2*u)/10))
```

Now choose basis functions and create tensor product space. Notice that one has to use complex Fourier space and not the real, because the integral measure is a function of u.


```python
N = 20
#B0 = FunctionSpace(N, 'C', bc=(0, 0), domain=(0, 2*np.pi))
B0 = FunctionSpace(N, 'F', dtype='D', domain=(0, 2*np.pi))
B1 = FunctionSpace(N, 'L', bc=(0, 0), domain=(-1, 1))

T = TensorProductSpace(comm, (B0, B1), dtype='D', coordinates=(psi, rv, sp.Q.negative(sp.sin(2*u)-10) & sp.Q.negative(sp.sin(2*u)/10-1)))
p = TrialFunction(T)
q = TestFunction(T)
b = T.coors.get_covariant_basis()
T.coors.sg
```


```python
sp.Matrix(T.coors.get_covariant_metric_tensor())
```


```python
Math(T.coors.latex_basis_vectors(symbol_names={u: 'u', v: 'v'}))
```

Plot the mesh to see the domain.


```python
mesh = T.local_cartesian_mesh()
x, y = mesh
plt.figure(figsize=(10, 4))
for i, (xi, yi) in enumerate(zip(x, y)):
    plt.plot(xi, yi, 'b')
    plt.plot(x[:, i], y[:, i], 'r')
```

Print the Laplace operator in curvilinear coordinates. We use `replace` to simplify the expression.


```python
dp = div(grad(p))
g = sp.Symbol('g', real=True, positive=True)
replace = [(1-sp.sin(2*u)/10, sp.sqrt(g)), (sp.sin(2*u)-10, -10*sp.sqrt(g)), (5*sp.sin(2*u)-50, -50*sp.sqrt(g))]
Math((dp*T.coors.sg**2).tolatex(funcname='p', symbol_names={u: 'u', v: 'v'}, replace=replace))
```

Solve Poisson's equation. First define a manufactured solution and assemble the right hand side


```python
ue = sp.sin(2*u)*(1-v**2)
f = (div(grad(p))).tosympy(basis=ue, psi=psi)
fj = Array(T, buffer=f*T.coors.sg)
f_hat = Function(T)
f_hat = inner(q, fj, output_array=f_hat)
```

Then assemble the left hand side and solve using a generic 2D solver


```python
M = inner(q, div(grad(p))*T.coors.sg)
#M = inner(grad(q*T.coors.sg), -grad(p))
u_hat = Function(T)
Sol1 = Solver2D(M)
u_hat = Sol1(f_hat, u_hat)
uj = u_hat.backward()
uq = Array(T, buffer=ue)
print('Error =', np.linalg.norm(uj-uq))
```


```python
for i in range(len(M)):
    print(len(M[i].mats[0].keys()), len(M[i].mats[1].keys()), M[i].mats[0].measure, M[i].mats[1].measure)
```

Plot the solution in the wavy Cartesian domain


```python
xx, yy = T.local_cartesian_mesh()
plt.figure(figsize=(12, 4))
plt.contourf(xx, yy, uj.real)
plt.colorbar()
```

Inspect the sparsity pattern of the generated matrix on the left hand side


```python
from matplotlib.pyplot import spy
plt.figure()
spy(Sol1.mat, markersize=0.1)
```


```python
from scipy.sparse.linalg import eigs
mats = inner(q*T.coors.sg, -div(grad(p)))
Sol1 = Solver2D(mats)
BB = inner(p, q*T.coors.sg)
Sol2 = Solver2D([BB])
f = eigs(Sol1.mat, k=20, M=Sol2.mat, which='LM', sigma=0)
```


```python
mats
```


```python
l = 10
u_hat = Function(T)
u_hat[:, :-2] = f[1][:, l].reshape(T.dims())
plt.contourf(xx, yy, u_hat.backward().real, 100)
plt.colorbar()
```

90e73c039f75d0592e3a9acdc7e4fd075d1b5153
ipynb
Jupyter Notebook
binder/non-orthogonal-poisson.ipynb
spectralDNS/shenfun
956633aa0f1638db5ebdc497ff68a438aa22b932
BSD-2-Clause
2017-06-17T13:30:27.000Z
2022-03-20T02:33:47.000Z
binder/non-orthogonal-poisson.ipynb
spectralDNS/shenfun
956633aa0f1638db5ebdc497ff68a438aa22b932
BSD-2-Clause
2017-05-16T06:53:04.000Z
2022-02-04T10:40:44.000Z
binder/non-orthogonal-poisson.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