import numpy as np
from pysit import RectangularDomain, CartesianMesh, PML
from pysit.util.implicit_surfaces import *
from pysit.gallery.gallery_base import GeneratedGalleryModel
__all__ = ['SonarModel', 'sonar_model', 'Submarine']
[docs]class Submarine(object):
""" Class which defines a submarine as the union of two spheres and a cyllinder. """
def __init__(self, c, scale=1.0, velocity=0.75):
self.c = np.array(c)
self.scale = scale
self.velocity = velocity
tube_length = 0.2
tube_radius = 0.05
end1_c = np.array(c)
end2_c = np.array(c)
end1_c[0] -= scale*tube_length/2.0
end2_c[0] += scale*tube_length/2.0
tube = ImplicitXAlignedCylinder(c=c, r=tube_radius*scale, length=tube_length*scale)
end1 = ImplicitSphere(c=end1_c, r=tube_radius*scale)
end2 = ImplicitSphere(c=end2_c, r=tube_radius*scale)
self.implicit_surface = ImplicitUnion(tube, end1, end2)
[docs] def __call__(self, grid, asarray=False):
return self.implicit_surface.interior(grid, asarray)
default_submarine_2D = Submarine(c=(1.25,0.75), scale=0.5, velocity=7.5)
default_submarine_3D = Submarine(c=(1.25,0.75,0.75), scale=0.5, velocity=7.5)
[docs]class SonarModel(GeneratedGalleryModel):
""" Gallery model for subsurface plus ssubmarine "sonar" model. """
model_name = "Sonar"
valid_dimensions = (2,3)
@property #read only
def dimension(self):
return self.domain.dim
supported_physics = ('acoustic',)
def __init__(self, n_pixels=(250,150),
submarine=None,
air_velocity=0.1,
water_velocity=1.0,
rock_velocity=10.0,
**kwargs):
""" Constructor for the sonar model.
Parameters
----------
n_pixels : tuple
The size, in pixels, of the model
submarine : none or PySIT Submarine
Specifies presence of a submarine in the model.
air_velocity : float
Velocity in air.
water_velocity : float
Velocity in water.
rock_velocity : float
Velocity in rock.
"""
GeneratedGalleryModel.__init__(self)
self.air_velocity = air_velocity
self.water_velocity = water_velocity
self.rock_velocity = rock_velocity
if len(n_pixels) not in [2,3]:
raise ValueError('Submarine-sonar model only works for dimensions greater than 1.')
if submarine is None:
if len(n_pixels) == 2:
submarine = default_submarine_2D
else: # len(n_pixels) == 3
submarine = default_submarine_3D
self.submarine = submarine
config_list = list()
# Configure X direction
x_lbc = kwargs['x_lbc'] if ('x_lbc' in list(kwargs.keys())) else PML(0.1, 100.0)
x_rbc = kwargs['x_rbc'] if ('x_rbc' in list(kwargs.keys())) else PML(0.1, 100.0)
xmin, xmax = 0.0, 2.5
x_config = (xmin, xmax, x_lbc, x_rbc)
config_list.append(x_config)
if len(n_pixels) == 3:
# If it is there, configure Y direction
y_lbc = kwargs['y_lbc'] if ('y_lbc' in list(kwargs.keys())) else PML(0.1, 100.0)
y_rbc = kwargs['y_rbc'] if ('y_rbc' in list(kwargs.keys())) else PML(0.1, 100.0)
ymin, ymax = 0.0, 2.5
y_config = (ymin, ymax, y_lbc, y_rbc)
config_list.append(y_config)
# Configure Z direction
z_lbc = kwargs['z_lbc'] if ('z_lbc' in list(kwargs.keys())) else PML(0.1, 100.0)
z_rbc = kwargs['z_rbc'] if ('z_rbc' in list(kwargs.keys())) else PML(0.1, 100.0)
zmin, zmax = 0.0, 1.5
z_config = (zmin, zmax, z_lbc, z_rbc)
config_list.append(z_config)
domain = RectangularDomain(*config_list)
mesh_args = [domain] + list(n_pixels)
mesh = CartesianMesh(*mesh_args)
self._mesh = mesh
self._domain = mesh.domain
# Set _initial_model and _true_model
self.rebuild_models()
[docs] def rebuild_models(self):
""" Rebuild the true and initial models based on the current configuration."""
zmin = self._domain.z.lbound
zmax = self._domain.z.rbound
xmin = self._domain.x.lbound
xmax = self._domain.x.rbound
grid = self.mesh.mesh_coords()
# the small number is added to prevent undesireable numerical effects
air_depth = (1e-8 + 2.0/15.0) * (zmax - zmin) + zmin
rock_bottom = 13.0/15.0 * (zmax - zmin) + zmin
coast_left = 3.0/25.0 * (xmax - xmin) + xmin
coast_right = 13.0/25.0 * (xmax - xmin) + xmin
max_depth = zmax
# Set up air layer
if self._domain.dim == 2:
n = (0., 1.)
p = (coast_right, air_depth)
else: # domain.dim == 3
n = (0.0, 0.0, 1.0)
p = (coast_right, coast_right, air_depth)
air_plane = ImplicitPlane(p,n)
air = air_plane
# Set up rock layer
if self._domain.dim == 2:
n = (coast_right - coast_left, -(1.0 - air_depth))
p = (coast_right, max_depth)
n2 = (0., -1.)
p2 = (0., rock_bottom)
else: # domain.dim == 3
n = (coast_right - coast_left, 0.0, -(1.0 - air_depth))
p = (coast_right, 0.0, max_depth)
n2 = (0., 0., -1.)
p2 = (0., 0., rock_bottom)
rock_plane = ImplicitPlane(p,n)
rock_plane2 = ImplicitPlane(p2,n2)
rock = ImplicitDifference(ImplicitUnion(rock_plane, rock_plane2), air_plane)
C0 = air.interior(grid, True) * self.air_velocity + \
rock.interior(grid, True) * self.rock_velocity
C0[np.where(C0 == 0.0)] = self.water_velocity
submarine = self.submarine
if submarine is not None:
sub = submarine.implicit_surface
C = air.interior(grid, True) * self.air_velocity + \
rock.interior(grid, True) * self.rock_velocity + \
sub.interior(grid, True) * submarine.velocity
C[np.where(C == 0.0)] = self.water_velocity
else:
C = C0.copy()
C.shape = self._mesh.shape()
C0.shape = self._mesh.shape()
self._true_model = C
self._initial_model = C0
[docs]def sonar_model( **kwargs ):
""" Friendly wrapper for instantiating the sonar model. """
# Setup the defaults
model_config = dict(n_pixels=(250,150),
submarine=None,
air_velocity=0.1,
water_velocity=1.0,
rock_velocity=10.0)
# Make any changes
model_config.update(kwargs)
return SonarModel(**model_config).get_setup()
if __name__ == '__main__':
from pysit import *
# Generate true wave speed
C, C0, m, d = sonar_model()
# C, C0, m, d = sonar_model(n_pixels=(250,250,150))
import matplotlib.pyplot as plt
fig = plt.figure()
fig.add_subplot(2,1,1)
vis.plot(C, m)
fig.add_subplot(2,1,2)
vis.plot(C0, m)
plt.show()