Skip to content

Add class TestSphereRemapEdgeIntegral to sphere_hodge.py #307

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add class TestSphereRemapEdgeIntegral to sphere_hodge.py #307

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
348 changes: 347 additions & 1 deletion pycompadre/examples/sphere_hodge.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ def run(level):
#gmls_obj=pycompadre.GMLS(pycompadre.ReconstructionSpace.ScalarTaylorPolynomial,
#gmls_obj=pycompadre.GMLS(pycompadre.ReconstructionSpace.VectorTaylorPolynomial,
gmls_obj=pycompadre.GMLS(pycompadre.ReconstructionSpace.ScalarTaylorPolynomial,
pycompadre.SamplingFunctionals['ScalarFaceIntegralSample'],
pycompadre.SamplingFunctionals['CellIntegralSample'],
pycompadre.SamplingFunctionals['PointSample'],
p_order,
DIM,
Expand Down Expand Up @@ -375,6 +375,352 @@ def run(level):
rel_l2_rates.append(np.log(rel_l2_errors[i]/rel_l2_errors[i-1])/np.log(0.5))
print('cell integrals :: rel_l2_rates', rel_l2_rates)

class TestSphereRemapEdgeIntegral(KokkosTestCase):

def test_edge_integrated_remap_2d_manifold(self):

def run(level):

dataset = Dataset('../../../test_data/grids/sphere_{0}.nc'.format(str(level)), "r", format="NETCDF4")

TWO = 2 # number of vertices on an edge

DIM = 3
assert DIM==3, "Only created for 3D problem with 2D manifolds (DIM==3)"
Q_DIM_EDGE = 1 # local manifold dimension for quadrature on edges
Q_DIM_CELL = 2 # local manifold dimension for quadrature on cells

remap_type = 1
edge_normal_remap_types = (pycompadre.SamplingFunctionals['FaceNormalAverageSample'],
pycompadre.SamplingFunctionals['FaceNormalIntegralSample'])
edge_tangent_remap_types = (pycompadre.SamplingFunctionals['EdgeTangentAverageSample'],
pycompadre.SamplingFunctionals['EdgeTangentIntegralSample'])
sampling_functional_edge_normal = edge_normal_remap_types[remap_type]
sampling_functional_edge_tangent = edge_tangent_remap_types[remap_type]

radius = 6371220.00
dimensions = dataset.dimensions
variables = dataset.variables

p_order = 2
q_order_edge = 4
q_order_cell = 4
if remap_type==0:
q_order_edge = 1 # point sample can be though of as 1 pt quadrature

qp_edge = pycompadre.Quadrature(q_order_edge, Q_DIM_EDGE, "LINE")
qp_cell = pycompadre.Quadrature(q_order_cell, Q_DIM_CELL, "TRI")

# make a triangle from this midpoint, previous midpoint, and centroid
qpoints_edge = qp_edge.getSites()
qweights_edge = qp_edge.getWeights()


# fields needed from nc file
vOnE = dataset['verticesOnEdge']
xE = dataset['xEdge']
yE = dataset['yEdge']
zE = dataset['zEdge']
xV = dataset['xVertex']
yV = dataset['yVertex']
zV = dataset['zVertex']
latE = dataset['latEdge']
lonE = dataset['lonEdge']
nEdges = dataset.dimensions['nEdges'].size

# two vertex coordinates plus a normal vector and optional tangent vector (we don't use)
extra_data = np.zeros(shape=(nEdges, DIM*TWO + 2*DIM), dtype='f8')

edge_points = np.zeros(shape=(nEdges, DIM), dtype='f8')
err = np.zeros(shape=(nEdges, DIM), dtype='f8')

# need to produce data for the remap (get edge normal integrated data)
# NOTE: This must be tangent to the sphere!!!
# NOTE: lon from 0 to 2*pi can exist in a neighborhood, creating a discontinuity
f_zonal = lambda lat, lon: np.cos(lat)
f_meridional = lambda lat, lon: lat**4-(np.pi/2.0)**4
f = lambda lat, lon, zonal, meridional: f_zonal(lat,lon)*zonal + f_meridional(lat,lon)*meridional

def get_sphere_basis(lat, lon):
sin_lat = np.sin(lat)
cos_lat = np.cos(lat)
sin_lon = np.sin(lon)
cos_lon = np.cos(lon)

zonalUnitVector = np.zeros(shape=(3), dtype='f8')
zonalUnitVector[0] = - sin_lon
zonalUnitVector[1] = cos_lon
zonalUnitVector[2] = 0

meridionalUnitVector = np.zeros(shape=(3), dtype='f8')
meridionalUnitVector[0] = - sin_lat * cos_lon
meridionalUnitVector[1] = - sin_lat * sin_lon
meridionalUnitVector[2] = cos_lat

verticalUnitVector = np.zeros(shape=(3), dtype='f8')
verticalUnitVector[0] = cos_lat * cos_lon
verticalUnitVector[1] = cos_lat * sin_lon
verticalUnitVector[2] = sin_lat

return (zonalUnitVector, meridionalUnitVector, verticalUnitVector)

def atan4(y, x):
result = 0.0
PI = np.pi
if ( x == 0.0 ):
if ( y > 0.0 ) :
result = 0.5 * PI
elif ( y < 0.0 ):
result = 1.5 * PI
elif ( y == 0.0 ):
result = 0.0
elif ( y == 0 ):
if ( x > 0.0 ):
result = 0.0
elif ( x < 0.0 ):
result = PI
else:
theta = np.arctan2( abs(y), abs(x) )
if ( x > 0.0 and y > 0.0 ):
result = theta
elif ( x < 0.0 and y > 0.0 ):
result = PI - theta
elif ( x < 0.0 and y < 0.0 ):
result = PI + theta
elif ( x > 0.0 and y < 0.0 ):
result = 2.0 * PI - theta
return result

def get_lat_lon(x):
lat = np.arctan2(x[2], np.sqrt(x[0]*x[0] + x[1]*x[1]))
lon = atan4(x[1], x[0])
return (lat, lon)

# directions
T = np.zeros(shape=(nEdges,DIM,DIM), dtype='f8')
v = np.zeros(shape=(DIM,TWO), dtype='f8') # 2 is from vertices on edge
exact_in_field_edge_normal = np.zeros(shape=(nEdges), dtype='f8')
exact_in_field_edge_tangent = np.zeros(shape=(nEdges), dtype='f8')

rot_rad = np.pi/2.0
cos_r = np.cos(rot_rad)
sin_r = np.sin(rot_rad)
#R = np.array([[1, 0, 0], [0, cos_r, -sin_r], [0, sin_r, cos_r]], dtype='f8')
R = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype='f8')

computed_total_length = 0.0
ref_total_length = 0.0
# get edges and put integral quantities on them
for edge in range(nEdges):
edge_points[edge,:] = np.asarray([xE[edge], yE[edge], zE[edge]], dtype='f8')
edge_length = 0.0
ref_total_length += dataset['dvEdge'][edge]
for local_vertex in range(TWO):
vertex_num = vOnE[edge, local_vertex] - 1
extra_data[edge,local_vertex*DIM:(local_vertex+1)*DIM] = \
np.asarray([xV[vertex_num], yV[vertex_num], zV[vertex_num]], dtype='f8')
extra_data[edge,local_vertex*DIM:(local_vertex+1)*DIM] = radius * extra_data[edge,local_vertex*DIM:(local_vertex+1)*DIM] / np.linalg.norm(extra_data[edge,local_vertex*DIM:(local_vertex+1)*DIM])

v[:,0] = extra_data[edge,0:DIM].copy()
v[:,1] = extra_data[edge,DIM:2*DIM].copy() - v[:,0]

a = extra_data[edge,0:DIM]
b = extra_data[edge,DIM:2*DIM]

#theta = np.arccos(np.dot(a,b)/(np.linalg.norm(a)*np.linalg.norm(b)))
# arctan less sensitive near the poles
theta = np.arctan(np.linalg.norm(np.cross(a,b)) / np.dot(a,b))
arclength = theta * radius

#lat = dataset['latEdge'][edge]
#lon = dataset['lonEdge'][edge]
lat, lon = get_lat_lon(edge_points[edge,:])
zonalUnitVector, meridionalUnitVector, verticalUnitVector = get_sphere_basis(lat, lon)
# can add assert that latEdge, lonEdge is same lat as calculated from midpoint
#print(lat, dataset['latEdge'][edge], (lat-dataset['latEdge'][edge]) / abs(dataset['latEdge'][edge]))
#print(lon, dataset['lonEdge'][edge], (lon-dataset['lonEdge'][edge]) / abs(dataset['lonEdge'][edge]))

# is verticalUnitVector close to xyz?
#print(radius*verticalUnitVector, edge_points[edge,:], np.linalg.norm(radius*verticalUnitVector-edge_points[edge,:])/radius)
#print(np.linalg.norm(radius*verticalUnitVector-edge_points[edge,:])/radius)

T[edge,0,:] = zonalUnitVector
T[edge,1,:] = meridionalUnitVector
T[edge,2,:] = verticalUnitVector

xyz = verticalUnitVector * radius
rot_xyz = np.matmul(R, xyz)

## not constant over all quadrature, just midpoint
#angle = dataset['angleEdge'][edge] # angle at centroid eastwards

## normal at the midpoint of edge (use for orientation)
## but constant on a great circle so valid at all quadrature as well
## not sufficiently accurate near poles
#norm_vec = np.cos(angle)*zonalUnitVector + np.sin(angle)*meridionalUnitVector

# replace with cross product of endpoint vertices
norm_vec = np.cross(a,b) / np.linalg.norm(np.cross(a,b))
extra_data[edge,TWO*DIM:(TWO+1)*DIM] = norm_vec

result_normal = 0.0
result_tangent = 0.0
for i in range(len(qpoints_edge)):

unscaled_transformed_qp = v[:,0].copy()
unscaled_transformed_qp += qpoints_edge[i][0] * v[:,1]
transformed_qp_norm = np.linalg.norm(unscaled_transformed_qp)
scaled_transformed_qp = unscaled_transformed_qp.copy() * radius / transformed_qp_norm

if remap_type!=0:
lat_qp, lon_qp = get_lat_lon(scaled_transformed_qp)
zonalUnitVector_qp, meridionalUnitVector_qp, verticalUnitVector_qp = get_sphere_basis(lat_qp, lon_qp)

scaling = arclength

## u_qp = v_1 + r_qp[1]*(v_2 - v_1)
## s_qp = u_qp * radius / norm(u_qp) = radius * u_qp / norm(u_qp)
##
## so G(:,1) is \partial{s_qp}/ \partial{r_qp[1]}
## where r_qp is reference quadrature point (R^1 in 1D manifold in R^3)
##
## G(:,i) = radius * ( \partial{u_qp}/\partial{r_qp[i]} * (\sum_m u_qp[k]^2)^{-1/2}
## + u_qp * \partial{(\sum_m u_qp[k]^2)^{-1/2}}/\partial{r_qp[i]} )
##
## = radius * ( T(:,i)/norm(u_qp) + u_qp*(-1/2)*(\sum_m u_qp[k]^2)^{-3/2}
## *2*(\sum_k u_qp[k]*\partial{u_qp[k]}/\partial{r_qp[i]}) )
##
## = radius * ( T(:,i)/norm(u_qp) + u_qp*(-1/2)*(\sum_m u_qp[k]^2)^{-3/2}
## *2*(\sum_k u_qp[k]*T(k,i)) )
##
#qp_norm_sq = transformed_qp_norm**2
#G = v.copy() / transformed_qp_norm
#for k in range(DIM):
# G[:,1] += unscaled_transformed_qp*(-0.5)*pow(qp_norm_sq,-1.5) \
# *2*(unscaled_transformed_qp[k]*v[k,1]);
#scaling = radius * np.linalg.norm(G[:,1])

if remap_type==0:
# good
#result += np.dot(f(lat, lon, zonalUnitVector, meridionalUnitVector), norm_vec)

# rotate the vector 90
alt_lat, alt_lon = get_lat_lon(rot_xyz/np.linalg.norm(rot_xyz))
alt_zonalUnitVector, alt_meridionalUnitVector, alt_verticalUnitVector = get_sphere_basis(alt_lat, alt_lon)
ans_f = f(alt_lat, alt_lon, alt_zonalUnitVector, alt_meridionalUnitVector)
alt_f = np.linalg.solve(R, ans_f)
result_normal += np.dot(alt_f, norm_vec)
else:
# norm_vec is constant for all qp on a great circle
#result += np.dot(f(lat_qp, lon_qp, zonalUnitVector_qp, meridionalUnitVector_qp), norm_vec)*qweights[i]*scaling

# rotate the vector 90
xyz_qp = verticalUnitVector_qp.copy() * radius
rot_xyz_qp = np.matmul(R, xyz_qp)
d_rot = np.linalg.norm(rot_xyz_qp-rot_xyz)
assert d_rot < dataset['dcEdge'][edge], "Quadrature and midpoint too far apart"
alt_lat, alt_lon = get_lat_lon(rot_xyz_qp/np.linalg.norm(rot_xyz_qp))
alt_zonalUnitVector, alt_meridionalUnitVector, alt_verticalUnitVector = get_sphere_basis(alt_lat, alt_lon)
ans_f = f(alt_lat, alt_lon, alt_zonalUnitVector, alt_meridionalUnitVector)
alt_f = np.linalg.solve(R, ans_f)
result_normal += np.dot(alt_f, norm_vec)*qweights_edge[i]*scaling

edge_length += qweights_edge[i]*scaling

exact_in_field_edge_normal[edge] += result_normal
computed_total_length += edge_length
err[edge] = abs(edge_length-arclength)/abs(edge_length)

#print('arc norm:',np.linalg.norm(err))

if remap_type!=0: # doesn't make sense to calculate arc length when not integrating
print("LENGTH", computed_total_length, " vs ", ref_total_length)

gmls_obj_from_edge_normal=pycompadre.GMLS(pycompadre.ReconstructionSpace.VectorTaylorPolynomial,
sampling_functional_edge_normal,
pycompadre.SamplingFunctionals['PointSample'],
p_order,
DIM,
"QR",
"MANIFOLD",
"NO_CONSTRAINT",
p_order)

gmls_obj_from_edge_tangent=pycompadre.GMLS(pycompadre.ReconstructionSpace.VectorTaylorPolynomial,
sampling_functional_edge_tangent,
pycompadre.SamplingFunctionals['PointSample'],
p_order,
DIM,
"QR",
"MANIFOLD",
"NO_CONSTRAINT",
p_order)

gmls_obj_from_edge_normal.setOrderOfQuadraturePoints(q_order_edge)
gmls_obj_from_edge_normal.setDimensionOfQuadraturePoints(Q_DIM_EDGE)
gmls_obj_from_edge_normal.setQuadratureType("LINE")

gmls_obj_from_edge_normal.setWeightingParameter(2)
gmls_obj_from_edge_normal.setWeightingType("power")

gmls_helper_from_edge_normal = pycompadre.ParticleHelper(gmls_obj_from_edge_normal)
gmls_helper_from_edge_normal.generateKDTree(edge_points)
gmls_obj_from_edge_normal.addTargets(pycompadre.TargetOperation.VectorPointEvaluation)
gmls_obj_from_edge_normal.setTargetExtraData(extra_data)
gmls_obj_from_edge_normal.setSourceExtraData(extra_data)

gmls_helper_from_edge_normal.generateNeighborListsFromKNNSearchAndSet(edge_points, p_order, DIM-1, 3.8)
gmls_helper_from_edge_normal.setTangentBundle(T)
gmls_obj_from_edge_normal.generateAlphas(1, True)

# get exact solution
exact_out_field_from_edge_normal_to_edge_tangent = np.zeros(shape=(nEdges,3), dtype='f8')
for edge in range(nEdges):
lat = dataset['latEdge'][edge]
lon = dataset['lonEdge'][edge]
zonalUnitVector, meridionalUnitVector, verticalUnitVector = get_sphere_basis(lat, lon)
#exact_out_field[edge,:] = f(lat, lon, zonalUnitVector, meridionalUnitVector)

# rotate the vector 90
xyz = verticalUnitVector * radius
rot_xyz = np.matmul(R, xyz)

alt_lat, alt_lon = get_lat_lon(rot_xyz/np.linalg.norm(rot_xyz))
alt_zonalUnitVector, alt_meridionalUnitVector, alt_verticalUnitVector = get_sphere_basis(alt_lat, alt_lon)
ans_f = f(alt_lat, alt_lon, alt_zonalUnitVector, alt_meridionalUnitVector)
alt_f = np.linalg.solve(R, ans_f)
exact_out_field[edge,:] = alt_f

# get computed solution
out_field = gmls_helper.applyStencil(exact_in_field,
pycompadre.TargetOperation.VectorPointEvaluation,
sampling_functional_edge_normal)
# remove extreme lats
for edge in range(nEdges):
lat = dataset['latEdge'][edge]
if abs(lat)>np.pi/3.:
out_field[edge,:] = 0.0
exact_out_field[edge,:] = 0.0

print('computed out:',out_field)
print('exact out:',exact_out_field)
print('diff:',out_field - exact_out_field)
#with np.printoptions(threshold=np.inf):
# print('diff:',out_field - exact_out_field)
print(np.max(out_field-exact_out_field, axis=0))
# calculate error norm (l2)
rel_l2_error = np.linalg.norm(out_field-exact_out_field, axis=0)/np.sqrt(out_field.shape[0])

rel_l2_errors = []
for i in range(0,args.grids):
rel_l2_errors.append(run(i))
print('edge integrals :: rel_l2_errors', rel_l2_errors)

rel_l2_rates = []
for i in range(1,len(rel_l2_errors)):
rel_l2_rates.append(np.log(rel_l2_errors[i]/rel_l2_errors[i-1])/np.log(0.5))
print('edge integrals :: rel_l2_rates', rel_l2_rates)

if __name__ == '__main__':

import argparse
Expand Down
Loading