Skip to content

Commit

Permalink
Merge 00a607e into 40a7249
Browse files Browse the repository at this point in the history
  • Loading branch information
ketch committed Dec 18, 2014
2 parents 40a7249 + 00a607e commit 845bc64
Show file tree
Hide file tree
Showing 3 changed files with 323 additions and 309 deletions.
17 changes: 6 additions & 11 deletions examples/advection_2d_annulus/advection_annulus.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
"""
import numpy as np

def mapc2p_annulus(grid,c_centers):
def mapc2p_annulus(xc, yc):
"""
Specifies the mapping to curvilinear coordinates.
Expand All @@ -31,8 +31,8 @@ def mapc2p_annulus(grid,c_centers):
p_centers = []

# Polar coordinates (first coordinate = radius, second coordinate = theta)
p_centers.append(c_centers[0][:]*np.cos(c_centers[1][:]))
p_centers.append(c_centers[0][:]*np.sin(c_centers[1][:]))
p_centers.append(xc[:]*np.cos(yc[:]))
p_centers.append(xc[:]*np.sin(yc[:]))

return p_centers

Expand Down Expand Up @@ -63,13 +63,10 @@ def ghost_velocities_upper(state,dim,t,qbc,auxbc,num_ghost):
Set the velocities for the ghost cells outside the outer radius of the annulus.
In the computational domain, these are the cells at the top of the grid.
"""
from mapc2p import mapc2p

grid=state.grid
if dim == grid.dimensions[0]:
dx, dy = grid.delta
X_edges, Y_edges = grid.c_edges_with_ghost(num_ghost=2)
R_edges,Theta_edges = mapc2p(X_edges,Y_edges) # Compute edge coordinates in physical domain
R_edges,Theta_edges = grid.p_edges_with_ghost(num_ghost=2)

auxbc[:,-num_ghost:,:] = edge_velocities_and_area(R_edges[-num_ghost-1:,:],Theta_edges[-num_ghost-1:,:],dx,dy)

Expand All @@ -82,13 +79,10 @@ def ghost_velocities_lower(state,dim,t,qbc,auxbc,num_ghost):
Set the velocities for the ghost cells outside the inner radius of the annulus.
In the computational domain, these are the cells at the bottom of the grid.
"""
from mapc2p import mapc2p

grid=state.grid
if dim == grid.dimensions[0]:
dx, dy = grid.delta
X_edges, Y_edges = grid.c_edges_with_ghost(num_ghost=2)
R_edges,Theta_edges = mapc2p(X_edges,Y_edges)
R_edges,Theta_edges = grid.p_edges_with_ghost(num_ghost=2)

auxbc[:,0:num_ghost,:] = edge_velocities_and_area(R_edges[0:num_ghost+1,:],Theta_edges[0:num_ghost+1,:],dx,dy)

Expand Down Expand Up @@ -193,6 +187,7 @@ def setup(use_petsc=False,outdir='./_output',solver_type='classic'):
theta = pyclaw.Dimension('theta',theta_lower,theta_upper,m_theta)
domain = pyclaw.Domain([r,theta])
domain.grid.mapc2p = mapc2p_annulus
domain.grid.num_ghost = solver.num_ghost

num_eqn = 1
state = pyclaw.State(domain,num_eqn)
Expand Down
40 changes: 12 additions & 28 deletions examples/shallow_sphere/Rossby_wave.py
Original file line number Diff line number Diff line change
Expand Up @@ -69,35 +69,27 @@ def fortran_src_wrapper(solver,state,dt):
state.q = problem.src2(mx,my,num_ghost,xlower,ylower,dx,dy,q,aux,t,dt,Rsphere)


def mapc2p_sphere_nonvectorized(grid,mC):
def mapc2p_sphere_nonvectorized(X,Y):
"""
Maps to points on a sphere of radius Rsphere. Nonvectorized version (slow).
Takes as input: array_list made by x_coordinates, y_ccordinates in the map
space.
Inputs: x-coordinates, y-coordinates in the computational space.
Returns as output: array_list made by x_coordinates, y_ccordinates in the
physical space.
Inputs: mC = list composed by two arrays
[array ([xc1, xc2, ...]), array([yc1, yc2, ...])]
Output: pC = list composed by three arrays
[array ([xp1, xp2, ...]), array([yp1, yp2, ...]), array([zp1, zp2, ...])]
Output: list of x-, y- and z-coordinates in the physical space.
NOTE: this function is not used in the standard script.
"""

# Get number of cells in both directions
mx, my = grid.num_cells[0], grid.num_cells[1]
mx, my = X.shape

# Define new list of numpy array, pC = physical coordinates
pC = []

for i in range(mx):
for j in range(my):
xc = mC[0][i][j]
yc = mC[1][i][j]
xc = X[i][j]
yc = Y[i][j]

# Ghost cell values outside of [-3,1]x[-1,1] get mapped to other
# hemisphere:
Expand Down Expand Up @@ -150,34 +142,26 @@ def mapc2p_sphere_nonvectorized(grid,mC):
return pC


def mapc2p_sphere_vectorized(grid,mC):
def mapc2p_sphere_vectorized(X,Y):
"""
Maps to points on a sphere of radius Rsphere. Vectorized version (fast).
Takes as input: array_list made by x_coordinates, y_ccordinates in the map
space.
Inputs: x-coordinates, y-coordinates in the computational space.
Returns as output: array_list made by x_coordinates, y_ccordinates in the
physical space.
Inputs: mC = list composed by two arrays
[array ([xc1, xc2, ...]), array([yc1, yc2, ...])]
Output: pC = list composed by three arrays
[array ([xp1, xp2, ...]), array([yp1, yp2, ...]), array([zp1, zp2, ...])]
Output: list of x-, y- and z-coordinates in the physical space.
NOTE: this function is used in the standard script.
"""

# Get number of cells in both directions
mx, my = grid.num_cells[0], grid.num_cells[1]
mx, my = X.shape

# 2D array useful for the vectorization of the function
sgnz = np.ones((mx,my))

# 2D coordinates in the computational domain
xc = mC[0][:][:]
yc = mC[1][:][:]
xc = X[:][:]
yc = Y[:][:]

# Compute 3D coordinates in the physical domain
# =============================================
Expand Down
Loading

0 comments on commit 845bc64

Please sign in to comment.