I’m currently trying to model the pressure field generated by a focussed ultrasound transducer in small a cavity. To that end, I’ve been looking into ways to solve the wave equation numerically. To quickly experiment with the concepts involved, I’ve set up a dummy problem which consists in modelling propagation of a pressure wave imposed at the boundary of a simple 2d domain.
According to the principle demonstrated in the Diffusion of a Gaussian function example, the general approach I’ve been following consists in approximating the time derivative using a finite difference scheme and utilizing the finite element method handled by FEniCSx for the spatial discretization.
Substituting the approximation of the derivative in the wave equation
\partial_{tt}^2 p = c_0^2 \nabla^2 p
yields
\frac{p_{n+1} - 2p_n + p_{n-1}}{\Delta t^2} = c_0^2 \nabla^2 p_n
where p_{n-1}, p_n and p_{n+1} correspond to the pressures at time t-\Delta t , t and t+\Delta t respectively.
To be interpreted by FEniCS the previous expression has to be rewritten under the weak form which leads to
a = \int (p_{n+1} v + c_0^2 \Delta t^2 \nabla p \cdot \nabla v) dx
L = \int ( 2p_n - p_{n-1}) v dx
After having built the mesh and implementing the equtation in FEniCSx, I impose a Dirichlet boundary condition at the transducer surface such that p(\bf{x}, t)\bigr\rvert_{x = S_{tx}} = \sin(\omega t). I then obtain the final solution by iterating over the time steps and repeatadly substituting p_{n} in p_{n-1} and p_{n+1} in p_n. Dispite choosing small spatial and temporal steps (\Delta x \approx \frac{\lambda}{10} and \Delta t = \frac{T}{100}) the pressure field I get shows a significant amount of attenuation. I assume that this effect comes from the temporal finite difference scheme as increasing the spatial resolution doesn’t seem to help, I’m quite surprised however to see such a visible effect with a time step in the order of a hundredth of the wave period. Does anybody know how I could fix this?
Here is my current version of the code:
import dolfinx, gmsh, gmsh_helpers, ufl, time
from petsc4py import PETSc
from mpi4py import MPI
from tqdm import tqdm
import numpy as np
# Physical quantities
tx_freq = 200e3 # Hz
c0 = 1481 # m.s-1
omega = 2*np.pi*tx_freq
k0 = omega / c0 # wavenumber
wave_len = c0 / tx_freq
# Transducer properties
tx_radius = 15e-2
tx_aperture_radius = 15e-2
alpha_aperture = np.arcsin(tx_aperture_radius / (2*tx_radius))
tx_marker, boundary_marker, ac_domain_marker = 1, 2, 3
# Domain parameters
topology_dim = 2
dims = (15e-2, 20e-2)
# Spatial discretisation parameters
degree = 3
n_wave = 10
h_elem = wave_len / n_wave
# Temporal discretisation parameters
dt = 1 / (100 * tx_freq)
simulation_duration = 15e-5
t_mesh = np.arange(0, simulation_duration, dt)
Mesh generation
gmsh.initialize()
model = gmsh.model
model.add("DummyDomain")
# Points defining the geometry of the domain
points = []
# Limits of the transducer surface
tx_center = model.geo.addPoint(tx_radius*np.cos(alpha_aperture),0,0, meshSize=h_elem, tag=10) # Transducer center coordinates
points.append(model.geo.addPoint(tx_radius*np.cos(-alpha_aperture+np.pi) + tx_radius*np.cos(alpha_aperture), tx_radius*np.sin(-alpha_aperture+np.pi), 0, meshSize=h_elem, tag=0))
points.append(model.geo.addPoint(tx_radius*np.cos(alpha_aperture+np.pi) + tx_radius*np.cos(alpha_aperture), tx_radius*np.sin(alpha_aperture+np.pi), 0, meshSize=h_elem, tag=1))
# Corners of the acoustic domain
rect_domain_geom = [[0., -1/2], [1., -1/2], [1., 1/2], [0., 1/2]]
for rect_geom in rect_domain_geom:
points.append(model.geo.addPoint(*[gg*dd for gg, dd in zip(rect_geom, dims)],0, meshSize=h_elem, tag=len(points)))
# Lines defining the limits of the domain
lines = []
# Transducer surface (line)
lines.append(model.geo.addCircleArc(points[0], tx_center, points[1]))
# Acoustic domain boundaries
for pt in points[1:]:
lines.append(model.geo.addLine(points[pt], points[(pt+1)%6]))
# Junction of the individual boubaries and definition of the acoustic domain surface
curveloop = model.geo.addCurveLoop(lines)
domain = model.geo.addPlaneSurface([curveloop])
# This command is mandatory and synchronize CAD with GMSH Model. The less you launch it, the better it is for performance purpose
gmsh.model.geo.synchronize()
# Assigns the various geometry elements to physical groups
gmsh.model.addPhysicalGroup(1, [lines[0]], tx_marker)
gmsh.model.setPhysicalName(1, tx_marker, "Tx")
gmsh.model.addPhysicalGroup(1, lines[1:], boundary_marker)
gmsh.model.setPhysicalName(1, boundary_marker, "Domain boundary")
gmsh.model.addPhysicalGroup(2, [domain], ac_domain_marker)
gmsh.model.setPhysicalName(2, ac_domain_marker, "Acoustic domain")
# Mesh generation and
model.mesh.generate(topology_dim)
#gmsh.write("DummyDomain.msh")
# Gmsh to dolfinx mesh
mesh, cell_tags, facet_tags = gmsh_helpers.gmsh_model_to_mesh(model, cell_data=True, facet_data=True, gdim=topology_dim)
# Finalize GMSH
gmsh.finalize()
Variational problem
P = dolfinx.FunctionSpace(mesh, ("Lagrange", degree))
# Retreive the transducer facet indexes to apply the bc later on
tx_facets = facet_tags.indices[facet_tags.values == tx_marker]
tx_dofs = dolfinx.fem.locate_dofs_topological(P, topology_dim - 1, tx_facets)
p_source = dolfinx.Function(P)
# Initial state of the acoustic domain
p_n_1 = dolfinx.Function(P) # p(t - dt)
with p_n_1.vector.localForm() as loc:
loc.set(0)
dolfinx.cpp.la.scatter_forward(p_n_1.x)
p_n = dolfinx.Function(P) # p(t)
with p_n.vector.localForm() as loc:
loc.set(0)
dolfinx.cpp.la.scatter_forward(p_n.x)
p_n1 = dolfinx.Function(P) # p(t + dt)
with p_n1.vector.localForm() as loc:
loc.set(0)
dolfinx.cpp.la.scatter_forward(p_n1.x)
# Define variational problem
p_trial, v = ufl.TrialFunction(P), ufl.TestFunction(P)
a = p_trial*v*ufl.dx + dt*dt*c0*c0*ufl.inner(ufl.grad(p_trial), ufl.grad(v))*ufl.dx
L = 2*p_n*v*ufl.dx - p_n_1*v*ufl.dx
xdmf_file = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "tx_prop_100T10nwave.xdmf", "w")
xdmf_file.write_mesh(mesh)
for ii, tt in enumerate(tqdm(t_mesh)):
# Compute the pressure at the transducer surface
with p_source.vector.localForm() as loc:
loc.set(np.sin(omega*tt))
tx_bc = dolfinx.DirichletBC(p_source, tx_dofs)
if ii == 0:
print('Setting up the problem...')
# Convert a to the matrix form
A = dolfinx.fem.assemble_matrix(a, bcs=[tx_bc])
A.assemble()
b = dolfinx.fem.create_vector(L)
# Create the linear solver using PETSc
solver = PETSc.KSP().create(mesh.mpi_comm())
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
dolfinx.fem.assemble_vector(b, L)
# Apply Dirichlet boundary condition to the vector
dolfinx.fem.apply_lifting(b, [a], [[tx_bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(b, [tx_bc])
# Solve linear problem
solver.solve(b, p_n1.vector)
dolfinx.cpp.la.scatter_forward(p_n1.x)
# Update the two solutions of the fintie difference discretization scheme
with p_n.vector.localForm() as p_loc_n, p_n_1.vector.localForm() as p_loc_n_1:
p_loc_n.copy(p_loc_n_1)
# a.copy(b) a to b
with p_n1.vector.localForm() as p_loc_n1, p_n.vector.localForm() as p_loc_n:
p_loc_n1.copy(p_loc_n)
p_n1.name = "pressure"
xdmf_file.write_function(p_n1, tt)
xdmf_file.close()
By the way, any suggestions to improve the general quality of the code are welcome!