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

Tue Mar 10 13:56:43 EDT 2009

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
>
>
> -- 
>
> Jinlun Zhang
> Polar Science Center, Applied Physics Laboratory
> University of Washington, 1013 NE 40th St, Seattle, WA 98105-6698
>
> Phone: (206)-543-5569;  Fax: (206)-616-3142
> zhang at apl.washington.edu
> http://psc.apl.washington.edu/pscweb2002/Staff/zhang/zhang.html
>
>
>
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel