Solve 1d initial value problem

General
1

Hi all,

I would like to solve 1D Helmholtz equation with k>0. The analytic solution is \sin(kx)/k.

\frac{\partial^2 \phi}{\partial x^2} +k^2\phi = 0\\ \phi(x=0) = 0\\ \frac{\partial\phi}{\partial x}\bigg|_{x=0} = 1

Following this tutorial, I supply the initial condition by applying Dirichlet and Neumann condition at the origin. Here is my attempt:

import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem import FunctionSpace

r_start = 1e-6
r_end = 20
nx = 1000
k = np.pi

# define mesh and function space
domain = mesh.create_interval(MPI.COMM_WORLD, nx, [r_start, r_end])
V = FunctionSpace(domain, ("CG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)
F = - ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + k ** 2 * ufl.inner(u, v) * ufl.dx

# define boundary and labels
boundaries = {
    1 : lambda x : np.isclose(x[0], r_start),
    2 : lambda x : np.isclose(x[0], r_end),
}

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for marker, locator in boundaries.items():
    facets = mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == 'Dirichlet':
            u_D = fem.Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = fem.locate_dofs_topological(V, fdim, facets)
            self._bc = fem.dirichletbc(u_D, dofs)
        elif type == 'Neumann':
                self._bc = ufl.inner(values, v) * ds(marker)
        elif type == 'Robin':
            self._bc = values[0] * ufl.inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type

# define boundary conditions
n = ufl.FacetNormal(domain)
uD = lambda x : np.sin(k * x[0]) / k
uN = - ufl.cos(k * x[0])

boundary_conditions = [
    BoundaryCondition("Dirichlet", 1, uD),
    BoundaryCondition("Neumann", 1, uN),
]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

# solve the problem
a = ufl.lhs(F)
L = ufl.rhs(F)
problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

plt.plot(V.tabulate_dof_coordinates()[:, 0], uh.x.array)

Using this script I obtained the wrong answer. I wonder if there is a simpler approach to this type of problem? Thanks a lot !!

2

Isn’t your problem an initial value problem and not a boundary value problem? I.e. not a problem defined on a finite domain.

3

Yes. Is it not possible (or not appropriate) to solve this kind of problem in fenics? I have tried the following script for 1d poisson which agrees with the exact solution.

script 
import numpy as np
import ufl
from mpi4py import MPI
from dolfinx import fem, mesh
from dolfinx.fem import FunctionSpace
from petsc4py.PETSc import ScalarType

x_start = 0.
x_end = 1.
nx = 100

domain = mesh.create_interval(MPI.COMM_WORLD, nx, [x_start, x_end])
V = FunctionSpace(domain, ("CG", 1))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(domain)

f = fem.Constant(domain, ScalarType(1))
F = - ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx + f * v * ufl.dx

boundaries = {
    1 : lambda x : np.isclose(x[0], x_start),
    2 : lambda x : np.isclose(x[0], x_end),
}

facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for marker, locator in boundaries.items():
    facets = mesh.locate_entities(domain, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))

facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)

class BoundaryCondition():
    def __init__(self, type, marker, values):
        self._type = type
        if type == 'Dirichlet':
            u_D = fem.Function(V)
            u_D.interpolate(values)
            facets = facet_tag.find(marker)
            dofs = fem.locate_dofs_topological(V, fdim, facets)
            self._bc = fem.dirichletbc(u_D, dofs)
        elif type == 'Neumann':
                self._bc = ufl.inner(values, v) * ds(marker)
        elif type == 'Robin':
            self._bc = values[0] * ufl.inner(u-values[1], v)* ds(marker)
        else:
            raise TypeError("Unknown boundary condition: {0:s}".format(type))
    @property
    def bc(self):
        return self._bc

    @property
    def type(self):
        return self._type
    
n = ufl.FacetNormal(domain)
uD = lambda x : - x[0] ** 2 / 2. + x[0]
uN = ufl.inner(n, ufl.grad(uD(x)))

boundary_conditions = [
    BoundaryCondition("Dirichlet", 1, uD),
    BoundaryCondition("Neumann", 1, uN),
]

bcs = []
for condition in boundary_conditions:
    if condition.type == "Dirichlet":
        bcs.append(condition.bc)
    else:
        F += condition.bc

a = ufl.lhs(F)
L = ufl.rhs(F)
problem = fem.petsc.LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
x = np.linspace(x_start, x_end, 100)
ax.plot(x, - x ** 2 / 2. + x, label='exact')
ax.plot(V.tabulate_dof_coordinates()[:,0], uh.x.array, label='fenicsx', linestyle='--')
ax.legend()
ax.set_xlabel('$x$')
ax.set_ylabel('$u$')
plt.tight_layout()
plt.show()
4

This isn’t an issue with FEniCS, but with your numerical approach to the problem. Consider e.g. scipy.integrate.odeint — SciPy v1.11.1 Manual. There’s an example at the bottom of the page.

5

I understand this problem can be solved using scipy. For my purpose I am trying to see if this problem can be implemented in fenics since providing boundary conditions on one side seems to be equivalent to a initial condition. Thanks.

6

It might help if you write out the weak formulation of your problem, being careful to keep track of boundary conditions. Even if you truncate the infinite domain x \in (0, \infty), the imposition of BCs will not yield a well posed problem.

If you want to commit variational crimes, you could try strongly enforcing \phi(x=0) = 0 and introducing a weak imposition of \left.\frac{\partial \phi}{\partial x}\right|_{x=0} = 0; however, I think the error convergence decays to \mathcal{O}(h).