Greetings everyone,
I am trying to create a code to compute the acoustic specter of a cavity, I want it to compute the specter and the fields for given frequencies and to save it to an xdmf file.
I try to adapt a tutorial found here :
The first part of the code (not shown here) appears to work properly and gives me back a specter however when I open the xdmf file in paraview, the fields appears to be zero everywhere. This is surprising since I get a coherent specter and I can’t figure out why this is the case.
Here is my code :
import numpy as np
import ufl
import time
from dolfinx import geometry
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, Constant
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_unit_square
from ufl import dx, grad, inner
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import XDMFFile, gmshio, VTXWriter
import dolfinx as dfx
# approximation space polynomial degree
deg = 1
n_elem1 = 64
n_elem2 = 64
mesh = create_unit_square(MPI.COMM_SELF, n_elem1, n_elem2)
n = ufl.FacetNormal(mesh)
#frequency range definition
f_axis = np.arange(100, 2000, 100)
# Test and trial function space
V = FunctionSpace(mesh, ("CG", deg))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = Function(V)
# Source amplitude
Q = 0.001
#Source definition position = (Sx,Sy)
Sx = 0.0
Sy = 0.0
#Narrow normalized gauss distribution
alfa = 0.015
delta_tmp = Function(V)
delta_tmp.interpolate(lambda x : 1/(np.abs(alfa)*np.sqrt(np.pi))*np.exp(-(((x[0]-Sx)**2+(x[1]-Sy)**2)/(alfa**2))))
int_delta_tmp = assemble_scalar(form(delta_tmp*dx))
delta = delta_tmp/int_delta_tmp
# fluid definition
c0 = 340
rho_0 = 1.225
omega = Constant(mesh, PETSc.ScalarType(1))
k0 = Constant(mesh, PETSc.ScalarType(1))
#building the weak form
f = 1j*rho_0*omega*Q*delta
a = inner(grad(u), grad(v)) * dx - k0**2 * inner(u, v) * dx
L = inner(f, v) * dx
uh = Function(V)
#building the problem
problem = LinearProblem(a, L, bcs = [],u = uh, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
# frequency for loop
print(" ")
for nf in range(0,len(f_axis)):
#print ("\033[A \033[A")
freq = f_axis[nf]
#print("Computing frequency: " + str(freq) + " Hz...")
all_u_values = []
# Compute solution
omega.value =f_axis[nf]*2*np.pi
k0.value = 2*np.pi*f_axis[nf]/c0
problem.solve()
#if f_axis[nf]%100 == 0 : #create an XDMFFile every 100 Hz
with XDMFFile(mesh.comm, "Specter_out/air_solution" + str(f_axis[nf]) + "Hz.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(uh)
Thank you for reading, have a nice day !
