Dimension issue when calculating passage flow rates for DART mesh run

[JanStreffing]

I tried to calculate the Drake Passage flow rates for Tidos DARTN AWI-CM3 TCO319-DART simulation on aleph today. I had previously used this for the CORE2 mesh without issues.

I now use:

FESOM mesh:
path                  = /proj/awi/input/fesom2/dart
alpha, beta, gamma    = 0.0, 0.0, 0.0
number of 2d nodes    = 3160340
number of 2d elements = 6262485

That data loads but the computation runs into the following dimensional issue:

Loading the data into memory...
A lot of velocity data (312.67GB)... This will take some time...
[########################################] | 100% Completed |  5min 10.4s

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 transport_DP, section_DP = pf.cross_section_transport([-65,-65,-55,-65.8],                                                          # select a section from the presets or [lon1, lon2, lat1, lat2]
      2                                                 mesh=mesh,                                                 # mesh
      3                                                 data_path='/scratch/awiiccp5/DARTN/outdata/fesom/',   # directory of the u, v, (T, S, uice, vice, m_ice, a_ice) files
      4                                                 mesh_diag=mesh_diag,                                            # mesh_diag
      5                                                 years=np.arange(2050,2129),                                     # years to compute
      6                                                 use_great_circle=False,                                  # compute the section as a great circle
      7                                                 how='ori',                                               # 'ori' do not apply mean, 'mean' apply time mean
      8                                                 add_TS=False,                                            # add temperature and salinity to the section
      9                                                 add_extent=2,                                             # the extent to look for gridcells nerby the section, choose large for low resolutions
     10                                                 n_points=1000                                          # number of waypoints between the start and end of the section
     11                                                )

File /mnt/lustre/home/awiiccp2/software/pyfesom2/pyfesom2/transport.py:1095, in cross_section_transport(section, mesh, data_path, years, mesh_diag, how, add_extent, abg, add_TS, chunks, use_great_circle, n_points)
   1091 effective_dx, effective_dy = _ComputeBrokenLineSegments(f_lat, f_lon, s_lat, s_lon, c_lat, c_lon, section)
   1093 vertical_cell_area_dx, vertical_cell_area_dy = _CreateVerticalGrid(effective_dx, effective_dy, mesh_diag)
-> 1095 ds = _UnrotateLoadVelocity(how, files, elem_box_indices, elem_box_nods, vertical_cell_area_dx, vertical_cell_area_dy, c_lon, c_lat, effective_dx, effective_dy, elem_order, chunks, mesh, abg)
   1097 ds = _TransportAcross(ds)
   1099 if add_TS:

File /mnt/lustre/home/awiiccp2/software/pyfesom2/pyfesom2/transport.py:921, in _UnrotateLoadVelocity(how, files, elem_box_indices, elem_box_nods, vertical_cell_area_dx, vertical_cell_area_dy, c_lon, c_lat, effective_dx, effective_dy, elem_order, chunks, mesh, abg)
    919 lon_elem_center = np.mean(mesh.x2[ds.elem_nods], axis=1)
    920 lat_elem_center = np.mean(mesh.y2[ds.elem_nods], axis=1)
--> 921 u, v = vec_rotate_r2g(abg[0], abg[1], abg[2], lon_elem_center[np.newaxis, :, np.newaxis],
    922                          lat_elem_center[np.newaxis, :, np.newaxis], ds.u_rot.values, ds.v_rot.values, flag=1)
    924 ds['u'] = (('time', 'elem', 'nz1'), u)
    925 ds['v'] = (('time', 'elem', 'nz1'), v)

File /mnt/lustre/home/awiiccp2/software/pyfesom2/pyfesom2/ut.py:246, in vec_rotate_r2g(al, be, ga, lon, lat, urot, vrot, flag)
    243 lon = lon * rad
    245 #   vector in rotated Cartesian
--> 246 txg = -vrot * np.sin(rlat) * np.cos(rlon) - urot * np.sin(rlon)
    247 tyg = -vrot * np.sin(rlat) * np.sin(rlon) + urot * np.cos(rlon)
    248 tzg = vrot * np.cos(rlat)

ValueError: operands could not be broadcast together with shapes (79,79,274) (1,274,1)

This seems to be quite far down, so I thought I ask here first before I go on a wild goose chase. Is this a known error with known solution?

[koldunovn]

This is code from @FinnHeu , so maybe he can react quickly to this. At first glance it seems that the code assumes [time, elem, depth] position of the dimensions, and you try to use new data from refactoring, which are [time, depth, elem].

[FinnHeu]

Sorry for the delayed reply @JanStreffing , I was on holiday. As far as I know, this error was not encountered before, so there is currently no known solution to this. I just tried exactly the same section [-65,-65,-55,-65.8] with my FESOM2.1 standalone runs on fArc-mesh and it worked perfectly fine. As @koldunovn already pointed out it seems the dimensions of your velocity data are in an unexpected order so the broadcasting fails. In your case instead of (79, 79, 274) it should be (79, 274, 79) to allow broadcasting with (1, 274, 1). My FESOM2.1 velocity output is in order [time, elem, nz1]. Please check if your velocity data is in the same dimensional order. In case it is not, transposing the velocity data before starting the calculation should do the job. If this dimensional order is the common thing for FESOM2.1 or in AWI-CM3 now (is it?), I can add a boolean switch as input to the code to handle both cases.

[JanStreffing]

I checked with new AWI-CM3.1 simulations and those also work fine. It seems to be an issue with the DART simulations on aleph. The dimension order of my U and V data on aleph are: time, nz1, elem

    try:
        u, v = vec_rotate_r2g(abg[0], abg[1], abg[2], lon_elem_center[np.newaxis, :, np.newaxis],
                             lat_elem_center[np.newaxis, :, np.newaxis], ds.u_rot.values, ds.v_rot.values, flag=1)
    except:
        u, v = vec_rotate_r2g(abg[0], abg[1], abg[2], lon_elem_center[np.newaxis, :, np.newaxis],
                             lat_elem_center[np.newaxis, :, np.newaxis], ds.u_rot.values.swapaxes(1,2), ds.v_rot.values.swapaxes(1,2), flag=1)

Which no longer crashes, but also the result does not make sense. For net Drake Passage flow through i get 28.84 Sv. Does this mistake potentially appear elsewhere without caused a crash, just bad results?

[FinnHeu]

Hm, that's weird. I have never encountered such bad results but I only used the software on high-resolution grids in the Arctic so far. I checked the Drake Passage section with my JRA55 hindcast (fArc) and it's producing a reasonable net transport of ~109 Sv and the u-section also looks okay to me (not an Antarctica expert): To check if the transport calculation goes wrong somewhere with your output you can compute e.g. a ring around Antarctica and check if the mean net transport is close to zero, so volume is conserved in a closed volume (https://github.com/FinnHeu/pyfesom2/blob/master/notebooks/transports.ipynb at the bottom). Currently, the software can't handle a full ring so you need to end the section below 180°E.

ds, section = pf.cross_section_transport(section=[-180,179.2,-60,-60],... , add_extent=.7)

Depending on your mesh you can get closer to 180°E by decreasing the add_extent parameter. If you get too close an index error will show up, you need to try out what fits best. All distances are recomputed in the calculation as they are currently not available in the fesom output, this adds a small error as well that might be larger on low-resolution meshes. At least in my fArc case this also produces good results even in the low-resolution part of the grid:

To me it seems the software is doing okay in Drake Passage, are you sure your output is alright? Maybe check Fram Strait transport or another section to eliminate your model results as the source of this strange transport.

[JanStreffing]

Hey Finn, I plotted the section and i don't see anything wrong with it:

The continental schelf backflow seems very strong, resulting is low total eastward flow:

I'm not sure if this really explains the low net value completely.