Boundary operators on physical group of mesh

[Felix Maiworm]

Dear Bempp team,

thank you for all your effort you put in this cool piece of software!

I have a question regarding the implementation of boundary operators concerning only one physical group of the mesh for my helmholtz/scattering problem. The mesh (see attachment) consists of 2 physical groups which i want to use to model sound hard (physical group 1) and sound absorbing faces (physical group 2). I tried using a blocked operator, but receive a value error ("ValueError: A must be a BoundaryOperator or BlockedBoundaryOperator").

Do you have a hint what i am doing wrong here? 

Thanks in advance,

best regards from Berlin

Felix

Ps. i am using P1 spaces because i want to implement a Burton-Miller formulation later

import bempp.api

import numpy as np

 

bempp.api.global_parameters.hmat.eps = 1e-3 # 1.5*1e-2

gmres_tolerance = 1e-2

f = 300 # frequency

c = 343 # speed of sound

rho = 1.2 # air density

k = 2 * np.pi * f / c # wave number

admittance_1 = 1e-10

admittance_2 = 0.5

grid = bempp.api.import_grid("GMSH/glatte_platte_segmente.msh")

physical_group_1 = [1]

physical_group_2 = [2]

space_general = bempp.api.function_space(grid, "P", 1)

space_1 = bempp.api.function_space(grid, "P", 1, domains=physical_group_1,element_on_segment=True, closed=True)

space_2 = bempp.api.function_space(grid, "P", 1, domains=physical_group_2,element_on_segment=True, closed=False)

identity_general = bempp.api.operators.boundary.sparse.identity(space_general, space_general, space_general)

identity_1 = bempp.api.operators.boundary.sparse.identity(space_general, space_1, space_1,)

identity_2 = bempp.api.operators.boundary.sparse.identity(space_general, space_2, space_2,)

dlp_general = bempp.api.operators.boundary.helmholtz.double_layer(space_general, space_general, space_general, k)

dlp_1 = bempp.api.operators.boundary.helmholtz.double_layer(space_general, space_1, space_1, k)

dlp_2 = bempp.api.operators.boundary.helmholtz.double_layer(space_general, space_2, space_2, k)

slp_general = bempp.api.operators.boundary.helmholtz.single_layer(space_general, space_general, space_general, k)

slp_1 = bempp.api.operators.boundary.helmholtz.single_layer(space_general, space_1, space_1, k)

slp_2 = bempp.api.operators.boundary.helmholtz.single_layer(space_general, space_2, space_2, k)

# left hands side of equations system, A

blocked = bempp.api.BlockedOperator(2, 1)

blocked[0,0] = 0.5 * identity_1 - dlp_1 - 1j * k * admittance_1 * slp_1

blocked[1,0] = 0.5 * identity_2 - dlp_2 - 1j * k * admittance_2 * slp_2

lhs = blocked.weak_form()

 

# right hands side of equations system, b    

    

theta = np.radians( 90 ) # polar angle in degrees

phi   = np.radians( 0 ) # azimuth angle in degrees

a = [ np.sin(theta)*np.cos(phi) , np.sin(theta)*np.sin(phi) , np.cos(theta) ]  # plane wave directional vector    

    

def combined_data(x, n, domain_index, result):

    result[0] = np.exp(1*1j * k * np.dot(x,a))

grid_fun_general = bempp.api.GridFunction(space_general, fun=combined_data) 

rhs = grid_fun_general

# solving

from bempp.api.linalg import gmres

surface_pressures, info = gmres(lhs, rhs,tol=gmres_tolerance)

[Timo Betcke]

Hi Felix,

thanks for your nice comments. You have the line

lhs = blocked.weak_form(). If you take it out the code should work.

Internally, the gmres function already calls weak_form and therefore expects

a boundary or blocked boundary operator.

Best wishes

Timo

[Felix Maiworm]

Hi Timo

thank you for answering.

Now with lhs = blocked i get the following error, see code below. I wonder if i use the spaces for the operators correctly. 

I want to calculate the scattering from partially absorbing surfaces. Is this possible together with using the h compression with bempp at all? The next unknown for me is how to set up the potential operators to calculate the field later.. So many questions ;)

Best regards

Felix

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-56-4ece98f228cc> in <module>()
     59 
     60 from bempp.api.linalg import gmres
---> 61 surface_pressures, info = gmres(lhs, rhs,tol=gmres_tolerance)

/usr/lib/python2.7/dist-packages/bempp/api/linalg/iterative_solvers.py in gmres(A, b, tol, restart, maxiter, use_strong_form, return_residuals, return_iteration_count)
     64         return _gmres_block_op_imp(
     65             A, b, tol, restart, maxiter, use_strong_form, return_residuals,
---> 66             return_iteration_count)
     67 
     68     raise ValueError(

/usr/lib/python2.7/dist-packages/bempp/api/linalg/iterative_solvers.py in _gmres_block_op_imp(A, b, tol, restart, maxiter, use_strong_form, return_residuals, return_iteration_count)
    199         A_op = A.weak_form()
    200         b_vec = projections_of_grid_function_list(
--> 201             b, A.dual_to_range_spaces)
    202 
    203     callback = _it_counter(return_residuals)

/usr/lib/python2.7/dist-packages/bempp/api/assembly/blocked_operator.py in projections_of_grid_function_list(grid_funs, projection_spaces)
    912     """
    913     projections = []
--> 914     for item, proj_space in zip(grid_funs, projection_spaces):
    915         projections.append(item.projections(proj_space))
    916     vec_len = 0

TypeError: zip argument #1 must support iteration

[felixhag...@gmail.com]

Hi Felix,

I could be wrong, but i think your right hand side of the equation has the wrong dimension. You set up a  [2x1]-blocked operator on the left hand side, but the right hand side is just a single GridFunction, but has to be an array

consisting of two GridFunctions.

Best regards,

Felix

[Felix Maiworm]

Hi Felix,

thanks for your hint! I changed the following part of the code for the right hand side and now gmres can solve the system:

grid_fun_1 = bempp.api.GridFunction(space_1, fun=combined_data) 

grid_fun_2 = bempp.api.GridFunction(space_2, fun=combined_data) 

rhs = grid_fun_1,grid_fun_2

Anyway the results are not what i wanted to achieve.. Using admittance_1 = admittance_2 and then comparing the results to results of an evenly distributed admittance over all surfaces, calulated by lhs = 0.5 * identity - dlp - 1j * k * admittance * slp, this is getting clear.

Do you have a hint how to formulate the left hand side for operators incorporating an additional coefficient for each face, like the 1j * k * admittance in my case? I am using the admittance to describe the sound absorbing behavior of each face of the mesh.

When i use the spaces the other way round, e.g. slp_1 = bempp.api.operators.boundary.helmholtz.single_layer(space_1, space_general, space_general, k), the results also make no sense.

Best regards,

Felix

[Felix Maiworm]

Hi Timo and Felix,

i solved the problem by constructing a Scipy LinearOperator that contains the diagonal matrix of my admittance values. This can be multiplied by the single layer boundary operator then. Anyway thanks for help!

Best regards

Felix

[mrban...@gmail.com]

Hello Timo and Felixs,

I have been trying to solve the interior acoustic problem with admittance boundary conditions, but I am still getting weird results for the cases where admittance≠0. Could you elaborate more on how you managed to solve the problem with the Scipy LinearOperator? 

Thank you so much, me and the research group I am in at UFSM, Santa Maria, Brazil would really appreciate some insights!

Best regards,

Luiz Augusto T. Ferraz Alvim

[Felix Maiworm]

Hi Luiz,

i created a admittance operator by

from scipy.sparse.linalg import LinearOperator

admittance_operator = scipy.sparse.linalg.aslinearoperator(scipy.sparse.diags( admittance_vector ,0))

then the left hand side is, in my case with the weak form operators

lhs_weak = 0.5 * i_wf - dlp_wf - 1j*k*(slp_wf * admittance_operator)

after solving i get the pressure field by

suface_pressures_admittance = bempp.api.GridFunction(space, coefficients=np.multiply(surface_pressures.coefficients,admittance[domain_indices]))

pressure_field =  u_inc(points_grid) + dlp_pot.evaluate(surface_pressures +  1j*k*slp_pot.evaluate(suface_pressures_admittance)

i hope this helps!

best regards from berlin,

Felix

[Luiz Augusto Alvim]

Hello Felix,

thank you so much for your response, with the new Bempp-cl I was able to run the admittance problem with the Multiplication Operator. I have had good results with real values of admittance, nonetheless, when using complex values I have found incorrect results. Have you encountered such problems? Once again, really appreciate the answers.

Best regards from Brazil!

Cheers!

Luiz Augusto T. Ferraz Alvim

[Timo Betcke]

Hi Luiz,

I would not yet recommend to use Bempp-cl. We are in the last steps of preparing the first release of it and are currently porting over all the tutorials. By doing so we have already found and fixed a number of bugs, but not yet merged into the master branch.

Best wishes

Timo

--
You received this message because you are subscribed to the Google Groups "bempp" group.
To unsubscribe from this group and stop receiving emails from it, send an email to bempp+un...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/bempp/5c31a707-b0e8-4a3f-ad54-6ffe7b3ccd3f%40googlegroups.com.

--

Timo Betcke

Professor of Computational Mathematics

University College London

Department of Mathematics

E-Mail: t.be...@ucl.ac.uk

Tel.: +44 (0) 20-3108-4068

[Luiz Augusto Alvim]

Thanks for your response Timo, really appreciate the effort put into bempp! I hope the bugs fixed could solve my problem.

Best regards,

Luiz Augusto T. Ferraz Alvim

Em quinta-feira, 7 de novembro de 2019 21:06:50 UTC-3, Timo Betcke escreveu:

Hi Luiz,

I would not yet recommend to use Bempp-cl. We are in the last steps of preparing the first release of it and are currently porting over all the tutorials. By doing so we have already found and fixed a number of bugs, but not yet merged into the master branch.

Best wishes

Timo

To unsubscribe from this group and stop receiving emails from it, send an email to be...@googlegroups.com.

To view this discussion on the web, visit https://groups.google.com/d/msgid/bempp/5c31a707-b0e8-4a3f-ad54-6ffe7b3ccd3f%40googlegroups.com.

--

Timo Betcke

Professor of Computational Mathematics

University College London

Department of Mathematics

E-Mail: t.b...@ucl.ac.uk

Tel.: +44 (0) 20-3108-4068