Skip to content

Added optional parameter to OSRC operators #295

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
evantwout wants to merge 1 commit into
base: development
Choose a base branch
from
Open

Added optional parameter to OSRC operators #295

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 27 additions & 19 deletions python/bempp/api/operators/boundary/helmholtz.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -320,7 +320,7 @@ def exterior_calderon_projector(grid, wave_number, parameters=None):
return .5 * multitrace_identity(grid, parameters) - multitrace_operator(grid, wave_number, parameters)


def osrc_dtn(space, wave_number, npade=2, theta=_np.pi / 3,
def osrc_dtn(space, wave_number, npade=2, theta=_np.pi / 3, damped_wavenumber=None,
parameters=None, label="osrc_dtn"):
"""Return the OSRC approximation to the DtN map.

Expand All @@ -334,6 +334,8 @@ def osrc_dtn(space, wave_number, npade=2, theta=_np.pi / 3,
Number of Pade terms to be used in the approximation
theta : float64
Angle of the branch cut of the square root operator.
damped_wavenumber : complex
Damped wavenumber of the OSRC operator.
parameters : bempp.api.common.ParameterList
Parameters for the operator. If none given
the default global parameter object
Expand All @@ -343,19 +345,20 @@ def osrc_dtn(space, wave_number, npade=2, theta=_np.pi / 3,

"""

return _OsrcDtn(space, wave_number, npade, theta, label=label)
return _OsrcDtn(space, wave_number, npade, theta, damped_wavenumber, label=label)


class _OsrcDtn(_BoundaryOperator):

def __init__(self, space, wave_number, npade=2, theta=_np.pi / 3.,
def __init__(self, space, wave_number, npade=2, theta=_np.pi / 3., damped_wavenumber=None,
parameters=None, label="osrc_dtn"):
super(_OsrcDtn, self).__init__(space, space, space, label=label)

self._space = space
self._wave_number = wave_number
self._npade = npade
self._theta = theta
self._damped_wavenumber = damped_wavenumber
self._parameters = parameters

def _weak_form_impl(self):
Expand All @@ -376,10 +379,11 @@ def _weak_form_impl(self):

# Compute damped wavenumber

bbox = self._space.grid.bounding_box
rad = np.linalg.norm(bbox[1, :] - bbox[0, :]) / 2
dk = self._wave_number + 1j * 0.4 * \
self._wave_number ** (1.0 / 3.0) * rad ** (-2.0 / 3.0)
if self._damped_wavenumber is None:
bbox = self._space.grid.bounding_box
rad = np.linalg.norm(bbox[1, :] - bbox[0, :]) / 2
self._damped_wavenumber = self._wave_number + 1j * 0.4 * \
self._wave_number ** (1.0 / 3.0) * rad ** (-2.0 / 3.0)

# Get the Pade coefficients

Expand All @@ -390,15 +394,15 @@ def _weak_form_impl(self):
term = ZeroBoundaryOperator(space, space, space).weak_form()

for i in range(self._npade):
term -= (alpha[i] / (dk ** 2) * stiff *
term -= (alpha[i] / (self._damped_wavenumber ** 2) * stiff *
InverseSparseDiscreteBoundaryOperator(
mass - beta[i] / (dk ** 2) * stiff))
mass - beta[i] / (self._damped_wavenumber ** 2) * stiff))
result = 1j * self._wave_number * (c0 * mass + term * mass)

return result


def osrc_ntd(space, wave_number, npade=2, theta=_np.pi / 3, parameters=None, label="osrc_ntd"):
def osrc_ntd(space, wave_number, npade=2, theta=_np.pi / 3, damped_wavenumber=None, parameters=None, label="osrc_ntd"):
"""Return the OSRC approximation to the NtD map.

Parameters
Expand All @@ -411,6 +415,8 @@ def osrc_ntd(space, wave_number, npade=2, theta=_np.pi / 3, parameters=None, lab
Number of Pade terms to be used in the approximation
theta : float64
Angle of the branch cut of the square root operator.
damped_wavenumber : complex
Damped wavenumber of the OSRC operator.
parameters : bempp.api.common.ParameterList
Parameters for the operator. If none given
the default global parameter object
Expand All @@ -419,19 +425,20 @@ def osrc_ntd(space, wave_number, npade=2, theta=_np.pi / 3, parameters=None, lab
Label for the operator.

"""
return _OsrcNtd(space, wave_number, npade, theta, label=label)
return _OsrcNtd(space, wave_number, npade, theta, damped_wavenumber, label=label)


class _OsrcNtd(_BoundaryOperator):

def __init__(self, space, wave_number, npade=2, theta=_np.pi / 3.,
def __init__(self, space, wave_number, npade=2, theta=_np.pi / 3., damped_wavenumber=None,
parameters=None, label="osrc_ntd"):
super(_OsrcNtd, self).__init__(space, space, space, label=label)

self._space = space
self._wave_number = wave_number
self._npade = npade
self._theta = theta
self._damped_wavenumber = damped_wavenumber
self._parameters = parameters

def _weak_form_impl(self):
Expand All @@ -452,10 +459,11 @@ def _weak_form_impl(self):

# Compute damped wavenumber

bbox = self._space.grid.bounding_box
rad = np.linalg.norm(bbox[1, :] - bbox[0, :]) / 2
dk = self._wave_number + 1j * 0.4 * \
self._wave_number ** (1.0 / 3.0) * rad ** (-2.0 / 3.0)
if self._damped_wavenumber is None:
bbox = self._space.grid.bounding_box
rad = np.linalg.norm(bbox[1, :] - bbox[0, :]) / 2
self._damped_wavenumber = self._wave_number + 1j * 0.4 * \
self._wave_number ** (1.0 / 3.0) * rad ** (-2.0 / 3.0)

# Get the Pade coefficients

Expand All @@ -466,11 +474,11 @@ def _weak_form_impl(self):
term = ZeroBoundaryOperator(space, space, space).weak_form()

for i in range(self._npade):
term -= (alpha[i] / (dk ** 2) * stiff *
term -= (alpha[i] / (self._damped_wavenumber ** 2) * stiff *
InverseSparseDiscreteBoundaryOperator(
mass - beta[i] / (dk ** 2) * stiff))
mass - beta[i] / (self._damped_wavenumber ** 2) * stiff))
result = 1. / (1j * self._wave_number) * (
mass * InverseSparseDiscreteBoundaryOperator(mass - 1. / (dk ** 2) * stiff) * (c0 * mass + term * mass))
mass * InverseSparseDiscreteBoundaryOperator(mass - 1. / (self._damped_wavenumber ** 2) * stiff) * (c0 * mass + term * mass))

return result

Expand Down