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: https://github.com/CICE-Consortium/CICE/blob/9808b5195b2af77d2b70f5fec3e417b0a6020c5f/cicecore/cicedynB/dynamics/ice_transport_driver.F90#L544 and https://github.com/CICE-Consortium/CICE/blob/9808b5195b2af77d2b70f5fec3e417b0a6020c5f/cicecore/cicedynB/dynamics/ice_transport_remap.F90#L667

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 also tested with avg_strength instead of avg_zeta. It does not fix the problem.

[JFLemieux73]

The only available option to test the corner problem is uniform_northeast. I modified the code to test the other corners with uniform "45 degree" winds. The problem is the same in the four corners (now shown).

[JFLemieux73]

For the other problem there is only eastblock that is available. I will come back later to this problem and focus on the corner problem.

[JFLemieux73]

I confirm what Till wrote above. uvelE and vvelN are not directly used as l_fixed_area = .false. If the remap is called as for the B-grid (uvelE and vvelN are not sent) the answer is the same (BFB) with obviously the same checkerboard pattern.

[eclare108213]

Reviewing and catching up on this discussion. Some thoughts:

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).

This is consistent with what Martin said, "that a C-grid is accurate for (flux) divergence operations, but can have some noise in the staggered velocities."

In your eastblock case above, presumably the only difference from the cyclic channel with the wind blowing north, other than the direction, is that there is an ice edge along the wall, i.e. a corner formed from the ice and the wall itself. Is that correct? I wonder if the CD grid velocities, with appropriate averaging, would be enough to damp the checkerboard instability.

The CICE documentation for incremental remapping includes this intriguing statement:

We made one other change in the scheme of [10] for locating triangles. In their paper, departure points are defined by projecting cell corner velocities directly backward. ... This approximation is only first-order accurate. Accuracy can be improved by estimating the velocity at the midpoint of the trajectory.

It would be fairly complicated to implement, but could using those velocities interpolated directly from the calculated uvelE, vvelN (without first interpolating them to the corners) bypass the checkerboard pattern associated with the corner velocities? If the issue stems from the momentum calculation, as Martin asserts, then I suspect this would not help. The increased accuracy might even make it worse.

Here's another thought: what if you did allow a small amount of free slip along coastlines, as @dabail10 suggested? ~puny. I guess it could be implemented by changing the inland value (currently set so that velocity is strictly zero at the coast) slightly, so that the effective coastline moves just a little. Maybe that's a crazy idea! Just throwing things out there...

[JFLemieux73]

Thanks Elizabeth. I am glad you are throwing ideas as I am a bit confused about what I should do next! :)

"In your eastblock case above, presumably the only difference from the cyclic channel with the wind blowing north, other than the direction, is that there is an ice edge along the wall, i.e. a corner formed from the ice and the wall itself. Is that correct?"

yes it is correct.

"I wonder if the CD grid velocities, with appropriate averaging, would be enough to damp the checkerboard instability."

I am not ready to go back to the CD grid yet...but let's keep that in mind.

"Here's another thought: what if you did allow a small amount of free slip along coastlines, as @dabail10 suggested? ~puny. I guess it could be implemented by changing the inland value (currently set so that velocity is strictly zero at the coast) slightly, so that the effective coastline moves just a little. Maybe that's a crazy idea! Just throwing things out there..."

Yes I could try something like that. I would like also to try to set u and/or v to zero close to land (kind of extend the BCs).

[JFLemieux73]

Another thing I would like to try is to use upwind close to land and remap elsewhere...

[JFLemieux73]

Finally for some tests the eastblock is better. The last ocean T cell on the right is T(79,j). In stepv_C, I added:

if (i .eq. 79) then vvel(i,j)=0.0 endif

This kind of extends the no slip condition. The checkerboard is still there (not shown).

Commenting the above and now adding in stepu_C:

if (i .eq. 78) then uvel(i,j)=0.0 endif

This kind of extends the no outflow condition. The checkerboard is still there (not shown).

Including both conditions at the same time (extend no slip, no outflow conditions) also leads to the checkerboard pattern...this suggests that the problem is not related to the boundary conditions.

[eclare108213]

This paper has some potentially useful background info and references. The POP manual for CESM also mentions the checkerboard (p 19).

[JFLemieux73]

Thanks Elizabeth. I will read it. I also wanted to test upwind at the beginning and then switch to remap. It does not help...the checkerboard develops when remap is switched on. I tested 2 or 5 days of upwind (remap after) with the eastblock test. Results after 15 days show the checkerboard...but less pronounced close to the wall (not shown).

[eclare108213]

See also Giacomo's paper on the CD approach he presented to us early on. He doesn't mention checkerboards here but does say that this method doesn't need stabilization as in other CD formulations for unstructured meshes. I'll ask him whether he encountered a checkerboard null space solution in any test cases on regular meshes. NB this paper does not show any results with coastlines. Implementing this CD approach would be a large effort, but some of his analysis methods might be helpful for understanding the current discretization.

Giacomo Capodaglio, Mark R. Petersen, Adrian K. Turner, Andrew F. Roberts, An unstructured CD-grid variational formulation for sea ice dynamics, Journal of Computational Physics, Volume 473, 2023, 111742, ISSN 0021-9991, https://doi.org/10.1016/j.jcp.2022.111742. (https://www.sciencedirect.com/science/article/pii/S0021999122008051) Abstract: We present a new unstructured CD-grid formulation of the elastic-viscous-plastic rheology, where the velocity unknowns are co-located at the edges. Our framework of choice is the Model for Prediction Across Scales (MPAS) within E3SM, the climate model of the U.S. Department of Energy, although our approach is general and could be applied to other models as well. The mesh cells in our analysis have n sides, with n greater than or equal to three. Given an initial mesh, we compute the divergence of the stress using a sub-grid, which introduces two separate sets of cells that may differ in shape or be equal depending on the initial topology of the initial mesh. The proposed approach is of higher order compared to the existing B-grid discretization of MPAS-Seaice, and does not need stabilization as other CD-grid formulations on unstructured grids. Our new methodology will contribute to transitioning the dynamics of MPAS-Seaice to a CD-grid mesh, facilitating improved coupling with MPAS-Ocean and reducing numerical errors associated with interpolation.

[JFLemieux73]

So I configured the McGill (C-grid) model to have similar test cases. I guess I don't have exactly the same parameters (e.g. drag coefficients) but here is a quick test. Here is the corner test case after 15 days with upwind:

As with CICE (C-grid) with the upwind we don't see the checkerboard pattern. Maybe to be expected as upwind is diffusive.

[JFLemieux73]

I recently coded a semi-lagrangian advection scheme in the McGill model. This is similar to the remapping scheme and a lot less diffusive than upwind. Note that close to land (2 cells) the order of the scheme has to be reduced. Close to land upwind is therefore applied. Here is the concentration after 15 days with the semi-lag scheme:

[JFLemieux73]

Apart from slightly different parameters, I have the impression that there are two things that are different between the McGill runs and the CICE runs: JFNK versus and EVP and the treatment of the BCs. I really don't think it is a solver issue (JFNK versus EVP).

Both models have no-slip, no outflow BCs but they are implemented differently: ghost cells for CICE and Taylor series expansions (complicated!!!) for McGill. Even though I wrote above that BCS are not the issue maybe it needs to be reconsidered.

So at this point I suspect that the checkerboard comes from one or the two following: 1) the implementation of BCs 2) the remapping scheme itself when used for the C-grid.

Elizabeth any ideas?

[JFLemieux73]

I have also tested the eastblock. There is no checkerboard in both upwind and semilag runs (not shown).

[JFLemieux73]

I would like to modify CICE to match as much as possible what is done in the McGill model. I though a bit more about the possible differences:

1) BCs (as mentioned above) 2) ice edge velocities: McGill uses freedrift values while CICE-C currently has zero velocities (see issue #713). 3) tauair and tauwater are not multiplied by aice in the McGill model.

I don't think 3) is important here. I think the main differences are 1) and 2). I will modify CICE for the eastblock test. It will be easier than the corner test case to test the BCs of the McGill model.

[eclare108213]

I agree that 3) is likely not important for this issue, but I do think it is wrong in the McGill model - see W. M. Connolley, J. M. Gregory, E. C. Hunke and A. J. McLaren (2004). On the consistent scaling of terms in the sea ice dynamics equation. J. Phys. Oceanogr., 34, 1776–1780.

2) should be fixed, but it's likely not impacting the northeast corner case.

That leaves 1)...

[JFLemieux73]

I agree. Ok I will start with 1)...

[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 - 3d0vvelN(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.