C-CD TEST: box2001

[dabail10]

No description provided.

[dabail10]

./cice.setup --case box2001b --mach cheyenne --env intel -s gridcd,debug,box2001 -p 40x1 --grid gbox80

[dabail10]

Hmm. There is a block issue with the winds in this case.

[dabail10]

I guess this is intended to be this way in box2001_data_atm. I am going to rerun with uniform_east.

[dabail10]

I simplified things a bit by doing:

atm_data_type = 'uniform_east'
ice_data_type = 'uniformp5'

So this means aice = 0.8 (I changed this in the code) and hi = 1.6 to start. Advection is still on in these cases. The strength field is constant at 805.888 in both B and CD. The ice cannot move I guess and so stays constant after 5 days.

Here are some plots of sigP (at gridcell centers):

CD-grid:

B-grid:

The patterns are qualitatively similar, but some interesting differences.

The values at the sides are very similar magnitude. Top =~ 108.4 versus 108.5, Left =~ 56.4 versus 56.3, Right =~ 749.6 versus 749.8, and Bottom =~ 697.5 for both.

The upper right corner is very different: 650. versus 0. The lower right corner is 180 versus 120. The upper left is 33 versus 103. Stresspu is only available for the CD grid:

The upper right corner here is -1545.

The uvelE versus uvel is very close to 0.0642 in the interior of both. uvel (B) is zero all along the north and east boundary. uvelE (CD) is only zero on the eastern boundary. In the lower left corner, uvel (B) is 0.03 m/s and uvelE (C) is 0.0013 m/s.

Going by the Hunke 2001 paper:

1. In most of the domain du/dx =~ 0.
2. On the left boundary du/dx > 0.
3. On the right boundary du / dx < 0.

[dabail10]

Sorry for the confusing things I am writing. False alarm. The stresspU value in the upper right corner is indeed zero. I have updated the earlier comment. The values of stresspU on the eastern boundary are big negative numbers. This is because du/dx < 0. The northern boundary has smaller negative values, but also because du/dx < 0. Maybe I should turn off advection and coriolis?

[dabail10]

Advection is already turned off by kstransport, but I set coriolis to 'zero'. Here are a couple other thoughts and then I am going to sign off.

1. Interpolating the strength for the viscous coefficients is probably not the way to go. For the B-grid we simply use the strength field at the grid cell centers. I think we should do the same here.

2. There is an interesting disconnect between uvelN and uvelE near the boundaries. On the B-grid, uvel/vvel will be zero at the eastern and northern boundaries. uvelE will be zero in the last column of cells against the eastern boundary, but nonzero along the northern boundary. uvelN will be zero in the last row along the northern boundary, but nonzero in the last column along the eastern boundary. So, I am thinking there are implications for the stresspU, stressmU, and stress12U terms. I believe this disconnect will then translate into the interior and cause the differences in the sig1, sig2, and sigP terms through the strain_rates_U calculation. Thinking of this in terms of uNip1j, uNij ...

On the last column in the interior say:

uNip1j = 0.
uNij = 0.003 (for example)
uip1j = uij = 0.
uEip1j = uEij = 0.

Forgetting about the ratios as the gridcells are all the same size from the code:

  uNip1j = uvelN(i+1,j) * npm(i+1,j) &
         +(npm(i,j)-npm(i+1,j)) * npm(i,j)   * ratiodxN(i,j)  * uvelN(i,j)
  uNij   = uvelN(i,j) * npm(i,j) &
         +(npm(i+1,j)-npm(i,j)) * npm(i+1,j) * ratiodxNr(i,j) * uvelN(i+1,j)

Actually uNip1j is equal to uNij here. So, divU is zero. However, don't you want divU to be nonzero? I thought the idea was to have uNip1j = -uNij? Or is only the v component reflected near the boundary?

What about if you go from ice to no ice. All of the npms are 1. Say, there was ice next to no ice. In this case, uNip1j = 0. and uNij = 0.003.

[dabail10]

Out of curiousity I did following:

  uNip1j = uvelN(i+1,j)
  uNij   = uvelN(i,j)
  vEijp1 = vvelE(i,j+1)
  vEij   = vvelE(i,j)

  uEijp1 = uvelE(i,j+1)
  uEij   = uvelE(i,j)
  vNip1j = vvelN(i+1,j)
  vNij   = vvelN(i,j)

and got pretty much the same answers. So, it doesn't appear to be related to the boundaries of the C/CD grid.

[dabail10]

I turned off Coriolis. I'm thinking this is related to the variational approach or not.

B-grid:

CD-grid:

The solutions are quite different for sigP. Maybe still boundary related. I am thinking I will try to code up a CD-grid variational version.

[JFLemieux73]

Thanks Dave. Coding a variational version is a lot of work. Other models do not use this variational approach (e.g. mitgcm) so I am not convinced that's the way to go.

[JFLemieux73]

We are not even sure that the B grid solution is right...Should we expect more symmetry? Can you use the implicit solver and crank up the number of nonlinear iterations?

[JFLemieux73]

Let<s talk about this at our meeting Thursday. We could then decide what are the next steps.

[dabail10]

Sounds good. I will trying setting ndte to 1200 in both.

[eclare108213]

Should we expect more symmetry?

Yes, ultimately it should be symmetric for symmetrically defined problems. It used to be -- not sure how long ago it became nonsymmetric. We can sample in time doing runs from older versions, to try to pinpoint when that happened, and target the revisions to the dynamics routines. A big change to the code was the addition of rEVP, but it might have happened before that. Or later.

[dabail10]

This is using ndte = 1200 in both. The B solution is more symmetric, but the CD solution is not. Note that the B-grid has a zero value in the upper right corner, while the CD-solution does not.

B-grid

CD-grid

[dabail10]

Would it be more symmetric with a cyclic boundary condition on the EW sides?

Dave

[eclare108213]

Possibly. But if the configuration (including any type of boundary conditions) and forcing are symmetric, the solution should be too.

[eclare108213]

In your B-grid plot above, what is the processor count and block size? The little box in the lower right interior looks suspicious. If you change the pe layout, do you get the same answer?

[dabail10]

I will try that. The blocks are 8x8 on 40 processors. I am trying a different layout.

[dabail10]

The solution is exactly the same with a 20x8x1 layout.

[dabail10]

Out of curiosity, I tried using umask instead of umaskCD in the CD grid solution. This does not change sigP. So, it does not appear to be boundary related.

[JFLemieux73]

Dave can you try as well with the implicit solver (kdyn=3)? You could increase the number of nonlinear iterations to 100 for example.

[dabail10]

I will give it a whirl.

[JFLemieux73]

I just finished a first draft on overleaf. The document describes the calculation of the strain rates, the viscous coefficients, the replacement pressure, the stresses, the divergence of the stress tensor and finally the solution of u and v for the C and CD grid. Please have a look here:
https://www.overleaf.com/project/620295ed9e158f02d60771f7

[dabail10]

Here is the kdyn=3 B grid solution. So which is correct above? :)

[JFLemieux73]

Thanks Dave! Hum...I will investigate this. Like I was telling you we are in a sense comparing apples and oranges as deltamin are not the same for kdyn=1 and kdyn=3. Also kdyn=1 has capping=1 while kdyn=3 has capping=0. I would like to add deltamin in the namelist so that we can have the same. I will also try to find out why we don't get symmetry (in all three...).

[apcraig]

We looked at symmetry briefly during the code camp in November. As @eclare108213 mentioned, it's clear that symmetry is broken at this point, that was clear in November too with the B grid cases we looked at. We don't know whether this is a roundoff issue in the numerics or a more devious issue. It's something we should try to understand. I recommend using the gbox12 case with 1 block and 1 pe for a simple boxsyme case to start with. You can manually write out any data you need for all points to try to understand when/where the bit-for-bit symmetry diverges. I am also willing to do that debugging if others don't. I think that's the fastest way to understand the issue.

[dabail10]

With closed boundaries and winds from the east at 5 m/s I wouldn't expect symmetry in the E-W direction. It is pretty symmetric in the N-S direction. The main differences are near the corners. Given that the sigP in the upper right is zero, but not in the other corners, perhaps this is why?

[apcraig]

@dabail10. I agree. Depending on the case you setup, different symmetry can be checked. You can use the open channel case too. In many ways, I created the open channel gbox12 case to verify identical answers at all points in the channel and then to make it even easier to debug symmetry. You only need to write out one column of data, not the entire box. If you switch to boxsymne, we should find symmetry across the diagonal. That was the original case we hoped to find symmetry, but I think we have to start with a simpler case.

[eclare108213]

When I fixed the discretization years ago (B grid) to make it symmetric, the order of ALL operations had to be identical as you rotate around the 4 corners, i.e. if a calculation is nw+ne+se for the ne corner, then the order has to be ne+se+sw for the se corner, etc. (That example is just for illustration -- I don't remember what the actual order was.)

[apcraig]

@eclare108213, that would be the roundoff error I'm thinking of too. Maybe we can have a cpp or runtime option to turn that on that code if we can get it to work. It would require special logic that checks the indices to decide what order to do operations in depending where you are on the grid and the case, and it would only work for very special cases that we would setup. The more devious issue I referenced would probably be a bug.

[JFLemieux73]

@apcraig @dabail10 I need your advice here. Right now tinyarea is computed in ice_grid.F90 as
tinyarea(i,j,iblk) = puny*tarea(i,j,iblk)
puny above is in fact deltamin for the calculation of the viscous coefficients. I would like to be able to define it in ice_in and have the possibility to define different values for EVP (kdyn=1) and VP (kdyn=3). I would like to do something like this in ice_grid:

 if (kdyn == 1) then
   tinyarea(i,j,iblk) = deltaminEVP*tarea(i,j,iblk)
 elseif (kdyn == 3) then
  tinyarea(i,j,iblk) = deltaminVP*tarea(i,j,iblk)
 endif

kdyn, deltamibnEVP and deltaminVP are part of the ice_dyn_shared module so I added in ice_grid:

use ice_dyn_shared, only: kdyn, deltaminEVP, deltaminVP

but it does not compile (a kind of circular issue I guess...).

Any suggestions?

[dabail10]

Yes, the problem is that ice_dyn_shared.F90 uses ice_grid.F90 for certain quantities, and you can't then use kdyn from ice_dyn_shared.F90 in ice_grid.F90. I would probably do this in ice_init.F90 when kdyn is read in from the namelist.

So, why is tinyarea two orders of magnitude larger in the VP? Are you proposing to do this for the B-grid as well?

[JFLemieux73]

deltamin is needed to regularize the viscous coeff so that they don't go to infinity. deltamin=2e-9 is the original value proposed by Hibler. In her papers Elizabeth wrote that for the EVP, the regularization is done by elastic waves and that deltamin is less important. I guess that;s why it was set to a different value (1e-11).

For small values of deltamin (I think 2e-9 is small...) the solution is not really affected by deltamin. However, implicit solvers struggle when deltamin is too small. This is why I don<t think it is a good idea to have 1e-11 for kdyn=3.

[eclare108213]

deltamin is the parameter that makes VP viscous. The E in EVP is a different way of regularizing the singularity that the V is addressing (viscosities become infinite as the ice becomes locked up), but we kept the V as an option for the most extreme cases and to maintain some similarity with the original VP approach. That means that deltamin doesn't need to be as large in EVP as in VP, to address the singularity problem, so we cranked it down by a couple orders of magnitude.

[JFLemieux73]

I was just a bit faster than you Elizabeth (9 seconds!). :)

[JFLemieux73]

Dave in the examples above is it possible the ocean currents are not zero (ocn_data_type = 'box2001')? (which might explain the lack of symmetry).

[apcraig]

A couple thoughts. First, why don't we have one parameter, deltamin, for all kdyns. It would be set in namelist with the default value being Hibler's or the EVP value. We don't need multiple variables. We just need one variable that depends on kdyn (and maybe other things later). Second, definition of tinyarea maybe should be in ice_grid anymore. grepping, it seems tinyarea really is just a dynamics thing. I would move it to init_dyn in ice_dyn_shared if that's possible and makes sense.

[dabail10]

Nope. I set it to 'calm'.

[JFLemieux73]

Ok thanks Dave and Tony. It seems to make sense to move it in the dynamics.

[apcraig]

You could even rename tinyarea to a new name. tinyarea=puny*tarea suggests it's a grid/roundoff thing. The new tinyarea is more like "dyn_area_limit", maybe?

[JFLemieux73]

tinyarea is only in the dynamics except here:
cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90

What should I do about this?

[dabail10]

That is the CESM driver. I'll take care of that one. Dave

…

On Thu, Mar 3, 2022 at 7:42 AM JFLemieux73 ***@***.***> wrote: tinyarea is only in the dynamics except here: cicecore/drivers/nuopc/cmeps/ice_mesh_mod.F90 What should I do about this? — Reply to this email directly, view it on GitHub <#694 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AC4ULUZSR5FJLLXXN2QBQALU6DFUFANCNFSM5PK2VWGA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you were mentioned.Message ID: ***@***.***>

-- David A Bailey email: dbailey_at_ucar.edu National Center for Atmospheric Research phone: 303-497-1737 PO Box 3000 fax : 303-497-1700 Boulder, CO 80307-3000 http://www.cgd.ucar.edu/staff/dbailey

[dabail10]

I removed the last comment. The issue was that I was using umask instead of umaskCD. Now everything is looking consistent at the T and U points. The sigP is interesting. Neither the B or CD look like the VP, but they are looking more similar.

I am now backing off to aice=0.5 in the three cases.

[apcraig]

FYI, box2001 forcing was both incorrect and not able to run on more than 1 block. This is fixed in apcraig#65

[apcraig]

Taking a step back, comparing standard box2001 on main and cgridDEV to make sure answers are not qualitatively different. Four plots are variable on main, variable on cgridDEV, absolute difference, relative difference. AICE is identical in both results.

UVEL:

main:

cgridDEV:

absolute difference:

relative difference:

VVEL:

main:

cgridDEV:

absolute difference:

relative difference:

[apcraig]

I think we should close this. Any disagree?

[apcraig]

Closing, active discussions of C-grid testing in #792. Feel free to reopen is needed.