Example problem: The Helmholtz equation – scattering problems

In this document we discuss the finite-element-based solution of the Helmholtz equation, an elliptic PDE that describes time-harmonic wave propagation problems. We start by reviewing the relevant theory and then present the solution of a simple model problem – the scattering of a planar wave from a circular cylinder.

Acknowledgement: This tutorial and the associated driver codes were developed jointly with Tarak Kharrat (EnstaParisTech, Paris). |

The Helmholtz equation governs time-harmonic solutions of problems governed by the linear wave equation

where is the wavespeed. Assuming that is time-harmonic, with frequency , we write the real function as

where is complex-valued. This transforms (1) into the Helmholtz equation

where

is the wave number. Like other elliptic PDEs the Helmholtz equation admits Dirichlet, Neumann (flux) and Robin boundary conditions.

If the equation is solved in an infinite domain (e.g. in scattering problems) the solution must satisfy the so-called Sommerfeld radiation condition which in 2D has the form

Mathematically, this conditions is required to ensure the uniqueness of the solution (and hence the well-posedness of the problem). In a physical context, such as a scattering problem, the condition ensures that scattering of an incoming wave only produces outgoing not incoming waves from infinity.

The discretisation of the Helmholtz equation itself only requires a trivial modification of `oomph-lib's`

Poisson elements – we simply add the term to the residual. Since most practical applications of the Helmholtz equation involve complex-valued solutions, we provide separate storage for the real and imaginary parts of the solution – each `Node`

therefore stores two unknowns values. By default, the real and imaginary parts are stored as values 0 and 1, respectively; see the section The enumeration of the unknowns for details.

The application of Dirichlet and Neumann boundary conditions is straightforward and follows the pattern employed for the solution of the Poisson equation:

- Dirichlet conditions are imposed by pinning the relevant nodal values and setting them to the appropriate prescribed values.
- Neumann (flux) boundary conditions are imposed via
`FaceElements`

(here the`HelmholtzFluxElements`

). As usual we attach these to the faces of the "bulk" elements that are subject to the Neumann boundary conditions.

The imposition of the Sommerfeld radiation condition for problems in infinite domains is slightly more complicated. In the following discussion we will restrict ourselves to two dimensions and assume that the infinite domain is truncated at a circular artificial boundary of radius [This assumption is also made in the implementation of `oomph-lib's`

`FaceElements`

that allow the (approximate) imposition of the Sommerfeld radiation condition. The methodology can easily be modified to deal with other geometries but this has not been done yet – any volunteers?] All methods exploit the fact that the relevant solution of the Helmholtz equation can be written in polar coordinates as

where the are suitable coefficients and is the -th-order Hankel function of the first kind.

It is possible to derive approximate versions of the Sommerfeld radiation condition in which the normal derivative of the solution on the artificial boundary is related to its value and possibly its tangential derivatives. Such boundary conditions (sometimes referred to as approximate or absorbing boundary conditions – ABCs) are typically derived from asymptotic expansions of the solution at large distances from the origin and become more accurate the larger the radius of the artificial boundary is. Higher accuracy can therefore only be achieved by increasing the size of the computational domain, with an associated increase in computational cost.

`oomph-lib`

provides an implementation of the following three boundary conditions (all taken from J. J. Shirron & I. Babuska's paper "A comparison of approximate boundary conditions and infinite element methods for exterior Helmholtz problems", Computer Methods in Applied Mechanics and Engineering **164** 121-139 (1998), in which the authors compare the accuracy of these and many other approximate boundary conditions).

**Feng's first order ABC:****Feng's second order ABC:****Feng's third order ABC:**

All three boundary conditions are implemented in the class `HelmholtzAbsorbingBCElement`

. The order of the approximation can be set via the member function `HelmholtzAbsorbingBCElement::abc_order()`

. All three boundary conditions are local (relating the function to its normal derivative) and do therefore not change the sparsity of the resulting finite element equations.

Using (4), it is easy to show (see, e.g., J. Jin "The Finite Element Method in Electromagnetics (second
edition)", Wiley (2002) p. 501ff – but note that Jin assumes that the potential varies like rather than as assumed here) that the normal (radial) derivative, on the artificial boundary

is given by

where

Equation (5) again provides a condition on the normal derivative of the solution along the artificial boundary and is implemented in the `HelmholtzDtNBoundaryElement`

class. Since depends on the solution everywhere along the artificial boundary (see (6)), the application of the boundary condition (5) introduces a non-local coupling between all the degrees of freedom located on that boundary. This is handled by classifying the unknowns that affect but are not associated with the element's own nodes as external `Data`

.

To facilitate the setup of the interaction between the `HelmholtzDtNBoundaryElements`

, `oomph-lib`

provides the class `HelmholtzDtNMesh`

which provides storage for (the pointers to) the `HelmholtzDtNBoundaryElements`

that discretise the artificial boundary. The member function `HelmholtzDtNMesh::setup_gamma()`

pre-computes the values required for the imposition of equation (5). The radius of the artificial boundary and the (finite) number of (Fourier) terms used in the sum in (6) are specified as arguments to the constructor of the `HelmholtzDtNMesh`

.

**NOTE:** Since depends on the solution, it must be recomputed whenever the unknowns are updated during the Newton iteration. This is best done by adding a call to `HelmholtzDtNMesh::setup_gamma()`

to `Problem::actions_before_newton_convergence_check()`

. [If Helmholtz's equation is solved in isolation (or within a coupled, but linear problem), Newton's method will converge in one iteration. In such cases the unnecessary recomputation of after the one-and-only Newton iteration can be suppressed by setting `Problem::Problem_is_nonlinear`

to `false`

.]

We will now demonstrate the methodology for a specific example: the scattering of sound waves in an acoustic medium of density and bulk modulus . Assuming that an incoming sound wave impacts a rigid, impermeable obstacle as shown in this sketch,

we wish to find the wave field that is scattered from the body.

For this purpose we denote the time-dependent displacement of the fluid particle in the acoustic medium by and introduce a displacement potential such that

(As usual we employ asterisks to distinguish dimensional quantities from their non-dimensional equivalents, to be introduced below.) It is easy to show that satisfies the linear wave equation (1) with wave speed .

Since the surface of the scatterer is impenetrable, the normal displacement of the fluid has to vanish on and the boundary condition for the displacement potential becomes

We non-dimensionalise all lengths and displacements on some problem-dependent lengthscale (e.g. the radius of the scatterer), non-dimensionalise the potential as and scale time on the period of the oscillation, The governing equation then becomes

where the square of the wavenumber is given by

Assuming that the incoming wave (already satisfying (8)) is described by a (known) non-dimensional displacement potential of the form

we write the total potential as

where represents the displacement potential associated with the scattered field which must satisfy (2). The boundary condition (7) then becomes a Neumann (flux) boundary condition for the scattered field,

For the special case of the incoming wave being a planar wave, propagating along the x-axis, the incoming field can be written in polar coordinates as

where is the Bessel function of the first kind of order . The exact solution for the scattering of such a wave from a circular disk is given by the series

where we have chosen the disk's radius, , as the lengthscale by setting . In the above expression, denotes the Hankel function of the first kind of order and the prime denotes differentiation with respect to the function's argument.

A quantity that is of particular interest in wave propagation problems is the time-average of the power radiated by the scatterer,

In the context of an acoustic wave, the total instantaneous power, radiated over a closed boundary is

where the pressure is related to the displacement potential via

The non-dimensional time-averaged radiated power can be expressed in terms of the complex potential as

The figure below shows an animation of the displacement potential for scattering from a circular disk for a non-dimensional wavenumber of over one period of the oscillation. The simulation was performed in an annular computational domain, bounded by the outer surface the (unit) disk and an artificial outer boundary of non-dimensional radius The Sommerfeld radiation condition was imposed using the DtN mapping and the simulation was performed with spatial adaptivity (note the non-uniform refinement).

The "carpet plot" compares the exact (green) and computed (red) solutions for the displacement potential. The colours in the contour plot at the bottom of the figure provide an alternative visualisation of the magnitude of the scattered field.

As usual, we define the problem parameters in a global namespace. The main physical parameter is the (square of the) wave number, . `N_fourier`

is the number of (Fourier) terms to be used in evaluation of the series in equations (6) and (10). The remaining parameters determine how the Sommerfeld radiation condition is applied.

//===== start_of_namespace=============================================

The function `get_exact_u`

returns the exact solution for the scattering problem. We will use this function for the validation of our results.

/// Exact solution for scattered field

/// (vector returns real and impaginary parts).

{

// Switch to polar coordinates

double r;

r=sqrt(x[0]*x[0]+x[1]*x[1]);

double theta;

theta=atan2(x[1],x[0]);

// Argument for Bessel/Hankel functions

double rr=sqrt(K_squared)*r;

// Evaluate Bessel/Hankel functions

complex <double > u_ex(0.0,0.0);

// We want to compute N_fourier terms but the function

// may return fewer than that.

int n_actual=0;

&jn_a[0],&yn_a[0],

&jnp_a[0],&ynp_a[0]);

// Shout if things went wrong

#ifdef PARANOID

{

std::ostringstream error_stream;

error_stream << "CRBond_Bessel::bessjyna() only computed "

<< n_actual << " rather than " << N_fourier

<< " Bessel functions.\n";

throw OomphLibError(error_stream.str(),

OOMPH_CURRENT_FUNCTION,

OOMPH_EXCEPTION_LOCATION);

}

#endif

// Evaluate Hankel at actual radius

Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier,rr,h,hp);

// Evaluate Hankel at inner (unit) radius

Hankel_functions_for_helmholtz_problem::Hankel_first(N_fourier

,sqrt(K_squared),

h_a,hp_a);

// Compute the sum: Separate the computation of the negative

// and positive terms

{

}

{

}

// Get the real & imaginary part of the result

u[0]=real(u_ex);

u[1]=imag(u_ex);

}// end of get_exact_u

void get_exact_u(const Vector< double > &x, Vector< double > &u)

Exact solution for scattered field (vector returns real and impaginary parts).

unsigned N_fourier

Number of terms used in the computation of the exact solution.

Next we provide a function that computes the prescribed flux (normal derivative) of the solution, , evaluated on the surface of the unit disk.

/// Flux (normal derivative) on the unit disk

/// for a planar incoming wave

complex<double>& flux)

{

// Switch to polar coordinates

double r;

r=sqrt(x[0]*x[0]+x[1]*x[1]);

double theta;

theta=atan2(x[1],x[0]);

// Argument of the Bessel/Hankel fcts

double rr=sqrt(K_squared)*r;

// Compute Bessel/Hankel functions

// We want to compute N_fourier terms but the function

// may return fewer than that.

int n_actual=0;

CRBond_Bessel::bessjyna(N_fourier,rr,n_actual,&jn[0],&yn[0],

&jnp[0],&ynp[0]);

// Shout if things went wrong...

#ifdef PARANOID

{

std::ostringstream error_stream;

error_stream << "CRBond_Bessel::bessjyna() only computed "

<< n_actual << " rather than " << N_fourier

<< " Bessel functions.\n";

throw OomphLibError(error_stream.str(),

OOMPH_CURRENT_FUNCTION,

OOMPH_EXCEPTION_LOCATION);

}

#endif

// Compute the sum: Separate the computation of the negative and

// positive terms

flux=std::complex<double>(0.0,0.0);

{

}

{

}

}// end of prescribed_incoming_flux

} // end of namespace

void prescribed_incoming_flux(const Vector< double > &x, complex< double > &flux)

Flux (normal derivative) on the unit disk for a planar incoming wave.

The driver code is very straightforward. We parse the command line to determine which boundary condition to use and set the flags in the global namespace accordingly.

//==========start_of_main=================================================

/// Solve 2D Helmholtz problem for scattering of a planar wave from a

/// unit disk

//========================================================================

{

// Store command line arguments

CommandLineArgs::setup(argc,argv);

// Define case to be run

unsigned i_case=0;

CommandLineArgs::specify_command_line_flag("--case",&i_case);

// Parse command line

CommandLineArgs::parse_and_assign();

// Doc what has actually been specified on the command line

CommandLineArgs::doc_specified_flags();

// Now set flags accordingly

switch(i_case)

{

case 0:

GlobalParameters::DtN_BC=true;

break;

case 1:

GlobalParameters::DtN_BC=false;

break;

case 2:

GlobalParameters::DtN_BC=false;

break;

case 3:

GlobalParameters::DtN_BC=false;

break;

}

int main(int argc, char **argv)

Solve 2D Helmholtz problem for scattering of a planar wave from a unit disk.

Next we build the problem, either with or without enabling spatial adaptivity and define the output directory.

//Set up the problem

//------------------

#ifdef ADAPTIVE

//Set up the problem with 2D nine-node elements from the

//QHelmholtzElement family.

problem;

#else

//Set up the problem with 2D nine-node elements from the

//QHelmholtzElement family.

problem;

#endif

// Create label for output

//------------------------

DocInfo doc_info;

// Set output directory

doc_info.set_directory("RESLT");

////////////////////////////////////////////////////////////////// //////////////////////////////////...

Finally, we solve the problem and document the results.

#ifdef ADAPTIVE

// Max. number of adaptations

unsigned max_adapt=1;

// Solve the problem with Newton's method, allowing

// up to max_adapt mesh adaptations after every solve.

problem.newton_solve(max_adapt);

#else

// Solve the problem with Newton's method

problem.newton_solve();

#endif

//Output solution

problem.doc_solution(doc_info);

} //end of main

void doc_solution(DocInfo &doc_info)

Doc the solution. DocInfo object stores flags/labels for where the output gets written to.

The problem class is very similar to that employed for the adaptive solution of the 2D Poisson equation with flux boundary conditions. The only difference is that we provide two separate meshes of `FaceElements:`

one for the inner boundary where the `HelmholtzFluxElements`

apply the Neumann condition (9), and one for the outer boundary where we apply the (approximate) Sommerfeld radiation condition. As discussed in section The Dirichlet-to-Neumann mapping (DtN) , we use the function `actions_before_newton_convergence_check()`

to recompute the integral whenever the unknowns are updated during the Newton iteration.

//========= start_of_problem_class=====================================

/// Problem class to compute scattering of planar wave from unit disk

//=====================================================================

template<class ELEMENT>

{

public:

/// Constructor

/// Destructor (empty)

/// Doc the solution. DocInfo object stores flags/labels for where the

/// output gets written to

void doc_solution(DocInfo& doc_info);

/// Update the problem specs before solve (empty)

void actions_before_newton_solve(){}

/// Update the problem specs after solve (empty)

void actions_after_newton_solve(){}

/// Recompute gamma integral before checking Newton residuals

{

{

Helmholtz_outer_boundary_mesh_pt->setup_gamma();

}

}

/// Actions before adapt: Wipe the mesh of prescribed flux elements

void actions_before_adapt();

/// Actions after adapt: Rebuild the mesh of prescribed flux elements

void actions_after_adapt();

/// Create BC elements on boundary b of the Mesh pointed

/// to by bulk_mesh_pt and add them to the specified survace Mesh

void create_outer_bc_elements(

const unsigned &b, Mesh* const &bulk_mesh_pt,

Mesh* const & helmholtz_outer_boundary_mesh_pt);

/// Create Helmholtz flux elements on boundary b of the Mesh pointed

/// to by bulk_mesh_pt and add them to the specified surface Mesh

Mesh* const & helmholtz_inner_boundary_mesh_pt);

/// Delete boundary face elements and wipe the surface mesh

/// Set pointer to prescribed-flux function for all

/// elements in the surface mesh on the surface of the unit disk

void set_prescribed_incoming_flux_pt();

/// Set up boundary condition elements on outer boundary

void setup_outer_boundary();

#ifdef ADAPTIVE

/// Pointer to the "bulk" mesh

RefineableTwoDAnnularMesh<ELEMENT>* Bulk_mesh_pt;

#else

/// Pointer to the "bulk" mesh

TwoDAnnularMesh<ELEMENT>* Bulk_mesh_pt;

#endif

/// Pointer to mesh containing the DtN (or ABC) boundary

/// condition elements

HelmholtzDtNMesh<ELEMENT>* Helmholtz_outer_boundary_mesh_pt;

/// Pointer to the mesh containing

/// the Helmholtz inner boundary condition elements

}; // end of problem class

void actions_after_adapt()

Actions after adapt: Rebuild the mesh of prescribed flux elements.

void create_outer_bc_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_outer_boundary_mesh_pt)

Create BC elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to the specified...

Mesh * Helmholtz_inner_boundary_mesh_pt

Pointer to the mesh containing the Helmholtz inner boundary condition elements.

void setup_outer_boundary()

Set up boundary condition elements on outer boundary.

RefineableTwoDAnnularMesh< ELEMENT > * Bulk_mesh_pt

Pointer to the "bulk" mesh.

void actions_before_newton_solve()

Update the problem specs before solve (empty)

HelmholtzDtNMesh< ELEMENT > * Helmholtz_outer_boundary_mesh_pt

Pointer to mesh containing the DtN (or ABC) boundary condition elements.

void delete_face_elements(Mesh *const &boundary_mesh_pt)

Delete boundary face elements and wipe the surface mesh.

void actions_after_newton_solve()

Update the problem specs after solve (empty)

void actions_before_adapt()

Actions before adapt: Wipe the mesh of prescribed flux elements.

void actions_before_newton_convergence_check()

Recompute gamma integral before checking Newton residuals.

void set_prescribed_incoming_flux_pt()

Set pointer to prescribed-flux function for all elements in the surface mesh on the surface of the un...

void create_flux_elements(const unsigned &b, Mesh *const &bulk_mesh_pt, Mesh *const &helmholtz_inner_boundary_mesh_pt)

Create Helmholtz flux elements on boundary b of the Mesh pointed to by bulk_mesh_pt and add them to t...

We start by building the bulk mesh, using the refineable or non-refineable version of the `TwoDAnnularMesh`

, depending on the macro `ADAPTIVE`

. (The error tolerances for the adaptive version are chosen such that the mesh is refined non-uniformly – with the default tolerances, `oomph-lib's`

automatic mesh adaptation procedure refine the mesh uniformly.)

//=======start_of_constructor=============================================

/// Constructor for Helmholtz problem

//========================================================================

template<class ELEMENT>

{

// Setup "bulk" mesh

// # of elements in theta

unsigned n_theta=15;

// # of elements in radius

unsigned n_r=5;

// Inner radius

double a=1.0;

// Thickness of annular computational domain

double h=0.5;

// Set outer radius

// Mesh is periodic

bool periodic=true;

// Full circle

double azimuthal_fraction=1.0;

#ifdef ADAPTIVE

// Build "bulk" mesh

Bulk_mesh_pt=

new RefineableTwoDAnnularMesh<ELEMENT>(periodic,

azimuthal_fraction,n_theta,n_r,a,h);

// Create/set error estimator

Bulk_mesh_pt->spatial_error_estimator_pt()=new Z2ErrorEstimator;

// Choose error tolerances to force some uniform refinement

Bulk_mesh_pt->min_permitted_error()=0.004;

Bulk_mesh_pt->max_permitted_error()=0.01;

#else

// Build "bulk" mesh

Bulk_mesh_pt=

new TwoDAnnularMesh<ELEMENT>(periodic,

azimuthal_fraction,n_theta,n_r,a,h);

#endif

Next we create the two (empty) meshes for the `FaceElements`

,

// Pointer to mesh containing the Helmholtz outer boundary condition

// elements. Specify outer radius and number of Fourier terms to be

// used in gamma integral

Helmholtz_outer_boundary_mesh_pt =

new HelmholtzDtNMesh<ELEMENT>(a+h,GlobalParameters::N_fourier);

// Pointer to mesh containing the Helmholtz inner boundary condition

// elements. Specify outer radius

Helmholtz_inner_boundary_mesh_pt = new Mesh;

and populate them using the functions `create_flux_elements(...)`

and `create_outer_bc_elements(...)`

.

// Create prescribed-flux elements from all elements that are

// adjacent to the inner boundary , but add them to a separate mesh.

create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);

// Create outer boundary elements from all elements that are

// adjacent to the outer boundary , but add them to a separate mesh.

create_outer_bc_elements(2,Bulk_mesh_pt,Helmholtz_outer_boundary_mesh_pt);

We add the various (sub-)meshes to the problem and build the global mesh

// Add the several sub meshes to the problem

add_sub_mesh(Bulk_mesh_pt);

add_sub_mesh(Helmholtz_outer_boundary_mesh_pt);

add_sub_mesh(Helmholtz_inner_boundary_mesh_pt);

// Build the Problem's global mesh from its various sub-meshes

build_global_mesh();

Finally, we complete the build of the various elements by by passing pointers to the relevant quantities to them, and assign the equation numbers.

// Complete the build of all elements so they are fully functional

// Loop over the Helmholtz bulk elements to set up element-specific

// things that cannot be handled by constructor: Pass pointer to

// wave number squared

unsigned n_element = Bulk_mesh_pt->nelement();

for(unsigned e=0;e<n_element;e++)

{

// Upcast from GeneralisedElement to Helmholtz bulk element

ELEMENT *el_pt = dynamic_cast<ELEMENT*>(Bulk_mesh_pt->element_pt(e));

//Set the k_squared pointer

el_pt->k_squared_pt() = &GlobalParameters::K_squared;

}

// Set up elements on outer boundary

setup_outer_boundary();

// Set pointer to prescribed flux function for flux elements

set_prescribed_incoming_flux_pt();

// Setup equation numbering scheme

cout <<"Number of equations: " << assign_eqn_numbers() << std::endl;

} // end of constructor

The problem is now ready to be solved.

The mesh adaptation is driven by the error estimates for the bulk elements. The various `FaceElements`

must therefore be removed from the global mesh before the adaptation takes place. We do this by calling the function `delete_flux_elements(...)`

for the two face meshes, before rebuilding the Problem's global mesh.

//=====================start_of_actions_before_adapt======================

/// Actions before adapt: Wipe the mesh of face elements

//========================================================================

template<class ELEMENT>

{

// Kill the flux elements and wipe the boundary meshs

delete_face_elements(Helmholtz_outer_boundary_mesh_pt);

delete_face_elements(Helmholtz_inner_boundary_mesh_pt);

// Rebuild the Problem's global mesh from its various sub-meshes

rebuild_global_mesh();

}// end of actions_before_adapt

After the (bulk-)mesh has been adapted, the flux elements must be re-attached. This is done by calling the functions `create_flux_elements(...)`

and `create_outer_bc_elements`

, followed by rebuilding the Problem's global mesh. Finally, we complete the build of the `FaceElements`

by calling the functions `setup_outer_boundary()`

and `set_prescribed_incoming_flux_pt()`

.

//=====================start_of_actions_after_adapt=======================

/// Actions after adapt: Rebuild the face element meshes

//========================================================================

template<class ELEMENT>

{

// Create prescribed-flux elements and BC elements

// from all elements that are adjacent to the boundaries and add them to

// Helmholtz_boundary_meshes

create_outer_bc_elements(2,Bulk_mesh_pt,Helmholtz_outer_boundary_mesh_pt);

create_flux_elements(0,Bulk_mesh_pt,Helmholtz_inner_boundary_mesh_pt);

// Rebuild the Problem's global mesh from its various sub-meshes

rebuild_global_mesh();

// Set pointer to prescribed flux function and DtN mesh

setup_outer_boundary();

set_prescribed_incoming_flux_pt();

}// end of actions_after_adapt

The helper function `delete_face_elements()`

is used to delete all `FaceElements`

in a given surface mesh before the mesh adaptation.

//============start_of_delete_face_elements================

/// Delete face elements and wipe the boundary mesh

//==========================================================

template<class ELEMENT>

delete_face_elements(Mesh* const & boundary_mesh_pt)

{

// Loop over the surface elements

unsigned n_element = boundary_mesh_pt->nelement();

for(unsigned e=0;e<n_element;e++)

{

// Kill surface element

delete boundary_mesh_pt->element_pt(e);

}

// Wipe the mesh

boundary_mesh_pt->flush_element_and_node_storage();

} // end of delete_outer_face_elements

The functions `create_flux_elements(...)`

and `create_outer_bc_elements(...)`

create the `FaceElements`

required to apply the boundary conditions on the inner and outer boundaries of the annular computational domain. They both loop over the bulk elements that are adjacent to the appropriate mesh boundary and attach the required `FaceElements`

to their faces. The newly created `FaceElements`

are then added to the appropriate mesh.

//============start_of_create_outer_bc_elements==============================

/// Create outer BC elements on the b-th boundary of

/// the Mesh object pointed to by bulk_mesh_pt and add the elements

/// to the Mesh object pointed to by helmholtz_outer_boundary_mesh_pt.

//===========================================================================

template<class ELEMENT>

Mesh* const & helmholtz_outer_boundary_mesh_pt)

{

// Loop over the bulk elements adjacent to boundary b?

unsigned n_element = bulk_mesh_pt->nboundary_element(b);

for(unsigned e=0;e<n_element;e++)

{

// Get pointer to the bulk element that is adjacent to boundary b

ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(

bulk_mesh_pt->boundary_element_pt(b,e));

//Find the index of the face of element e along boundary b

int face_index = bulk_mesh_pt->face_index_at_boundary(b,e);

// Build the corresponding outer flux element

// Dirichlet to Neumann boundary conditon

{

HelmholtzDtNBoundaryElement<ELEMENT>* flux_element_pt = new

HelmholtzDtNBoundaryElement<ELEMENT>(bulk_elem_pt,face_index);

//Add the flux boundary element to the helmholtz_outer_boundary_mesh

helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);

}

// ABCs BC

else

{

HelmholtzAbsorbingBCElement<ELEMENT>* flux_element_pt = new

HelmholtzAbsorbingBCElement<ELEMENT>(bulk_elem_pt,face_index);

//Add the flux boundary element to the helmholtz_outer_boundary_mesh

helmholtz_outer_boundary_mesh_pt->add_element_pt(flux_element_pt);

}

} //end of loop over bulk elements adjacent to boundary b

} // end of create_outer_bc_elements

(We omit the listing of the function `create_flux_elements(...)`

because it is very similar. Feel free to inspect in the source code.)

The post-processing function `doc_solution(...)`

computes and outputs the total radiated power, and plots the computed and exact solutions (real and complex parts).

//=====================start_of_doc=======================================

/// Doc the solution: doc_info contains labels/output directory etc.

//========================================================================

template<class ELEMENT>

void ScatteringProblem<ELEMENT>::doc_solution(DocInfo&

doc_info)

{

ofstream some_file,some_file2;

char filename[100];

// Number of plot points

unsigned npts;

npts=5;

// Compute/output the radiated power

//----------------------------------

sprintf(filename,"%s/power%i.dat",doc_info.directory().c_str(),

doc_info.number());

some_file.open(filename);

// Accumulate contribution from elements

double power=0.0;

unsigned nn_element=Helmholtz_outer_boundary_mesh_pt->nelement();

for(unsigned e=0;e<nn_element;e++)

{

HelmholtzBCElementBase<ELEMENT> *el_pt =

dynamic_cast< HelmholtzBCElementBase<ELEMENT>*>(

Helmholtz_outer_boundary_mesh_pt->element_pt(e));

power += el_pt->global_power_contribution(some_file);

}

some_file.close();

oomph_info << "Total radiated power: " << power << std::endl;

// Output solution

//-----------------

sprintf(filename,"%s/soln%i.dat",doc_info.directory().c_str(),

doc_info.number());

some_file.open(filename);

Bulk_mesh_pt->output(some_file,npts);

some_file.close();

// Output exact solution

//----------------------

sprintf(filename,"%s/exact_soln%i.dat",doc_info.directory().c_str(),

doc_info.number());

some_file.open(filename);

Bulk_mesh_pt->output_fct(some_file,npts,GlobalParameters::get_exact_u);

some_file.close();

double error,norm;

sprintf(filename,"%s/error%i.dat",doc_info.directory().c_str(),

doc_info.number());

some_file.open(filename);

Bulk_mesh_pt->compute_error(some_file,GlobalParameters::get_exact_u,

error,norm);

some_file.close();

// Doc L2 error and norm of solution

oomph_info << "\nNorm of error : " << sqrt(error) << std::endl;

oomph_info << "Norm of solution: " << sqrt(norm) << std::endl << std::endl;

Finally, we create the data required to produce an animation of the actual (real) potential at 40 instants during a period of the oscillation.

// Do animation of Helmholtz solution

//-----------------------------------

unsigned nstep=40;

for (unsigned i=0;i<nstep;i++)

{

sprintf(filename,"%s/helmholtz_animation%i_frame%i.dat",

doc_info.directory().c_str(),

doc_info.number(),i);

some_file.open(filename);

sprintf(filename,"%s/exact_helmholtz_animation%i_frame%i.dat",

doc_info.directory().c_str(),

doc_info.number(),i);

some_file2.open(filename);

double phi=2.0*MathematicalConstants::Pi*double(i)/double(nstep-1);

unsigned nelem=Bulk_mesh_pt->nelement();

for (unsigned e=0;e<nelem;e++)

{

ELEMENT* el_pt=dynamic_cast<ELEMENT*>(

Bulk_mesh_pt->element_pt(e));

el_pt->output_real(some_file,phi,npts);

el_pt->output_real_fct(some_file2,phi,npts,

}

some_file.close();

some_file2.close();

}

} // end of doc

As discussed in the introduction, most practically relevant solutions of the Helmholtz equation are complex valued. Since `oomph-lib's`

solvers only deal with real (double precision) unknowns, the equations are separated into their real and imaginary parts. In the implementation of the Helmholtz elements, we store the real and imaginary parts of the solution as two separate values at each node. By default, the real and imaginary parts are accessible via `Node::value(0)`

and `Node::value(1)`

. However, to facilitate the use of the elements in multi-physics problems we avoid accessing the unknowns directly in this manner but provide the virtual function

std::complex<unsigned> HelmholtzEquations<DIM>::u_index_helmholtz()

which returns a complex number made of the two unsigneds that indicate which nodal value represents the real and imaginary parts of the solution. This function may be overloaded in combined multi-physics elements in which a Helmholtz element is combined (by multiple inheritance) with another element, using the strategy described in the Boussinesq convection tutorial.

Confirm that the (costly) re-computation of the integral in `actions_before_newton_convergence_check()`

after the first (and only) linear solve in the Newton iteration can be avoided by declaring the problem to be linear.

Explore the accuracy (and computational cost) of the various `FaceElements`

that apply the Sommmerfeld radiation condition. In particular, confirm that the accuracy of the DtN boundary condition is (nearly) independent of the radius of the artificial outer boundary, whereas the accuracy of the ABC boundary condition can only be improved by increasing the size of the computational domain.

- The source files for this tutorial are located in the directory:
demo_drivers/helmholtz/scattering/ - The driver code is:
demo_drivers/helmholtz/scattering/scattering.cc

A pdf version of this document is available.