Wrong solution to partial differential equation in polar coordinate

1

Dear FEniCS community,

I would like to solve the following partial differential equation (PDE)
v_0 \cos(\theta-\phi) \frac{\partial}{\partial r} T(r, \phi, \theta)+ v_0 \sin(\theta-\phi)\frac{1}{r}\frac{\partial}{\partial \phi} T(r, \phi, \theta) \\\quad+ \gamma v_0 \sin(\theta-\phi)\frac{1}{r}\frac{\partial}{\partial \theta}T(r, \phi, \theta)+ \frac{g}{2} \frac{\partial^2}{\partial \theta^2} T(r,\phi, \theta) = -1

with domain D_{r_0}\times [0, 2\pi], where D_{r_0} is disk with radius r_0.
r and \phi are the polar coordinate of the disk. \theta\in [0,2\pi].

The boundary conditions are Dirichlet boundary conidtion on the disk edge T(r,\phi,\theta)|_{r=r_0} = 0 and periodic boundary conditions with respect to
\phi and \theta, T(r, \phi+2\pi, \theta) = T(r,\phi, \theta+2\pi) = T(r, \phi, \theta).

Below is the code.

import dolfin as df
import ufl
from math import pi
import matplotlib.pyplot as plt

v0 = 1.
gamma = 7.
g = 10.**3
r0 = 1.0
angle_period = 2*pi

class PeriodicBoundary(df.SubDomain):
    def inside(self, x, on_boundary):
        return bool((df.near(x[1], 0) or df.near(x[2], angle_period)) and
                    (not ((df.near(x[1], 0) and df.near(x[2], angle_period)) or
                     (df.near(x[1], angle_period) and df.near(x[2], 0)))) and
                    on_boundary)

    def map(self, x, y):
        if df.near(x[1], angle_period) and df.near(x[2], angle_period):
            y[0] = x[0]
            y[1] = x[1] - angle_period
            y[2] = x[2] - angle_period
        elif df.near(x[1], angle_period):
            y[0] = x[0]
            y[1] = x[1] - angle_period
            y[2] = x[2]
        else:   # near(x[2], 1)
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - angle_period


# set mesh
nn = 10
nx = nn
ny = nn
nz = nn
# Create mesh and define function spaces
mesh = df.BoxMesh(
        df.Point(0., 0., 0.), df.Point(r0, angle_period, angle_period),
        nx, ny, nz)
pbc = PeriodicBoundary()
V = df.FunctionSpace(mesh, 'P', 2, constrained_domain=pbc)

# Define boundaries
def outer_boundary(x, on_boundary):
    tol = 1.E-14
    return df.near(x[0], r0, tol) and on_boundary


bc = df.DirichletBC(V, df.Constant(0), outer_boundary)

u = df.TrialFunction(V)
v = df.TestFunction(V)
xs = df.SpatialCoordinate(mesh)

F = v0*ufl.cos(xs[2]-xs[1])*u.dx(0)*v*xs[0]*df.dx \
    + v0*ufl.sin(xs[2]-xs[1])*u.dx(1)*v*df.dx \
    + v0*gamma*ufl.sin(xs[2]-xs[1])*u.dx(2)*v*df.dx \
    - g/2*u.dx(2)*v.dx(2)*xs[0]*df.dx \
    + 1.*v*xs[0]*df.dx

a, L = df.lhs(F), df.rhs(F)
u = df.Function(V)
df.solve(a == L, u, bc)

# verify u is indeed the solution of the pde
w = df.project(
        v0*ufl.cos(xs[2]-xs[1])*u.dx(0)
        + v0/xs[0]*ufl.sin(xs[2]-xs[1])*u.dx(1)
        + v0*gamma/xs[0]*ufl.sin(xs[2]-xs[1])*u.dx(2)
        + g/2*u.dx(2).dx(2) + 1.,
        V
        )

# expecting w should be near 0.
# however the command below gives the value of w
#   [93.26907903 363.77849385 297.36876578 ...  -168.24578545 -420.92284536 -1026.43193945]
print(w.compute_vertex_values())

The problem is that FEniCS seems to give incorrect solution to the PDE.
When I substitute the solution back into the PDE, it does not give zero as expected.

How can I solve this problem?
Thank you for your time and consideration.

I am using FEniCS 2019.1.0 installed via conda-forge on Manjaro Linux.

2

Hi James,

I think you may have an error in the definition of your coordinates. Spherical coordinates are usually defined as:

\{r,\phi,\theta\}

with
r \in [0,\infty)\\ \phi \in [0, 2\pi]\\ \theta \in [0, \pi]

Which makes the periodic boundary condition:
T(r,\phi + 2\pi,\theta) = T(r,\phi,\theta +\pi) = T(r,\phi,\theta)

3

Hi Diego. Thanks for the response.
It is not spherical coordinates I am using.
r and \phi depict the position of a particle in polar coordinate.
\theta actually represents the orientation of the particle.
The function T(r, \phi, \theta) depends on both position r and \phi and orientation \theta.
It is actually more like a cylinder with the upper and lower surfaces identified (periodic boundary condition for \theta).
Sorry that this notation resembling the spherical coordinates causes confusion.