Hi,
I am trying to solve 1D Helmholtz equation using DG elements. I am getting correct eigenvalues but eigenvectors are discontinuos even if I implement the stabilization terms. I built this script upon this tutorial. Here is the MWE;
from dolfinx.fem import FunctionSpace, form, Function,dirichletbc,locate_dofs_topological
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.mesh import meshtags,locate_entities, create_interval
from ufl import TestFunction,TrialFunction,Measure,inner,grad,jump,avg,Circumradius, FacetNormal
import numpy as np
from slepc4py import SLEPc
from mpi4py import MPI
def OneDimensionalSetup(n_elem, x_MP = 0.2, r_outer = 0.25):
# 1 ------------------------------------ 2
# x=r_outer
""" This function builds one dimensional setup.
For boundaries, Tag 1 specifies left end and Tag 2 specifies right end.
Returns:
mesh, subdomains, facet_tags
"""
mesh = create_interval(MPI.COMM_WORLD, n_elem, [0.0, r_outer])
def subdomain_func(x, x_MP=x_MP, eps=1e-16):
x = x[0]
return np.logical_and(x_MP - eps <= x, x <= x_MP + eps)
tdim = mesh.topology.dim
marked_cells = locate_entities(mesh, tdim, subdomain_func)
fl = 0
subdomains = meshtags(mesh, tdim, marked_cells, np.full(len(marked_cells), fl, dtype=np.int32))
boundaries = [ (0, lambda x: np.isclose(x[0], x_MP)),
(1, lambda x: np.isclose(x[0], 0)),
(2, lambda x: np.isclose(x[0], r_outer))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full(len(facets), marker))
facet_indices = np.array(np.hstack(facet_indices), dtype=np.int32)
facet_markers = np.array(np.hstack(facet_markers), dtype=np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
return mesh, subdomains, facet_tags
n_elem = 50
mesh, subdomains, facet_tags = OneDimensionalSetup(n_elem)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tags)
dS = Measure("dS", domain=mesh, subdomain_data=facet_tags)
V = FunctionSpace(mesh, ("DG", 1))
# Dirichlet boundary conditions for left and right ends (equals to zero)
fdim = mesh.topology.dim-1
u_bc_left = Function(V)
u_bc_right = Function(V)
facets_left = np.array(facet_tags.indices[facet_tags.values == 1])
facets_right = np.array(facet_tags.indices[facet_tags.values == 2])
dofs_left = locate_dofs_topological(V, fdim, facets_left)
dofs_right = locate_dofs_topological(V, fdim, facets_right)
bc_left = dirichletbc(u_bc_left, dofs_left)
bc_right = dirichletbc(u_bc_right, dofs_right)
bcs=[bc_left,bc_right]
u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(mesh)
h = 2 * Circumradius(mesh)
c = 347
alpha = 10
gamma = 10
h_avg = avg(h)
mass = -c**2 * (inner(grad(u), grad(v))*dx + inner(inner(grad(u), n), v) * ds)
# Add DG/IP terms
mass += c**2 * (inner(jump(u, n), avg(grad(v)) )*dS + inner(avg(grad(u)), jump(v, n))*dS)
mass += -c**2*gamma/h_avg*inner(jump(u, n), jump(v, n))*dS
# Add Nitsche terms
mass += c**2*(inner(u, inner(grad(v), n)) * ds - alpha / h * inner(u, v) * ds)
mass_form = form(mass)
A = assemble_matrix(mass_form,bcs=bcs)
A.assemble()
stiffness = inner(u , v) * dx
stiffnes_form = form(stiffness)
C = assemble_matrix(stiffnes_form)
C.assemble()
target = 694*2*np.pi
nev=6
# This function prints the results of the eps_solver1 function below
def results(E):
print("******************************")
print("*** SLEPc Solution Results ***")
print("******************************")
print("Number of iterations of the method: %d" % E.getIterationNumber())
print("Solution method: %s" % E.getType())
nev, ncv, mpd = E.getDimensions()
print("Number of requested eigenvalues: %d" % nev)
print("Stopping condition: tol=%.4g, maxit=%d" % (E.getTolerances()))
print("Number of converged eigenpairs %d" % E.getConverged())
if E.getConverged() > 0:
print()
for i in range(E.getConverged()):
k = E.getEigenvalue(i)
w = np.sqrt(k)
f = w/np.pi/2
print("%6f, %6f" % (f.real, f.imag))
print()
# This function solves eigensystem
def eps_solver1(A, C, target, nev, two_sided=False, print_results=False):
E = SLEPc.EPS().create(MPI.COMM_WORLD)
C = - C
E.setOperators(A, C)
st = E.getST()
st.setType('sinvert')
E.setTarget(target)
E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE) # TARGET_REAL or TARGET_IMAGINARY
E.setDimensions(nev, SLEPc.DECIDE)
E.setTolerances(1e-15)
E.setFromOptions()
E.solve()
if print_results:
results(E)
return E
E = eps_solver1(A, C, target**2, nev=nev, print_results= True)
vr = Function(V)
vi = Function(V)
ind = 0
E.getEigenvector(ind, vr.vector, vi.vector)
import matplotlib.pyplot as plt
fig, ax = plt.subplots(2, figsize=(12, 6))
ax[0].plot(vr.x.array.real)
ax[1].plot(vr.x.array.imag)
plt.show()
This script outputs;
******************************
*** SLEPc Solution Results ***
******************************
Number of iterations of the method: 2
Solution method: krylovschur
Number of requested eigenvalues: 6
Stopping condition: tol=1e-15, maxit=100
Number of converged eigenpairs 7
694.093164, -0.000000
1388.745568, 0.000000
2084.517709, -0.000000
2781.972505, -0.000000
3481.676277, 0.000000
4184.199466, 0.000000
4890.117030, -0.000000
But the eigenvector I get has discontinuities;
My question is, how should I manipulate my weak form further if necessary?
