Problem solving Harmonic Oscillator

General
,
1

Hello everyone,

I am trying to solve the following simple harmonic oscillator problem:

\begin{align} \frac{du(t)}{dt} &= v(t) \\ \frac{dv(t)}{dt} &= -\frac{k_1}{m_{1}}u(t) \\ u(0) &= u_{\text{init}}\\ v(0) &= v_{\text{init}} \end{align}

with this weak form:

\begin{align} a((u,u_\text{test}), \ (v, v_\text{test})) &= \int_{0}^{T} u'(t) \cdot u_\text{test} -v(t) \cdot u_\text{test} \ dt \ + \ ... \\ & \quad \int_{0}^{T} v'(t) \cdot v_\text{test} + \frac{k_{1}}{m_{1}} u(t) \cdot v_\text{test} \ dt \\ L(u_\text{test}, \ v_\text{test}) &= 0 \\ \\ \Rightarrow A\vec{x} = L \\ \end{align}

I think I got the weak form down correctly, however, I am not getting any sensible solution, which I assume is connected to the boundary values. This is my code:

# %%
import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh, default_scalar_type, mesh, plot
import ufl
from ufl import dot, dx
import basix
import matplotlib.pyplot as plt
from dolfinx.fem.petsc import LinearProblem

# Define parameters
k_1 = 1.0
m_1 = 1.0
v_0 = 0.0
u_0 = 1.0

# Create mesh and define function space
domain = mesh.create_interval(MPI.COMM_WORLD, 100, (0, 5))
element = basix.ufl.element("CG", domain.basix_cell(), 1)
V = fem.functionspace(domain, basix.ufl.mixed_element([element, element]))

# Define trial and test functions
u, v = ufl.TrialFunctions(V)
u_test, v_test = ufl.TestFunctions(V)

# Define variational problem
a = (dot(u.dx(0), u_test) - dot(v, u_test)) * dx             # Displacement
a += (dot(v.dx(0), v_test) + k_1/m_1 * dot(u, v_test)) * dx  # Velocity

f = fem.Constant(domain, 0.0) # No driving force
L = dot(f, u_test) * dx
L += dot(f, v_test) * dx

# Apply boundary conditions
ic_facets = mesh.locate_entities_boundary(
    domain, dim=0, marker=lambda x: np.isclose(x[0], 0.0)
)
dofs = fem.locate_dofs_topological(V=V.sub(0), entity_dim=0, entities=ic_facets)
bc_u = fem.dirichletbc(value=default_scalar_type(u_0), dofs=dofs, V=V.sub(0))
dofs = fem.locate_dofs_topological(V=V.sub(1), entity_dim=0, entities=ic_facets)
bc_v = fem.dirichletbc(value=default_scalar_type(v_0), dofs=dofs, V=V.sub(1))
bcs = [bc_u, bc_v]

# Solve problem
problem = LinearProblem(a, L, bcs=bcs)
solution = problem.solve()

print("Solved!")

# Split the solution to obtain the individual components
uh_u, uh_v = solution.split()
print(solution.x.array.real)

The code above results in an “inf” solution vector. Not prescribing any boundary conditions results in a zero vector solution, which makes sense if one assumes 0 initial velocity and displacement. This is why I assume that I am not prescribing the boundary conditions correctly.

Also, I am not sure if I understand the split function on the solution vector. I would assume that I would get the u and v vector respectively, however this does not seem to be the case, since the length of those arrays is the same as the solution vector (twice the number of mesh points, I would assume that it would be half of that, for u and v respectively)

Thanks in advance!

2

My concern here is that perhaps you’re using the wrong tool for the job. The problem you’ve stated is an initial value problem (rather than a boundary value problem, for which finite elements shine). Wouldn’t you be better suited with one of the many ODE solvers as interfaced with, e.g., SciPy?

https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#scipy.integrate.solve_ivp

3

Hi Nate,
thanks for your response. I am aware that FEM might not be the classical approach to solve a IVP, however, it should still be possible. There are a class of methods called Space-Time Finite Element methods, where you discretize the space and time with a finite element approximation. The pendulum is only one of the simplest mechanics examples I could think of, and I wanted to test fenicsx with that as an example.
I am pretty sure the issue is my implementation, not the method. Do you maybe see an issue with my code?

4

I think the problem comes from the fact that the system is non-symmetric. It seems to work with:

problem = LinearProblem(a, L, bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)
3 Likes
5

Hi @bleyerj,
thank you very much for your suggestion, this now gives me the expected solution!

However, I am having trouble post-processing the result, due to not quite understanding how the split function works.

This is my code for plotting:

# Solve problem
problem = LinearProblem(a, L, bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    })
solution = problem.solve()

print("Solved!")


# Split the solution to obtain the individual components
uh_u, uh_v = solution.split()
print(solution.x.array.real)

# Extract solution values at each mesh point
u_values = uh_u.x.array.real[::2]
v_values = uh_v.x.array.real[1::2]

# Plotting
plt.figure()
plt.plot(u_values, label="Displacement (u)")
plt.plot(v_values, label="Velocity (v)")
plt.xlabel("Position")
plt.ylabel("Values")
plt.legend()
plt.show()

image

Even though I use the split command to separate the displacement and velocity, both field are in each uh_u and uh_v. I tried to just take every second value, however this is not quite right for the first few values and far away from an easy automatic way.

Does someone have an idea what I might be doing wrong, and how I could fix it?

6

You may want to apply collapse() to uh_u and/or uh_v. After that, you will not need to slice the arrays.

7

Consider the following for post processing:

# Solve problem
problem = LinearProblem(a, L, bcs=bcs,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    })
solution = problem.solve()

print("Solved!")

from dolfinx import geometry
# Compute FE function at visualisation points
n_pts = 100
x = np.vstack((np.linspace(0.0, 5.0, n_pts), [np.zeros(n_pts)] * 2)).T
bbox = geometry.bb_tree(domain, domain.topology.dim)
candidate_cells = geometry.compute_collisions_points(bbox, x)
collided_cells = geometry.compute_colliding_cells(
    domain, candidate_cells, x)
collided_cells = collided_cells.array[collided_cells.offsets[:-1]]

# Extract solution values at each mesh point
u_values = solution.sub(0).collapse().eval(x, collided_cells)
v_values = solution.sub(1).collapse().eval(x, collided_cells)
# Plotting
plt.figure()
plt.plot(x[:, 0], u_values, label="Displacement (u)")
plt.plot(x[:, 0], v_values, label="Velocity (v)")
plt.xlabel("Position")
plt.ylabel("Values")
plt.legend()
plt.show()

image

1 Like