Hello everyone,
I am trying to solve the following problem in its variational form:
In the first step, I want to solve the homogeneous eigenvalue problem, where the right-hand side is set to 0. I am using SLEPc for this. In the second step, I plan to take the computed eigenvalue, \omega_0, and substitute it into the right-hand side to solve the inhomogeneous equation.
Now, to my issue: I am unsure how to evaluate \nabla \hat{p} at (x_{\text{ref}}) in the variational form. Is there a way to implement something like a delta function that equals 1 at x_{\text{ref}} and 0 elsewhere?
I also have another problem related to the shapes in the matrix definition for matrix D:
File "/home/hell/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/ufl/tensoralgebra.py", line 166, in __new__
raise ValueError(f"Shapes do not match: {ufl_err_str(a)} and {ufl_err_str(b)}")
ValueError: Shapes do not match: <ComponentTensor id=139710411729280> and <Argument id=139710413172800>
I’m working with FEniCSx 0.8
Thank you in advance. If something is unclear, I could provide more information.
Code:
from dolfinx.fem import locate_dofs_geometrical, functionspace, Function, dirichletbc
from mpi4py import MPI
from numpy.array_api import complex128
from petsc4py import PETSc
from ufl import TrialFunction, TestFunction, dx, inner, grad, SpatialCoordinate
from dolfinx import default_scalar_type, fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector
from dolfinx.mesh import CellType, create_rectangle, locate_entities
import numpy as np
# ---------- Domain definition ----------
w = 0.5
h = 0.1
nx = 100
ny = 10
cell_w = w / nx
msh = create_rectangle(
MPI.COMM_WORLD, np.array([[0, 0], [w, h]]), np.array([nx, ny]), CellType.quadrilateral
)
cell_w = w / nx
msh = create_rectangle(
MPI.COMM_WORLD, np.array([[0, 0], [w, h]]), np.array([nx, ny]), CellType.quadrilateral
)
# ---------- Define functions and area for c1, c2, rho1, rho2 and flame region ----------
Q = functionspace(msh, ("DG", 0))
def Omega_0(x):
return x[0] <= w/2
def Omega_1(x):
return x[0] >= w/2
def flame_region(x):
return np.logical_and(x[0] >= w/2 - cell_w, x[0] <= w/2 + cell_w)
def ref_point(x):
return np.logical_and(x[0] >= w/2-cell_w, x[0] <= w/2)
cells_1 = locate_entities(msh, msh.topology.dim, Omega_0)
cells_2 = locate_entities(msh, msh.topology.dim, Omega_1)
cells_flame_region = locate_entities(msh, msh.topology.dim, flame_region)
ref_point_region = locate_entities(msh,msh.topology.dim, ref_point)
c = Function(Q)
c.x.array[cells_1] = np.full_like(cells_1, 347.18, dtype=default_scalar_type)
c.x.array[cells_2] = np.full_like(cells_2, 694.36, dtype=default_scalar_type)
rho = Function(Q)
rho.x.array[cells_1] = np.full_like(cells_1, 4, dtype=default_scalar_type)
rho.x.array[cells_2] = np.full_like(cells_2, 1, dtype=default_scalar_type)
# ---------- Define functionspace, functions and BCs ----------
degree = 1
V = functionspace(msh, ('Lagrange', degree))
u = TrialFunction(V)
v = TestFunction(V)
e = Function(V)
o = Function(V)
x = SpatialCoordinate(msh)
dofs = locate_dofs_geometrical(V, lambda x: np.isclose(x[0], w))
bc = [dirichletbc(default_scalar_type(0), dofs, V)]
# Eigenvalues from SLEPc
omega0 = [(1709.4992392760605+3.1592511086743223e-17j), (4363.241033918917-3.6488298900508794e-17j), (7018.447703170593-2.4364863230542956e-15j), (10443.252332619693-3.7277152271679734e-15j)]
rho_ref = 4
n_ = 0.01
delta_f = cell_w
tau = 10e-4
k = 1 / rho_ref * n_ / delta_f * np.exp(1j * tau * omega0[0])
e.x.array[cells_flame_region] = np.full_like(cells_flame_region, k, dtype=complex128)
o.x.array[ref_point_region] = np.full_like(ref_point_region,k, dtype=complex128)
# ---------- Define D matrix for flame model ----------
A = inner((1 / rho) * grad(u), grad(v)) * dx
M = inner(1 / (c ** 2 * rho) * u, v) * dx
A = assemble_matrix(fem.form(A), bcs=bc)
A.assemble()
M = assemble_matrix(fem.form(M), bcs=bc)
M.assemble()
D = inner(grad(u)*e, u) * dx
D = assemble_vector(fem.form(D))
D.assemble()
P = A + D
# ---------- Solve eigenvalue problem with SLEPc ----------