uvel instability with CD grid

[JFLemieux73]

No description provided.

[JFLemieux73]

I use gx3 (with kstrength=0, Pstar=2.75e4). The large scale patterns of uvel,vvel are similar to the ones on the B grid (good!). Close to land however there could be large values of u,v (e.g. 3 m/s) flipping sign from one cell to the next.
-is it due to the fact that we have not added diffusion yet (as proposed by Carolin)?
-is it due to boundary conditions?

i=75,j=109 is one grid cell with the problem. The mask in the region is

0 0 0
1 1 0
1 1 0
1 0 0

The lowest left T cell is 75,108.

[JFLemieux73]

Comparing the code with the design document, there is a mistake in strain_rates_U
shearU = dxU(i,j) * ( uEijp1 - uEij ) &
- uvelU(i,j) * uvm(i,j) * ( dxE(i,j+1) - dxE(i,j) ) &
+ dyU(i,j) * ( vNip1j - vNij ) &
+ vvelU(i,j) * uvm(i,j) * ( dyN(i+1,j) - dyN(i,j) )

The last line should start with a minus sign.

[JFLemieux73]

I should also modify div_stress to use earear and narear.

[eclare108213]

Weren't there some 2-point stencils (C-grid) used that could instead be 4-point stencils (CD grid)? Or maybe those were only for shifting/averaging grid quantities?

[JFLemieux73]

uvelU above is one example. It could be calc from uvelN and uvelE around the U point. I am not sure how it is calculated for the moment. @dabail10 can you answer this question?

[dabail10]

Aha. It looks like uvelU,vvelU are only computed at the end of the evp subroutine. There should be another call at the beginning of the substep. Looks like that fell through the cracks. I will commit something today.

[JFLemieux73]

@dabail10 I think your last PR (uvelU, vvelU) corrected the high velocity problem. I just did a quick gx3 run and the velocities seem to be ok. Good job!!!

I wrote too quickly...there are still high values of velocities (e.g. 3 m/s)...need to investigate more.

[JFLemieux73]

The minus sign is corrected in the code:
apcraig#45

[dabail10]

The boxisland test is still blowing up as well. I will do some digging to see if I can find the points where it is blowing up.

Dave

[JFLemieux73]

Going back to the gx3 test case. uvelE: -0.41 to 0.46 (ok?) vvelE: -1.1 to 5.6.

1. capping=0 (like implicit solver) instead of 1: same issue.
2. purely viscous (tmpcalc = strength/tinyarea in viscous_coeffs_and_rep_pressure_Tviscous_coeffs_and_rep_pressure_T): blows up...this is interesting.

[JFLemieux73]

3. replacement = false: uvelN: -5.9 to 6
4. replacement = .false., viscous = .true., blows up
5. coriolis = zero, replacement = .false., viscous = .true., blows up
6. kstrength =0, coriolis = zero, replacement = .false., viscous = .true., blows up

I will use the settings of 6 for debugging (simpler problem).

[JFLemieux73]

Just to confirm that it runs fine (with realistic velocities) if I use the settings of 6 but with Pstar=0. I go back to the standard Pstar=2.75e4.

[JFLemieux73]

Some obs:

It blows up quickly during the first call of evp.

ksub=1 : u=v=0 everywhere (initialization)
ksub=2: |u|,|v| < 1 m/s everywhere
ksub=3: |u| or |v| larger than 5 m/s at some locations.

Is it possible the model struggles with an initialization at rest and wind blowing directly at full strength?

[dabail10]

Try adding the following to ice_restart_driver.F90 so that the E and N velocities are initialized from the B grid on a restart:

call grid_average_X2Y('U2ES',uvel,uvelE)
call grid_average_X2Y('U2ES',vvel,vvelE)
call grid_average_X2Y('U2NS',uvel,uvelN)
call grid_average_X2Y('U2NS',vvel,vvelN)

[JFLemieux73]

Thanks @dabail10. In your box test have you tried with a ramping up of the wind or it was directly at full strength?

[dabail10]

I have not tried ramping up the winds. I could do that easily. Or reduce the wind speed.

[JFLemieux73]

I think I should use your box test instead. How do I set it up?

[dabail10]

I reduced the wind speed to 1 m/s eastward and this still blows up. It takes more substeps though. If you chose -s boxislandse that is what I am using. I'm noticing that all of the stress terms are still zero at the point I am looking at near the western boundary. I did close the boundaries, but this doesn't seem to help. I suspect my issue is related to land/ocean mask issues. There are a number of box options in configuration/scripts/options.

[dabail10]

By the way, this is not my test. Elizabeth and David Hebert set this up.

[JFLemieux73]

7. kstrength =0, coriolis = zero, replacement = .false., viscous = .true., specialBC = .true., blows up

I have coded a simpler BC...it still blows up.

[dabail10]

So, I set ndte to 3 and histfreq to '1' so I could at least get a history file. Here is uvelN from the box domain. I have set the colorbar limits to -10 to 10 m/s, but there are values here of -60 m/s. This is definitely related to land/ocean masking. The wind is uniform_east at 10 cm/s. The sig1, sig2, and sigP look very similar to this plot. This is kind of a chicken or egg though. I am not sure what blows up first.

ncview.uvelN.pdf

[dabail10]

I did another run with ndte = 1. Everything looks ok at step 1. By step 2, the velocities have already gone over 10 m/s and we are getting +/- values near land. The sig1, sig2, and sigP are looking pretty big. Values in each are a max of 0.48, 0.61, and 37920 respectively. I think the sigP values are unreasonable? The divu and shear values look ok to me. For the CD case we send stresspT, stressmT and stress12T to the subroutine principle_stress. sigP == stresspT. So, I will focus on the disctretization of stresspT and stresspU for now.

[JFLemieux73]

Do you keep the same evp time step when you change ndte like this?

sig1 and sig2 are normalized. They should be between 0 and -1. Max values of sigP scale with the ice strength. If you set kstrength=0, the max value at a point should be Pstar*h. h=mean thickness in a grid cell. With Pstar=27500, a value of 37920 is realistic.

[dabail10]

Good points. I should try reducing the timestep as well.

[JFLemieux73]

Hum...maybe my simplified test is not great. I tested the following for the B-grid

kstrength =0, coriolis = zero, replacement = .false., viscous = .true. and it blows up...
kstrength =0, coriolis = zero, replacement = .false. works...it is viscous=.true. that makes it blow up.

[dabail10]

Just reran with dt = 30s. So, here is something strange. This plot is vvelE on the very first step. The blue values are -8.66e-10. So, pretty much everywhere there is a southward velocity. The green values are identically 0. So, where we have land on the east, V = 0. However, the one gridpoint E-W channels have vvelE values of -8.66e-10. I would certainly expect a uvelE here, but does it make sense to have a vvelE?

ncview.vvelE.pdf

[dabail10]

I think I found another issue. In stress_T we compute zetax2T and etax2T. These are then averaged from adjoining T-cells to compute zetax2U and etax2U. If I'm not mistaken, I think this requires a halo update? Also, stresspT, stressmT, and stress12T (plus U counterparts) are used to compute strintxE/N, and strintyE/N. These also use i+1, i-1, j+1, j-1 values and should also have halo updates? I think everything else is updated in the halos ... but I might be missing other things.

[apcraig]

As noted before, if you are calling grid_average_X2Y, the input fields should have valid halos values before the average is called. We should definitely review all the fields where averaging is called.

[dabail10]

We don't use grid_average for these quantities.

[apcraig]

I see, but the halo's are used in parts of the dynamics calculations then. Obviously, the halos need to be valid before used. Great!

[JFLemieux73]

Going back to the boxsyme test. I want to make sure I get a converged solution with the B grid. I have set Pstar=2.75e3 (10 times smaller than default) and ndte=3600. I look at the fields after one time step (after one hour).

Here is the Bgrid uvel solution:
uvelB.pdf

The uvelE CD solution looks ok but is noisy:
uvelE_CD.pdf

However the uvelN CD solution is quite different than the ref (B grid) solution:
uvelN_CD.pdf

[JFLemieux73]

Let's assume that the CD grid requires a stabilization scheme (which we don't have) but that a regular C grid would be ok. I have coded a quick and dirty C grid for the boxsyme test. After each subcycle, uvelN and vvelE are calc from their neighbors:
uvelN(i,j)=( uvelE(i,j)*epm(i,j) + uvelE(i-1,j)*epm(i-1,j) + &
uvelE(i-1,j+1)*epm(i-1,j+1) + uvelE(i,j+1)*epm(i,j+1) )/4d0
Not sure if BCs are ok. Nevertheless, here is the uvelE field I get with this (called uvel_C):
uvelE_C.pdf

uvelN is very similar (not shown). It also works for the vvel field (not shown).

[dabail10]

What a great idea! I did a quick C-grid implementation as you have and the kmtislands test remains stable for 5 days! There are still some weird speeds in excess of 300 m/s, but otherwise it keeps going. So, it is definitely related to the tangential velocity components and perhaps the stabilization is the key.

[JFLemieux73]

Thanks Dave. I also tested what I call a rotated C grid. uvelN and vvelE are calculated the regular way while uvelE and vvelN are obtained from their neighbors:
uvelN_rotC.pdf

This is noisier than what I got for the regular C grid...The noise is even more obvious for vvelE:
vvelE_rotC.pdf

This suggest that the code is ok for uvelE and vvelN but that something is wrong for uvelN and/or vvelE.

[JFLemieux73]

uvelE and vvelN are calculated from stresspT, stressmT and stress12U

while

uvelN and vvelE are calculated from stresspU, stressmU and stress12T

Is there something wrong in the calculation of stresspU and/or stressmU and/or stress12T? I would suspect stresspU and/or stressmU. Maybe we could do a test by setting stress12=0 (e.g. e_ratio very large)?

[dabail10]

While the velocities are not as large, there is still an issue with uvelE and vvelN in the hacked C grid. So, I think there is something going on with all of these. I did not see any issues with the code compared to the document, but I could have missed something.

[apcraig]

The following two cases are not bit-for-bit on the cgridDEV branch,

smoke gbox80 1x1x80x80x1 boxsyme,diag1,gridcd,reprosum
smoke gbox80 1x1x40x40x4 boxsyme,diag1,gridcd,reprosum

The default is run with ndte=240. I'll set to 1 or 2 and then try to debug that a bit. Do I need it to be set to at least 2? I can't remember what was recommended. Thanks.

[JFLemieux73]

@apcraig I would like to set up similar test cases as boxsym but with a block of ice in the middle instead of ice everywhere. With the following code in ice_init.F90:

     else if (trim(ice_data_type) == 'uniform') then

        icells = 0
        do j = jlo+celloffset, jhi-celloffset
    do i = ilo+celloffset, ihi-celloffset
       if (tmask(i,j)) then
          icells = icells + 1
              indxi(icells) = i
          indxj(icells) = j
           endif
        enddo
        enddo

it could work for the regular boxsym (celloffset=0) or for the new test (e.g. celloffset=10). Is it ok? Should I define a new variable (celloffset) in the namelist? Do you have a better name for celloffset? Thanks!

[apcraig]

I'm trying the same sort of test. What I'm doing is creating a new ice_data_type called "smallblock". That is going to create a 2x2 block of ice in the middle on the domain. I'm also adding a 12x12 grid, gbox12. I can PR those mods soon if that works for you. If you prefer your method, we could do that too, but I think if we want something permanent, we should have a new ice_data_type. Also, no new namelist would be required with my approach. I'm testing it now.

[apcraig]

See apcraig#51

Following cases are not bit-for-bit on first timestep. Now debugging

smoke gbox12 1x1x12x12x1 boxsyme,diag1,gridcd,reprosum
smoke gbox12 1x1x6x6x4 boxsyme,diag1,gridcd,reprosum

with manual changes in ice_in for ice_data_type='smallblock', ndte=2

[apcraig]

If I merge @dabail10 dynhalo PR apcraig#49 to my sandbox, then the above tests are bit-for-bit. @dabail10, we should merge your PR, is there any reason we shouldn't do that? It's a draft PR right now.

[dabail10]

Excellent. No reason. I just wanted to make sure this was fine.

[apcraig]

With @dabail10's PR, the following tests are bit-for-bit,

smoke gbox12 1x1x12x12x1 boxsyme,diag1,gridcd,reprosum
smoke gbox12 1x1x6x6x4 boxsyme,diag1,gridcd,reprosum

and

smoke gbox80 1x1x80x80x1 boxsyme,diag1,gridcd,reprosum
smoke gbox80 1x1x40x40x4 boxsyme,diag1,gridcd,reprosum

Doesn't mean everything is OK, but we know this fixes a problem that may have affected the solution.

Separately, does it make sense to add some stupidly simple smoothing on the D velocities? @JFLemieux73 showed that averaging the 4 gridcells to override the D velocities smooths out the solutions nicely. Should we try to do something like

vel_D = vel_D_computed * (1.-eta) + vel_D_average * (eta)

again, just to see how much eta we need to keep things stable. We know 0 is too small, 1 is big enough. Does it work with 0.01, .1, .5, 0.95, 0.99? Just an idea. This would not be a long term stabilizer, but might provide an ability to at least evaluate the instability as well as the D solution a bit.

[apcraig]

I will merge apcraig#49, apcraig#50, and apcraig#51 once @JFLemieux73 approves. @JFLemieux73, please review apcraig#49 and apcraig#51 when you get a chance.

[JFLemieux73]

The smoothing of the D velocities helps to get a solution closer to the B grid one but I agree with @dabail10 that there is still something wrong. These new tests will be really helpful for debugging.

[apcraig]

I just ran a handful of tests with gbox12, 1x1x12x12x1 (one block), boxsyme, gridcd.

smallblock + ndte=2 or ndte=240 produces totally stable and uninteresting solution for 5 days.
uniform + dt=10 + ndte=2 produces what seem to be crazy solutions on the first timestep.

This looks to me like a boundary condition issue. When I set dt=10 and ndte=2, should I expect a stable solution if things are working fine? The subcycling is just decreasing the timestep? It's not an algorithm that needs to converge to a stable solution? dt=10 + ndte=2 should produce the same results as dt=240 and ndte=480, more or less?

[apcraig]

Is the boundary condition no normal and no slip? That's what the velocities suggest on the boundary. So maybe those are OK, but something is happening to the interior solution.

[JFLemieux73]

Hi Tony,

if you set ndte=2 then elastic waves do not damp and the evp time step is too large. I would not expect a realistic solution.

When I want to compare the B and CD grid solutions after a few subcycles (to investigate when they start to diverge), what I do is that I set npt=1 and I keep ndte=240. Then I change the code in the evp:

Let's say I want to look at the solutions after 2 subcycles I do:

do ksub=1,2 (instead of 1,ndte

Yes the BCs are no slip no normal flow.

[apcraig]

Thanks for the clarification on ndte @JFLemieux73.

I setup another case which I think is even more enlightening. If I change ew_boundary to 'cyclic', I get a channel. This should be an easier solution to debug. That produces a really interesting and wrong result. So that's gbox12, 1x1x12x12x1, boxsyme, gridcd then set ew_boundary=cyclic. close_boundaries is ignored in this case. A good thing is that the solution is identical at each i location (as it should be), so it should be easier to debug. I ran with dt=3600 and ndte=240.

[apcraig]

I attach some output from a 1 hour run with some output added.
cice.runlog.220125-130443.txt
This writes out the 4 velocities at i=6 across the channel at each subcycle. You can see how over the first 10 subcycles, the v velocity slowly comes into the solution, the N and E solutions eventually decouple from each other, then the v velocity grows rapidly in vvelE. It's pretty clear from all the runs that the N and E solutions are too independent. Are we doing 2 point averaging when we should be doing 4 pt averaging for some of the dynamics variables (or something like that)?

--- a/cicecore/cicedynB/dynamics/ice_dyn_evp.F90
+++ b/cicecore/cicedynB/dynamics/ice_dyn_evp.F90
@@ -883,6 +883,11 @@ subroutine evp (dt)
                uvel(:,:,:) = uvel(:,:,:)*uvm(:,:,:)
                vvel(:,:,:) = vvel(:,:,:)*uvm(:,:,:)
 
+               i = 6
+               do j = 2,ny_block-1
+                  write(nu_diag,*) 'tcx1E',ksub,j,uvelE(i,j,1),vvelE(i,j,1)
+                  write(nu_diag,*) 'tcx1N',ksub,j,uvelN(i,j,1),vvelN(i,j,1)
+               enddo
             end select

[JFLemieux73]

I dont think the issue is the two point averaging versus four. I think all these terms disappear when we run on a cartesian grid with dx=dy=constant (like we do in the box tests). Look at the design document. For example in the section strain rates at the U point, divUij = ...uNu(dyNi+1j - dyNij). uNu is u from the N points interpolated to the U point. But it should not contribute as dyNi+1j = dyNij.

But I agree with you that we get two solutions with the CD grid...

[JFLemieux73]

Closing this issue: see C-CD TEST issues for follow up.