Calls to outfld in atm that create temp arrays hit invalid float (when pcols != ncols) Intel/debug/cori-knl

[ndkeen]

A test in acme_developer SMS_D_Ln1.ne30_oEC.F1850C5AV1C-02.

While adjusting PE layouts, I hit the error below. I assumed it was something I did wrong, but then I happened to notice the code stopped in a similar call on edison with a different problem.

Note issue: #1241

480: forrtl: error (65): floating invalid
480: Image              PC                Routine            Line        Source
480: acme.exe           000000000ABEA8D1  Unknown               Unknown  Unknown
480: acme.exe           000000000ABE8A0B  Unknown               Unknown  Unknown
480: acme.exe           000000000AB97114  Unknown               Unknown  Unknown
480: acme.exe           000000000AB96F26  Unknown               Unknown  Unknown
480: acme.exe           000000000AAFF789  Unknown               Unknown  Unknown
480: acme.exe           000000000AB0B8DC  Unknown               Unknown  Unknown
480: acme.exe           000000000A6CB910  Unknown               Unknown  Unknown
480: acme.exe           0000000004236D4D  hetfrz_classnuc_c         965  hetfrz_classnuc_cam.F90
480: acme.exe           00000000030B15C3  microp_aero_mp_mi         852  microp_aero.F90
480: acme.exe           0000000000CE8883  physpkg_mp_tphysb        2407  physpkg.F90
480: acme.exe           0000000000CC0681  physpkg_mp_phys_r        1008  physpkg.F90

This happens before any time stepping.

Below is my attempt to summarize the issue in the code components/cam/src/physics/cam/hetfrz_classnuc_cam.F90. @gold2718 suggested I try to change the temp array calculation to only work over ncol. As I'm not that familiar with this code, is it OK to change just those calls to outfld() to work over ncol and leave the others as pcol? Or are we hiding another issue?

 ! ndk -- just noting how these arrays are initialized -- if I simply change this to niimm_bc=0, the problem goes away
   niimm_bc(:ncol,:) = 0._r8
   nicnt_bc(:ncol,:) = 0._r8
   nidep_bc(:ncol,:) = 0._r8
...
  ! ndk just noting there are several calls to outfld
   call outfld('BCFREZIMM', nnuccc_bc, pcols, lchnk)
   call outfld('BCFREZCNT', nnucct_bc, pcols, lchnk)
   call outfld('BCFREZDEP', nnudep_bc, pcols, lchnk)

  ! ndk however these are the lines that cause problems
   !call outfld('NIMIX_IMM', niimm_bc+niimm_dst, pcols, lchnk)  !ndk                                                                                                                                  
   !call outfld('NIMIX_CNT', nicnt_bc+nicnt_dst, pcols, lchnk)
   !call outfld('NIMIX_DEP', nidep_bc+nidep_dst, pcols, lchnk)

 !ndk -- goldy suggested making this change:
   call outfld('NIMIX_IMM', niimm_bc(:ncol,:)+niimm_dst(:ncol,:), ncol, lchnk)  !ndk                                                                                                                                  
   call outfld('NIMIX_CNT', nicnt_bc(:ncol,:)+nicnt_dst(:ncol,:), ncol, lchnk)
   call outfld('NIMIX_DEP', nidep_bc(:ncol,:)+nidep_dst(:ncol,:), ncol, lchnk)

case dir:

/global/cscratch1/sd/ndk/acme_scratch/m26f/SMS_D_Ln1.ne30_oEC.F1850C5AV1C-02.cori-knl_intel.20170209_142749

[gold2718]

My suggested change does not hide another issue because those array locations should not be used by anyone. Almost all of the physics code is loops from 1 to ncol which is correct and efficient.
Physics columns are always packed into the first ncol slots of each chunk (physics columns are broken into chunks for threading). pcols (a compile time constant) is the maximum size of chunks.
To summarize, I think my suggested change is the correct way to write outfld calls. You will find that usage in many modules in the ACME atmosphere code but (unfortunately), its usage is not consistent. Memory testing (compiler options that initializes all memory to NaN) is a useful method to find code like this where neither correct option is used (the other 'correct' option is to initialize all values of arrays, even if they are not used in calculations).

[ndkeen]

OK, so just to be clear, it should be OK (and suggested) to change all outfld calls to use ncol?

[gold2718]

Yes for any call that is currently using pcols. There are a few that are different (e.g., micro_mg_cam.F90) and you should not change those.

[worleyph]

is it OK to change just those calls to outfld() to work over ncol and leave the others as pcol?

The way that I read this, niimm_bc+niimm_dst is doing arithmetic on arrays that contain uninitialized values. The other calls are not doing any calculation.

I am concerned with the proposed fix:

call outfld('NIMIX_CNT', nicnt_bc(:ncol,:)+nicnt_dst(:ncol,:), ncol, lchnk)

as I am not sure that the declared first dimention of the generated temporary array would be. The ncol argument in outfld represents the declared first dimension in the 2D array.

[worleyph]

So, I disagree with @gold2718 - the array declarations in outfld are

real(r8), intent(in) :: field(idim,*) ! Array containing field values

If the field is declared as (pcols, pver), then the parameter idim needs to be pcols, doesn't it?

[gold2718]

Yes which is why it is important to pass the field as field(:ncol,:), not just a naked field.
What you are saying is that the size of the first dimension of the field array must match the idim input parameter.
I claim that field(:ncol,:) meets that requirement.

[worleyph]

Thanks. So we should always be creating temporary arrays in these calls, not just passing pointers? Guess for chunks this won't cause memory issues.

[worleyph]

Note that this does cause performance problems in threaded regions for the PGI compiler, and is less efficient even for the Intel compiler. Outfld calls typically are not in inner loops (at least, I hope not), so probably irrelevant.

[worleyph]

Too bad outfld does not have both actual and declared first dimensions in the call. Would eliminate the need to generate temporary arrays, though this would not have affected @ndkeen 's error since it was adding these fields before the call to outfld.

[gold2718]

outfld calls are never supposed to be in the inner loop but only at the parameterization interface level or higher.
Still, I see your point that we should use native arrays where convenient. How does this sound for advice? Assume we have two fields:

real(r8) :: pfield(pcols, pver)
real(r8) :: field(ncol, pver)
...
call outfld ('name1', pfield(:ncol,:), ncol, c)
call outfld('name2', field, ncol, c)

This should always work, no?

[gold2718]

Yes, one of these days, I will rewrite the history processing to use real buffers so that the outfld interface will take proper Fortran arrays (not FORTRAN 77 as it does currently).
One of these days . . .

[worleyph]

call outfld ('name1', pfield(:ncol,:), pcols, c)

I thought that this is what you just pointed out to me would not work? A temporary array is created with the dimensions (:ncol,:) that is then passed in, so the declared dimension should be ncols?

[gold2718]

Gaa, I shouldn't type when I'm tired. I fixed the example and your comment is correct. Where possible, work arrays should be defined with dimension ncol (e.g., the chemistry code does this in most places), so no temporary is required.
However, given that many outfld calls are using fields already defined as pcols, those will have to be chopped into a temporary array (e.g., pfield(:ncol,:)) when passed to outfld.
@ndkeen is this clear enough for you to proceed?

[worleyph]

I'm late to the game here, but I wonder if there are two separate issues - a potential refactor and fixing this particular bug.

My preference is not to do array arithmetic as part of a subroutine call parameter. While "uglier", doing the computation outside the call with an explicitly declared work array, and then passing in this work array does not depend on the compiler to be doing the right thing (functionally and performance-wise).

This array could be declared with ncols, and then the choice of outfld array dimension parameters becomes obvious. Or it could be declared with pcols, and initialized appropriately, and so be consistent with all of the other calls in this section of code (until a refactor happens).

[gold2718]

I second Pat's recommendation.

[ndkeen]

I was going to reply to this with a fix to see if both of you agreed.

But I'm also seeing a different problem in a call to oufld(). It's still an acme_dev test, and on cori-knl with DEBUG Intel, but it's a floating overflow and in clubb_intr.F90.

   call outfld( 'RCMINLAYER_CLUBB', rcm_in_layer*1000._r8,   pcols, lchnk )

000: MCT::m_Router::initp_: GSMap indices not increasing...Will correct
174: forrtl: error (73): floating divide by zero
174: Image              PC                Routine            Line        Source
174: acme.exe           000000000ABE4BD1  Unknown               Unknown  Unknown
174: acme.exe           000000000ABE2D0B  Unknown               Unknown  Unknown
174: acme.exe           000000000AB91414  Unknown               Unknown  Unknown
174: acme.exe           000000000AB91226  Unknown               Unknown  Unknown
174: acme.exe           000000000AAF9A89  Unknown               Unknown  Unknown
174: acme.exe           000000000AB060B9  Unknown               Unknown  Unknown
174: acme.exe           000000000A6C5C10  Unknown               Unknown  Unknown
174: acme.exe           0000000002A84B5A  clubb_intr_mp_clu        2456  clubb_intr.F90
174: acme.exe           0000000000CE7ACC  physpkg_mp_tphysb        2362  physpkg.F90
174: acme.exe           0000000000CC0661  physpkg_mp_phys_r        1008  physpkg.F90
174: acme.exe           000000000A75C0A3  Unknown               Unknown  Unknown
174: acme.exe           000000000A714D10  Unkn

Is it possible that it's still a similar issue: there are invalid values in the (unused part of) array that are multiplied by 1000?

[gold2718]

@ndkeen, This is very likely the same problem with the same solution.

[worleyph]

@ndkeen, I agree with @gold2718 . This is also another example of array arithmetic in a subroutine call (as well as computing with uninitialized values).

[ndkeen]

OK. Have not made progress on this -- was planning on doing it today.

[worleyph]

In case it is not clear (lots of suggestions above), my suggestion is to declare and use a work array, either

(a)

 real(r8), allocatable :: tmp_array(:,:)
 integer :: istat
 character(len=*), parameter :: routine = 'hetfrz_classnuc_cam_calc'
 ...
 allocate(tmp_array(ncols,pver), stat=istat)
 call alloc_err(istat, routine, 'tmp_array', ncols*pver)
 tmp_array = 0._r8
 ...
 tmp_array = niimm_bc(:ncols,:)+niimm_dst(:ncols,:)
 call outfld('NIMIX_IMM', tmp_array, ncols, lchnk) 
 (etc.)
 ...
 deallocate(tmp_array)

or
(b)

 real(r8) :: tmp_array(pcols,pver)
 ...
 tmp_array = 0._r8
 ...
 tmp_array = niimm_bc(:ncols,:)+niimm_dst(:ncols,:)
 call outfld('NIMIX_IMM', tmp_array, pcols, lchnk) 
 (etc.)

for each of the problematic outfld calls (any that invole array arithmetic in the call). (b) is simpler, but retains the style that @gold2718 doesn't like.

[singhbalwinder]

I am not sure how frequently these outfld calls are executed. If they are executed every time step (or often), won't the solution prescribed in (a) above fragment the memory with repeated "allocate" "deallocate" calls?

[mt5555]

I also vote for (b). for temp variables inside subroutines, it's much more efficient to create them on the stack than it is to allocate them. If the allocation is inside a threaded region (not sure if that is true here), allocate is especially slow.

[worleyph]

Only reason (a) requires allocate/deallocate is that ncols isn't available to be used for declarations. Probably not worth adding this as a parameter to the subroutine call just to enable this option.

[gold2718]

outfld calls are typically called from a threaded region
threading is over chunks and ncol is a property of the chunk
ncol is typically available as state%ncol (state is an input in almost every routine that calls outfld).

[worleyph]

@gold2718 - thanks. My mistake. I didn't look at the code closely enough.

[worleyph]

So, (a) should be replaced by (c)

 real(r8) :: tmp_array(state%ncol,pver)
  ...
  tmp_array = niimm_bc(:ncols,:)+niimm_dst(:ncols,:)
  call outfld('NIMIX_IMM', tmp_array, ncols, lchnk) 
  (etc.)

[ndkeen]

I started a branch ndk/cam/add-tmp_arrays-for-outfld and used tmp_array(:,:) to do the temporary calcs for several outfld() calls. I verified that at least 2 tests that were failing are now OK -- one on cori-knl, one on edison. I changed several outfld() calls in surrounding code, but certainly not all instances of outfld() throughout the code. I will do some more testing with this branch.

[ndkeen]

@gold2718 or @worleyph, can you look at those changes in this branch to see if this is what you had in mind? Obv I will remove the comments. I made the "tmp_array fix" to several outfld calls, even if they weren't causing an immediate problem.

Or I could make a PR and someone could review that instead. ?

[worleyph]

@ndkeen , I think that you want to use

 real(r8) :: tmp_array(state%ncol,pverp)

if you are going to use ncol in

 call outfld( 'RTP2_CLUBB',       tmp_array,               ncol, lchnk )

Otherwise you still have a pcol-dimensioned array but are telling outfld that this is an ncol-dimensioned array.

[worleyph]

and

 real(r8) :: tmp_array(state%ncol,pver) !ndk

in the other routine.

[ndkeen]

got it -- changed both of those and retesting

[ndkeen]

Update: I got caught up with other issues, but wanted to follow-through with this. I will try it again with master, make same changes, etc.

[ndkeen]

And after module & PE layout adjustments, one of the tests (SMS_D_Ld1.ne16_ne16.FC5ATMMOD) hit the same error on cori-knl. So another reason to get a PR going.

[ndkeen]

I had
real(r8) :: tmp_array(pcols,pverp)
instead of
real(r8) :: tmp_array(state%ncol,pverp)

in both places

trying to see if that fixes issues we are now seeing