Asked by Antoine Lefebvre

I am trying to implement the Boudary Layer Impedance model to approximate thermo-viscous losses at solid boundaries for pressure acoustics frequency domain problems, such as described by Bossart et al. (2003) and also by Kampinga (2010).

The idea is to use an impedance boundary condition at solid walls which is first estimated and then updated based on the current solution iteratively until convergence.

In order to calculate the value of the impedance at the boundaries at each iteration, I have to calculate the tangential wavenumber, which is obtained from the tangential laplacian of the pressure at the boundaries.

First question, I tried to calculate the laplacian on the boundaries using:

sol = pde.getSolution()
g_bnd = grad(sol, where = FunctionOnBoundary(domain))
d_bnd = div(g_bnd)

and I get :

RuntimeError: FinleyAdapterException: Finley_Assemble_gradient: Cannot calculate gradient of data because of unsuitable input data representation.

Then, assuming I found a solution to that, this is not the value I am looking for because I need the divergence of the components of the gradient parallel to the boundary only. Do I have access to tangential vectors at boundaries?

Many thanks,



Bossart, R., Joly, N., & Bruneau, M. (2003). Hybrid numerical and analytical solutions for acoustic boundary problems in thermo-viscous fluids. Journal of Sound and Vibration, 263(1), 69-84.

Kampinga, R. (2010). Viscothermal acoustics using finite elements. Ph.D. thesis. University of Twente, Enschede, The Netherlands.

Question information

English Edit question
esys-escript Edit question
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Lutz Gross (l-gross) said :

In general escript-finley does not allow 2nd order derivatives due its internal representation of
gradients. you could calcute s=Laplace(sol) by solving the PDE

s = - (- grad(sol)_i)_{,i}

which makes X= - grad(sol) and y= grad(sol)_i * n_i in the escript PDE . This then can be
interpolated to FunctionOnBoundary(domain)

The calculation of the tangential Laplace will not work with this approach.

Revision history for this message
Antoine Lefebvre (antoine-lefebvre2) said :

Thanks Lutz,

I found out a way to reformulate the problem using only the gradient and the normal vectors on the boundary.

After performing the calculations, I have a Data object Y, with values for all points on the boundaries but I only need to update the value of the points tagged with a specific value, so I do:

 for i in xrange( BC.getNumberOfDataPoints() ):
    if BC.getTagNumber(i) == 2:
 Yi = Y.getTupleForDataPoint(i)
        tmp = wavenumber*Z0*Yi[0]
        BC.setValueOfDataPoint(i, [[-tmp,-tmp],[tmp,-tmp]])


Is this the proper way to update the BC point by point? For a large model, that loop can become expansive, no?

Revision history for this message
Joel Fenwick (j-fenwick1) said :

Yes the loop would be expensive and it also uses
which we recommend avoiding _if possible_.

I haven't actually tried this but perhaps something like the following would do what you want:

# create a function for all points in Y
update=tmp*[[-1, -1], [1, -1]]

# create a function which is 1 for tag 2 and zero for everything else
mask=Data(0, BC.getFunctionSpace())

#selectively copy values into BC where the mask is >0
BC.copyWithMask(update, mask)

Revision history for this message
Antoine Lefebvre (antoine-lefebvre2) said :

Thanks Joel,
The copyWithMask method appears to work properly.

I am still having some issues with that boundary condition I am trying to implement.

I am currently doing:
g_bnd = grad(sol, where = FunctionOnBoundary(domain))

to get the gradient on the boundaries but the component of the gradient normal to the boundary is always 0. Is this the intended behavior?

In my case, I need the normal component. Should I calculate the gradient on Function(domain) instead and use a Projector to get the values on the Face elements?

Thanks again for the help.

Revision history for this message
Antoine Lefebvre (antoine-lefebvre2) said :

I am trying to calculate the laplacian as suggested by Lutz in answer #1 above but I am not getting correct results. Here is my script with a simple 2d test case. The equation is:

laplacian(p) + k**2*p = 0

with n_i*grad(p)_i = 0 on all boundaries but one where it is set to 1.

I am looking forward for your comments.


cs = 343. # speed of sound (m/sec )
rho = 1.25 # density (kg/m**3)
frequency = 100. # Hz

L = 0.3

wavenumber = 2*pi*frequency/cs

domain = Rectangle(l0=L,l1=0.02,n0=30,n1=10, order=2)

pde = LinearPDE(domain, numEquations=1)

A = pde.createCoefficient("A")
for i in range(domain.getDim()):

pde.setValue(A=A, D=wavenumber**2)

# normal acceleration BC on input plane
anBC = pde.createCoefficient("y")
anBC.setTaggedValue(1, 1)

sol = pde.getSolution()

laplacian = LinearPDE(domain, numEquations=1)
laplacian.setValue(D=1, X=-grad(sol) , y=inner(grad(sol,where=FunctionOnBoundary(domain)), domain.getNormal()))

l = laplacian.getSolution()

saveVTK("p.%04d.vtu"%frequency, p=sol, gp=grad(sol), l=l)

Revision history for this message
Launchpad Janitor (janitor) said :

This question was expired because it remained in the 'Open' state without activity for the last 15 days.

Revision history for this message
Lutz Gross (l-gross) said :

Will take a look at this problem in the next few days.

Revision history for this message
Antoine Lefebvre (antoine-lefebvre2) said :

Ok. Thanks.

Revision history for this message
Lutz Gross (l-gross) said :

What is the result for sol you are expecting?

Revision history for this message
Antoine Lefebvre (antoine-lefebvre2) said :

Do you refer to my original question or to the comment #5 above?

The equation I solved is simply laplacian(p)+k**2*p=0, so when calculating the gradient I should obtain k**2*p as the answer. This is not very useful but I should be able to make it work. This is not what I currently get.

Then, if it works, I'll try to reformulate the method that you proposed in #1 and do a system of three equations to retrieve the 3 spacial components of the vectorial gradient of grad(sol), from which I should be able to calculate the laplacian parallel to the surface and finally go on with implementing the boundary layer impedance model.

Thanks again for your support.

Revision history for this message
Lutz Gross (l-gross) said :

You need to use rich face element when calling "grad(sol,where=FunctionOnBoundary(domain)" and using the component normal to the boundary, see .
Otherwise only the tangential components are calculated properly while the normal component is set to zero - which in you case means that y=0 in the laplacian PDE.

Unfortunately, Rectangle does not use rich face elements by default. To use rich face element you need to call

    Rectangle(l0=L,l1=0.02,n0=60,n1=20, order=2, useElementsOnFace=True)

Hope this helps.

Can you help with this problem?

Provide an answer of your own, or ask Antoine Lefebvre for more information if necessary.

To post a message you must log in.