<jupyter_start><jupyter_text>Scattering of a plane wave by a homogeneous circleThis demo illustrates how to: * Compute the scattering of a plane wave by a homogeneous dielectric obstable * Solve the volume integral equation using an iterative method * Postprocess the solution to evaluate the total field * Check the accuracy by comparing to the analytical solution * Make a nice plot of the solution in the domain<jupyter_code># Import packages
import os
import sys
# FIXME: avoid this sys.path stuff
sys.path.append(os.path.join(os.path.abspath(''), '../../'))
import numpy as np
from scipy.sparse.linalg import LinearOperator, gmres
import time
from vines.geometry.geometry import shape_2d, generatedomain2d
from vines.operators.acoustic_operators import get_operator_2d, circulant_embedding
from vines.fields.plane_wave import PlaneWave_2d
from vines.operators.acoustic_matvecs import mvp_2d
from vines.precondition.circulant_acoustic import mvp_circ_2d, circulant_preconditioner
import matplotlib.pyplot as plt
# Set problem parameters and discretization resolution
ko = 20          # wavenumber
refInd = 1.2     # refractive index
shape = 'circle'    # choose shape (hex, circle, ellipse)
radius = 1       # radius of shape
n_per_lam = 10   # number of points per wavelength
angle = 0  # Incident wave angle to x-axis
d_inc = np.array([np.cos(angle), np.sin(angle)])
lambda_ext = 2 * np.pi / ko  # wavelength
# Generate grid points (r), indices of interior points (idx), pixel size (dx), 
# shape vertices (verts), interior wavelength (lambda_int)

r, idx, dx, verts, lambda_int = shape_2d(shape, refInd, lambda_ext, radius, n_per_lam)

M, N, _ = r.shape  # number of voxels in x,y directions (M, N, respectively)
# Plot the "idx" variable, which indicates which pixels lie within the circular scatterer
plt.imshow(idx)
plt.show()
# Get Toeplitz operator
a = np.sqrt(dx**2 / np.pi)  # radius of equivalent area circle
toep = get_operator_2d(dx**2, ko, r, a)
# Circulant embedding of Toeplitz matrix (required for FFT matvec)
opCirc = circulant_embedding(toep, M ,N)
# Set up the mu^2-1 matrix, call it MR
mu_sq = np.ones((M, N))
mu_sq[idx] = refInd ** 2
MR = mu_sq - 1
# Define matrix-vector product and corresponding linear operator
mvp = lambda x: mvp_2d(x, opCirc, idx, MR)

A = LinearOperator((M*N, M*N), matvec=mvp)
# Construct circulant approximation of Toeplitz matrix in x-direction
start = time.time()
circ_inv = circulant_preconditioner(toep, M, N, refInd)
end = time.time()
print('Preconditioner assembly time = ', end - start)
# Set up matrix-vector product with circulant preconditioner and establish preconditioner operator
mvp_prec = lambda x: mvp_circ_2d(x, circ_inv, M, N, idx)
prec = LinearOperator((M*N, M*N), matvec=mvp_prec)
# Assemble right-hand side (u_inc). Use a plane wave.
u_inc = PlaneWave_2d(1, ko, d_inc, r)

# Create an array that equals the incident field inside the scatterer and is zero outside
rhs = np.zeros((M, N), dtype=np.complex128)
rhs[idx] = u_inc[idx]
rhs_vec = rhs.reshape((M*N, 1), order='F')
# Perform iterative solve
it_count = 0

def iteration_counter(x):
    global it_count
    it_count += 1

start = time.time()
solp, info = gmres(A, rhs_vec, M=prec, tol=1e-5, callback=iteration_counter)
end = time.time()
print("The linear system was solved in {0} iterations".format(it_count))
print("Solve time {0} seconds".format(end-start))
from vines.operators.acoustic_matvecs import scattered_field
# Scattered field
u_sca = scattered_field(solp, opCirc, M, N, MR)
# Total field
u = u_inc + u_sca
# Plot the field
%matplotlib inline
import matplotlib
from matplotlib import pyplot as plt
matplotlib.rcParams.update({'font.size': 20})
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
fig = plt.figure(figsize=(8, 5))
ax = fig.gca()
plt.imshow(np.real(u.T), extent=[r[0, 0, 0], r[-1, 0, 0], r[0, 0, 1], r[0, -1, 1]],
           cmap=plt.cm.get_cmap('seismic'), interpolation='spline16')

circle = plt.Circle((0., 0.), radius, color='black', fill=False)
ax.add_artist(circle)

plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar()
plt.show()
# Compare with analytical solution
from analytical.analytical import penetrable_circle
u_exact = penetrable_circle(ko, ko*refInd, radius, r.T)
# Compute the error (transposes are required due to different defintions of coordinates shape)
error_l2 = np.linalg.norm(u_exact.T - u) / np.linalg.norm(u_exact.T)
print('error = ', error_l2)<jupyter_output>error =  0.020115714560849308
{'fraction_non_alphanumeric': 0.11309789897183728, 'fraction_numerical': 0.01743406347787215, 'mean_word_length': 5.0967302452316074, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '639471', 'license': 'MIT', 'n_tokens_mistral': 1598, 'n_tokens_neox': 1521, 'n_words': 494, 'path': 'demo/acoustics/demo0_circle.ipynb', 'repo_name': 'AndrewGibbs/vines', 'size': '4474'}