[MITgcm-devel] (not so) funny things happen in seaice_lsr and pickups

[Jinlun Zhang]

I fully agree that discretizations of metric terms may break symmetry, 
even not mathematically. Good luck, Martin :-) .
Jinlun

Samar Khatiwala wrote:
> Hi Martin et al
>
> Its been fun following this conversation passively, but the last 
> couple of emails make me want to say something.
>
> I think Martin has it right that improper treatment of the metric 
> terms can break the symmetry. More generally its important to
> keep in mind that symmetry or self-adjointness depends on the choice 
> of inner product. This is where the metric terms are
> critical. I don't pretend to know what your sea ice equations are, but 
> it would be a mistake to think that the nice properties of
> the continuum operators will somehow carry over to their discrete 
> matrix representations without a lot of special care. As an
> example think of the shallow water equations. These are self-adjoint 
> (in a particular inner product), but most standard discretizations
> will break this (and there are many papers in which you find this). 
> The correct discretization (derived by Platzman) is not at all
> obvious (at least to me). (Martin: since I know you will enjoy 
> 'wasting' your time on this :-), let me know if you want my matlab code
> to play with.)
>
> Samar
>
> On Mar 10, 2009, at 1:36 PM, Jinlun Zhang wrote:
>
>> Martin Losch wrote:
>>> The question marathon continues ...
>>> Jinlun, my observation is, that the metric terms cause all the 
>>> trouble (no problems so far in cubed sphere integration of any 
>>> resolution, nor on cartesian grid).  I have the suspicion that the 
>>> symmetry of the laplacian operator is broken, because the metric 
>>> terms are distributed unevenly between left and right hand side of 
>>> the u and v equations. In order to ensure the diagonal dominance of 
>>> the coefficient matrix, would it be possible to move all metric 
>>> terms to the rhs? Just wondering ...
>> Hi Martin, equations governing sea ice physics would be incorrect if 
>> metric terms are moved to rhs. I doubt metric terms would break  
>> symmetry under any conditions (I have not seen that). This is because 
>> a rotation of a grid (or a change of coordinate system) would not 
>> change the properties of a matrix.
>>>
>>> BTW, the system seems to be stable when I have a non-zero ZMIN.
>> Non-zero ZMIN is not right and should not be used. With non-zero 
>> ZMIN, the system is heavily damped. I doubt non-symmetry caused 
>> trouble would be fixed by non-zero ZMIN, it would eventually blow up 
>> no matter what ZMIN is.
>>
>> Jinlun
>>>
>>> Martin
>>>
>>> On Mar 9, 2009, at 5:54 PM, Jinlun Zhang wrote:
>>>
>>>> Hi Martin,
>>>>
>>>> As you said, for B-grid, the value does not matter because of the 
>>>> non-slip bc's. I am not sure about C-grid non-slip bc's. If you 
>>>> need that value at the boundary for C-grid, it seems that your 
>>>> approach probably makes better sense.
>>>>
>>>> Jinlun
>>>>
>>>> Martin Losch wrote:
>>>>> Hi again,
>>>>>
>>>>> I have another question for Jinlun. I am now quite confident, that 
>>>>> I have found all (serious) problems in the metric terms. However, 
>>>>> my code still blows up. Could the discretization of \eta be a 
>>>>> problem? On Z-points ("bottom left" corner of C-grid cell, on a 
>>>>> B-grid this is where the velocities are), I average eta like this:
>>>>>       DO j=1,sNy+1
>>>>>        DO i=1,sNx+1
>>>>>         etaMeanZ (I,J,bi,bj) =
>>>>>    &         (           ETA (I,J  ,bi,bj)  + ETA (I-1,J  ,bi,bj)
>>>>>    &         +           ETA (I,J-1,bi,bj)  + ETA (I-1,J-1,bi,bj) )
>>>>>    &         / MAX(1.D0,maskC(I,J,  k,bi,bj)+maskC(I-1,J,  k,bi,bj)
>>>>>    &         +          maskC(I,J-1,k,bi,bj)+maskC(I-1,J-1,k,bi,bj) )
>>>>>        ENDDO
>>>>>       ENDDO
>>>>> In the B-grid code I find this:
>>>>>       DO j=1-Oly+1,sNy+Oly
>>>>>        DO i=1-Olx+1,sNx+Olx
>>>>>         ETAMEAN(I,J,bi,bj) =QUART*(
>>>>>    &          ETA(I,J-1,bi,bj) + ETA(I-1,J-1,bi,bj)
>>>>>    &         +ETA(I,J  ,bi,bj) + ETA(I-1,J  ,bi,bj))
>>>>>         ZETAMEAN(I,J,bi,bj)=QUART*(
>>>>>    &          ZETA(I,J-1,bi,bj) + ZETA(I-1,J-1,bi,bj)
>>>>>    &         +ZETA(I,J  ,bi,bj) + ZETA(I-1,J  ,bi,bj))
>>>>>        ENDDO
>>>>>       ENDDO
>>>>> This means that \eta/\zeta on the boundaries are only half as 
>>>>> large as in my C-grid discretization (however, these values may 
>>>>> actually never be used because of the no-slip BCs on the B-grid). 
>>>>> So I am wondering, what the "correct" value of \eta on the 
>>>>> boundary is. I observe, that code with
>>>>>       DO j=1,sNy+1
>>>>>        DO i=1,sNx+1
>>>>>         etaMeanZ (I,J,bi,bj) = 0.25 _d 0 *
>>>>>    &         (           ETA (I,J  ,bi,bj)  + ETA (I-1,J  ,bi,bj)
>>>>>    &         +           ETA (I,J-1,bi,bj)  + ETA (I-1,J-1,bi,bj) )
>>>>>        ENDDO
>>>>>       ENDDO
>>>>> tends to be stable (at least in the shorter test that I made). 
>>>>> Jinlun, what's your opinion on this?
>>>>>
>>>>> Martin
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>>