Repository: UXARRAY/uxarray

Issue:
Optimize Face Bounds Construction
### Proposed new feature or change:

Optimize the performance when constructing `Grid.face_bounds`

Discussion:
github-actions[bot]: 
## ASV Benchmarking

<details>
  <summary>Benchmark Comparison Results</summary>

  
Benchmarks that have improved:

| Change   | Before [1790e0e3]    | After [21239b9d]    |   Ratio | Benchmark (Parameter)                           |
|----------|----------------------|---------------------|---------|-------------------------------------------------|
| -        | 401M                 | 348M                |    0.87 | mpas_ocean.Integrate.peakmem_integrate('480km') |

Benchmarks that have stayed the same:

| Change   | Before [1790e0e3]    | After [21239b9d]    |   Ratio | Benchmark (Parameter)                                                                                                                 |
|----------|----------------------|---------------------|---------|---------------------------------------------------------------------------------------------------------------------------------------|
|          | 243M                 | 242M                |    1    | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))    |
|          | 303M                 | 303M                |    1    | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))  |
|          | 321M                 | 308M                |    0.96 | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc')) |
|          | 319M                 | 304M                |    0.95 | face_bounds.FaceBounds.peakmem_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))  |
|          | 14.0±0.03s           | 13.9±0.05s          |    0.99 | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/mpas/QU/oQU480.231010.nc'))       |
|          | 2.07±0.03s           | 2.07±0.01s          |    1    | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/scrip/outCSne8/outCSne8.nc'))     |
|          | 19.6±0.01s           | 19.4±0.07s          |    0.99 | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/geoflow-small/grid.nc'))    |
|          | 74.4±0.9ms           | 72.8±0.9ms          |    0.98 | face_bounds.FaceBounds.time_face_bounds(PosixPath('/home/runner/work/uxarray/uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc'))     |
|          | 1.65±0.01s           | 1.66±0.01s          |    1.01 | import.Imports.timeraw_import_uxarray                                                                                                 |
|          | 641±20ms             | 626±6ms             |    0.98 | mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('120km')                                                              |
|          | 39.3±0.7ms           | 39.2±0.3ms          |    1    | mpas_ocean.ConnectivityConstruction.time_face_face_connectivity('480km')                                                              |
|          | 489±20μs             | 482±6μs             |    0.99 | mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('480km')                                                                    |
|          | 1.22±0μs             | 1.24±0.01μs         |    1.02 | mpas_ocean.ConstructTreeStructures.time_ball_tree('120km')                                                                            |
|          | 296±1ns              | 314±7ns             |    1.06 | mpas_ocean.ConstructTreeStructures.time_ball_tree('480km')                                                                            |
|          | 815±2ns              | 854±10ns            |    1.05 | mpas_ocean.ConstructTreeStructures.time_kd_tree('120km')                                                                              |
|          | 285±3ns              | 305±10ns            |    1.07 | mpas_ocean.ConstructTreeStructures.time_kd_tree('480km')                                                                              |
|          | 396M                 | 394M                |    1    | mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('120km', False)                                                                       |
|          | 383M                 | 383M                |    1    | mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('120km', True)                                                                        |
|          | 355M                 | 355M                |    1    | mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('480km', False)                                                                       |
|          | 354M                 | 354M                |    1    | mpas_ocean.GeoDataFrame.peakmem_to_geodataframe('480km', True)                                                                        |
|          | 1.21±0s              | 1.21±0.01s          |    1    | mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', False)                                                                          |
|          | 58.8±0.6ms           | 59.7±1ms            |    1.02 | mpas_ocean.GeoDataFrame.time_to_geodataframe('120km', True)                                                                           |
|          | 98.0±0.5ms           | 95.3±0.6ms          |    0.97 | mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', False)                                                                          |
|          | 5.52±0.08ms          | 5.64±0.1ms          |    1.02 | mpas_ocean.GeoDataFrame.time_to_geodataframe('480km', True)                                                                           |
|          | 270M                 | 266M                |    0.98 | mpas_ocean.Gradient.peakmem_gradient('120km')                                                                                         |
|          | 244M                 | 243M                |    1    | mpas_ocean.Gradient.peakmem_gradient('480km')                                                                                         |
|          | 2.66±0.01ms          | 2.67±0.01ms         |    1    | mpas_ocean.Gradient.time_gradient('120km')                                                                                            |
|          | 290±5μs              | 290±5μs             |    1    | mpas_ocean.Gradient.time_gradient('480km')                                                                                            |
|          | 365M                 | 365M                |    1    | mpas_ocean.Integrate.peakmem_integrate('120km')                                                                                       |
|          | 177±2ms              | 181±4ms             |    1.03 | mpas_ocean.Integrate.time_integrate('120km')                                                                                          |
|          | 12.4±0.08ms          | 12.0±0.04ms         |    0.97 | mpas_ocean.Integrate.time_integrate('480km')                                                                                          |
|          | 348±3ms              | 349±1ms             |    1    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'exclude')                                                  |
|          | 357±3ms              | 352±2ms             |    0.99 | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'include')                                                  |
|          | 350±4ms              | 345±1ms             |    0.99 | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('120km', 'split')                                                    |
|          | 22.6±0.1ms           | 22.6±0.3ms          |    1    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'exclude')                                                  |
|          | 22.5±0.4ms           | 22.5±0.2ms          |    1    | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'include')                                                  |
|          | 22.5±0.1ms           | 23.7±0.3ms          |    1.05 | mpas_ocean.MatplotlibConversion.time_dataarray_to_polycollection('480km', 'split')                                                    |
|          | 55.0±0.1ms           | 55.3±0.2ms          |    1    | mpas_ocean.RemapDownsample.time_inverse_distance_weighted_remapping                                                                   |
|          | 44.4±0.07ms          | 44.4±0.02ms         |    1    | mpas_ocean.RemapDownsample.time_nearest_neighbor_remapping                                                                            |
|          | 362±2ms              | 361±2ms             |    1    | mpas_ocean.RemapUpsample.time_inverse_distance_weighted_remapping                                                                     |
|          | 265±1ms              | 265±0.9ms           |    1    | mpas_ocean.RemapUpsample.time_nearest_neighbor_remapping                                                                              |
|          | 232M                 | 234M                |    1.01 | quad_hexagon.QuadHexagon.peakmem_open_dataset                                                                                         |
|          | 232M                 | 232M                |    1    | quad_hexagon.QuadHexagon.peakmem_open_grid                                                                                            |
|          | 6.98±0.06ms          | 6.95±0.1ms          |    0.99 | quad_hexagon.QuadHexagon.time_open_dataset                                                                                            |
|          | 6.08±0.03ms          | 5.93±0.04ms         |    0.98 | quad_hexagon.QuadHexagon.time_open_grid                                                                                               |

Benchmarks that have got worse:

| Change   | Before [1790e0e3]    | After [21239b9d]    |   Ratio | Benchmark (Parameter)                                              |
|----------|----------------------|---------------------|---------|--------------------------------------------------------------------|
| +        | 1.77±0.03ms          | 2.06±0.1ms          |    1.16 | mpas_ocean.ConnectivityConstruction.time_n_nodes_per_face('120km') |

</details>

philipc2: @hongyuchen1030 

For the first step in optimizing, I've looked at decorating some of the commonly called functions with `numba`, and used a `numba` implementation of `isclose`. Looks to be about a 15-20% performance improvement to begin with.

Next step will be to eliminate any unnecessary coordinate conversion calls, such as `_xyz_to_lonlat_rad`

This might look like us providing both the `latlon_rad` and `cartesian` points in the function call to avoid needing to convert. Something like:

```Python
 def _point_within_gca(pt_latlon,
                       pt_xyz,
                       gca_a_latlon,
                       gca_a_xyz,
                       gca_b_latlon,
                       gca_b_xyz,
                       is_directed=False):
```

Hopefully once we've moved these conversions out of each function call, we can see if the remaining logic can be decorated with `numba`.
philipc2: @hongyuchen1030 


https://github.com/UXARRAY/uxarray/blob/1790e0e30abf6b8869c35747d6393af03bccac88/uxarray/grid/arcs.py#L77

I had a question about this and other parts of the code where the Cartesian coordinates are converted to latlon in radians. Are we able to use the `node_lon` and `node_lat` that are part of the `Grid` (in degrees) instead of needing to run the conversion? Additionally, is radians here just a preference or can we use degrees? 
hongyuchen1030: > Next step will be to eliminate any unnecessary coordinate conversion calls, such as `_xyz_to_lonlat_rad`
>
> This might look like us providing both the `latlon_rad` and `cartesian` points in the function call to avoid needing to convert. Something like:
> 
> ```python
>  def _point_within_gca(pt_latlon,
>                        pt_xyz,
>                        gca_a_latlon,
>                        gca_a_xyz,
>                        gca_b_latlon,
>                        gca_b_xyz,
>                        is_directed=False):
> ```

I have also been considering this design today. It will be highly beneficial to call the coordinate conversion only once, both for performance and accuracy.

We can extract all functions that take node coordinates as input and ensure that both Cartesian and spherical coordinates are required as parameters.

> I had a question about this and other parts of the code where the Cartesian coordinates are converted to latlon in radians. Are we able to use the `node_lon` and `node_lat` that are part of the Grid (in degrees) instead of needing to run the conversion? Additionally, is radians here just a preference or can we use degrees?

That is another reason why I have been advocating for better documentation and conventions for the coordinate system.

Currently, we have the following coordinate systems:
1. Spherical coordinates in degrees for users.
2. Spherical coordinates in radians because the trigonometry functions in `numpy` use radians  (`np.sin`, `np.cos`, etc.).
3. Cartesian coordinate system with normalization.
4. Cartesian coordinate system without normalization (we still haven't confirmed if the node data in MPAS or other unnormalized coordinates are equal to the `normalized_coordinates * earth_radius` yet).

Before we attempt to rewrite everything, we must ensure which coordinate systems we are using across the entire project.

philipc2: > I have also been considering this design today. It will be highly beneficial to call the coordinate conversion only once, both for performance and accuracy.

Agreed, was thinking of this as well.

> That is another reason why I have been advocating for better documentation and conventions for the coordinate system.

We have some documentation here in our user guide. https://uxarray.readthedocs.io/en/latest/user-guide/representation.html

Any suggestions on how we could improve this? 

> Currently, we have the following coordinate systems:

Wondering if it would make sense to include a `node_lon_rad` and `node_lat_rad` extension to how we represent our coordinates. I've seen some other grids store the coordinates both in degrees and radians, so this might be a worth extension (and would help with some of the conversions we are doing).

Hopefully the normalization inconsistencies will be resolved with #878 




hongyuchen1030: > We have some documentation here in our user guide. https://uxarray.readthedocs.io/en/latest/user-guide/representation.html
> 
> Any suggestions on how we could improve this?
![image](https://github.com/user-attachments/assets/6b906b4c-4568-4c81-968f-8f8b465c70e3)

1. I understand we already specify the convention underneath; however, the `node_lat` has the description specifying the units while others do not. It would be great to keep them uniform. Additionally, placing the note before the detailed description would enhance clarity.

2. Just a note: I think we didn't use the convention of [-180, 180] until we updated the coordinates conversion. The old conversion function used the convention of [0, 360). All my algorithms are based on [0, 2π) as well. I highly recommend using [0, 2π) or [0, 360) since it provides an easier way to handle the antimeridian faces. This also necessitates a thorough re-check after any major functions are updated.

3. It might also be beneficial to specify that we have a unique definition for the antimeridian points and the pole points.

> Wondering if it would make sense to include a `node_lon_rad` and `node_lat_rad` extension to how we represent our coordinates. I've seen some other grids store the coordinates both in degrees and radians, so this might be a worth extension (and would help with some of the conversions we are doing).

I think what we've discovered so far is that the degrees are solely for the users, while all calculations use radians since the `numpy` trigonometry functions operate in radians. I'm always open to utilizing storage space to reduce time complexity.

philipc2: > I understand we already specify the convention underneath; however, the node_lat has the description specifying the units while others do not. It would be great to keep them uniform. Additionally, placing the note before the detailed description would enhance clarity.

That's a good catch! 

> Just a note: I think we didn't use the convention of [-180, 180] until we updated the coordinates conversion. The old conversion function used the convention of [0, 360). All my algorithms are based on [0, 2π) as well. I highly recommend using [0, 2π) or [0, 360) since it provides an easier way to handle the antimeridian faces. This also necessitates a thorough re-check after any major functions are updated.

We use -180 to 180 since it makes certain visualizations tasks (i.e. projections and the use of the `antimeridian` package) simpler. Most of our datasets are also already in the -180 to 180 range.

> I think what we've discovered so far is that the degrees are solely for the users, while all calculations use radians since the numpy trigonometry functions operate in radians. I'm always open to utilizing storage space to reduce time complexity.

This isn't fully the case. Many of our internal functions (i.e. BallTree) use degrees. 

 

> numpy trigonometry functions operate in radians

One option could be to use Scipy's degree implementations, but I'm not sure how much this implementation differs from Numpy's  

https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sindg.html#scipy.special.sindg
philipc2: @hongyuchen1030 

Would you mind explaining this code snippet for me?

https://github.com/UXARRAY/uxarray/blob/1790e0e30abf6b8869c35747d6393af03bccac88/uxarray/grid/arcs.py#L288-L292


hongyuchen1030: > @hongyuchen1030
> 
> Would you mind explaining this code snippet for me?
> 
> https://github.com/UXARRAY/uxarray/blob/1790e0e30abf6b8869c35747d6393af03bccac88/uxarray/grid/arcs.py#L288-L292

Sure, I wonder if you want to hop into a quick call?
hongyuchen1030: > We use -180 to 180 since it makes certain visualization tasks (i.e., projections and the use of the antimeridian package) simpler. Most of our datasets are also already in the -180 to 180 range.

I understand the benefits of using the -180 to 180 range for visualization tasks and consistency with existing datasets. My algorithms do not expose the longitude range to users, nor do they modify existing longitudes, so this shouldn't pose a significant issue. However, I wanted to bring this up because if we require functions to take `node_lat` and `node_lon` from the `grid`, we may need to quickly check if longitude conversion is necessary. This note could be helpful for troubleshooting future bugs.

> This isn't fully the case. Many of our internal functions (i.e., BallTree) use degrees.

Is this because we directly use the `grid.node_lat` and `grid.node_lon` without calling the `numpy` trigonometric functions?

> One option could be to use Scipy's degree implementations, but I'm not sure how much this implementation differs from Numpy's.

I compared the source code for Scipy's `sindg` and Numpy's `sin`. Here’s what I found:

1. Numpy’s `sin` is a wrapper for the compiled `c_sin`, with the source code for `s_sin` available [here](https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/x86_64/fpu/s_sinf.S;h=1529fb1f2452cefa1fdf482f8fa9165ffef4c197;hb=HEAD).
2. Scipy’s `sindg` was implemented recently in this [`.h` file](https://github.com/scipy/scipy/blob/94532e74b902b569bfad504866cb53720c5f4f31/scipy/special/xsf/cephes/sindg.h#L5), and it still converts some variables into radians.

I haven't delved deeply into the internal workings of both functions, but we could perform a simple benchmark to compare their running times and compatibility if needed.

philipc2: @hongyuchen1030 

After some further optimizations, this is the point that I've gotten to. Looks like the `pyfma` implementations still contribute to a pretty significant slowdown, taking up the majority of the execution time. 

## FMA Disabled 
~5.73s Execution Time 
![image](https://github.com/user-attachments/assets/5aab3d54-bd78-4d93-8cd5-3f84ccf7b194)


## FMA Enabled
~14.1s Execution Time 

![image](https://github.com/user-attachments/assets/0ced72c4-b9b0-4efb-80a5-ae878d632daa)

philipc2: A few more refinements and with `fma` disabled, able to get down for 3.74s

![image](https://github.com/user-attachments/assets/f1b059b5-feed-4e58-98a4-8eb53b41438e)

philipc2: ![image](https://github.com/user-attachments/assets/f21d1659-eb46-4528-a2fa-64e538a84321)

From the benchmarks run above, looks like almost a 8-10x improvement in performance across the board with these changes. 

I was also able to compute the bounds on a 30km MPAS atmosphere grid (655k Faces, 1,310,000 nodes) in about 40 seconds, which is a significant improvement and wouldn't of been possible without these optimizations. 
hongyuchen1030: > @hongyuchen1030
> 
> After some further optimizations, this is the point that I've gotten to. Looks like the `pyfma` implementations still contribute to a pretty significant slowdown, taking up the majority of the execution time.
> 
> ## FMA Disabled
> ~5.73s Execution Time ![image](https://private-user-images.githubusercontent.com/67855069/356027928-5aab3d54-bd78-4d93-8cd5-3f84ccf7b194.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjMwODk4NzQsIm5iZiI6MTcyMzA4OTU3NCwicGF0aCI6Ii82Nzg1NTA2OS8zNTYwMjc5MjgtNWFhYjNkNTQtYmQ3OC00ZDkzLThjZDUtM2Y4NGNjZjdiMTk0LnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MDglMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODA4VDAzNTkzNFomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWI0YmFkNTZmNmZlZGMzMmZjMTJmMWQxNzIyOWZjZjY2ZTg4NTg0N2UwMjZlM2ZkYmFhNmFhZDM1NTczZWM4NGMmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.2GF-LrJ7CynplPu-FyJugUijaGmYFKN7Lyfer5B_roM)
> 
> ## FMA Enabled
> ~14.1s Execution Time
> 
> ![image](https://private-user-images.githubusercontent.com/67855069/356030502-0ced72c4-b9b0-4efb-80a5-ae878d632daa.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MjMwODk4NzQsIm5iZiI6MTcyMzA4OTU3NCwicGF0aCI6Ii82Nzg1NTA2OS8zNTYwMzA1MDItMGNlZDcyYzQtYjliMC00ZWZiLTgwYTUtYWU4NzhkNjMyZGFhLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA4MDglMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwODA4VDAzNTkzNFomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWVhMWZiY2E0NzI3YjZkYmExYTcwMDE4ZmQ3YjZlNDk0M2NhYTZlMDVlOGMwYWE0MzIwYTlkNGYzMGNiNWJlYWUmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.oqe6rm5B4S1PUnEOiY60UBV_ip9dWkAA9_H_F2RhZ1U)
> 
> Do you think anything from here could possible act as an alternative to using `pyfma`?
> 
> https://docs.scipy.org/doc/scipy/reference/linalg.blas.html

I've been giving this case a lot of thought, considering my research findings and the intricacies of floating point operations.

Given the technical details involved with floating point operations, I won't delve into too much detail to avoid overwhelming you. Here's a summary of my thoughts:

1. **Maintain Current Approach**: We can keep the current implementation as it is, but set `is_FMA_enable` to `False` globally and use the basic `np.cross`. We'll document this choice, informing users that they can opt for a more accurate method by enabling FMA.
2. **Adjust Error Tolerance**: The immediate task would be to reset the error tolerance for intersection point snapping. This warrants a separate issue to be opened specifically for this. And since we are not allowed to use any compiled files in our project code base, we cannot do much about it now.

### Rationale:

1. **Algorithm Accuracy**: While CDO uses Kahan's Algorithm with FMA, my research indicates there's a more accurate method available. This method works with and without FMA, though the non-FMA version requires more flops. We can explore this choice further as we develop in Python.
2. **Floating Point Operations**: Standard floating point operations generally produce acceptable results. We can monitor their reliability by calculating the condition number. I'm working on an adaptive algorithm to address this in the future.
3. **Task Complexity**: As you might have noticed, all these floating point operations in Python inevitably lead to compiled C. Although all these packages use the same compiled C function, their binding and wrapping are different, making them difficult to connect properly.

### Conclusion:

Since FMA isn't currently ready for the Python package, we can focus on the general case first and handle the ill-conditioned cases later. The package should primarily use normal floating point operations, and when an ill-conditioned case arises, it's not overly time-consuming to apply the more accurate operations to this subset.

If you have any questions or need further clarification, please don't hesitate to ask. I'm more than happy to discuss and explain everything.

philipc2: > Maintain Current Approach: We can keep the current implementation as it is, but set is_FMA_enable to False globally and use the basic np.cross. We'll document this choice, informing users that they can opt for a more accurate method by enabling FMA.

I think this option makes the most sense to me and would fit the expectations of our users also. 

> And since we are not allowed to use any compiled files in our project code base, we cannot do much about it now.

I wouldn't say that we are **_not allowed_**, it's more something that we haven't explored yet. It might be worth considering. 

I appreciate the write up and glad that we're making some good improvements both to the performance and accuracy of our operators. 


philipc2: Merging this so that we can continue our work with #785 

Overall, was able to get about a 10-12x performance improvement by using `numba` without major rewrites to the algorithm. There's still some more optimizations we could do, but this would require significant rewrites, which should not be a priority right now. 

Patch:
diff --git a/benchmarks/face_bounds.py b/benchmarks/face_bounds.py
index c41602b6..b249e7b9 100644
--- a/benchmarks/face_bounds.py
+++ b/benchmarks/face_bounds.py
@@ -10,8 +10,6 @@ grid_geoflow = current_path / "test" / "meshfiles" / "ugrid" / "geoflow-small" /
 grid_scrip = current_path / "test" / "meshfiles" / "scrip" / "outCSne8" / "outCSne8.nc"
 grid_mpas= current_path / "test" / "meshfiles" / "mpas" / "QU" / "oQU480.231010.nc"
 
-
-
 class FaceBounds:
 
     params = [grid_quad_hex, grid_geoflow, grid_scrip, grid_mpas]
diff --git a/uxarray/__init__.py b/uxarray/__init__.py
index ed5c31ca..d15b3807 100644
--- a/uxarray/__init__.py
+++ b/uxarray/__init__.py
@@ -1,6 +1,8 @@
 # Sets the version of uxarray currently installeds
 # Attempt to import the needed modules
 
+import uxarray.constants
+
 from .core.api import open_grid, open_dataset, open_mfdataset
 
 from .core.dataset import UxDataset
@@ -21,6 +23,21 @@ except Exception:
     # Placeholder version incase an error occurs, such as the library isn't installed
     __version__ = "999"
 
+# Flag for enabling FMA instructions across the package
+
+
+def enable_fma():
+    """Enables Fused-Multiply-Add (FMA) instructions using the ``pyfma``
+    package."""
+    uxarray.constants.ENABLE_FMA = True
+
+
+def disable_fma():
+    """Disable Fused-Multiply-Add (FMA) instructions using the ``pyfma``
+    package."""
+    uxarray.constants.ENABLE_FMA = False
+
+
 __all__ = (
     "open_grid",
     "open_dataset",
@@ -30,4 +47,6 @@ __all__ = (
     "INT_DTYPE",
     "INT_FILL_VALUE",
     "Grid",
+    "enable_fma",
+    "disable_fma",
 )
diff --git a/uxarray/constants.py b/uxarray/constants.py
index d42befac..a8baa9c1 100644
--- a/uxarray/constants.py
+++ b/uxarray/constants.py
@@ -15,4 +15,6 @@ ERROR_TOLERANCE = 1.0e-8
 ENABLE_JIT_CACHE = True
 ENABLE_JIT = True
 
+ENABLE_FMA = False
+
 GRID_DIMS = ["n_node", "n_edge", "n_face"]
diff --git a/uxarray/grid/arcs.py b/uxarray/grid/arcs.py
index 381d02e4..63344473 100644
--- a/uxarray/grid/arcs.py
+++ b/uxarray/grid/arcs.py
@@ -2,9 +2,13 @@ import numpy as np
 
 # from uxarray.grid.coordinates import node_xyz_to_lonlat_rad, normalize_in_place
 
-from uxarray.grid.coordinates import _xyz_to_lonlat_rad, _normalize_xyz
+from uxarray.grid.coordinates import _xyz_to_lonlat_rad_no_norm, _normalize_xyz_scalar
 from uxarray.constants import ERROR_TOLERANCE
 
+from uxarray.utils.computing import isclose, cross, dot
+
+from numba import njit
+
 
 def _to_list(obj):
     if not isinstance(obj, list):
@@ -17,78 +21,29 @@ def _to_list(obj):
     return obj
 
 
-def point_within_gca(pt, gca_cart, is_directed=False):
-    """Check if a point lies on a given Great Circle Arc (GCA). The anti-
-    meridian case is also considered.
-
-    Parameters
-    ----------
-    pt : numpy.ndarray (float)
-        Cartesian coordinates of the point.
-    gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr)
-        Cartesian coordinates of the Great Circle Arc (GCR).
-    is_directed : bool, optional, default = False
-        If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected,
-        and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False.
-        For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1
-
-    Returns
-    -------
-    bool
-        True if the point lies between the two endpoints of the GCR, False otherwise.
-
-    Raises
-    ------
-    ValueError
-        If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes.
-        In such cases, consider breaking the GCR into two separate GCRs.
-
-    ValueError
-        If the input GCR spans more than 180 degrees (π radians).
-        In such cases, consider breaking the GCR into two separate GCRs.
-
-    Notes
-    -----
-    The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and
-    accounting for the anti-meridian case.
-
-    The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude).
-    In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation.
-
-    The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations.
-
-    Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons.
-    """
-    # Convert the cartesian coordinates to lonlat coordinates
-    pt_lonlat = np.array(_xyz_to_lonlat_rad(pt[0], pt[1], pt[2]))
-    GCRv0_lonlat = np.array(
-        _xyz_to_lonlat_rad(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2])
-    )
-    GCRv1_lonlat = np.array(
-        _xyz_to_lonlat_rad(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2])
-    )
-
-    # Convert the list to np.float64
-    gca_cart[0] = np.array(gca_cart[0], dtype=np.float64)
-    gca_cart[1] = np.array(gca_cart[1], dtype=np.float64)
-
-    # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes
+@njit
+def _point_within_gca_body(
+    angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed
+):
     angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1])
-    if np.allclose(angle, np.pi, rtol=0, atol=ERROR_TOLERANCE):
+    if isclose(angle, np.pi, rtol=0.0, atol=ERROR_TOLERANCE):
         raise ValueError(
             "The input Great Circle Arc is exactly 180 degree, this Great Circle Arc can have multiple planes. "
             "Consider breaking the Great Circle Arc"
             "into two Great Circle Arcs"
         )
 
-    if not np.allclose(
-        np.dot(np.cross(gca_cart[0], gca_cart[1]), pt), 0, rtol=0, atol=ERROR_TOLERANCE
+    if not isclose(
+        dot(cross(np.asarray(gca_cart[0]), np.asarray(gca_cart[1])), pt),
+        0,
+        rtol=0.0,
+        atol=ERROR_TOLERANCE,
     ):
         return False
 
-    if np.isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0, atol=ERROR_TOLERANCE):
+    if isclose(GCRv0_lonlat[0], GCRv1_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE):
         # If the pt and the GCA are on the same longitude (the y coordinates are the same)
-        if np.isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE):
+        if isclose(GCRv0_lonlat[0], pt_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE):
             # Now use the latitude to determine if the pt falls between the interval
             return in_between(GCRv0_lonlat[1], pt_lonlat[1], GCRv1_lonlat[1])
         else:
@@ -96,16 +51,16 @@ def point_within_gca(pt, gca_cart, is_directed=False):
             return False
 
     # If the longnitude span is exactly 180 degree, then the GCA goes through the pole point
-    if np.isclose(
-        abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0, atol=ERROR_TOLERANCE
+    if isclose(
+        abs(GCRv1_lonlat[0] - GCRv0_lonlat[0]), np.pi, rtol=0.0, atol=ERROR_TOLERANCE
     ):
         # Special case, if the pt is on the pole point, then set its longitude to the GCRv0_lonlat[0]
-        if np.isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0, atol=ERROR_TOLERANCE):
+        if isclose(abs(pt_lonlat[1]), np.pi / 2, rtol=0.0, atol=ERROR_TOLERANCE):
             pt_lonlat[0] = GCRv0_lonlat[0]
-        if not np.isclose(
-            GCRv0_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE
-        ) and not np.isclose(
-            GCRv1_lonlat[0], pt_lonlat[0], rtol=0, atol=ERROR_TOLERANCE
+        if not isclose(
+            GCRv0_lonlat[0], pt_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE
+        ) and not isclose(
+            GCRv1_lonlat[0], pt_lonlat[0], rtol=0.0, atol=ERROR_TOLERANCE
         ):
             return False
         else:
@@ -169,6 +124,69 @@ def point_within_gca(pt, gca_cart, is_directed=False):
             )
 
 
+def point_within_gca(pt, gca_cart, is_directed=False):
+    """Check if a point lies on a given Great Circle Arc (GCA). The anti-
+    meridian case is also considered.
+
+    Parameters
+    ----------
+    pt : numpy.ndarray (float)
+        Cartesian coordinates of the point.
+    gca_cart : numpy.ndarray of shape (2, 3), (np.float or gmpy2.mpfr)
+        Cartesian coordinates of the Great Circle Arc (GCR).
+    is_directed : bool, optional, default = False
+        If True, the GCA is considered to be directed, which means it can only from v0-->v1. If False, the GCA is undirected,
+        and we will always assume the small circle (The one less than 180 degree) side is the GCA. The default is False.
+        For the case of the anti-podal case, the direction is v_0--> the pole point that on the same hemisphere as v_0-->v_1
+
+    Returns
+    -------
+    bool
+        True if the point lies between the two endpoints of the GCR, False otherwise.
+
+    Raises
+    ------
+    ValueError
+        If the input GCR spans exactly 180 degrees (π radians), as this GCR can have multiple planes.
+        In such cases, consider breaking the GCR into two separate GCRs.
+
+    ValueError
+        If the input GCR spans more than 180 degrees (π radians).
+        In such cases, consider breaking the GCR into two separate GCRs.
+
+    Notes
+    -----
+    The function checks if the given point is on the Great Circle Arc by considering its cartesian coordinates and
+    accounting for the anti-meridian case.
+
+    The anti-meridian case occurs when the GCR crosses the anti-meridian (0 longitude).
+    In this case, the function handles scenarios where the GCA spans across more than 180 degrees, requiring specific operation.
+
+    The function relies on the `_angle_of_2_vectors` and `is_between` functions to perform the necessary calculations.
+
+    Please ensure that the input coordinates are in radians and adhere to the ERROR_TOLERANCE value for floating-point comparisons.
+    """
+    # Convert the cartesian coordinates to lonlat coordinates
+    pt_lonlat = np.array(_xyz_to_lonlat_rad_no_norm(pt[0], pt[1], pt[2]))
+    GCRv0_lonlat = np.array(
+        _xyz_to_lonlat_rad_no_norm(gca_cart[0][0], gca_cart[0][1], gca_cart[0][2])
+    )
+    GCRv1_lonlat = np.array(
+        _xyz_to_lonlat_rad_no_norm(gca_cart[1][0], gca_cart[1][1], gca_cart[1][2])
+    )
+    gca_cart = np.asarray(gca_cart)
+
+    # First if the input GCR is exactly 180 degree, we throw an exception, since this GCR can have multiple planes
+    angle = _angle_of_2_vectors(gca_cart[0], gca_cart[1])
+
+    out = _point_within_gca_body(
+        angle, gca_cart, pt, GCRv0_lonlat, GCRv1_lonlat, pt_lonlat, is_directed
+    )
+
+    return out
+
+
+@njit
 def in_between(p, q, r) -> bool:
     """Determines whether the number q is between p and r.
 
@@ -190,6 +208,7 @@ def in_between(p, q, r) -> bool:
     return p <= q <= r or r <= q <= p
 
 
+@njit
 def _decide_pole_latitude(lat1, lat2):
     """Determine the pole latitude based on the latitudes of two points on a
     Great Circle Arc (GCA).
@@ -229,6 +248,7 @@ def _decide_pole_latitude(lat1, lat2):
     return closest_pole
 
 
+@njit
 def _angle_of_2_vectors(u, v):
     """Calculate the angle between two 3D vectors u and v in radians. Can be
     used to calcualte the span of a GCR.
@@ -281,22 +301,24 @@ def extreme_gca_latitude(gca_cart, extreme_type):
         raise ValueError("extreme_type must be either 'max' or 'min'")
 
     n1, n2 = gca_cart
-    dot_n1_n2 = np.dot(n1, n2)
+
+    dot_n1_n2 = dot(np.asarray(n1), np.asarray(n2))
     denom = (n1[2] + n2[2]) * (dot_n1_n2 - 1.0)
     d_a_max = (n1[2] * dot_n1_n2 - n2[2]) / denom
 
     d_a_max = (
         np.clip(d_a_max, 0, 1)
-        if np.isclose(d_a_max, [0, 1], atol=ERROR_TOLERANCE).any()
+        if isclose(d_a_max, 0, atol=ERROR_TOLERANCE)
+        or isclose(d_a_max, 1, atol=ERROR_TOLERANCE)
         else d_a_max
     )
 
-    _, lat_n1 = _xyz_to_lonlat_rad(n1[0], n1[1], n1[2])
-    _, lat_n2 = _xyz_to_lonlat_rad(n2[0], n2[1], n2[2])
+    _, lat_n1 = _xyz_to_lonlat_rad_no_norm(n1[0], n1[1], n1[2])
+    _, lat_n2 = _xyz_to_lonlat_rad_no_norm(n2[0], n2[1], n2[2])
 
     if 0 < d_a_max < 1:
         node3 = (1 - d_a_max) * n1 + d_a_max * n2
-        node3 = np.array(_normalize_xyz(node3[0], node3[1], node3[2]))
+        node3 = np.array(_normalize_xyz_scalar(node3[0], node3[1], node3[2]))
         d_lat_rad = np.arcsin(np.clip(node3[2], -1, 1))
 
         return (
diff --git a/uxarray/grid/coordinates.py b/uxarray/grid/coordinates.py
index ceb8aceb..d34f0fc2 100644
--- a/uxarray/grid/coordinates.py
+++ b/uxarray/grid/coordinates.py
@@ -10,6 +10,8 @@ from typing import Union
 
 from numba import njit
 
+import math
+
 
 @njit(cache=True)
 def _lonlat_rad_to_xyz(
@@ -25,6 +27,47 @@ def _lonlat_rad_to_xyz(
     return x, y, z
 
 
+@njit
+def _xyz_to_lonlat_rad_no_norm(
+    x: Union[np.ndarray, float],
+    y: Union[np.ndarray, float],
+    z: Union[np.ndarray, float],
+):
+    """Converts a Cartesian x,y,z coordinates into Spherical latitude and
+    longitude without normalization, decorated with Numba.
+
+    Parameters
+    ----------
+    x : float
+        Cartesian x coordinate
+    y: float
+        Cartesiain y coordinate
+    z: float
+        Cartesian z coordinate
+
+
+    Returns
+    -------
+    lon : float
+        Longitude in radians
+    lat: float
+        Latitude in radians
+    """
+
+    lon = math.atan2(y, x)
+    lat = math.asin(z)
+
+    # set longitude range to [0, pi]
+    lon = np.mod(lon, 2 * np.pi)
+
+    z_mask = np.abs(z) > 1.0 - ERROR_TOLERANCE
+
+    lat = np.where(z_mask, np.sign(z) * np.pi / 2, lat)
+    lon = np.where(z_mask, 0.0, lon)
+
+    return lon, lat
+
+
 def _xyz_to_lonlat_rad(
     x: Union[np.ndarray, float],
     y: Union[np.ndarray, float],
@@ -126,6 +169,15 @@ def _normalize_xyz(
     return x_norm, y_norm, z_norm
 
 
+@njit
+def _normalize_xyz_scalar(x: float, y: float, z: float):
+    denom = np.linalg.norm(np.asarray(np.array([x, y, z]), dtype=np.float64), ord=2)
+    x_norm = x / denom
+    y_norm = y / denom
+    z_norm = z / denom
+    return x_norm, y_norm, z_norm
+
+
 def _populate_node_latlon(grid) -> None:
     """Populates the latitude and longitude coordinates of a Grid (`node_lon`,
     `node_lat`)"""
diff --git a/uxarray/grid/geometry.py b/uxarray/grid/geometry.py
index d2c43e8f..08229d1a 100644
--- a/uxarray/grid/geometry.py
+++ b/uxarray/grid/geometry.py
@@ -6,6 +6,8 @@ from uxarray.grid.utils import (
     _get_cartesian_face_edge_nodes,
     _get_lonlat_rad_face_edge_nodes,
 )
+
+from uxarray.utils.computing import isclose, allclose, all
 import warnings
 import pandas as pd
 import xarray as xr
@@ -467,7 +469,7 @@ def _check_intersection(ref_edge, edges):
         intersection_point = gca_gca_intersection(ref_edge, edge)
 
         if intersection_point.size != 0:
-            if np.allclose(intersection_point, pole_point, atol=ERROR_TOLERANCE):
+            if allclose(intersection_point, pole_point, atol=ERROR_TOLERANCE):
                 return True
             intersection_count += 1
 
@@ -488,9 +490,9 @@ def _classify_polygon_location(face_edge_cart):
         Returns either 'North', 'South' or 'Equator' based on the polygon's location.
     """
     z_coords = face_edge_cart[:, :, 2]
-    if np.all(z_coords > 0):
+    if all(z_coords > 0):
         return "North"
-    elif np.all(z_coords < 0):
+    elif all(z_coords < 0):
         return "South"
     else:
         return "Equator"
@@ -584,7 +586,7 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True):
     >>> _insert_pt_in_latlonbox(np.array([[1.0, 2.0], [3.0, 4.0]]),np.array([1.5, 3.5]))
     array([[1.0, 2.0], [3.0, 4.0]])
     """
-    if np.all(new_pt == INT_FILL_VALUE):
+    if all(new_pt == INT_FILL_VALUE):
         return old_box
 
     latlon_box = np.copy(old_box)  # Create a copy of the old box
@@ -612,18 +614,17 @@ def _insert_pt_in_latlonbox(old_box, new_pt, is_lon_periodic=True):
     # Check for pole points and update latitudes
     is_pole_point = (
         lon_pt == INT_FILL_VALUE
-        and np.isclose(
-            new_pt[0], [0.5 * np.pi, -0.5 * np.pi], atol=ERROR_TOLERANCE
-        ).any()
+        and isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE)
+        or isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE)
     )
 
     if is_pole_point:
         # Check if the new point is close to the North Pole
-        if np.isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE):
+        if isclose(new_pt[0], 0.5 * np.pi, atol=ERROR_TOLERANCE):
             latlon_box[0][1] = 0.5 * np.pi
 
         # Check if the new point is close to the South Pole
-        elif np.isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE):
+        elif isclose(new_pt[0], -0.5 * np.pi, atol=ERROR_TOLERANCE):
             latlon_box[0][0] = -0.5 * np.pi
 
         return latlon_box
@@ -769,9 +770,7 @@ def _populate_face_latlon_bound(
             )
 
             # Check if the node matches the pole point or if the pole point is within the edge_cart
-            if np.allclose(
-                n1_cart, pole_point, atol=ERROR_TOLERANCE
-            ) or point_within_gca(
+            if allclose(n1_cart, pole_point, atol=ERROR_TOLERANCE) or point_within_gca(
                 pole_point, np.array([n1_cart, n2_cart]), is_directed=False
             ):
                 is_center_pole = False
@@ -850,16 +849,16 @@ def _populate_face_latlon_bound(
             )
 
             # Insert extreme latitude points into the latlonbox if they differ from the node latitudes
-            if not np.isclose(
+            if not isclose(
                 node1_lat_rad, lat_max, atol=ERROR_TOLERANCE
-            ) and not np.isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE):
+            ) and not isclose(node2_lat_rad, lat_max, atol=ERROR_TOLERANCE):
                 # Insert the maximum latitude
                 face_latlon_array = _insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_max, node1_lon_rad])
                 )
-            elif not np.isclose(
+            elif not isclose(
                 node1_lat_rad, lat_min, atol=ERROR_TOLERANCE
-            ) and not np.isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE):
+            ) and not isclose(node2_lat_rad, lat_min, atol=ERROR_TOLERANCE):
                 # Insert the minimum latitude
                 face_latlon_array = _insert_pt_in_latlonbox(
                     face_latlon_array, np.array([lat_min, node1_lon_rad])
@@ -939,7 +938,7 @@ def _populate_bounds(
     intervals_tuple_list = []
     intervals_name_list = []
 
-    faces_edges_cartesian = _get_cartesian_face_edge_nodes(
+    face_edges_cartesian = _get_cartesian_face_edge_nodes(
         grid.face_node_connectivity.values,
         grid.n_face,
         grid.n_max_face_edges,
@@ -948,7 +947,7 @@ def _populate_bounds(
         grid.node_z.values,
     )
 
-    faces_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes(
+    face_edges_lonlat_rad = _get_lonlat_rad_face_edge_nodes(
         grid.face_node_connectivity.values,
         grid.n_face,
         grid.n_max_face_edges,
@@ -956,17 +955,26 @@ def _populate_bounds(
         grid.node_lat.values,
     )
 
-    for face_idx, face_nodes in enumerate(grid.face_node_connectivity):
-        face_edges_cartesian = faces_edges_cartesian[face_idx]
+    face_node_connectivity = grid.face_node_connectivity.values
 
-        # Skip processing if the face is a dummy face
-        if np.any(face_edges_cartesian == INT_FILL_VALUE):
-            continue
+    # # TODO: vectorize dummy face check
+    s = face_edges_cartesian.shape
+    dummy_face_face_edges_cart = np.any(
+        face_edges_cartesian.reshape((s[0], s[1] * s[2] * s[3])) == INT_FILL_VALUE,
+        axis=1,
+    )
+
+    s = face_edges_lonlat_rad.shape
+    dummy_face_face_edges_latlon = np.any(
+        face_edges_lonlat_rad.reshape((s[0], s[1] * s[2] * s[3])) == INT_FILL_VALUE,
+        axis=1,
+    )
 
-        face_edges_lonlat_rad = faces_edges_lonlat_rad[face_idx]
+    dummy_faces = dummy_face_face_edges_cart | dummy_face_face_edges_latlon
 
-        # Skip processing if the face is a dummy face
-        if np.any(face_edges_lonlat_rad == INT_FILL_VALUE):
+    for face_idx, face_nodes in enumerate(face_node_connectivity):
+        if dummy_faces[face_idx]:
+            # skip dummy faces
             continue
 
         is_GCA_list = (
@@ -974,14 +982,15 @@ def _populate_bounds(
         )
 
         temp_latlon_array[face_idx] = _populate_face_latlon_bound(
-            face_edges_cartesian,
-            face_edges_lonlat_rad,
+            face_edges_cartesian[face_idx],
+            face_edges_lonlat_rad[face_idx],
             is_latlonface=is_latlonface,
             is_GCA_list=is_GCA_list,
         )
 
-        assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1]
-        assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1]
+        # # do we need this ?
+        # assert temp_latlon_array[face_idx][0][0] != temp_latlon_array[face_idx][0][1]
+        # assert temp_latlon_array[face_idx][1][0] != temp_latlon_array[face_idx][1][1]
         lat_array = temp_latlon_array[face_idx][0]
 
         # Now store the latitude intervals in the tuples
@@ -1001,7 +1010,7 @@ def _populate_bounds(
         attrs={
             "cf_role": "face_latlon_bounds",
             "_FillValue": INT_FILL_VALUE,
-            "long_name": "Provides the latitude and longitude bounds for each face in radians.",
+            "long_name": "Latitude and Longitude bounds for each face in radians.",
             "start_index": INT_DTYPE(0),
             "latitude_intervalsIndex": intervalsIndex,
             "latitude_intervals_name_map": df_intervals_map,
diff --git a/uxarray/grid/grid.py b/uxarray/grid/grid.py
index 5b92c102..14e0245d 100644
--- a/uxarray/grid/grid.py
+++ b/uxarray/grid/grid.py
@@ -891,6 +891,7 @@ class Grid:
                 "Constructing of `Grid.bounds` has not been optimized, which may lead to a long execution time."
             )
             _populate_bounds(self)
+
         return self._ds["bounds"]
 
     @property
diff --git a/uxarray/grid/intersections.py b/uxarray/grid/intersections.py
index 3d3b03f4..d57969db 100644
--- a/uxarray/grid/intersections.py
+++ b/uxarray/grid/intersections.py
@@ -4,10 +4,12 @@ from uxarray.grid.utils import _newton_raphson_solver_for_gca_constLat
 from uxarray.grid.arcs import point_within_gca
 import platform
 import warnings
-from uxarray.utils.computing import cross_fma
+from uxarray.utils.computing import cross_fma, allclose, cross, dot
 
+import uxarray.constants
 
-def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False):
+
+def gca_gca_intersection(gca1_cart, gca2_cart):
     """Calculate the intersection point(s) of two Great Circle Arcs (GCAs) in a
     Cartesian coordinate system.
 
@@ -54,12 +56,13 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False):
     w0, w1 = gca1_cart
     v0, v1 = gca2_cart
 
-    # Compute normals and orthogonal bases using FMA
-    if fma_disabled:
-        w0w1_norm = np.cross(w0, w1)
-        v0v1_norm = np.cross(v0, v1)
-        cross_norms = np.cross(w0w1_norm, v0v1_norm)
+    if not uxarray.constants.ENABLE_FMA:
+        # Compute normals and orthogonal bases without FMA
+        w0w1_norm = cross(w0, w1)
+        v0v1_norm = cross(v0, v1)
+        cross_norms = cross(w0w1_norm, v0v1_norm)
     else:
+        # Compute normals and orthogonal bases using FMA
         w0w1_norm = cross_fma(w0, w1)
         v0v1_norm = cross_fma(v0, v1)
         cross_norms = cross_fma(w0w1_norm, v0v1_norm)
@@ -73,29 +76,29 @@ def gca_gca_intersection(gca1_cart, gca2_cart, fma_disabled=False):
             )
 
     # Check perpendicularity conditions and floating-point arithmetic limitations
-    if not np.allclose(
-        np.dot(w0w1_norm, w0), 0, atol=ERROR_TOLERANCE
-    ) or not np.allclose(np.dot(w0w1_norm, w1), 0, atol=ERROR_TOLERANCE):
+    if not allclose(dot(w0w1_norm, w0), 0, atol=ERROR_TOLERANCE) or not allclose(
+        dot(w0w1_norm, w1), 0, atol=ERROR_TOLERANCE
+    ):
         warnings.warn(
             "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution."
         )
 
-    if not np.allclose(
-        np.dot(v0v1_norm, v0), 0, atol=ERROR_TOLERANCE
-    ) or not np.allclose(np.dot(v0v1_norm, v1), 0, atol=ERROR_TOLERANCE):
+    if not allclose(dot(v0v1_norm, v0), 0, atol=ERROR_TOLERANCE) or not allclose(
+        dot(v0v1_norm, v1), 0, atol=ERROR_TOLERANCE
+    ):
         warnings.warn(
             "The current input data cannot be computed accurately using floating-point arithmetic.  Use with caution. "
         )
 
-    if not np.allclose(
-        np.dot(cross_norms, v0v1_norm), 0, atol=ERROR_TOLERANCE
-    ) or not np.allclose(np.dot(cross_norms, w0w1_norm), 0, atol=ERROR_TOLERANCE):
+    if not allclose(
+        dot(cross_norms, v0v1_norm), 0, atol=ERROR_TOLERANCE
+    ) or not allclose(dot(cross_norms, w0w1_norm), 0, atol=ERROR_TOLERANCE):
         warnings.warn(
             "The current input data cannot be computed accurately using floating-point arithmetic. Use with caution. "
         )
 
     # If the cross_norms is zero, the two GCAs are parallel
-    if np.allclose(cross_norms, 0, atol=ERROR_TOLERANCE):
+    if allclose(cross_norms, 0, atol=ERROR_TOLERANCE):
         return np.array([])
 
     # Normalize the cross_norms
@@ -156,7 +159,7 @@ def gca_constLat_intersection(
     x1, x2 = gca_cart
 
     if fma_disabled:
-        n = np.cross(x1, x2)
+        n = cross(x1, x2)
 
     else:
         # Raise a warning for Windows users
diff --git a/uxarray/utils/computing.py b/uxarray/utils/computing.py
index 3d801e92..f9929deb 100644
--- a/uxarray/utils/computing.py
+++ b/uxarray/utils/computing.py
@@ -1,6 +1,65 @@
 import numpy as np
 import sys
 
+from numba import njit
+
+
+@njit
+def all(a):
+    """Numba decorated implementation of ``np.all()``
+
+    See Also
+    --------
+    numpy.all
+    """
+
+    return np.all(a)
+
+
+@njit
+def isclose(a, b, rtol=1e-05, atol=1e-08):
+    """Numba decorated implementation of ``np.isclose()``
+
+    See Also
+    --------
+    numpy.isclose
+    """
+
+    return np.isclose(a, b, rtol=rtol, atol=atol)
+
+
+@njit
+def allclose(a, b, rtol=1e-05, atol=1e-08):
+    """Numba decorated implementation of ``np.allclose()``
+
+    See Also
+    --------
+    numpy.allclose
+    """
+    return np.allclose(a, b, rtol=rtol, atol=atol)
+
+
+@njit
+def cross(a, b):
+    """Numba decorated implementation of ``np.cross()``
+
+    See Also
+    --------
+    numpy.cross
+    """
+    return np.cross(a, b)
+
+
+@njit
+def dot(a, b):
+    """Numba decorated implementation of ``np.dot()``
+
+    See Also
+    --------
+    numpy.dot
+    """
+    return np.dot(a, b)
+
 
 def _fmms(a, b, c, d):
     """


Test patch:
diff --git a/test/test_arcs.py b/test/test_arcs.py
index c9859e72..17d04fe5 100644
--- a/test/test_arcs.py
+++ b/test/test_arcs.py
@@ -36,7 +36,7 @@ class TestIntersectionPoint(TestCase):
         ]
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
+            point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart))
 
         gcr_180degree_cart = [
             _lonlat_rad_to_xyz(0.0, np.pi / 2.0),
@@ -45,7 +45,7 @@ class TestIntersectionPoint(TestCase):
 
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(pt_same_lon_in, gcr_180degree_cart)
+            point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_180degree_cart))
 
         # Test when the point and the GCA all have the same longitude
         gcr_same_lon_cart = [
@@ -53,19 +53,19 @@ class TestIntersectionPoint(TestCase):
             _lonlat_rad_to_xyz(0.0, -1.5)
         ]
         pt_same_lon_in = _lonlat_rad_to_xyz(0.0, 0.0)
-        self.assertTrue(point_within_gca(pt_same_lon_in, gcr_same_lon_cart))
+        self.assertTrue(point_within_gca(np.asarray(pt_same_lon_in), np.asarray(gcr_same_lon_cart)))
 
         pt_same_lon_out = _lonlat_rad_to_xyz(0.0, 1.500000000000001)
-        res = point_within_gca(pt_same_lon_out, gcr_same_lon_cart)
+        res = point_within_gca(np.asarray(pt_same_lon_out), np.asarray(gcr_same_lon_cart))
         self.assertFalse(res)
 
         pt_same_lon_out_2 = _lonlat_rad_to_xyz(0.1, 1.0)
-        res = point_within_gca(pt_same_lon_out_2, gcr_same_lon_cart)
+        res = point_within_gca(np.asarray(pt_same_lon_out_2), np.asarray(gcr_same_lon_cart))
         self.assertFalse(res)
 
         # And if we increase the digital place by one, it should be true again
         pt_same_lon_out_add_one_place = _lonlat_rad_to_xyz(0.0, 1.5000000000000001)
-        res = point_within_gca(pt_same_lon_out_add_one_place, gcr_same_lon_cart)
+        res = point_within_gca(np.asarray(pt_same_lon_out_add_one_place), np.asarray(gcr_same_lon_cart))
         self.assertTrue(res)
 
         # Normal case
@@ -76,7 +76,7 @@ class TestIntersectionPoint(TestCase):
                                                         -0.997]])
         pt_cart_within = np.array(
             [0.25616109352676675, 0.9246590335292105, -0.010021496695000144])
-        self.assertTrue(point_within_gca(pt_cart_within, gcr_cart_2, True))
+        self.assertTrue(point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_2), True))
 
         # Test other more complicate cases : The anti-meridian case
 
@@ -87,16 +87,16 @@ class TestIntersectionPoint(TestCase):
         gcr_cart = np.array([[0.351, -0.724, 0.593], [0.617, 0.672, 0.410]])
         pt_cart = np.array(
             [0.9438777657502077, 0.1193199333436068, 0.922714737029319])
-        self.assertTrue(point_within_gca(pt_cart, gcr_cart, is_directed=True))
+        self.assertTrue(point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True))
         # If we swap the gcr, it should throw a value error since it's larger than 180 degree
         gcr_cart_flip = np.array([[0.617, 0.672, 0.410], [0.351, -0.724,
                                                           0.593]])
         with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart_flip, is_directed=True)
+            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=True)
 
         # If we flip the gcr in the undirected mode, it should still work
         self.assertTrue(
-            point_within_gca(pt_cart, gcr_cart_flip, is_directed=False))
+            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart_flip), is_directed=False))
 
         # 2nd anti-meridian case
         # GCR vertex0 in radian : [4.104711496596806, 0.5352983676533828],
@@ -107,9 +107,9 @@ class TestIntersectionPoint(TestCase):
         pt_cart_within = np.array(
             [0.6136726305712109, 0.28442243941920053, -0.365605190899831])
         self.assertFalse(
-            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=True))
+            point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=True))
         self.assertFalse(
-            point_within_gca(pt_cart_within, gcr_cart_1, is_directed=False))
+            point_within_gca(np.asarray(pt_cart_within), np.asarray(gcr_cart_1), is_directed=False))
 
         # The first case should not work and the second should work
         v1_rad = [0.1, 0.0]
@@ -119,10 +119,10 @@ class TestIntersectionPoint(TestCase):
         gcr_cart = np.array([v1_cart, v2_cart])
         pt_cart = _lonlat_rad_to_xyz(0.01, 0.0)
         with self.assertRaises(ValueError):
-            point_within_gca(pt_cart, gcr_cart, is_directed=True)
+            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_cart), is_directed=True)
         gcr_car_flipped = np.array([v2_cart, v1_cart])
         self.assertTrue(
-            point_within_gca(pt_cart, gcr_car_flipped, is_directed=True))
+            point_within_gca(np.asarray(pt_cart), np.asarray(gcr_car_flipped), is_directed=True))
 
     def test_pt_within_gcr_cross_pole(self):
         gcr_cart = np.array([[0.351, 0.0, 0.3], [-0.351, 0.0, 0.3]])
diff --git a/test/test_grid.py b/test/test_grid.py
index 10faba69..57eb3f3b 100644
--- a/test/test_grid.py
+++ b/test/test_grid.py
@@ -952,8 +952,7 @@ class TestLatlonBounds(TestCase):
         face_bounds = bounds_xarray.values
         nt.assert_allclose(grid.bounds.values, expected_bounds, atol=ERROR_TOLERANCE)
 
-    def test_populate_bounds_MPAS(self):
-        xrds = xr.open_dataset(self.gridfile_mpas)
-        uxgrid = ux.Grid.from_dataset(xrds, use_dual=True)
-        bounds_xarray = uxgrid.bounds
-        pass
+    def test_opti_bounds(self):
+        import uxarray
+        uxgrid = ux.open_grid(gridfile_CSne8)
+        bounds = uxgrid.bounds