[MITgcm-devel] more seaice blues for the adjoint

[Martin Losch]

Hi Jinlun,

I check with the b-code and at first I got zero motion!!! But why?  
the reason is not the grid, but this line:
zmin = 4e8
if I set this to
zmin = 0
as we have agree a few days/weeks ago, I do get motion.

So it's not the C-grid, but this lower limit for zeta, that makes  
zeta = constant if delta = 0, and thus P = zeta*2*delt2 = const, so  
no dP/dx -> no motion. if zmin = 0, then zeta is not constant,  
neither P and we have motion.
so if zmin = 0, then I need to do something else to ensure P =  
constant if delta = 0., eg. P = P * delta/delt2

Martin

On 8 Jun 2007, at 08:38, Jinlun Zhang wrote:

> Hi Martin,
> Before your fix, is it occurring with both LSR and EVP C-grid  
> codes? Also, would it be possible for you to try the B-grid that I  
> put in originally (if it is still one of the options)? It is OK for  
> not trying if you don't have time. I remember I tested this thing  
> with B-grid years ago without finding similar problem, but I am not  
> sure now after so many years without dealing with it again. I just  
> hope this is not one of the differences between B and C. Your  
> explanation certainly indicates that this thing does not depend on  
> grid. But a test on B-grid would clear any residual doubts.
> Thanks, Jinlun
>
>
> Martin Losch wrote:
>> Chris,
>> when there is no strain (in the absence of any forcing), the ice  
>> should not move, right? Unless you expect it to spread like a pile  
>> of sand with time, but I don't know what the time scale for that  
>> would be, very slow.
>> If delta = 0, we reset it to delta_min, because we have to devide  
>> by it zeta = press/(2*delta), this effictively makes zeta <=  
>> zeta_max = press/(2*delta_min), after capping zeta, we recompute  
>> press_replace = zeta*2*delta, if delta has been replaced by  
>> delta_min, this give non zero pressure and thus a forcing due to  
>> internal stress. That's the numerical issue. I am not shure that I  
>> solved it very well (the adjoint breaks), but this way the ice  
>> does not move in the absence of forcing.
>>
>> Martin
>> On 7 Jun 2007, at 18:05, chris hill wrote:
>>
>>> Martin,
>>>
>>>  Quick ice question.
>>>
>>> Is it just a numerical issue i.e. is it clear what the underlying  
>>> equations are, since ice is neither a classical fluid or a  
>>> classical solid. What time scale does it spread on?
>>>
>>> Chris
>>> Martin Losch wrote:
>>>> No, but it's a numerical issure. I did not pay any attention to  
>>>> this before and maybe this was no problem in the original code  
>>>> that you supplied but if you have zero strain, Delta = 0, and  
>>>> this is replaced by some Delta_min (SEAICE_EPS in the code) in  
>>>> order to avoid division by zero (and infinite viscosities). P  
>>>> however is still non-zero and divergence(stress) end up being  
>>>> (dP/dx, dP/dy) .NE. 0, so that you have a forcing of down the  
>>>> slope of P. This does not make sense (why would ice  
>>>> spontaneously spread?), so that the pressure is replaced by P =  
>>>> P*Delta/max(Delta,SEAICE_EPS), so that in the (rare) case of no  
>>>> strain P = zeta = eta = 0 and thus no forcing by stress. Makes  
>>>> sense to me (and is really only relevant in idealized test case,  
>>>> just like the spurious motion of sigma coordinates vs. z- 
>>>> coordinates in ocean models).
>>>> Hope I did not misunderstand anything here.
>>>> Martin
>>>> On 7 Jun 2007, at 13:40, Jinlun Zhang wrote:
>>>>> Hi Martin,
>>>>> Would this be due to some kind of finite differencing problem?
>>>>> Jinlun
>>>>>
>>>>> Martin Losch wrote:
>>>>>> Hi Patrick,
>>>>>>
>>>>>> I have found (or rather, was pointed to) a problem with the  
>>>>>> seaice solvers: The start to move spontaneously (in the  
>>>>>> absence of forcing), if the sea ice distribution is NOT uniform.
>>>>>> I have implemented a fix but this will cause problems with the  
>>>>>> adjoint: I need terms like
>>>>>> SQRT(deltaC), which used to be SQRT(MAX 
>>>>>> (deltaC,SEAICE_EPS_SQ)), so that the derivate code will be  
>>>>>> involve 1/sqrt(deltaC). Should I put this into #ifdefs?
>>>>>>
>>>>>> Martin
>>>>>
>>>>> _______________________________________________
>>>>> MITgcm-devel mailing list
>>>>> MITgcm-devel at mitgcm.org
>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>> _______________________________________________
>>>> MITgcm-devel mailing list
>>>> MITgcm-devel at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>
> _______________________________________________
> MITgcm-devel mailing list
> MITgcm-devel at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-devel