-
Notifications
You must be signed in to change notification settings - Fork 0
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
umaskCD and more #33
umaskCD and more #33
Conversation
…ransport files are changed in a way so that velocities are interpolated to b grid for depature point function and kept at E or N grid possible.
…nted out in vp and eap
…ransport files are changed in a way so that velocities are interpolated to b grid for depature point function and kept at E or N grid possible.
…nted out in vp and eap
Added function to fint zeta, eta and replacement pressure at U
I cloned a fresh sandbox from @apcraig's branch (through #37), pulled the changes from this branch onto it, fixed the conflict, and modified the code so that it compiles. A basic smoke test with |
Merging #38 didn't help my runs, but I might not have all of Till's intended changes -- I only fixed what had been committed so far. For the gx3 grid, is the dynamic forcing being interpolated from T to the N, E cell faces? Or is JRA coming in at U points and shouldn't be shifted? Some of the grid_sys suite runs complete but the velocities are not correct. Those will be good debugging platforms. |
Ok, thanks. |
See #39 |
Only quick suite has been run.
Without having merge your changes. I can now build and run this pull request for quick_suites. This passes the quick test |
I will compare the two and merge tomorrow. |
I believe the JRA forcing all comes in on T points just like the CESM coupling (for now). |
I should note that I only ran b grid solution not cd. |
Before finalizing PR, please run both B and CD if possible. B to check bit-for-bit. CD to check that it runs (maybe even with debug flags on). Neither is absolutely required, and whether they pass or not, I can also test/debug before (or after) the PR is merged and then created fixes and/or provide feedback. Thanks! |
emask(:,:,iblk) = .false. | ||
do j = jlo-nghost, jhi | ||
do i = ilo-nghost, ihi | ||
!!!!!WARNING THIS SUM is not calculating uvmCD mask on the northern and eastern halo |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is fixed in #38, and you can probably just replace this ice_grid.F90 file with that one.
I have merged 33 and 39. 33 is still binary identical. My new merge is not. |
I'm a little confused what's happening (which is OK), but might it be a good idea to create a new branch from the head of cgridDEV (picking up #38 and others), then merging in #33 and #39 to that branch, resolving conflicts, and debugging from there? Then you could cancel #33 and #39 and create a new PR. |
I did make a new based on my merge of 39. |
I think the current cgridDEV head is bit-for-bit with the cice-consortium main. But I can check again. That's with intel and gnu on cheyenne. With pgi on cheyenne, it isn't. So there is a sensitivity to compiler version and we know pgi is particular tricky. You could also try running a consortium main and the cgridDEV with mods with the debug flags. This reduces the optimization and should make it more likely to be bit-for-bit. Let me know if you have questions or want me to test/debug anything in particular. |
Thanks |
Sounds good. When you think it's ready, I can run some independent tests just to confirm. Let me know when you feel it's ready. |
Instead an additional if else gridsystem is introduced. subroutine evp should be refactored
uvel (:,:,iblk), vvel (:,:,iblk), & | ||
Tbu (:,:,iblk)) | ||
|
||
if (grid_system == trim('B')) then |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
why is the trim around 'B'?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
should have been around grid_system (if needed)
Currently not doing anything.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please move it to grid_system. We may clean all of these up later.
! COMING SOON!!! | ||
zetax2U = c0 | ||
etax2U = c0 | ||
rep_prsU = c0 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It looks like these initializations were lost in the new code modifications. Or were they moved somewhere else?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
they dont need to be initialized. They are calculated in the function after as intent out.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
These lines were added as a temporary fix when viscous_coeffs_and_rep_pressure_U was not called. Now it is. I think this change is correct in the big picture.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
zetax2T_00,zetax2T_10,zetax2T_11,zetax2T_01, & | ||
etax2T_00, etax2T_10, etax2T_11, etax2T_01, & ! 2 x visous coeffs, replacement pressure | ||
maskT_00, maskT_10, maskT_11, maskT_01, & | ||
tarea_00, tarea_10, tarea_11, tarea_01, & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
please comment the code to explain what 00, 10, 01, 11 are
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
will do. These are i+0,j+0,i+1 and j+1
@@ -1899,7 +1903,7 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & | |||
real (kind=dbl_kind), intent(in):: & | |||
Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner | |||
|
|||
logical , intent(in):: capping | |||
real(kind=dbl_kind) , intent(in):: capping |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It looks like the intent here is that capping is either c0 or c1 as if it were a logical but by setting it to a real it avoids an if check. That should be documented in the code and should there be a check added that capping must be either c0 or c1 valued?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
capping was introduced by @JFLemieux73 when he reorganized the EVP solver, This introduced an if statement within the viscous_coeffs_and_rep_pressure function. To avoid the if statement I replaced it with a c0/c1. It is c1 in the evp and c0 in the vp and set as parameters further up. Currently hard coded. The intent is to replace it with a namelist parameter both in main and in the cd grid code. This should be a string perhaps "smooth_cap" and "hard_cap" which is then turned into real (c0 and c1)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Will it always be either 0 or 1? And will it always be 1 for EVP (and presumably the variants) and 0 for VP? If so, could it be set based on kdyn at initialization, instead of namelist?
I ran a full test suite on cheyenne (all B configurations) with cgridDEV plus this PR and all results are bit-for-bit. So I think the only issue is that the CD cases don't complete (not surprising) and whether we want to do anything more to this PR (outside maybe some cleanup) before merging and continuing to debug the CD configuration. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This can be merged once Till is happy with it. I've requested a few changes.
uvel (:,:,iblk), vvel (:,:,iblk), & | ||
Tbu (:,:,iblk)) | ||
|
||
if (grid_system == trim('B')) then |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please move it to grid_system. We may clean all of these up later.
! COMING SOON!!! | ||
zetax2U = c0 | ||
etax2U = c0 | ||
rep_prsU = c0 |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
@@ -1899,7 +1903,7 @@ subroutine viscous_coeffs_and_rep_pressure (strength, tinyarea, & | |||
real (kind=dbl_kind), intent(in):: & | |||
Deltane, Deltanw, Deltasw, Deltase ! Delta at each corner | |||
|
|||
logical , intent(in):: capping | |||
real(kind=dbl_kind) , intent(in):: capping |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Will it always be either 0 or 1? And will it always be 1 for EVP (and presumably the variants) and 0 for VP? If so, could it be set based on kdyn at initialization, instead of namelist?
Currently yes it is the case that capping for evp=c1 and capping vp =c0 |
I would like to have it merged. I fixed the trim('B') but I did not add comments nor did change capping nor did I add comments. |
This bug is removed if the creation of umaskCD iincludes an icemask just as in dyn_prep2 for umask (B grid) |
nmassdti (:,:,iblk), fcorN_blk (:,:,iblk), & | ||
nmask (:,:,iblk), & | ||
nmassdti (:,:,iblk), fcorN_blk (:,:,iblk),& | ||
nmask (:,:,iblk), & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
can we avoid unnecessary whitespace changes ? this introduces false positive in 'git blame' output, if care is not taken..,.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
changed
|
||
call viscous_coeffs_and_rep_pressure (strength(i,j), tinyarea(i,j),& | ||
Deltane, Deltanw, & | ||
Deltasw, Deltase, & | ||
zetax2ne, zetax2nw, & | ||
zetax2sw, zetax2se, & | ||
etax2ne, etax2nw, & | ||
etax2sw, etax2se, & | ||
rep_prsne, rep_prsnw, & | ||
rep_prssw, rep_prsse, & | ||
capping) | ||
|
||
call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& | ||
Deltane, zetax2ne, & | ||
etax2ne, rep_prsne, & | ||
capping) | ||
|
||
call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& | ||
Deltanw, zetax2nw, & | ||
etax2nw, rep_prsnw, & | ||
capping) | ||
|
||
call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& | ||
Deltasw, zetax2sw, & | ||
etax2sw, rep_prssw, & | ||
capping) | ||
|
||
call viscous_coeffs_and_rep_pressure_T (strength(i,j), tinyarea(i,j),& | ||
Deltase, zetax2se, & | ||
etax2se, rep_prsse, & | ||
capping) | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That is not strictly necessary, right ? it's just cleaning up the B-grid implementation. So maybe it could be done in a later step ? If it's bit for bit, then all for the better and I guess it does not matter much...
epm, npm, uvm, & | ||
epm, npm, hm, uvm, & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it would be preferable to restrict ourselves to 2 arguments per line, for better readability.
! ice extent mask (U-cells) | ||
iceumask(i,j) = (umask(i,j)) .and. (aiu (i,j) > a_min) & | ||
.and. (umass(i,j) > m_min) | ||
else ! ice mask shpuld be applied to cd grid. For now it is not implemented. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we agreed on spelling out elseif ( trim(grid_sytem) == 'CD')
(or using a select
statement instead). It's also clearer when we read the code, I think.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Obsolete. see soon upcoming push that make cd grid run complete
@@ -155,7 +155,7 @@ module ice_grid | |||
real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: & | |||
hm , & ! land/boundary mask, thickness (T-cell) | |||
bm , & ! task/block id | |||
uvm , & ! land/boundary mask, velocity (U-cell) | |||
uvm , & ! land/boundary mask, velocity (U-cell) - land if one point is one |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure about that comment means exactly .
Also below in alloc_grid
you used a different wording: "water in case of all water point"
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For concistency. Comment is removed
if (uvm(i,j,iblk) > p5 ) umask (i,j,iblk) = .true. | ||
if (uvmCD(i,j,iblk) > c1p5) umaskCD(i,j,iblk) = .true. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
alignment:
if (uvm(i,j,iblk) > p5 ) umask (i,j,iblk) = .true. | |
if (uvmCD(i,j,iblk) > c1p5) umaskCD(i,j,iblk) = .true. | |
if (uvm(i,j,iblk) > p5 ) umask (i,j,iblk) = .true. | |
if (uvmCD(i,j,iblk) > c1p5) umaskCD(i,j,iblk) = .true. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
resolved
enddo | ||
!$OMP END PARALLEL DO | ||
|
||
!$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block) | ||
do iblk = 1, nblocks |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Any reason why we are creating a second loop here ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
i and j loops are slightly different, iblk loop could be merged.
indxui (:,iblk), indxuj (:,iblk), & | ||
aiu (:,:,iblk), umass (:,:,iblk), & | ||
umassdti (:,:,iblk), fcor_blk (:,:,iblk), & | ||
umaskCD (:,:,iblk), & |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's unclear to me why we need yet another call to dyn_prep2 here. I thought that all CD grid quantities were being dealt with in the two calls to dyn_prep2 below (for 'E' and 'N' quantities, line 440 and 473 in this version). Here your are adding a new call, identical ( as far as I can tell) to the B-grid case, apart from this new umaskCD
array instead of umask
... But on the CD grid, we are not using most of these arrays that are sent as arguments here...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
umaskCD is defined different, icellu will be different (stress_u). It can be that calling dyn_prep2 is too much and a call to a simpler function could be used.
@TillRasmussen, if you are making changes, they are not showing up in the PR. You will need to push to the cgrid-something branch. No rush, just wanted to point that out. Thanks! |
…ails on 2 test. 1 is pending. Comments from Phillipe also included
small correction to cd grid. One run failed to finish when including the divergence stress #34. Otherwize no more corrections from here |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think we should go ahead and merge this as is, so that we all have a common code to debug. There will be a lot of cleanup on the new code before we merge it into the main repo, but for now we can put these things on our cleanup list. From this PR, these are
- comment the code to explain that 00, 10, 01, 11 are i+0,j+0,i+1 and j+1
- possibly define capping based on kdyn
- join iblk loops in ice_grid.F90 (and review loop structure everywhere)
- white space, alignments, comments, etc
I would not define capping based on kdyn. There should be an additional input in the namelist. For example I would like to have the possibility to have capping=true or false for the implicit solver. Sorry I am a bit busy this week and cannot help for the debugging. I should be able to get back to this next week. |
ok, thanks @JFLemieux73 |
I am running the gridsys test suite quickly and then will merge. I believe the CD cases are now running to completion, so that's great. This should be an excellent place for us to start testing/debugging the results. I have also added @eclare108213's bullets to CICE-Consortium#660. |
The test suite looks OK. The only problem is that the boxislands cases are aborting when CD is set. A sample traceback is (abort_ice)ABORTED: I wonder if there are invalid values in that array or something like it. We should debug this case further as we move forward. The equivalent symmetry tests without boxislands run. The gridsys results for de8f6e9 are I will merge this and #40 now and then fire up a full test suite on cheyenne. |
You can assign me to it. |
I changed diagfreq to 1 and it now aborts on the first step in the diagnostics. The speed is blowing up here at least at one point. I'll keep looking. Dave total ice area (km^2) = 2.64780799999989523E+06 0.00000000000000000E+00 |
Quick update. Looks like the internal stresses are blowing up. This is a print inside step_vel. aiu 0.999999999999923 |
Thanks Dave. Is it strintx,y that make uvel,vvel blow up or the opposite? Could you please print strintx,y and uvelE,vvelE at the beginning of the 1st subcycle? |
I agree there is a bit of a chicken or egg here. Here are some print statements after stressT and stressU It does appear that stress12U goes pretty big on the third subcycle. Then the velocities at the fourth subcycle go wacky. I believe the stress12U components are going bad first ... or maybe not. Is a stress value of 1639176.08411872 too big? Dave |
Well...the fact that stresspT and stressmT are zero is weird. There is something wrong there. I suggest to use kstrength=0 for debugging (it is easier to interpret the values of stresses). With kstrength = 0, the stresses have values on the order of Pstar (default Pstar is 27500). We should expect values between 0 and 10*Pstar. I will help you as soon as I can with the debugging. What test case are you running? gx3? |
This is a boxislande case. It is a new case with a limited area mask. I am not that familiar with all the settings. |
Full test results on cheyenne with 3 compilers is here, The test suites with 232 cases are the standard full suites compared on the B grid with Consortium main. The test suites with 28 tests are the gridsys suite which consists of B and CD tests and those are compared to the last cgridDEV version (which is pretty meaningless at this point). Everything looks pretty good. B configurations are still bit-for-bit with Consortium main except pgi and some box tests that have changes. CD cases generally run except the case with boxislands. The pgi non-bit-for-bit is expected. The boxislands will require a bit more debugging, but that seems to be underway. And again, this does not suggest any of the CD results are reasonable, only that the short CD test runs completed. |
Here are some questions that might provide clues. Maybe you've already looked at these things.
|
I have started to investigate what is going on. I think there is an instability in compact regions close to land. I have created a new issue specifically for this: CICE-Consortium#666 It is issue 666...not a good sign! :) |
For detailed information about submitting Pull Requests (PRs) to the CICE-Consortium,
please refer to: https://github.com/CICE-Consortium/About-Us/wiki/Resource-Index#information-for-developers
Still a few bugs
PR checklist
ENTER INFORMATION HERE
ENTER INFORMATION HERE
ENTER INFORMATION HERE