Script to difference (i.e. Baseline compare) two MPAS meshes using L2 Error Norm.

diff.py

#!/usr/bin/env python3

import click
import numpy as np 
import xarray as xr 

def reindex(ds, var):
    
    dims = ds[var].dims
    
    _dict = None

    if "nCells" in dims:
        _dict = dict(nCells = ds.indexToCellID - 1)
    elif "nEdges" in dims:
        _dict = dict(nEdges = ds.indexToEdgeID - 1)
    elif "nVertices" in dims: 
        _dict = dict(nVertices = ds.indexToVertexID - 1)
    
    if _dict: 
        return ds[var].sel(_dict)
    else: 
        return ds[var]

def calc_diff(ds1, ds2, var): 

    da1 = reindex(ds1, var)
    da2 = reindex(ds2, var)

    return da2 - da1

@click.command()
@click.option(
    '-v', '--verbose', is_flag=True, default=False,
    help='Print variables that are different'
    )
@click.argument(
    'file_a', type=click.Path(exists=True), #help="Filepath to file_a"
    )
@click.argument(
    'file_b', type=click.Path(exists=True), #help="Filepath to file_b"
    )
def difference(file_a, file_b, verbose=False):

    ds1 = xr.open_dataset(file_a)
    ds2 = xr.open_dataset(file_b)

    different_vars = [] 

    for var in ds1: 
        
        assert var in ds2
    
        diff = calc_diff(ds1, ds2, var)
        
        L2_norm = np.linalg.norm(diff.values.ravel(), 2)
        Linf_norm = np.linalg.norm(diff.values.ravel(), np.inf)
        
        if (L2_norm > 0.0) | (Linf_norm > 0.0): 

            print(f"{var:<40} L2   Norm={L2_norm:1.3e}, Linf Norm={Linf_norm:1.3e}")
            different_vars.append(var)

    n = len(different_vars) 
    
    if n != 0: 
        click.echo(f"ERROR: {n} fields are not within tolerance")
    else: 
        click.echo(f"SUCCESS: all fields are within tolerance")

if __name__ == "__main__":
    difference()