[MITgcm-support] instability issue with sponge
Jody Klymak
jklymak at uvic.ca
Wed Feb 8 14:18:18 EST 2012
Hi John,
This all looks fine, though your diffusivities are largish. You might simply try a smaller time step, or look up the CFL criteria for diffusivities.
> While I've got you guys on the line I have another quick question about the local density. At one point it occurred to me that it would be cleaner if I forced the body of water in a legit m=1 mode instead of the sinusoidal profile used in
> internal_wave/code/obcs_calc.F
> which really only holds for no stratification.
I don't understand what you are trying to say here. The linear mode-1 response in constant-N ocean is a half-cosine in horizontal velocity and pressure, and a half sine wave in vertical velocity. Are you worried about non-linearities for some reason?
> To do this I need to weight by the local density and thought I could extract the density by invoking the routine
> FIND_RHO_SCALAR
> I defined a few local variables in code/obcs_calc.F just to be safe and to make it easy to print them out
>
> _RL origMom,vbar,nrmConst
> _RL plocJGP,rhoRefJGP(Nr),salJGP
>
> then added this call
>
> salJGP=35.
> DO k=1,Nr
> pLocJGP = -rhoConst*rC(k)*gravity
> CALL FIND_RHO_SCALAR(
> I tRef(k), salJGP, pLocJGP,
> O rhoRefJGP(k), myThid )
> ENDDO
Mine is in seawater.F, though I see it has moved to FindRho.F in more recent checkpoints, but:
rholoc = rhoNil*(
& sBeta *(sLoc-sRef(1))
& -tAlpha*(tLoc-tRef(1))
& ) + rhoNil
so, I'd expect your density to be:
rhoNil*(1 -0.2e-4*(tRef(k)-tRef(1))
rhoNil is usually 999.98. So I guess its a little strange that your deep densities are not greater than 1000, but at the surface it makes sense.
This is pretty confusing, because for find_Rho_2d, in find_rho.F, the definition is different:
dRho = rhoNil-rhoConst
DO j=jMin,jMax
DO i=iMin,iMax
rhoLoc(i,j)=rhoNil*(
& sBeta*(sFld(i,j)-refSalt)
& -tAlpha*(tFld(i,j)-refTemp) )
& + dRho
where refSalt and RefTemp are the sRef(k) and tRef(k).
Other than the density issue, I'd be interested to see a plot of V and T when you start having your instabilities (around day 29), particularly in your sponge, and just beyond.
I'd also be interested to see data.obcs
Cheers, Jody
>
> The tRef values are OK, the salinity is OK, the pressure values (pLocJGP) are OK, but the funny thing is that the densities returned by the subroutine are always less than 1000, as if the routine were refusing to use the salinity in the calculation.
>
>
>
>
> Thanks very much,
> John
>
>
>
>
> Begin forwarded message:
>
>> From: Jody Klymak <jklymak at uvic.ca>
>> Date: February 8, 2012 5:50:34 AM AKST
>> To: mitgcm-support at mitgcm.org
>> Subject: Re: [MITgcm-support] instability issue with sponge
>> Reply-To: mitgcm-support at mitgcm.org
>>
>> Hi John,
>>
>> On Feb 8, 2012, at 0:52 AM, Martin Losch wrote:
>>
>>>
>>> the sponge layer can indeed cause perturbation. After all you are restoring your model to prescribed boundary values (in your case OBNt, OBNs, OBNu, OBNv) over the sponge layer with a restoring time scale that is determined by the appropriate parameters. Unless thes OBNt/s/u/v are not exactly the same as the corresponding model values (i.e. in your case with no forcing with uniform salinity=35, temperature with N2=10^-6, u=v=0), the sponge layer will make model move. I have the impression that this is not the primary problem.
>>
>> Just to add to that, by default the sponge sponges to Sref and Tref, so if your temperature profile is set by TRef you should be fine. But if TRef is a constant, and your temperature profile set by setting hydrogThetaFile, then obviously you'll sponge your density differently, and you'll have to modify obcs_sponge.F appropriately.
>>
>> But this should be seen if you look at your T output.
>>
>> Otherwise, its hard to see how w/ no forcing you are getting instabilities in the sponge layer, though I have to admit I've rarely run for 30 days. re you sure you've set your timescales to something sensible? a couple of years ago Nicolas Grisouard suggested:
>>
>>> Say you have a forcing frequency w.
>>> - number of points: as long as the horizontal wavelength of mode-1 internal waves of frequency w,
>>> - [UV]relaxobcsinner = (2*pi/w)/1000,
>>> - [UV]relaxobcsbound = 2*pi/w.
>>>
>>
>>
>> I think this is a bit draconian, and usually get away with much shorter sponges w/o noticeable reflections, but its a starting point. You may also consider telescoping the sponge layer so dy gets larger w/in it. No need to waste a lot of computation cells on a-physical space. Just make sure your telescope is relatively gradual.
>>
>>> The stability issue is more difficult to answer without additional information. With your parameters the maximum allowed velocity is around .5*dy/dt ~ 3m/s before the CFL criterion is violated. I assume that's what is happing. Other parameters affect this velocity, i.e. the amount of viscosity and diffusivity you prescribe, and the choice of advection scheme.
>>
>>
>> These could be problems, though I'm not sure how the bare advection scheme could exceed CFL if there is no forcing. But maybe the viscosity and diffusivity.
>>
>> I guess even in the no-forcing case with a well-prescribed sponge, there will eventually be motion due to diffusivity at the surface and seafloor. You'll start to make mixed layers, but the sponge will be towards stratified conditions. That will drive some secondary circulation. But even at 30 days, your vertical diffusivity would have to be pretty high for this to be noticeable.
>>
>> Cheers, Jody
>>
>>
>>
>>> I suggest that you post your "code" directory and your "input" directory (the namelist "data*" files)
>>>
>>> Martin
>>>
>>> PS:
>>> - unless you have Nx>1 and a wall in y-direction, you are not dealing with a channel, but with a periodic domain. I am not sure if you want that.
>>>
>>> - For performance reasons, I suggest that you use Nx=300 (on an f-plane this makes no difference). Some compilers can optimize the innermost loops (which are i-loops), but if Nx=1 (or, 2) there isn't much to optimize.
>>>
>>> On Feb 8, 2012, at 4:25 AM, John Pender wrote:
>>>
>>>> Greetings, All -
>>>>
>>>> I have been tasked with crawling up the MITgcm learning curve using a 2D ocean channel as a vehicle. I have been running into a stability issue that I haven't managed to fix and could use some advice. The broad features of the channel are:
>>>>
>>>> 2D channel running north/south for 600 km
>>>> f0 and beta are both zero for now
>>>> uniform depth of 4000m
>>>> uniform salinity of 35
>>>> stratified temperature with N2 = 10^-6
>>>> Ny = 300
>>>> Nz = 20
>>>> deltaT = 300 sec
>>>>
>>>> The idea is to force the southern end in a specific way by enabling OBCs and modifying obcs_calc.F to suit. The forcing is baroclinic so it creates a progression of cells (m=1) which propagate northward in a slow, stately fashion. The time window of this experiment is very long - on the order of 100 days - so the cells have more than enough time to reach the northern end of the channel. I don't want the waves reflecting back so I have installed a sponge at the northern end. I've used many permutations of the sponge thickness and the inner/outer parameters and it is indeed crushing the baroclinic wave train as advertised. For a while.
>>>>
>>>> The problem is stability. I need this experiment to run for, as I say, 100 days but the model is self terminating at around 30 days because of wildly erratic large numbers. What's interesting is that I have turned the forcing off and the experiment still self terminates in about the same number of days. I would have expected very little (if any) motion in the unforced case, and that is indeed what happens when I also turn off the sponge. The model runs endlessly when the sponge is gone and there is no forcing. When the sponge is on, however, motion is generated for no discernible reason at the inner edge of the sponge. The motion is very small at first (it appears after about 3 time steps) and if it stayed small I wouldn't worry about it. The problem is that the motion gets larger and larger over time (ultimately crashing the model) so there must be an additive feature or perhaps even a nonlinear feedback effect that I obviously don't understand.
>>>>
>>>> It has been suggested to me that turning on the sponge somehow changes the density of the water within the sponge's extent (though obviously this shouldn't happen), which then forces motion. I don't actually know if this is the case because the model doesn't much *like* keeping me appraised of the local density, but this explanation is definitely consistent with my observations of v, w, and eta. I would think this would produce a finite amount of motion (assuming energy's conserved) so it doesn't really explain the steady growth in mechanical energy over the 30 days or so it takes to run the numbers into the stratosphere.
>>>>
>>>> Thanks very much,
>>>> John
>>>>
>>>>
>>>> _______________________________________________
>>>> MITgcm-support mailing list
>>>> MITgcm-support at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>>>
>>>
>>> _______________________________________________
>>> MITgcm-support mailing list
>>> MITgcm-support at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>>
>> --
>> Jody Klymak
>> http://web.uvic.ca/~jklymak/
>>
>>
>>
>>
>>
>> _______________________________________________
>> MITgcm-support mailing list
>> MITgcm-support at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-support
>
> _______________________________________________
> MITgcm-support mailing list
> MITgcm-support at mitgcm.org
> http://mitgcm.org/mailman/listinfo/mitgcm-support
--
Jody Klymak
http://web.uvic.ca/~jklymak/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20120208/dcc8bbd8/attachment-0001.htm>
More information about the MITgcm-support
mailing list