Corners for SCRIPR/CONSERV 1st order conservative remapping

[JanStreffing]

Currently we are using global conservative flux residual redistribution for AWICM (1,2 & 3). Never versions of OASIS3 (MCT3 & MCT4) are offering a locally conservative remapping. I had a short discussion with the EC-Earth community and their idea is to switch from GAUSWGT remapping mit conservative redistribution to CONSERV remapping. @helgegoessling was also interested

1. I think this would in principle also be of interest to us. Maybe @dsidoren, can comment.
2. Similar to how we already write out the center lon/lat (coord_nod2D(1, i) & coord_nod2D(2, i)) for the current remapping schemes we would need corner lon/lats for CONSERV.

Would this from oce_mesh be the right vector? mesh%x_corners(n,mesh%nod_in_elem2D_num(n))

For reference, this is what the OASIS manual states:

If the SCRIPR/CONSERV remapping is specified, longitudes and latitudes for the source and target grid corners must also be available in the grids.nc file as double precision REAL arrays dimensioned(nx,ny,4) or (nbrpts,1,4) where 4 is the number of corners (in the counterclockwize sense, starting by any corner). The names of the arrays must be composed of the grid prefix and the suffix “.clo” or “.cla” for respectively the grid corner longitudes or latitudes. As for the other grid information, the corners can be provided in grids.nc before the run by the user or directly by the component code through specific calls (see section 2.2.4). Longitudes must be given in degrees East in the interval -360.0 to 720.0. Latitudes must be given in degrees North in the interval -90.0 to 90.0. Note that if some grid points overlap, it is recommended to define those points with the same number (e.g. 360.0 for both, not 450.0 for one and 90.0 for the other) to ensure automatic detection of overlap by OASIS3-MCT. The corners of a cell cannot be defined modulo 360 degrees. For example, a cell located over Greenwich will have to be defined with corners at -1.0 deg and 1.0 deg but not with corners at 359.0deg and 1.0 deg.

[helgegoessling]

One needs to distinguish between (1) "cell" quantities (u,v(,w); are they used in the coupling?), centered at triangle centers, which have just the three triangle nodes as boundaries, and (2) "node" quantities which have a more complex polygon composed of the adjacent triangle centers and edge medians (on average 12 points) around them. Both of these are contained in the separate "CDO grid description" files that are provided with the FESOM meshes.

The OASIS manual states those must be 4 points, which would mean it can be done only with regular meshes. Similarly, the "SCRIP" library seems to have (or had) issues with unstructured meshes, which is why extra remapping functions have been introduced in CDOs a couple of years ago (e.g., "remapcon" based on SCRIP didn't work, whereas "remapycon" ("y" comes from "YAC") works also with our crazy polygons (which have concave segment junctions at places)). Are you sure that the new OASIS versions can only use SCRIP for the remapping? Then we might need to compute the weights offline with CDOs.

[JanStreffing]

Interesting, yes. I had been looking at the description of the SCRIPR/CONSERV section. It thats that it supports LR, D and U grids, meaning linear rectangular, gaussian and unstructured. Yet for the definition of corners they ask for 4 exactly, which sort of contradicts this.

CONSERV performs 1st or 2nd order conservative remapping, which means that the weight of a source cell is proportional to area intersected by the target cell (plus some other terms proportional to the gradient of the field in the longitudinal and latitudinal directions for the second order).The configuring line is:# SCRIPR (for CONSERV) $CMETH $CGRS $CFTYP $REST $NBIN $NORM $ORDER where: $CMETH = CONSERV. $CGRS is the source grid type: LR, D and U. Note that the grid corners have to given by the user in the grid data file grids.nc or by the code itself in the initialisation phase by calling routine oasiswritecorner

I'll write an email asking for clarification.

[koldunovn]

I am not sure if it's super relative here, but at some point we discuss with @trackow , that coupling with velocities can be done with velocities interpolated to nodes that are internally available in FESOM2 and that might have positive implications for computational speed.

[dsidoren]

conservative remapping is only needed for A2O. currently we send only the scalar fields to the ocean and the allocation mesh%x_corners(n,mesh%nod_in_elem2D_num(n)) is correct. I will check how it is computed. Also as I remember in the earlier versions of OASIS the corners must have been ordered (clockwise or so) and do not know how it is now.

[JanStreffing]

Hello everyone, I got a message back from the oasis team (a while ago already)

Dear Jan,

Thanks for you mail.
Indeed there is an error in the User Guide. The number of corners/vertices can be anything but must be the same for all cells. If your cells have a varying number of corners/vertices, you should used the largest number and just repeat the last one as needed.
Finally, please pay attention to the fact that for second order conservative remapping, it is the model that must provide the gradients of the field. This is explained in details in section 2.2.7, on p. 19. Please note that we made a correction to the user guide in this paragraph about two months ago so you should make sure to use an up-to-date version.

I hope this helps and please do not hesitate if you have additional questions,

regards,
Sophie

So it can be done with unstructured meshes.

However the calculation of the gradient fields would require some attention still, as this would have to be done in the source code of the respective atmospheric model. I hope this can be done without extra MPI communication at the edges.

Since the EC-Earth community is also interested and would have to calculate these gradients as well, I may have a chat with some time with e.g. @uwefladrich

[JanStreffing]

I have found some time to revive this topic. In https://github.com/FESOM/fesom2/tree/feat/oasis_corners_2nd_try I have implemented the changes needed for FESOM2 to write out its 'corner lat lon', in this case using only the points C1-6 in: https://fesom2.readthedocs.io/en/latest/geometry.html#

We could make it better by using the vertices center points as well, but that's for later. While I can generate 1st order conservative remapping files with the resulting grids.nc, the unit test of those remapping files fail. And when plugged into AWI-CM3 we get a blowup. I have contacted the oasis team about this, hopefully they can point out what's wrong with the clo & cla fields. If any you want to have a look: https://swiftbrowser.dkrz.de/public/dkrz_b94ae58d-50e8-4ad6-b453-afff299daabf/oasis_remappings/

[koldunovn]

Do you implement it in a way that we can output them in fesom.mesh.diag.nc?

[JanStreffing]

I calculate the values distributed and communicate them to the localroot. I the output them in grids.nc via the oasis_write_corner function. I'm unsure how much work it would be to call the fesom.mesh.diag.nc from after here. What do you think?

fesom2/src/cpl_driver.F90

Lines 414 to 415 in ced6233

	print *, 'FESOM before write corner'
	CALL oasis_write_corner (grid_name, number_of_all_points, 1, rmax, all_x_corners(:,:,:), all_y_corners(:,:,:))

[koldunovn]

Oh, well, if it involves oasis it's a bit too much trouble :) I was thinking of having the corners outputted in the stand alone :)

[JanStreffing]

Perhaps, yes. You would have to copy past the computation to a place that does not sit behind ifdef oasis. It would work there, but it's a bit of work. Not 100% for free.

[patrickscholz]

@JanStreffing Stupid Question from my side: You need the element centroids as the corner points of the scalar cell, right? To compute the area or whatever you do for your conservative remapping. But are you aware that the scalar cell is not just defined by the element centroids but also by the edge mid points? So we connect element centroid --> edge mid points --> element centroid --> edge mid point ... . So our scalar cell is at the end some weirdly shaped polygon.

[JanStreffing]

@patrickscholz Yes, I am aware that I'm makeing an error by neglecting half the polygon points at the moment. I tried to write this here:

FESOM2 to write out its 'corner lat lon', in this case using only the points C1-6 in:

We could make it better by using the vertices center points as well, but that's for later.
Originally posted by @JanStreffing in #21 (comment)

but I probably used wrong vocabulary. Oasis sounds flexible enough to deal with the strange shapes that we feed it.

[JanStreffing]

So far I implemented the element centroid corners on #432
@patrickscholz If I want to avoid the error made by assuming no grid distortion, I need the edge mid points as you mentioned. I had a look through oce_mesh if they would be available already like the element centroids. I didn't find them, correct? I could calculate them, eg. based on element edge length + the nodal center coords +element edge angle. The latter I could also compute from having the two nodal center coords at the two ends of the element edge. Can you tell me what I have to work with and what I would need to calculate myself?

[patrickscholz]

Richtig die edge midpkt coordinaten werden im moment nicht berechnet/gespeichert. Ist aber relativ einfach zu berechnen: du brauchst eigentlich nur das edge(1:2, mydim_edge) array darin stehen die beiden knoten indices die die edge bilden . Mit den knoten indices kannst du die xy coordinaten der beiden edge punkte bestimmen. Und dann must du nur über die coordinaten der beiden edge punkte mitteln und solltest den edge mittel punkt erhalten.

[JanStreffing]

Right.. much easier! I'll do that.

[JanStreffing]

Acutally, I found:

fesom2/src/oce_mesh.F90

Lines 2080 to 2096 in f96cf3d

	subroutine edge_center(n1, n2, x, y, mesh)
	USE MOD_MESH
	USE o_PARAM
	USE g_CONFIG
	implicit none
	integer :: n1, n2 ! nodes of the edge
	real(kind=WP), intent(inout) :: x, y
	type(t_mesh), intent(inout), target :: mesh
	real(kind=WP) :: a(2), b(2)
	a=mesh%coord_nod2D(:,n1)
	b=mesh%coord_nod2D(:,n2)
	if(a(1)-b(1)> cyclic_length/2.0_WP) a(1)=a(1)-cyclic_length
	if(a(1)-b(1)<-cyclic_length/2.0_WP) b(1)=b(1)-cyclic_length
	x=0.5_WP*(a(1)+b(1))
	y=0.5_WP*(a(2)+b(2))

I don't quite understand lines 93 and 94, but before we take two nodes, and after average their coordinates and call it edge center. So while we don't have a vector storing all of them, the routine of calculation is already there.

[patrickscholz]

Line 93 and 94 consider edges that span across the periodic boundary close to lon = -180-->+180 if you would take there simply the arithmetic mean you would end up at lon~0° for the edge midpoint which would be wrong therefore you correct one of the edge points by the length of the cyclic boundary. Like you have it there looks ok !

[patrickscholz]

We dont store these coordinates in the model because we dont need them for the computation itself, what we therefore store are dx & dy the vectors from the edge midpoints towards both of the element centroids that participate on that edge this is stored i think in edge_cross_dxdy (something like that) this is what we actually need to compute the fluxes in the tracer advection. I didnt remember that there was already a routine to compute the edge centers ... than just use that!

[helgegoessling]

That code to compute edge midpoints by averaging lon-lat type coordinates will still be quite inaccurate close to the poles. It's thus safer to transform to Cartesian (x,y,z) coordinates, average those, and then transform back to lon,lat. However, if this is all done in rotated coordinates where the poles are in Greenland and Antarctica, then the inaccuracy without the transformation should be small and thus OK.

[koldunovn]

@pgierz isn't it something you did for #354 ?

[JanStreffing]

That code to compute edge midpoints by averaging lon-lat type coordinates will still be quite inaccurate close to the poles. It's thus safer to transform to Cartesian (x,y,z) coordinates, average those, and then transform back to lon,lat. However, if this is all done in rotated coordinates where the poles are in Greenland and Antarctica, then the inaccuracy without the transformation should be small and thus OK.

I am first calculating the mid point and then unrotating. So perhaps it's ok?

fesom2/src/cpl_driver.F90

Lines 341 to 346 in 9ac6d19

	do i = 1, my_number_of_points
	do n = 1, nn_num(i)
	call edge_center(i, nn_pos(n+1,i), this_x_coord, this_y_coord, mesh)
	call r2g(coord_e_edge_center(1,i,n), coord_e_edge_center(2,i,n), this_x_coord, this_y_coord)
	end do
	end do

The resulting coords look ok to me, but the simulation stops with: srun: error: l30549: tasks 1-4: Exited with exit code 152 and I have not yet found the real error in the lengthy logfile. I think it has something to do with the lines above. I'm skipping the first neighbour, because that's the center point of the node again. But probably I should skip by starting from 2 instead of going to n+1

[pgierz]

@koldunovn, if I understand @JanStreffing correctly: no. I just re-shuffled already existing information with "more beautiful" names to get it to fit into the U-GRID standards, but did not add any new information to the output.

Code 152 is (normally) a "request acquisition error", unless you are using an exotic compiler with different error codes (I severely doubt that case). Probably an MPI resource isn't being cleanly released...? But that is admittedly poking in the dark. It could however also be that you're accessing memory where there is just array garbage. I could give the error a skeptical look until it explains itself, let me know where to look if you want that help...

[pgierz]

The intel handbook reference is here: https://www.intel.com/content/www/us/en/docs/fortran-compiler/developer-guide-reference/2023-0/list-of-runtime-error-messages.html

[dsidoren]

That code to compute edge midpoints by averaging lon-lat type coordinates will still be quite inaccurate close to the poles. It's thus safer to transform to Cartesian (x,y,z) coordinates, average those, and then transform back to lon,lat. However, if this is all done in rotated coordinates where the poles are in Greenland and Antarctica, then the inaccuracy without the transformation should be small and thus OK.

I fully agree. (1) even if you do the rotation differently the lon/lat representation is by itself inaccurate close to the poles and (2) inside the model we operate on rotated meshes only

[JanStreffing]

I was able to write the edge and element based corners. I did not yet go and upgrade the edge center calculation to Cathesian coordinates, because it currently looks like SCRIP can't deal with our convex control volume geometry anyway. Once again the issue is on hold.

[JanStreffing]

Superseded by #550