# Lighthill Acoustic´s Analogy

Hi community! (Firstly, sorry for mi English) I´m new in Fenics, and I´m trying to simulate an aerodynamic noise (from a Known velocity field (w) previously calculated with NavierStokes), but I have a problem with the Lightill Tensor
Here my code:

``````import meshio
import numpy as np
from dolfin import *

mesh_NS = Mesh()
with XDMFFile("mesh_NS.xdmf") as infile:

mvc = MeshValueCollection("size_t", mesh_NS, 2)
with XDMFFile("mf_NS.xdmf") as infile:
mf_NS = cpp.mesh.MeshFunctionSizet(mesh_NS, mvc)

T =70   # final time
num_steps =8000  # number of time steps
dt = T / num_steps # time step size
c = 340         # undisturbed speed of sound (m/s)
rho = 1.20          # density

W = VectorFunctionSpace(mesh_NS,'P', 2)

p=TrialFunction(W)
v=TestFunction(W)
w=Function(W)
w_1,w_2=as_vector((w,w))
w_11=inner(w_1,w_1)
w_12=inner(w_1,w_2)
w_22=inner(w_2,w_2)
def f(w):
return as_vector([[rho*(w_11.dx(0)+w_12.dx(1))],[rho*(w_12.dx(0)+w_22.dx(1))]])

p_0=Function(W)
p_1=Function(W)

n=FacetNormal(mesh_NS)
k=Constant(dt)
c=Constant(c)

L=((2*p_1) - p_0 + c*c*k*k*inner(f(w),v))*dx

timeseries_w = TimeSeries('NS_SC_200/velocitySC_200_series')

# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
timeseries_w.retrieve(w.vector(),t)

# Compute solution
solve(a == L, p, bcp)

#Update previous solution
p_0.assign(p_1)
p_1.assign(p)
``````

etc.

My problem is in the inner product(f(w),v), in relation with de shape of f(w) function. Here is my error:

``````    L=((2*p_1) - p_0 + c*c*k*k*inner(f(w),v))*dx
File "/.../anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/operators.py", line 169, in inner
return Inner(a, b)
File ".../anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/tensoralgebra.py", line 161, in __new__
error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
File "/.../anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/log.py", line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: <ListTensor id=140540322902224> and <Argument id=140540322720352>.
``````

Please!!! Does anyone how to solve it??

Please format your code using 3x` encapsulation, to format it properly.

`````````
def f(x):
return x
```
``````

I don´t know if I understand what you are appointing
I tried to do it in this way:

def f(w_11,w_12,w_22):
return as_vector([[rho*(w_11.dx(0)+w_12.dx(1))],[rho*(w_12.dx(0)+w_22.dx(1))]])

But it results in the same error:

File “/–/Acustica_lighthill.py”, line 75, in
L=((2p_1) - p_0 + cckk*inner(f(w_11,w_12,w_22),v))*dx
File “–/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/operators.py”, line 169, in inner
return Inner(a, b)
File “/–/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/tensoralgebra.py”, line 161, in new
error(“Shapes do not match: %s and %s.” % (ufl_err_str(a), ufl_err_str(b)))
File “/–/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Shapes do not match: and .

I´m lost…

I mean that you should format your code in the post such that it is readable and copyable for others. Ive modified your post with the suggested changes.

Ups, I didn´t notice ¡¡Many thanks!!

I can’t run your code, but it looks like `f()` is 2D and `v` is possible 3d. It would explain your mismatched shape error.

I can’t remember the legacy dolfinx syntax to check, it’s something like `print(mesh.geometry().dim())`

Sorry! My mesh it´s really tedious and I can´t write such a long code
If it helps I could try to upload de .xdmf files.
And the data for w is stored in a .h5 file (0,75 GB)

I tried to do what you said:

def f(w_11,w_12,w_22):
return as_vector([[rho*(w_11.dx(0)+w_12.dx(1))],[rho*(w_12.dx(0)+w_22.dx(1))], ])

And I got the same error Well, consider looking at `f().ufl_shape` and `v.ufl_shape` to make sure they are appropriate.

The way you’ve used `ufl.as_vector` is a little odd, consider the following:

``````In : import ufl

In : ufl.as_vector([, ]).ufl_shape
Out: (2, 1)

In : ufl.as_vector([1, 2]).ufl_shape
Out: (2,)
``````

which if you try to compute in an inner product throws a similar error:

``````In : ufl.inner(ufl.as_vector([, ]), ufl.as_vector([1, 2]))

Shapes do not match: <ListTensor id=140128052927232> and <ListTensor id=140127900825984>.
---------------------------------------------------------------------------
UFLException                              Traceback (most recent call last)
----> 1 ufl.inner(ufl.as_vector([, ]), ufl.as_vector([1, 2]))

``````
1 Like

Yes!! I´ve done what you said, and effectively, it was the problem in the inner product I changed the shape of f(), and now it seems that shapes match (both, f() and v, are(2,)).
I wrote

``````def f(w_11,w_12,w_22):
return as_vector(([rho*(w_11.dx(0)+w_12.dx(1)),rho*(w_12.dx(0)+w_22.dx(1))]))
``````

But…the problem remains with the sum of the inner product and the error says:

Can’t add expressions with different shapes.
Traceback (most recent call last):
File “/–.py”, line 85, in
L=((2p_1) - p_0 + cckk*inner(f(w_11,w_12,w_22),v))*dx
File “/–/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/exproperators.py”, line 212, in _add
return Sum(self, o)
File “/–/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/algebra.py”, line 53, in new
error(“Can’t add expressions with different shapes.”)
File “/–/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can’t add expressions with different shapes.

I checked the shapes, and I´m trying to sum shapes(2,) and ().

`print(((2*p_1) - p_0).ufl_shape)`
The answer is (2,)

`print((c*c*k*k*inner(f(w_11,w_12,w_22),v)).ufl_shape)`
The answer is ()

Hello,
I don’t understand the general shape of your formulation and the RHS of your equation, since seems that f(w) is a vector. Lighhill’s analogy returns a scalar, through the double divergence of Lighthill Tensor. What you have is basically something that resembles a specific volume velocity distribution that is, indeed, a scalar.
What you should put in the wave equation has to be scalar and you should work always with scalar fields, derived from the vector fields.

I agree with you, I´m sure I need to work with scalars and I´m sure the error is in the variational formulation, but…

That´s the basis of my f() expression as the source term

Anyway, I think my mistake lies in the variational formulation because I didn´t apply the Green theorem in the integral of this term (f()). I´m trying to…

I´m very grateful for your help! Being able to solve this is very important to me.

oook , I’ve just noticed…and is not equal to  = So, I need to recalculate my FV.

Many many thanks!!

1 Like

Sorry for my late response but yes, the mistake is that that is not a divergence but a double divergence, that makes the lightill tensor source a scalar field.
Anyway, my advise is to compute the source field inside the cfd code and import it in fenics, that solves the standard non-homogeneous wave Equation

1 Like

You don´t have to apologize! I´m very grateful for all your responses. I resolved the problem with the variational formulation, but, what I was trying to do doesn´t work. I mean: I tried to import mi flux solution (w) and split its components (w_1, w_2) to compute the Lighthill tensor, but it seems that when you export a timeseries .h5 file, you don’t have access to the components of the vector of the flux field (in this case).

So, I took your advice and I included the computation of the Lighthill Tensor into the Navier Stokes code in fenics. But, I´m lost with the postprocesing.
Here it´s my code:

``````from dolfin import *
import matplotlib.pyplot as plt
set_log_level(LogLevel.ERROR)

mesh_NS = Mesh()
with XDMFFile("mesh_NS.xdmf") as infile:

mvc = MeshValueCollection("size_t", mesh_NS, 2)
with XDMFFile("mf_NS.xdmf") as infile:
mf_NS = cpp.mesh.MeshFunctionSizet(mesh_NS, mvc)

# Set parameter values (air)

T =70   # final time
num_steps =8000  # number of time steps
dt = T / num_steps # time step size
mu = 0.000018         # dynamic viscosity
rho = 1.20          # density

# Define function spaces
V = VectorFunctionSpace(mesh_NS, 'P', 2)
Q = FunctionSpace(mesh_NS, 'P', 1)

#Boundaries
#Already defined in the mesh
#12:'fronteras superior e inferior'.
#14:'salida flujo'

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define boundary conditions
U_x=0.044

bcu_inflow = DirichletBC(V,  Constant((U_x, 0)),mf_NS, 13)
bcu_cylinder = DirichletBC(V, Constant((0, 0)),mf_NS, 11)
bcu_walls = DirichletBC(V.sub(1), Constant(0),mf_NS, 12) #symmetry conditions

bcp_outflow = DirichletBC(Q, Constant(0),mf_NS, 14)

bcu = [bcu_inflow,bcu_cylinder,bcu_walls]
bcp = [bcp_outflow]

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh_NS)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):

# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('NS_SC_200/velocitySC_200.xdmf')
xdmffile_p = XDMFFile('NS_SC_200/pressureSC_200.xdmf')
xdmffile_TL = XDMFFile('NS_SC_200/Lighthill_tensor_200.xdmf')

# Create time series (for use in other .py. Mirar https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft09_reaction_system.py)
timeseries_u = TimeSeries('NS_SC_200/velocitySC_200_series')
timeseries_p = TimeSeries('NS_SC_200/pressureSC_200_series')
timeseries_TL = TimeSeries('NS_SC_200/Lighthill_tensor_200_series')

# Create progress bar
progress = Progress('Time-stepping', num_steps)

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
t += dt

# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')

#Ligthill Tensor
def TL(u):
V = u.function_space()
mesh_NS=V.mesh
(u_1,u_2)=u.split(True)
u_11=inner(u_1,u_1)
u_12=inner(u_1,u_2)
u_22=inner(u_2,u_2)
u_11=Function(V)
u_12=Function(V)
u_22=Function(V)
degree=V.ufl_element().degree()
L=VectorFunctionSpace(mesh_NS,'P',degree)
TL=project(([rho*(u_11.dx(0)+u_12.dx(1)),rho*(u_12.dx(0)+u_22.dx(1))]),L)
return TL

# Save solution to file (XDMF/HDF5) Para visualización
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
xdmffile_TL.write(TL, t)
# Save nodal values to file (para importar en otro .py)
timeseries_u.store(u_.vector(), t)
timeseries_p.store(p_.vector(), t)
timeseries_TL.store(TL.vector(), t)

# Update previous solution
u_n.assign(u_)
p_n.assign(p_)

# Update progress bar
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
print('u max:', u_.vector().get_local().max())
print('p max:', p_.vector().get_local().max())

``````

…and it always results in the same error:

Traceback (most recent call last):
File “/mnt/—.py”, line 220, in
xdmffile_TL.write(TL, t)
AttributeError: ‘function’ object has no attribute ‘_cpp_object’

Thanks!

I think the question is…¿How to compute a function when it depends on .dx()(as an Expression do) and the results for the trial function? May be it deserves another topic…

The function TL(u) that you define wants an argument u and seems that when you are writing it you are not providing any “u”.

1 Like

I tried as follow:

``````    def TL(u_11,u_12,u_22):
"Return de Lighthill Tensor projected into same space as u"
V = VectorFunctionSpace(mesh_NS, 'P', 2)
u=TrialFunction(V)
mesh_NS = V.mesh()
u_1,u_2=split(u)
u_1=Function(V)
u_2=Function(V)
u_11=inner(u_1,u_1)
u_12=inner(u_1,u_2)
u_22=inner(u_2,u_2)
u_11=Function(V)
u_12=Function(V)
u_22=Function(V)
degree=V.ufl_element().degree()
L=VectorFunctionSpace(mesh_NS,'P',degree)
TL=project(as_vector([rho*(u_11.dx(0)+u_12.dx(1)),rho*(u_12.dx(0)+u_22.dx(1))]),L)
return TL
``````

But it results in the same error:

File “/mnt/–.py”, line 222, in
xdmffile_TL.write(TL, t)
AttributeError: ‘function’ object has no attribute ‘_cpp_object’

You are defining `TL` twice, once as

And then as the return type

You should use unique names for each of these definitions to avoid conflicts

1 Like

Yes!!  I think I did it…Now the code is running, and it’s something new!
(It will finish in several hours, I´ll tell you if it definitely works)
Here is my code:

``````from dolfin import *
import matplotlib.pyplot as plt
set_log_level(LogLevel.ERROR)

mesh_NS = Mesh()
with XDMFFile("mesh_NS.xdmf") as infile:

mvc = MeshValueCollection("size_t", mesh_NS, 2)
with XDMFFile("mf_NS.xdmf") as infile:
mf_NS = cpp.mesh.MeshFunctionSizet(mesh_NS, mvc)

T =70   # final time
num_steps =8000  # number of time steps
dt = T / num_steps # time step size
mu = 0.000018         # dynamic viscosity
rho = 1.20          # density

# Define function spaces
V = VectorFunctionSpace(mesh_NS, 'P', 2)#para la velocidad
Q = FunctionSpace(mesh_NS, 'P', 1)#para la presión

#Boundaries
#Los tenemos definidos ya en los grupos físicos de la malla:
#12:'fronteras superior e inferior'
#14:'salida flujo'

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

u_1,u_2=split(u)
u_1=Function(V)
u_2=Function(V)
u_11=Function(V)
u_12=Function(V)
u_22=Function(V)
u_11=inner(u_1,u_1)
u_12=inner(u_1,u_2)
u_22=inner(u_2,u_2)

# Define boundary conditions
U_x=0.044

bcu_inflow = DirichletBC(V,  Constant((U_x, 0)),mf_NS, 13)
bcu_cylinder = DirichletBC(V, Constant((0, 0)),mf_NS, 11)
bcu_walls = DirichletBC(V.sub(1), Constant(0),mf_NS, 12) #symmetry conditions

bcp_outflow = DirichletBC(Q, Constant(0),mf_NS, 14)

bcu = [bcu_inflow,bcu_cylinder,bcu_walls]
bcp = [bcp_outflow]

# Define expressions used in variational forms
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh_NS)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):

# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile('NS_and_TL/velocitySC_200.xdmf')
xdmffile_p = XDMFFile('NS_and_TL/pressureSC_200.xdmf')
xdmffile_TL = XDMFFile('NS_and_TL/Lighthill_tensor_200.xdmf')

# Create time series (for use in acoustic .py)
timeseries_u = TimeSeries('NS_and_TL/velocitySC_200_series')
timeseries_p = TimeSeries('NS_and_TL/pressureSC_200_series')
timeseries_TL = TimeSeries('NS_and_TL/Lighthill_tensor_200_series')

# Create progress bar
progress = Progress('Time-stepping', num_steps)

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
t += dt

# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1, 'bicgstab', 'hypre_amg')

# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2, 'bicgstab', 'hypre_amg')

# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, 'cg', 'sor')

#Resuelvo tensor de Lighthill
# Tensor_Lighthill
def TL(u_11,u_12,u_22):
"Return de Lighthill Tensor projected into same space as u"
u_1,u_2=split(u)
u_1=Function(V)
u_2=Function(V)
u_11=Function(V)
u_12=Function(V)
u_22=Function(V)
u_11=inner(u_1,u_1)
u_12=inner(u_1,u_2)
u_22=inner(u_2,u_2)
degree=V.ufl_element().degree()
L=VectorFunctionSpace(mesh_NS,'P',degree)
TL=project(as_vector([rho*(u_11.dx(0)+u_12.dx(1)),rho*(u_12.dx(0)+u_22.dx(1))]),L)
return TL

# Save solution to file (XDMF/HDF5) Para visualización
xdmffile_u.write(u_, t)
xdmffile_p.write(p_, t)
xdmffile_TL.write(TL(u_11,u_12,u_22), t)
# Save nodal values to file (para importar en otro .py)
timeseries_u.store(u_.vector(), t)
timeseries_p.store(p_.vector(), t)
timeseries_TL.store(TL(u_11,u_12,u_22).vector(), t)

# Update previous solution
u_n.assign(u_)
p_n.assign(p_)

# Update progress bar
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
print('u max:', u_.vector().get_local().max())
print('p max:', p_.vector().get_local().max())
``````

I did some changes, but the most important, appart from the definition of TL() as a function of u_11_u_12 and u_22 is that I had to define u_11,u_12,u_22 in the beginning too, and when a I tried to write the TL to the xdmffile and timeseries, I had to write TL(u_11,u_12,u_22) (I agree with you dokken, I should avoid different names for the same function).

Many thanks to all of you!

1 Like