NaNs in tracer fields in the first timestep of simulations using meshes with cavities

[svenharig]

This problem has been reported every now and again, I remember it first by @JanStreffing towards the end of last year. Especially in the CORE3 development context while trying out different configurations and distributions it occurred frequently. In the very first timestep there might be NaNs occurring in temperature and/or salinity (probably also other tracers) resulting in the model to blow up. It only happens in meshes with cavities, and NaNs occur only in mesh partitions including a fraction of a cavity. It may depend on the actual partition used - running with a distribution for one node might work perfectly while using several nodes and thus partitioning for that size may fail.

Since this affects the test simulations in the CORE3 mesh development I tried to track it down in FESOM2.7.1 and in the meantime I have a suspicion it might be related to tracer advection routines. One of the first NaN occurrences I could find, happens in the routine oce_tra_adv_fct in the file oce_adv_tra_fct.F90.

In that routine a large auxiliary field AUX is used and seemingly filled and evaluated inconsistently in some cases.

It is filled in a loop over elements

    do elem=1, myDim_elem2D
        enodes=elem2D_nodes(:,elem)
        nu1 = ulevels(elem)
        nl1 = nlevels(elem)

        !$ACC LOOP VECTOR
        do nz=nu1, nl1-1
            AUX(1,nz,elem)=max(fct_ttf_max(nz,enodes(1)), fct_ttf_max(nz,enodes(2)), fct_ttf_max(nz,enodes(3)))
            AUX(2,nz,elem)=min(fct_ttf_min(nz,enodes(1)), fct_ttf_min(nz,enodes(2)), fct_ttf_min(nz,enodes(3)))
        end do
       . . .
    end do

and later on evaluated in a loop over nodes

    do n=1, myDim_nod2D
       nu1 = ulevels_nod2D(n)
       nl1 = nlevels_nod2D(n)
       !___________________________________________________________________
       !$ACC LOOP VECTOR
       do nz=nu1,nl1-1
          ...
          do elem=2,nod_in_elem2D_num(n)
              tvert_max(nz, n) = dmax1(tvert_max(nz, n), AUX(1,nz, nod_in_elem2D(elem,n)))
              tvert_min(nz, n) = dmin1(tvert_min(nz, n), AUX(2,nz, nod_in_elem2D(elem,n)))
         end do
      end do

    ....

and in the case with NaN trouble that I investigate, the level range covered by ulevels(elem) and ulevels_nod2D(n) is quite different, so in this particular case AUX in the problematic mesh node is set in the level range 19 to 25 but evaluated in the range 7 to 25, so here the result in the nonoverlapping part depends on prior settings in AUX which may randomly put NaNs into tvert_max or tvert_min and from here the trouble spreads in my case. This could explain the dependency on the partition and the restriction of the problem to cavity regions where the level range varies.

My workaround for the moment is to fill the proper range by replacing the ulevels(elem) by the max level range over the nodes in the elem:

    do elem=1, myDim_elem2D
        enodes=elem2D_nodes(:,elem)
        !sh nu1 = ulevels(elem)
        !sh nl1 = nlevels(elem)
        ! ====== WORKAROUND ======
        nu1=minval(ulevels_nod2D(enodes(:)))
        nl1=maxval(nlevels_nod2D(enodes(:)))
        !=========================
       ...

    end do

which avoids the random entries but of course I cannot be sure if the resulting values in AUX are meaningful. In any case with the workaround this problem does no longer occur, the test simulations run so far or blow up for other reasons ;) Maybe I was looking here at some weird exceptional case but please take a look if this discrepancy might happen in other setups with cavities as well. My ongoing experiments are located on albedo in the folder

/albedo/work/user/sharig/FESOM271_CORE3_test/run.CORE3.cav.jig.dev/

[mandresm]

We've done a git blame on those lines and many of the changes are from the NatESM people. Also some changes come from @patrickscholz and others from @dsidoren (related to OpenMP implementation).

@PengyangSong, have you or any of your collaborators experienced this issue before? If so did you found a solution for it?

[PengyangSong]

We've done a git blame on those lines and many of the changes are from the NatESM people. Also some changes come from @patrickscholz and others from @dsidoren (related to OpenMP implementation).

@PengyangSong, have you or any of your collaborators experienced this issue before? If so did you found a solution for it?

Yes, I will soon push a solution

[mandresm]

Thanks @PengyangSong !

[PengyangSong]

Hi @svenharig @JanStreffing @GMHuettner
See commit 7ed3c49 and try out branch fix/cavity_surf_nan
Make sure it work before merging it back to the main.
Pengyang

[patrickscholz]

@svenharig , @PengyangSong : Although Pengyangs solution will work it puts quite some unnecessary if conditions. I think what should work better is to do the same as in ...

fesom2/src/oce_adv_tra_fct.F90

Lines 157 to 165 in 787d046

	!$ACC END LOOP
	if (nl1<=nl-1) then
	!$ACC LOOP VECTOR
	do nz=nl1,nl-1
	AUX(1,nz,elem)=-bignumber
	AUX(2,nz,elem)= bignumber
	end do
	!$ACC END LOOP
	endif

but for the range over the cavity ...

if (nu1>1) then
        !$ACC LOOP VECTOR
        do nz=1,nu1-1
            AUX(1,nz,elem)=-bignumber
            AUX(2,nz,elem)= bignumber
        end do
        !$ACC END LOOP
endif

... than the rest of the code should work as it is?!