Multiple Scattering Beam Challenge

Below describes the coding challenge for the INI Winter school on Multiple Scattering (https://www.newton.ac.uk/event/mwsw05/)

To do the challenge, you will need to

Install Julia (https://julialang.org/downloads/). Usually the long-term support version, currently v1.6.7, is best
Get familiar with creating plots (https://docs.juliaplots.org/latest/install/)
Add the package MultipleScattering.jl (https://github.com/JuliaWaveScattering/MultipleScattering.jl)

With everything setup, please go through MultipleScattering's documentation (https://juliawavescattering.github.io/MultipleScattering.jl/dev/manual/intro/) paying specific attention to how to place particles, and how to create plots.

On Friday morning, each team will present their solution. Each team will be judged on three aspects, all of equal weighting:

Mathematical insight that was used
The presentation itself, including images and videos made to illustrate the solution
The specific goal given below. That is, to focus the wave in the region_of_interest for ω = 1.0 and not for ω = 1.5.

Incident beam

Before we get into the details, let's setup a simulation in MultipleScattering.jl

using MultipleScattering, LinearAlgebra, Statistics

# choose a background medium with 2 spatial dimensions with soundspeed c = 1
medium = Acoustic(2; c = 1.0, ρ = 1.0)

# the incident beam. Do not change this source
λ = 2pi / real(medium.c)

# The incident beam. Do not change this source
source = plane_source(medium;
    amplitude = 1.0,
    direction = [1.0, 1.0],
    position = [0.0, 0.0],
    causal = true,
    beam_width = λ
);

NOTE: This incident wave is a plane wave of the form

$$ \displaystyle u_\text{in}(x,y) = \exp(\mathrm i k (x + y) / \sqrt{2} - \mathrm i \omega t), $$

if casuality is satified: $x + y > 0$, and if we are within the beam width: $W < $beam_width$/2$, where $W = \lVert [x,y] - [1,1] (x+y)/ 2 \rVert$. Otherwise, $u_\text{in}(x,y) = 0$.

As a consequence, $u_\text{in}(x,y)$ is not exactly a solution to the wave equation. Note it is possible to define a beam which is a solution to the wave equation. To do this, you just add a line of point sources to imitate a finite transducer as shown in RegularSources.

Ploting

Let's plot this source. Noting that configuring Plots for the first time can take a while.

using Plots

plot_origin = [2.0,8.0]
plot_dimensions = [20.0, 20.0]
plot_domain = Box(plot_origin, plot_dimensions);

ω = 1.0 # angular frequency
resolution = 40
plot(source, ω;
    seriestype = :heatmap,
    region_shape = plot_domain,
    res = resolution
)
# savefig("beam-heatmap.png") #this would save the figure

The default type of plot is a contour plot. Because aesthetics is important, you should also look at the contour plot by replacing :heatmap with :contour.

Next we define a region of interest

region_of_interest = Circle([0.0,10.0], 0.5);

plot!(region_of_interest, linewidth = 2.0, linecolor = :green)

The Challenge

Your goal is to place circular particles that achieve two things

for ω = 1.0 the beam will be focused into the region_of_interest
for ω = 1.5 the region_of_interest will not have any waves inside it.

You can alter the size, position, and physical properties of the particles.

To be more specific, we will show how we will evalute the field in region_of_interest by using

Xs, is = points_in_shape(region_of_interest; res = 20)
Xs = Xs[is] # Xs are the points in a square grid and Xs[is] are the points in the region.

To evaluate the field in this region we use

ω = 1.0
result1 = run(source, Xs, ω)

ω = 2.0
result2 = run(source, Xs, ω)

# We take the mean of the absolute value of the field in this region
region_measure(result) = mean(abs.(field(result)))

# For example, for both ω = 1.0 and ω = 2.0 we, unsurprisingly, find that this region_measure is zero
region_measure(result1)
region_measure(result2)

To be specfic, your goal is to use the same set of particles such that

region_measure(result1) is as large as possilbe
region_measure(result2) is as small as possible

Placing particles

Let's place a wall made of particles and see how the region_measure changes

particle_medium = Acoustic(2; ρ=0.1, c=0.1);

# circles are the only shape you should use
shapes = [Circle([5.0,y], 0.2) for y = 1.0:0.401:9.5];
particles = [Particle(particle_medium, s) for s in shapes];

ω = 1.0 # choose angular frequncy
# for larger simulations it is better to split the step that calculates the solution from the step that plots the solution
result = run(particles, source, plot_domain, ω; res = resolution)
plot(result, ω;
    seriestype = :heatmap
)
plot!(particles)
plot!(region_of_interest, linewidth = 2.0, linecolor = :green)

We can see that this wall of particles has increased the measure in the region for both frequencies. However, your goal is to increase measure only for ω = 1.0.

ω = 1.0
result1 = run(particles, source, Xs, ω)
region_measure(result1)

ω = 1.5
result2 = run(particles, source, Xs, ω)
region_measure(result2)

Particles will only scatter the beam if their centre is within the beam width. One nice trick is to store only the scattered field using

result = run(particles, source, plot_domain, ω;
    res = resolution,
    only_scattered_waves = true
)

then you can plot the absolute value of the field with

plot(result, ω;
    seriestype = :heatmap,
    field_apply = abs
)

Fine tuning

Like any numerical simulation you should check convergence. In MultipleScattering.jl the simulation can be refined by increasing the basis_order, for example

ω = 2.0

measures = map(0:6) do basis_order
    result2 = run(particles, source, Xs, ω; basis_order = basis_order)
    region_measure(result2)
end
errors = abs.(measures[1:end-1] ./ measures[end] .- 1.0)

plot(log.(errors))

The default is basis_order = 5 which works well for particles with radius no large than 1.0. For larger particles you will need to increase this.

Things to try

You could try to position the particles to create a:

filter for a specific frequency by spacing out particles in layers
frequency sensitive reflector
waveguide, where the walls are just particles side-by-side
resonator, such as a Helmholtz resonator
random media and then optimise for particle positions
several of the above, and many more, all working together!

arturgower/beam_challenge.md