C grid noise in corners

[JFLemieux73]

See also issue #791.

[JFLemieux73]

There is noise in the thickness and concentration fields in corners with the C grid. Here is again the figure from issue #791:

Note that this is obtained with the uniform grid.

[JFLemieux73]

Reducing the time step from 60 min to 15 min does not change the pattern:

[JFLemieux73]

Hum...the noise disappears with upwind:

[dabail10]

I do think that the incremental remapping on the B-grid is a problem with the C-grid. We could workaround this by introducing a slip condition along the land boundaries for velocity.

[eclare108213]

It looks like a null-space solution is being energized. The checkerboard pattern is in the null space for some operators on B grids, and C grids probably have similar issues. Upwind advection is very diffusive, enough to keep null-space solutions from growing, but incremental remapping is much less diffusive. Initially, you could try averaging the C-grid velocities to the B-grid corners and running the remapping in the standard way, without the C-grid changes we made. That defeats the purpose of developing the C-grid discretization, but it's just a test. Allowing a slip condition with remapping could help, or just use upwind for the C-grid, for now. The best solution to this problem might be to rewrite incremental remapping for a C grid. That's been started but needs a funding source.

[JFLemieux73]

Thanks both of you for your inputs. Elizabeth, is the null-space issue related to the C grid (momentum) discretization, the remapping (using C grid velocities interpolated to the U points) or a combination of both?

I don't understand your sentence "Initially, you could try averaging the C-grid velocities to the B-grid corners and running the remapping in the standard way, without the C-grid changes we made. "...isn't it what we are doing right now?

[JFLemieux73]

Maybe a solution would be to apply upwind only close to land and remapping for the rest of the domain...

[eclare108213]

is the null-space issue related to the C grid (momentum) discretization, the remapping (using C grid velocities interpolated to the U points) or a combination of both?

It's probably a combination, but note that the remapping is not using U velocities - see below.

It's definitely related to the C-grid discretization -- avoiding the B-grid checkerboard instability is the main reason I do 4 calculations for each grid cell (ne, se, sw, nw) in EVP. I also think it's important to ensure that energy is dissipated in the discretized equations in the same manner as in the continuous equations, which the 4 calculations do, when combined for F1 and F2; I'm not sure whether your discretization guarantees that, but others have been successful with it so there's hope.

"Initially, you could try averaging the C-grid velocities to the B-grid corners and running the remapping in the standard way, without the C-grid changes we made. "...isn't it what we are doing right now?

No, incremental remapping needs the edge velocities, so for the B grid, the corner velocities are interpolated to the edges. For the C grid, we took that out and use the edge velocities directly. Not doing that (i.e. reverting to the B-grid process) would add some diffusion into the solution associated with the averaging and also probably change the behavior near land, due to the masking.

[JFLemieux73]

I am confused. I think remapping uses the U velocities directly. See section 2.4 in the doc:

"as the remapping scheme was originally developed for B grid applications, uvel and vvel are directly used for the advection."

I could check the code but my impression is that for the C grid we interpolate uvelE and vvelN to the U point in order to use the remapping.

[eclare108213]

Look at the code here:

CICE/cicecore/cicedynB/dynamics/ice_transport_driver.F90

Line 544 in 9808b51

if (grid_ice == 'CD' .or. grid_ice == 'C') then

and

CICE/cicecore/cicedynB/dynamics/ice_transport_remap.F90

Line 667 in 9808b51

if (grid_ice == 'CD' .or. grid_ice == 'C') then ! velocities are already on the center

I think the doc needs to be corrected by adding E, N, U, etc to uvel and vvel, to not be so confusing.

[eclare108213]

I looked at the doc, and I don't think it's right. But maybe I'm wrong. @TillRasmussen please straighten us out on this!

Quote from https://cice-consortium-cice.readthedocs.io/en/main/science_guide/sg_horiztrans.html

Two transport schemes are available: upwind and the incremental remapping scheme of [10] as modified for sea ice by [37]. The upwind scheme is naturally suited for a C grid discretization. As such, the C grid velocity components (i.e. uvelE=u at the E point and vvelN=v at the N point) are directly passed to the upwind transport scheme. On the other hand, if the B grid is used, uvel and vvel (respectively u and v at the U point) are interpolated to the E and N points such that the upwind advection can be performed. Conversely, as the remapping scheme was originally developed for B grid applications, uvel and vvel are directly used for the advection. If the remapping scheme is used for the C grid, uvelE and vvelN are first interpolated to the U points before performing the advection.

This is what I think it should say, starting with the 'Conversely' sentence above:

... Because the remapping scheme was originally developed for B grid applications, uvelU and vvelU are are interpolated to the E and N points internally in the advection code, when running on B grids. If the remapping scheme is used for the C grid, then this interpolation step is skipped and uvelE and vvelN are used directly.

[TillRasmussen]

For the incremental remapping then there are two sides of the story. For the departure_points the corner points (uvel and vvel) are used. In order to calculate the edgearea uvelE and uvelN are used. Edge areas are currently not calculated as l_fixed_area = .false. in transport_driver

For consistency with the rest of the code uvel and vvel should probably be renamed to uvelU and vvelU

[JFLemieux73]

Thanks Till.

Elizabeth, you wrote above: "I'm not sure whether your discretization guarantees that, but others have been successful with it so there's hope."

Looking at the overleaf document again (see section 5), our implementation is not exactly what others have done. Again, this is due to the fact that we started with a CD grid. Maybe we should test exactly what Kimmritz had? Basically, we would need a different approach for calculating Ds at the T point.

[eclare108213]

Well, I'm also not sure that theirs guarantees it, or that yours is somehow worse than theirs for any reason, or on the other hand that either is bad. We would need to do some actual numerical analysis on the discretization. But if it's straightforward to test Kimmritz's approach, then you could try that. What do the others use for advection?

[JFLemieux73]

Ii is not written what they used in Kimmritz et al 2016.

I wrote an email to Martin Losch yesterday showing him the checkerboard pattern. He replied:

Hi JF,

hard to say. I learned that a C-grid is accurate for (flux) divergence operations, but can have some noise in the staggered velocities.

Whenever I find noise in scalar fields like concentration, it’s related to noise in the velocities, but the advection schemes should not generate this. Then again, I have no experience with the remapping advection scheme (something I should try at some point), and I usally use “stable” advection schemes that have a little bit of diffusion in them. 1st order upwind is extreme, but I usually use a 3rd order direct space time DST method or a second order scheme with flux limiting (to avoid overshoots).

Bottom line: In my experience noise is almost always related to the dynamics solver (on a c-grid).

[JFLemieux73]

What I wrote above is wrong...we already calculate DsT^2 as done by Kimmritz 2016.

comment in stressC_T:

!-----------------------------------------------------------------
! Square of shear strain rate at T obtained from interpolation of
! U point values (Bouillon et al., 2013, Kimmritz et al., 2016
!-----------------------------------------------------------------

I will update the overleaf document.

[JFLemieux73]

Is it possible Kimmritz et al. had the same issue but did not see it because their advection scheme was more diffusive (like upwind) than remapping?

Or the problem comes from our formulation of boundary conditions for corners (maybe different than in Kimmritz et al....not clear what they did)?

[JFLemieux73]

I have looked at the C grid divu field (at different times) when using remapping. It is indeed noisy. If the noise comes from the solution of the momentum equation, I guess I should be able to see it in instantaneous divu fields even using upwind. I have looked at divu at many different times and don't see any noise. The C and B grid divu fields are very similar when using upwind (not shown).

[eclare108213]

Does it only happen with corners? i.e. if you set up a box case with a wall on one side and periodic in the other direction, would it show the checkerboarding?

[JFLemieux73]

Good idea. I will test it.

[eclare108213]

You could try with the forcing perpendicular to the wall and also at an angle, e.g. northeastward if the wall is on the east side.

[JFLemieux73]

Here is a channel test case with winds blowing north. The north and south boundaries have walls and cyclic condition is used along the x axis (west-east). We don't see the checkerboard pattern (2005-01-15) with C grid and remap:

[JFLemieux73]

As suggested by Elizabeth, I kept the same test case but used atm_data_type = 'uniform_northeast' instead. The conclusion is the same: there is no checkboard pattern (not shown).

[JFLemieux73]

Hum...I found another case where we see the checkerboard. I set atm_data_type = 'uniform_east' and ice_data_type = 'eastblock'. This is aice after 15 days with the C grid and remap:

and with upwind:

[JFLemieux73]

Tony I think we have seen this before...I thought we had fixed that problem but I guess I was wrong.

[JFLemieux73]

I think with avg_zeta the only thing I need to change for the BCs is dvdx in strain_rates_U. dudy is zero at the wall (eastern wall). The U point at the wall is (i,j). dvdx at Uij can be found with Taylor series expansions. It is expressed as:

dvdx ~ vvelN(i-1,j) / (3 dx) - 3 vvelN(i,j) / dx

Indeed the stencil is larger with this than with our current (ghost cell) method.

[JFLemieux73]

Ok I have added this in strain_rates_U after the calculation of shearU:

if (npm(i,j) == 1 .and. npm(i+1,j) == 0) then
shearU(i,j) = dyU(i,j) * (vvelN(i-1,j) / 3d0 - 3d0*vvelN(i,j))
endif

I am not available for the rest of the day...we will know tomorrow morning if this works.

[JFLemieux73]

It is disappointing but the new BC does not change anything; the checkerboard is still there (not shown).

[JFLemieux73]

Elizabeth, I am a bit confused about what I should try next...should I try the ice edge velocities 2)?

[dabail10]

Did you already change the subcycles in the EVP? I was looking back at the paper led by Bill Lipscomb:

https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2005JC003355

This is more about high resolution ... However, is there a B/C grid issue with the ridging?

[JFLemieux73]

Thanks Dave. I don't think it is an EVP issue. I use elasticDamp = 0.1d0 and ndte=2400 (10 times bigger than default). These values of elasticDamp and ndte should increase the numerical convergence. I know that paper led by Bill. This problem should be less of an issue with kstrength=0 (Hibler...this is what I use) and as dt is reduced. I use 60 min right now and I think I tested 15 min before. Maybe I could give it a try with 5 min. Stay tuned.

[eclare108213]

I think implementing the ice edge velocities would be interesting and useful for other purposes, but I do not think that it will help this issue. I want to dig back into my evp notes to see what the discretization might look like in that approach (bigger stencil). I've rederived it since the original derivation, so I should be able to reproduce the steps needed, I'm just not sure what it would look like for a C grid.

[JFLemieux73]

If we want to implement another discretization we should consider what Bill proposed;

#714

This could be our next coding camp! 😎

[JFLemieux73]

Dave,
I ran it with a dt=120 s...the checkerboard is still there.

[JFLemieux73]

I read again Till's comment above about l_fixed_area = .false. in the transport_driver. So I decided to set l_fixed_area = .true. I am not exactly sure what this does and if I am allowed to use it but interestingly when I set it to true the checkerboard disappears for both test cases as shown below:

Any comments about this?

[eclare108213]

Wow! I'm not familiar with how l_fixed_area is supposed to work, not having ever used it, but my impression is that you need to externally calculate some things to send into the transport routine in addition to changing the flag. So this might not be actually transporting anything? @whlipscomb might be able to shed some light on what's needed here.

[dabail10]

Are you sure you had incremental remapping on?

[JFLemieux73]

Yes. advection = 'remap'

[dabail10]

Cool.

[JFLemieux73]

Here are some comments taken from ice_transport_remap.F90:

! Finally, a few words about the edgearea fields:
!
! In earlier versions of this scheme, the divergence of the velocity
! field implied by the remapping was, in general, different from the
! value of del*u computed in the dynamics. For energetic consistency
! (in CICE as well as in layered ocean models such as HYPOP),
! these two values should agree. This can be ensured by setting
! l_fixed_area = T and specifying the area transported across each grid
! cell edge in the arrays edgearea_e and edgearea_n. The departure
! regions are then tweaked, following an idea by Mats Bentsen, such
! that they have the desired area. If l_fixed_area = F, these regions
! are not tweaked, and the edgearea arrays are output variables.

[TillRasmussen]

Would it be worthwhile to reach out to Mats Bentsen?

…

On Thu, Dec 22, 2022 at 7:34 PM JFLemieux73 ***@***.***> wrote: Here are some comments taken from ice_transport_remap.F90: ! Finally, a few words about the edgearea fields: ! ! In earlier versions of this scheme, the divergence of the velocity ! field implied by the remapping was, in general, different from the ! value of del*u computed in the dynamics. For energetic consistency ! (in CICE as well as in layered ocean models such as HYPOP), ! these two values should agree. This can be ensured by setting ! l_fixed_area = T and specifying the area transported across each grid ! cell edge in the arrays edgearea_e and edgearea_n. The departure ! regions are then tweaked, following an idea by Mats Bentsen, such ! that they have the desired area. If l_fixed_area = F, these regions ! are not tweaked, and the edgearea arrays are output variables. — Reply to this email directly, view it on GitHub <#792 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH7N7DQRGIVIHIDOGSFIKZLWOSNMDANCNFSM6AAAAAASBKJUVQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[JFLemieux73]

We could start with Bill's opinion and then maybe we could contact Mats. In the meantime I am trying to better understand the remapping...

[eclare108213]

I am working up a discretization using the variational approach, as in the current B-grid EVP. I'm not sure it will work, since C-grid doesn't have enough velocity nodes to use a fully bilinear approximation -- this means that the discretization will be lower order than in a bilinear case. It might be less likely to checkerboard, but there's no guarantee that it won't, and we also might get some insight into how the boundary conditions play in.

[JFLemieux73]

Happy new year! I still have to investigate but results above suggest that the checkerboard comes from the remapping not the C grid discretization...do we need the variational approach?

[eclare108213]

I think it's important, but I'm willing to entertain other ideas. My understanding is that the advection scheme is not creating the checkerboard so much as allowing it, because a feature of incremental remapping is its very low diffusion -- a good thing as far as advection goes. I interpret your results above as advection schemes (or maybe energetic inconsistencies) in the other codes covering up the problems, rather than fixing the source. On the B-grid, the variational approach eventually led to a discretization that did not checkerboard while using the same advection scheme being used here, and it's not a discretization that could simply be written down based on the form of the continuous equations. In deriving the current B-grid EVP scheme, I went through several iterations that did checkerboard, even while using the variational approach for the derivation -- those schemes were much simpler (lower order). The key was to use a higher order basis function (bilinear) for the velocity AND stresses, which is what produces the 4 values per grid cell (nw, ne, sw, se) for all of the components.

The problem I'm having with the C-grid is that it does not have enough degrees of freedom (which Sergei also pointed out, I think). I can use a bilinear basis function for the stresses but I haven't found a basis that works for velocities located on the edge mid-points. The discretization I just worked out uses a bilinear stress approximation (on corners) but only linear for the velocities (i.e. assuming u varies in the x direction and is constant in y, similarly for v). I just reached the end of the derivation and it is not correct for sigma12, which needs du/dy etc. (I should have realized that earlier -- lesson learned). I could write down what the sigma12 discretization would be based on the forms of sigma1 and sigma2, but then it wouldn't be consistent with the divergence, tension and stress terms. (Unless maybe we somehow change the boundary conditions?) The variational approach ensures that all of these things are consistent with each other, as they must be to guarantee that the discretized equations dissipate energy in the same manner as the continuous equations. That's why I think it's important.

[whlipscomb]

@JFLemieux73, I'm having to revisit my long-ago youth. But fortunately, I wrote some comments on l_fixed_area in ice_transport_remap.F90.

According to those comments, this tweak ensures that the area fluxed across each cell edge is equal to the area implied by del*u in the dynamics. On a C-grid, u is simply the cell-edge velocity, instead of an average from neighboring corners as on the B-grid.

It would be easier to show with a picture, but I'll try in words. The idea is to tweak the size of the departure region without too much changing its basic shape. Try drawing a north cell edge with a departure region above. The departure region is a quadrilateral defined by the two edge vertices and the departure points of the two back trajectories. The velocity in this case would be mainly in the negative-y direction. For simplicity, assume that both departure points lie inside the cell above the edge (i.e., they don't lie in the adjacent cells to the left and right).

On the segment connecting the two departure points, mark the midpoint. From the midpoint, draw a segment to each edge vertex. This divides the departure region into three triangles.

Now slide the midpoint up a little bit. This slightly increases the area of the three triangles. If you slide the midpoint down, the three triangles get smaller.

The way the code works is that it figures out how much to move the midpoint. Otherwise, the midpoint stays put, and the logic is the same as for l_fixed_area = F.

It doesn't surprise me that setting l_fixed_area = T suppresses noise, but I'm not sure I'd have predicted it either. It's been too long since I thought about it.

But if it works, I'd say to go ahead and use this setting with the C grid. At some later time, we could write a simpler remapping scheme that's more natural on a C-grid. But I'm guessing that for now, you're mainly looking for a fix that works.

Please let me know if this description isn't clear.

[JFLemieux73]

The C grid noise has been eliminated with l_fixed_area=T (see issues #809 and #811). The code does not crash anymore for gx simulations (see PR #833).