-
Notifications
You must be signed in to change notification settings - Fork 96
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
11 changed files
with
1,103 additions
and
215 deletions.
There are no files selected for viewing
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,297 @@ | ||
#!/usr/bin/env python | ||
# encoding: utf-8 | ||
r""" | ||
Two-dimensional variable-coefficient acoustics on a mapped grid | ||
=============================================================== | ||
Solve the variable-coefficient acoustics equations in 2D: | ||
.. math:: | ||
p_t + K(x,y) (u_x + v_y) & = 0 \\ | ||
u_t + p_x / \rho(x,y) & = 0 \\ | ||
v_t + p_y / \rho(x,y) & = 0. | ||
Here p is the pressure, (u,v) is the velocity, :math:`K(x,y)` is the bulk modulus, | ||
and :math:`\rho(x,y)` is the density. | ||
This example shows how to solve a problem with variable coefficients on a mapped grid. | ||
The domain contains circular inclusions with different acoustic properties. | ||
""" | ||
from __future__ import absolute_import | ||
import numpy as np | ||
from six.moves import range | ||
|
||
# Circle radius, square radius, circle center: | ||
# ((r1, r2), (x0, y0)) | ||
circles = ( ((0.15,0.204),( 0.45,0.3)), | ||
((0.15,0.204),( 0.7, 0.75))) | ||
|
||
impedance = (12.,10.) | ||
sound_speed = (0.3,1.5) | ||
|
||
def inclusion_mapping(xc, yc): | ||
"""Apply mapping from square to circle for each circle specified. | ||
Leave rest of grid cartesian. | ||
""" | ||
xp = xc+0. | ||
yp = yc+0. | ||
|
||
for circle in circles: | ||
circle_center = circle[1] | ||
r1 = circle[0][0] | ||
r2 = circle[0][1] | ||
|
||
x0, y0 = circle_center | ||
xdist = np.abs(xc-x0) | ||
ydist = np.abs(yc-y0) | ||
insquare = np.where(np.maximum(xdist,ydist)<=r2) # Square section of grid to be deformed | ||
xc0 = (xc[insquare]-x0)/r2 | ||
yc0 = (yc[insquare]-y0)/r2 | ||
|
||
xc1 = np.abs(xc0) | ||
yc1 = np.abs(yc0) | ||
d = np.maximum(xc1, yc1) | ||
d = np.maximum(d,1.e-10) | ||
d = np.minimum(d, 0.99999) | ||
d1 = d*r2/np.sqrt(2) | ||
|
||
# Pick mapping | ||
#R = np.sqrt(2) * d1 | ||
R = r1*np.ones(d1.shape) | ||
R = r1**2 / (r2*d) | ||
|
||
# Modify d1 and R ouside circle to morph back to square: | ||
ij = np.where(d>(r1/r2)) | ||
d1[ij] = r1/np.sqrt(2) + (d[ij]-r1/r2)*(r2-r1/np.sqrt(2))/(1.-r1/r2) | ||
R[ij] = r1 * ((1.-r1/r2) / (1.-d[ij]))**(r2/r1 + 0.5) | ||
|
||
xp2 = d1/d * xc1 | ||
yp2 = d1/d * yc1 | ||
center = d1 - np.sqrt(R**2 - d1**2) | ||
|
||
ij = np.where(xc1>=yc1) | ||
xp2[ij] = center[ij] + np.sqrt(R[ij]**2 - yp2[ij]**2) | ||
|
||
ij = np.where(yc1>=xc1) | ||
yp2[ij] = center[ij] + np.sqrt(R[ij]**2 - xp2[ij]**2) | ||
|
||
xp2 = np.sign(xc0) * xp2 | ||
yp2 = np.sign(yc0) * yp2 | ||
|
||
xp[insquare] = x0 + xp2 | ||
yp[insquare] = y0 + yp2 | ||
|
||
return xp, yp | ||
|
||
def compute_geometry(grid): | ||
r"""Computes | ||
a_x | ||
a_y | ||
length_ratio_left | ||
b_x | ||
b_y | ||
length_ratio_bottom | ||
cell_area | ||
""" | ||
|
||
dx, dy = grid.delta | ||
area_min = 1.e6 | ||
area_max = 0.0 | ||
|
||
x_corners, y_corners = grid.p_nodes | ||
|
||
lower_left_y, lower_left_x = y_corners[:-1,:-1], x_corners[:-1,:-1] | ||
upper_left_y, upper_left_x = y_corners[:-1,1: ], x_corners[:-1,1: ] | ||
lower_right_y, lower_right_x = y_corners[1:,:-1], x_corners[1:,:-1] | ||
upper_right_y, upper_right_x = y_corners[1:,1: ], x_corners[1:,1: ] | ||
|
||
a_x = upper_left_y - lower_left_y #upper left and lower left | ||
a_y = -(upper_left_x - lower_left_x) | ||
anorm = np.sqrt(a_x**2 + a_y**2) | ||
a_x, a_y = a_x/anorm, a_y/anorm | ||
length_ratio_left = anorm/dy | ||
|
||
b_x = -(lower_right_y - lower_left_y) #lower right and lower left | ||
b_y = lower_right_x - lower_left_x | ||
bnorm = np.sqrt(b_x**2 + b_y**2) | ||
b_x, b_y = b_x/bnorm, b_y/bnorm | ||
length_ratio_bottom = bnorm/dx | ||
|
||
area = 0*grid.c_centers[0] | ||
area += 0.5 * (lower_left_y+upper_left_y)*(upper_left_x-lower_left_x) | ||
area += 0.5 * (upper_left_y+upper_right_y)*(upper_right_x-upper_left_x) | ||
area += 0.5 * (upper_right_y+lower_right_y)*(lower_right_x-upper_right_x) | ||
area += 0.5 * (lower_right_y+lower_left_y)*(lower_left_x-lower_right_x) | ||
area = area/(dx*dy) | ||
area_min = min(area_min, np.min(area)) | ||
area_max = max(area_max, np.max(area)) | ||
|
||
return a_x, a_y, length_ratio_left, b_x, b_y, length_ratio_bottom, area | ||
|
||
def incoming_square_wave(state,dim,t,qbc,auxbc,num_ghost): | ||
""" | ||
Incoming square wave at left boundary. | ||
""" | ||
if t<0.05: | ||
s = 1.0 | ||
else: | ||
s = 0.0 | ||
|
||
for i in range(num_ghost): | ||
reflect_ind = 2*num_ghost - i - 1 | ||
alpha = auxbc[0,i,:] | ||
beta = auxbc[1,i,:] | ||
u_normal = alpha*qbc[1,reflect_ind,:] + beta*qbc[2,reflect_ind,:] | ||
u_tangential = -beta*qbc[1,reflect_ind,:] + alpha*qbc[2,reflect_ind,:] | ||
u_normal = 2.0*s - u_normal | ||
qbc[0,i,:] = qbc[0,reflect_ind,:] | ||
qbc[1,i,:] = alpha*u_normal - beta*u_tangential | ||
qbc[2,i,:] = beta*u_normal + alpha*u_tangential | ||
|
||
|
||
def setup(kernel_language='Fortran', use_petsc=False, outdir='./_output', | ||
solver_type='classic', time_integrator='SSP104', lim_type=2, | ||
num_output_times=20, disable_output=False, num_cells=200): | ||
from clawpack import riemann | ||
|
||
if use_petsc: | ||
import clawpack.petclaw as pyclaw | ||
else: | ||
from clawpack import pyclaw | ||
|
||
riemann_solver = riemann.acoustics_mapped_2D | ||
|
||
if solver_type=='classic': | ||
solver=pyclaw.ClawSolver2D(riemann_solver) | ||
solver.dimensional_split=False | ||
solver.limiters = pyclaw.limiters.tvd.MC | ||
elif solver_type=='sharpclaw': | ||
solver=pyclaw.SharpClawSolver2D(riemann_solver) | ||
solver.time_integrator=time_integrator | ||
|
||
solver.bc_lower[0]=pyclaw.BC.custom | ||
solver.user_bc_lower = incoming_square_wave | ||
solver.bc_upper[0]=pyclaw.BC.extrap | ||
solver.bc_lower[1]=pyclaw.BC.extrap | ||
solver.bc_upper[1]=pyclaw.BC.extrap | ||
|
||
solver.aux_bc_lower[0]=pyclaw.BC.extrap | ||
solver.aux_bc_upper[0]=pyclaw.BC.extrap | ||
solver.aux_bc_lower[1]=pyclaw.BC.extrap | ||
solver.aux_bc_upper[1]=pyclaw.BC.extrap | ||
|
||
x = pyclaw.Dimension(0.,1.0,num_cells,name='x') | ||
y = pyclaw.Dimension(0.,1.0,num_cells,name='y') | ||
domain = pyclaw.Domain([x,y]) | ||
|
||
num_eqn = 3 | ||
num_aux = 9 # geometry (7), impedance, sound speed | ||
state = pyclaw.State(domain,num_eqn,num_aux) | ||
state.grid.mapc2p = inclusion_mapping | ||
|
||
a_x, a_y, length_left, b_x, b_y, length_bottom, area = compute_geometry(state.grid) | ||
|
||
state.aux[0,:,:] = a_x | ||
state.aux[1,:,:] = a_y | ||
state.aux[2,:,:] = length_left | ||
state.aux[3,:,:] = b_x | ||
state.aux[4,:,:] = b_y | ||
state.aux[5,:,:] = length_bottom | ||
state.aux[6,:,:] = area | ||
state.index_capa = 6 # aux[6,:,:] holds the capacity function | ||
|
||
grid = state.grid | ||
xp, yp = grid.p_centers | ||
state.aux[7,:,:] = 1.0 # Impedance | ||
state.aux[8,:,:] = 1.0 # Sound speed | ||
|
||
for i, circle in enumerate(circles): | ||
# Set impedance and sound speed in each inclusion | ||
radius = circle[0][0] | ||
x0, y0 = circle[1] | ||
distance = np.sqrt( (xp-x0)**2 + (yp-y0)**2 ) | ||
in_circle = np.where(distance <= radius) | ||
state.aux[7][in_circle] = impedance[i] | ||
state.aux[8][in_circle] = sound_speed[i] | ||
|
||
# Set initial condition | ||
state.q[0,:,:] = 0. | ||
state.q[1,:,:] = 0. | ||
state.q[2,:,:] = 0. | ||
|
||
claw = pyclaw.Controller() | ||
claw.keep_copy = True | ||
if disable_output: | ||
claw.output_format = None | ||
claw.solution = pyclaw.Solution(state,domain) | ||
claw.solver = solver | ||
claw.outdir = outdir | ||
claw.tfinal = 0.9 | ||
claw.num_output_times = num_output_times | ||
claw.write_aux_init = True | ||
claw.setplot = setplot | ||
if use_petsc: | ||
claw.output_options = {'format':'binary'} | ||
|
||
return claw | ||
|
||
|
||
def setplot(plotdata): | ||
""" | ||
Plot solution using VisClaw. | ||
This example shows how to mark an internal boundary on a 2D plot. | ||
""" | ||
|
||
from clawpack.visclaw import colormaps | ||
|
||
plotdata.clearfigures() # clear any old figures,axes,items data | ||
plotdata.mapc2p = inclusion_mapping | ||
|
||
# Figure for pressure | ||
plotfigure = plotdata.new_plotfigure(name='Pressure', figno=0) | ||
|
||
# Set up for axes in this figure: | ||
plotaxes = plotfigure.new_plotaxes() | ||
plotaxes.title = 'Pressure' | ||
plotaxes.scaled = True # so aspect ratio is 1 | ||
plotaxes.afteraxes = plot_circles | ||
|
||
# Set up for item on these axes: | ||
plotitem = plotaxes.new_plotitem(plot_type='2d_contour') | ||
plotitem.contour_nlevels = 100 | ||
plotitem.contour_min = -2.51 | ||
plotitem.contour_max = 2.51 | ||
plotitem.plot_var = 0 | ||
|
||
# Figure for x-velocity plot | ||
plotfigure = plotdata.new_plotfigure(name='x-Velocity', figno=1) | ||
|
||
# Set up for axes in this figure: | ||
plotaxes = plotfigure.new_plotaxes() | ||
plotaxes.title = 'u' | ||
plotaxes.afteraxes = plot_circles | ||
|
||
plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') | ||
plotitem.plot_var = 1 | ||
plotitem.pcolor_cmap = 'RdBu' | ||
plotitem.add_colorbar = True | ||
plotitem.pcolor_cmin = -1.0 | ||
plotitem.pcolor_cmax= 1.0 | ||
|
||
return plotdata | ||
|
||
def plot_circles(current_data): | ||
import matplotlib.pyplot as plt | ||
ax = plt.gca() | ||
for circle in circles: | ||
x0, y0 = circle[1] | ||
radius = circle[0][0] | ||
circle1=plt.Circle((x0,y0),radius,color='k',fill=False,lw=3) | ||
ax.add_artist(circle1) | ||
|
||
|
||
if __name__=="__main__": | ||
import sys | ||
from clawpack.pyclaw.util import run_app_from_main | ||
output = run_app_from_main(setup,setplot) |
Oops, something went wrong.