Closed JFLemieux73 closed 2 years ago
Components of rheology term strintx.pdf strintxE.pdf strinty.pdf strintyN.pdf
To setup the channel case, do something like the following.
./cice.setup -m cheyenne -e intel -g gbox12 -p 1x1x12x12x1 -s boxsyme,diag1,run1day,gridb --test smoke --testid gD01
in ice_in, set
histfreq = ‘d’,‘1’,‘x’,‘x’,‘x’
close_boundaries = .false. ! not required
ice_data_type = ‘channel’
ew_boundary_type = ‘cyclic’
f_aice = ‘d1’
Results in table form as follows. This is output at i=5 for all j (1:12). Results at each i location are identical. The channel case here is 12x12 with land at T points j=1,2,11,12 and aice=0.5 at j=4,5,6,7,8. No ice at j=3,9,10. Column 1 = B grid uvel Column 2 = C grid uvelE Column 3 = CD grid uvelE Column 4 = CD grid uvelN Column 5 = CD grid uvelE with pstar=0. Column 6 = CD grid uvelN with pstar=0.
UVEL_1 UVELE_1 UVELE_1 UVELN_1 UVELE_1 UVELN_1
1 / 1: .... .... .... .... .... ....
2 / 2: .... .... .... .... .... ....
3 / 3: 0.08421 0.00000 0.000000 0.08421 0.00000 0.08421
4 / 4: 0.08421 0.08420 0.001599 0.08421 0.08421 0.08421
5 / 5: 0.08421 0.08420 0.002078 0.08421 0.08421 0.08421
6 / 6: 0.08421 0.08420 0.002238 0.08421 0.08421 0.08421
7 / 7: 0.08421 0.08420 0.002078 0.08421 0.08421 0.08421
8 / 8: 0.08421 0.08420 0.001599 0.08421 0.08421 0.08421
9 / 9: 0.00000 0.00000 0.000000 0.00000 0.00000 0.00000
10 / 10: 0.00000 0.00000 0.000000 0.00000 0.00000 0.00000
11 / 11: .... .... .... .... .... ....
12 / 12: .... .... .... .... .... ....
Same for v velocities,
VVEL_1 VVELN_1 VVELE_1 VVELN_1 VVELE_1 VVELN_1
1 / 1: .... .... .... .... .... ....
2 / 2: .... .... .... .... .... ....
3 / 3: -1.420E-17 -6.481E-04 0.000E+00 -2.033E-19 0.0000 0.0000
4 / 4: -8.501E-18 1.780E-07 1.677E-04 -4.794E-19 0.0000 0.0000
5 / 5: -2.831E-18 9.961E-08 1.198E-04 -6.700E-20 0.0000 0.0000
6 / 6: 2.831E-18 -1.093E-07 -1.290E-18 1.572E-19 0.0000 0.0000
7 / 7: 8.501E-18 -1.574E-07 -1.198E-04 9.004E-19 0.0000 0.0000
8 / 8: 1.420E-17 6.481E-04 -1.677E-04 4.571E-19 0.0000 0.0000
9 / 9: 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.0000 0.0000
10 / 10: 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.0000 0.0000
11 / 11: .... .... .... .... .... ....
12 / 12: .... .... .... .... .... ....
The obvious issue is the decoupled CD solution associated with Column 3 in U. The other thing to note is the relatively poor performance of the V velocity in the C solution, Column 2 in V. Pstar=0 fixes the problems with the CD solution.
Here are the latest tables for the same channel cases as above with ndte=1200 and a new column for VP results. Results are unchanged. VP "v" is not quite as zero as some other cases.
Column 1 = B grid uvel Column 2 = C grid uvelE Column 3 = CD grid uvelE Column 4 = CD grid uvelN Column 5 = CD grid uvelE with pstar=0. Column 6 = CD grid uvelN with pstar=0. Column 7 = B grid uvel with kdyn=3 (vp)
UVEL_1 UVELE_1 UVELE_1 UVELN_1 UVELE_1 UVELN_1 UVEL_1
1 / 1: .... .... .... .... .... .... ....
2 / 2: .... .... .... .... .... .... ....
3 / 3: 0.08421 0.00000 0.000000 0.08421 0.00000 0.08421 0.08421
4 / 4: 0.08421 0.08420 0.001599 0.08421 0.08421 0.08421 0.08421
5 / 5: 0.08421 0.08420 0.002078 0.08421 0.08421 0.08421 0.08421
6 / 6: 0.08421 0.08420 0.002238 0.08421 0.08421 0.08421 0.08421
7 / 7: 0.08421 0.08420 0.002078 0.08421 0.08421 0.08421 0.08421
8 / 8: 0.08421 0.08420 0.001599 0.08421 0.08421 0.08421 0.08421
9 / 9: 0.00000 0.00000 0.000000 0.00000 0.00000 0.00000 0.00000
10 / 10: 0.00000 0.00000 0.000000 0.00000 0.00000 0.00000 0.00000
11 / 11: .... .... .... .... .... .... ....
12 / 12: .... .... .... .... .... .... ....
VVEL_1 VVELE_1 VVELE_1 VVELN_1 VVELE_1 VVELN_1 VVEL_1
1 / 1: .... .... .... .... .... .... ....
2 / 2: .... .... .... .... .... .... ....
3 / 3: -1.999E-17 -3.240E-04 0.000E+00 -4.366E-19 0.0000 0.0000 -1.007E-11
4 / 4: -1.196E-17 -3.240E-04 1.677E-04 -1.903E-19 0.0000 0.0000 -1.117E-12
5 / 5: -3.983E-18 3.929E-08 1.198E-04 -1.449E-19 0.0000 0.0000 -1.134E-13
6 / 6: 3.983E-18 -4.227E-09 4.461E-18 1.590E-21 0.0000 0.0000 1.134E-13
7 / 7: 1.196E-17 -4.222E-08 -1.198E-04 2.518E-19 0.0000 0.0000 1.117E-12
8 / 8: 1.999E-17 3.240E-04 -1.677E-04 5.885E-19 0.0000 0.0000 1.007E-11
9 / 9: 0.000E+00 3.240E-04 0.000E+00 0.000E+00 0.0000 0.0000 0.000E+00
10 / 10: 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.0000 0.0000 0.000E+00
11 / 11: .... .... .... .... .... .... ....
12 / 12: .... .... .... .... .... .... ....
If I force v=0 is step_vel for the channel CD case, v==0 in output and the u solutions are unchanged.
Baseline uE, vE, uN, vN:
UVELE_1 VVELE_1 UVELN_1 VVELN_1
1 / 1: .... .... .... ....
2 / 2: .... .... .... ....
3 / 3: 0.000000 0.000E+00 0.08421 -4.366E-19
4 / 4: 0.001599 1.677E-04 0.08421 -1.903E-19
5 / 5: 0.002078 1.198E-04 0.08421 -1.449E-19
6 / 6: 0.002238 4.461E-18 0.08421 1.590E-21
7 / 7: 0.002078 -1.198E-04 0.08421 2.518E-19
8 / 8: 0.001599 -1.677E-04 0.08421 5.885E-19
9 / 9: 0.000000 0.000E+00 0.00000 0.000E+00
10 / 10: 0.000000 0.000E+00 0.00000 0.000E+00
11 / 11: .... .... .... ....
12 / 12: .... .... .... ....
v set to zero:
UVELE_1 VVELE_1 UVELN_1 VVELN_1
1 / 1: .... .... .... ....
2 / 2: .... .... .... ....
3 / 3: 0.000000 0.0000 0.08421 0.0000
4 / 4: 0.001599 0.0000 0.08421 0.0000
5 / 5: 0.002078 0.0000 0.08421 0.0000
6 / 6: 0.002238 0.0000 0.08421 0.0000
7 / 7: 0.002078 0.0000 0.08421 0.0000
8 / 8: 0.001599 0.0000 0.08421 0.0000
9 / 9: 0.000000 0.0000 0.00000 0.0000
10 / 10: 0.000000 0.0000 0.00000 0.0000
11 / 11: .... .... .... ....
12 / 12: .... .... .... ....
A few other results before I forget as I continue testing,
Here is the first subcycle: stress12T 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 stress12U 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 strintxE 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 strintxN 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 uvelE 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 uvelN 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 3.078947368421053E-004 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
Second subcycle: stress12T 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 stress12U 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.863973461411587 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -0.863973461411587 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 strintxE 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -5.399834133822417E-005 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -5.399834133822417E-005 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 strintxN 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 uvelE 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 6.149286482003538E-004 6.157812421928050E-004 6.157812421928050E-004 6.157812421928050E-004 6.149286482003538E-004 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 uvelN 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 6.157812421928050E-004 6.157812421928050E-004 6.157812421928050E-004 6.157812421928050E-004 6.157812421928050E-004 6.157812421928050E-004 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
End of run: stress12T 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000
stress12U 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 779.587470870799 467.699878980278
155.889843963600 -155.889843963600 -467.699878980278
-779.587470870798 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000
strintxE 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 -1.949297449315755E-002
-1.948812718854237E-002 -1.948623049545003E-002 -1.948812718854232E-002
-1.949297449315756E-002 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000
strintxN 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000
uvelE 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
0.000000000000000E+000 1.598515885989748E-003 2.078016700778288E-003
2.237839921854726E-003 2.078016700778291E-003 1.598515885989755E-003
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
uvelN 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
8.421243560722592E-002 8.421243560722592E-002 8.421243560722592E-002
8.421243560722592E-002 8.421243560722592E-002 8.421243560722592E-002
0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000
Here are the zeta/eta
ksub = 1zetax2T 0.000000000000000E+000
etax2T 0.000000000000000E+000
zetax2U 243.847278997721
etax2U 60.9618197494304
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 243.847278997721
etax2U 60.9618197494304
ksub=2
zetax2T 0.000000000000000E+000
etax2T 0.000000000000000E+000
zetax2U 243.847278997721
etax2U 60.9618197494304
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 243.847278997721
etax2U 60.9618197494304
end of run:
zetax2T 0.000000000000000E+000
etax2T 0.000000000000000E+000
zetax2U 243.847278997721
etax2U 60.9618197494304
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 243.847278997721
etax2U 60.9618197494304
Fixed values.
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
zetax2T 487.694557995443
etax2T 121.923639498861
zetax2U 487.694557995443
etax2U 121.923639498861
With fixes to viscous coeff U calc using capping, CD grid
UVELE_1 VVELE_1 UVELN_1 VVELN_1
1 / 1: .... .... .... ....
2 / 2: .... .... .... ....
3 / 3: 0.00000 0.000E+00 0.08421 -1.879E-17
4 / 4: 0.08420 5.362E-08 0.08421 -1.699E-17
5 / 5: 0.08420 3.234E-08 0.08421 -7.977E-18
6 / 6: 0.08420 -4.004E-22 0.08421 7.977E-18
7 / 7: 0.08420 -3.234E-08 0.08421 1.699E-17
8 / 8: 0.08420 -5.362E-08 0.08421 1.879E-17
9 / 9: 0.00000 0.000E+00 0.00000 0.000E+00
10 / 10: 0.00000 0.000E+00 0.00000 0.000E+00
11 / 11: .... .... .... ....
12 / 12: .... .... .... ....
For the C grid with update viscous coeffs,
UVELE_1 VVELN_1
1 / 1: .... ....
2 / 2: .... ....
3 / 3: 0.00000 -6.481E-04
4 / 4: 0.08419 -1.158E-07
5 / 5: 0.08421 3.232E-08
6 / 6: 0.08421 -3.232E-08
7 / 7: 0.08421 1.158E-07
8 / 8: 0.08419 6.481E-04
9 / 9: 0.00000 0.000E+00
10 / 10: 0.00000 0.000E+00
11 / 11: .... ....
12 / 12: .... ....
I have started to look at the boxsyme test for the B/C/CD cases. The default setup has advection and thermodynamics off just like our channel flow. It also has aice=1.0. The ice is basically "stuck". This produces a "stuck" solution and while the velocities in the three cases have similar properties, they are also quite different. I believe we are just looking at different numerics. If I change aice to 0.5, then I get free drift in the interior, zero velocity against the wall, and the solutions are a lot more similar except in the corners. I don't see how we can ever have identical results near the boundaries given the different staggers. I am running for 6 days and have increased ndte to 1200. I want to turn on advection next.
Should the boxsym tests be using aice=1.0 in general. Is that really the test we want? Should I continue to test aice=1.0 or is aice=0.5 better? Or do we really need to do both?
I attach images of uvel and vvel for the B/C/CD results for the boxsyme test with aice=0.5 after 5 days with ndte=1200. The scales are the consistent. I am not plotting u=v=0 at the walls. Outside of the corners, you can see that the "change in color" is slightly different for each solution, but those differences are small. I think these look pretty good. Looking at some raw numbers, the E and N solutions seems to be quite consistent for the CD solution which is good. The one thing that does jump out is the solution on the left and right edge for the C grid especially noticeable in the v solution in the corners. The maximum magnitude is not right at the edge as it is for other solutions.
Should the boxsym tests be using aice=1.0 in general. Is that really the test we want? Should I continue to test aice=1.0 or is aice=0.5 better? Or do we really need to do both?
aice=1.0 is too stiff of a problem. aice=0.5 doesn't really test the stress terms, although it is useful for checking symmetry and similar kinds of debugging. aice=0.9 or a little higher could be okay. This is why the box2001 test is set up with a ramp in aice from 0 to 1 -- you get all the various kinds of behavior.
Some other test results from today
-boxwallblockp5 (an aice=0.5 box on the east wall with east forcing) with advection after 3 months is as follows. Now we have an interesting case that is starting to show noise in both the C and CD cases. The block is being push against the wall and spreading. Maybe this is the next thing to debug?
Wow, very interesting. How small are the values left behind (in purple)? Wondering if I can use this to debug residual ice issues.
I believe the purple values are less than 0.01. I tried to enhance the shading there too, but it got smaller than I thought we cared about. That purple basically represents where the box started. We have eastward forcing and it mostly piles up and spreads on the east wall, but a tiny bit is left behind apparently.
I redid the kmtislands (CD) case and printed some values. On the second sub-cycle, I get:
uvelE 2.472165410902706E-004 (i=50,j=33)
uvelE 2.472165410902706E-004 (i=49,j=33)
uvelE 2.472165410902706E-004 (i=51,j=33)
uvelE 0.000000000000000E+000 (i=50,j=34)
uvelE 0.000000000000000E+000 (i=50,j=32)
stresspU -397.468827643619
stresspU -397.468827643619
stresspU -397.468827643619
stresspU 0.000000000000000E+000
stresspU -397.468827643619
stress12U -99.3672069109048
stress12U -99.3672069109048
stress12U -99.3672069109048
stress12U 0.000000000000000E+000
stress12U 99.3672069109048
So, doing some back of the envelope calculations, assuming vvel = 0.
`` divU = 0. tensionU = 0. shearU = 16000. (0. - 2.47x10^-4) = 3.952 deltaU = 3.952 tmpcalc =~ 20000./3.952 =~ 5060. zetax2U = 5060. etax2U = 0.255060. = 1265. rep_prsU =~ 20000.
stresspU =~ -20000. stress12U =~ 0.51265.3.952 `` This assumes a strength of 20000. It is probably less than this. Anyhow, do these numbers make sense for a one-gridpoint wide channel? Should there be a stresspU / stress12U value at the land gridcell below, but not above?
Hi Dave,
I think we should get nonzero values for stresspU and stress12U both above and below.
I am working on a special test case for one-gridpoint wide channels. I think we can get an analytical solution. It will be easier to understand what is going on.
There are analytical solutions for one grid cell channels with cyclic BCs. Have a look at the pdf below. For strong winds there is a plastic solution. For weak winds we can also find a viscous solution. It also allows us to calculate the shear stresses at the walls.
Hi Tony,
About the figures above...Have you tried to reduce the time step to see if you still get the noise in the C (and CD) solutions?
I have not tried reduced timestep. I assume that means ndte=1200? I can give that a try.
Sorry I meant the advective time step. It is probably set to one hour. I would try 15min.
Hi Tony,
As you can see we have decided to create one issue per test case. You could go ahead and create issues for the other test cases you are testing (e.g. boxwallblockp5). It would be nice if the first comment of each issue was the command for the setup (and possibly things that need to be modified in ice_in). Finally maybe we could close this issue. Thanks!
Closing, see #684 #686 #687 #688 #689
ncview.uvel.pdf ncview.uvelE.pdf ncview.vvel.pdf ncview.vvelN.pdf