Merge pull request #127 from cbegeman/ocn-standardize-great-circle-di…

…stance

Use MPAS-Tools for great circle distance in cosine bell tests

docs/developers_guide/ocean/api.md

@@ -98,6 +98,7 @@
 
    init.Init
    init.Init.run
+   init.cosine_bell
 
    forward.Forward
    forward.Forward.compute_cell_count

docs/users_guide/ocean/tasks/cosine_bell.md

@@ -10,7 +10,7 @@ Bell test case as first described in
 but using the variant from Sec. 3a of
 [Skamarock and Gassmann](https://doi.org/10.1175/MWR-D-10-05056.1).  A flow
 field representing solid-body rotation transports a bell-shaped perturbation
-in a tracer $psi$ once around the sphere, returning to its initial location.
+in a tracer $\psi$ once around the sphere, returning to its initial location.
 
 The task is a convergence test with time step varying proportionately to grid
 size. The result of the `analysis` step of the task is a plot like the
@@ -133,6 +133,11 @@ sphere. `psi_0` and `radius`, $R$, are given as config options and may be
 altered by the user. In the `init` step we assign `debug_tracers_1`
 to $\psi$.
 
+```{image} images/cosine_bell_init.png
+:align: center
+:width: 500 px
+```
+
 The initial velocity is equatorial:
 
 $$
@@ -172,8 +177,8 @@ time_integrator = RK4
 rk4_dt_per_km = 3.0
 ```
 
-The convergence_eval_time, run duration and output interval are the period for 
-advection to make a full rotation around the globe, 24 days:
+The `convergence_eval_time`, `run_duration` and `output_interval` are the
+period for advection to make a full rotation around the globe, 24 days:
 
 ```cfg
 # config options for spherical convergence tests

polaris/ocean/tasks/cosine_bell/analysis.py

@@ -1,7 +1,10 @@
 import numpy as np
 import xarray as xr
+from mpas_tools.transects import lon_lat_to_cartesian
+from mpas_tools.vector import Vector
 
 from polaris.ocean.convergence.spherical import SphericalConvergenceAnalysis
+from polaris.ocean.tasks.cosine_bell.init import cosine_bell
 
 
 class Analysis(SphericalConvergenceAnalysis):
@@ -73,30 +76,39 @@ def exact_solution(self, mesh_name, field_name, time, zidx=None):
                   'solution for the cosine_bell test case')
 
         config = self.config
-        latCent = config.getfloat('cosine_bell', 'lat_center')
-        lonCent = config.getfloat('cosine_bell', 'lon_center')
+        lat_center = config.getfloat('cosine_bell', 'lat_center')
+        lon_center = config.getfloat('cosine_bell', 'lon_center')
         radius = config.getfloat('cosine_bell', 'radius')
         psi0 = config.getfloat('cosine_bell', 'psi0')
         vel_pd = config.getfloat('cosine_bell', 'vel_pd')
 
         ds_mesh = xr.open_dataset(f'{mesh_name}_mesh.nc')
+        sphere_radius = ds_mesh.sphere_radius
+
         ds_init = xr.open_dataset(f'{mesh_name}_init.nc')
-        # find new location of blob center
-        # center is based on equatorial velocity
-        R = ds_mesh.sphere_radius
-        # distance in radians is
-        distRad = 2.0 * np.pi * time / (86400.0 * vel_pd)
-        newLon = lonCent + distRad
-        if newLon > 2.0 * np.pi:
-            newLon -= 2.0 * np.pi
-        tracer = np.zeros_like(ds_init.tracer1[0, :, 0].values)
-        latC = ds_init.latCell.values
-        lonC = ds_init.lonCell.values
-        # TODO replace this with great circle distance
-        temp = R * np.arccos(np.sin(latCent) * np.sin(latC) +
-                             np.cos(latCent) * np.cos(latC) * np.cos(
-            lonC - newLon))
-        mask = temp < radius
-        tracer[mask] = (psi0 / 2.0 *
-                        (1.0 + np.cos(3.1415926 * temp[mask] / radius)))
-        return xr.DataArray(data=tracer, dims=('nCells',))
+        latCell = ds_init.latCell.values
+        lonCell = ds_init.lonCell.values
+
+        # distance that the cosine bell center traveled in radians
+        # based on equatorial velocity
+        distance = 2.0 * np.pi * time / (86400.0 * vel_pd)
+
+        # new location of blob center
+        lon_new = lon_center + distance
+        if lon_new > 2.0 * np.pi:
+            lon_new -= 2.0 * np.pi
+
+        x_center, y_center, z_center = lon_lat_to_cartesian(
+            lon_new, lat_center, sphere_radius, degrees=False)
+        x_cells, y_cells, z_cells = lon_lat_to_cartesian(
+            lonCell, latCell, sphere_radius, degrees=False)
+        xyz_center = Vector(x_center, y_center, z_center)
+        xyz_cells = Vector(x_cells, y_cells, z_cells)
+        ang_dist_from_center = xyz_cells.angular_distance(xyz_center)
+        distance_from_center = ang_dist_from_center * sphere_radius
+
+        bell_value = cosine_bell(psi0, distance_from_center, radius)
+        tracer1 = np.where(distance_from_center < radius,
+                           bell_value,
+                           0.0)
+        return xr.DataArray(data=tracer1, dims=('nCells',))

polaris/ocean/tasks/cosine_bell/init.py

@@ -1,6 +1,8 @@
 import numpy as np
 import xarray as xr
 from mpas_tools.io import write_netcdf
+from mpas_tools.transects import lon_lat_to_cartesian
+from mpas_tools.vector import Vector
 
 from polaris import Step
 from polaris.ocean.vertical import init_vertical_coord
@@ -77,13 +79,15 @@ def run(self):
         ds['temperature'] = temperature_array.expand_dims(dim='Time', axis=0)
         ds['salinity'] = salinity * xr.ones_like(ds.temperature)
 
-        distance_from_center = sphere_radius * np.arccos(
-            np.sin(lat_center) * np.sin(latCell) +
-            np.cos(lat_center) * np.cos(latCell) *
-            np.cos(lonCell - lon_center))
-        bell_value = psi0 / 2.0 * (
-            1.0 + np.cos(np.pi *
-                         np.divide(distance_from_center, radius)))
+        x_center, y_center, z_center = lon_lat_to_cartesian(
+            lon_center, lat_center, sphere_radius, degrees=False)
+        x_cells, y_cells, z_cells = lon_lat_to_cartesian(
+            lonCell, latCell, sphere_radius, degrees=False)
+        xyz_center = Vector(x_center, y_center, z_center)
+        xyz_cells = Vector(x_cells, y_cells, z_cells)
+        ang_dist_from_center = xyz_cells.angular_distance(xyz_center)
+        distance_from_center = ang_dist_from_center * sphere_radius
+        bell_value = cosine_bell(psi0, distance_from_center, radius)
         debug_tracers = xr.where(distance_from_center < radius,
                                  bell_value,
                                  0.0)
@@ -104,3 +108,21 @@ def run(self):
         ds['fVertex'] = xr.zeros_like(ds_mesh.xVertex)
 
         write_netcdf(ds, 'initial_state.nc')
+
+
+def cosine_bell(max_value, ri, r):
+    """
+    Compute values according to cosine bell function
+
+    Parameters
+    ----------
+    max_value : float
+        Maximum value of the cosine bell function
+
+    ri : np.ndarray of type float
+        Distance from the center of the cosine bell in meters
+
+    r : float
+        Radius of the cosine bell in meters
+    """
+    return max_value / 2.0 * (1.0 + np.cos(np.pi * np.divide(ri, r)))