[MITgcm-support] instability issue with sponge
John Pender
jgpender at alaska.edu
Wed Feb 8 13:56:12 EST 2012
Thanks very much Martin and Jody for your prompt replies.
It occurred to me a week or so ago that there might be a discrepancy between the reference temperature values and the initial temperature values so
1) I use gendata.m to create a file of the reference temperature data
2) I read this file in to the model rather than typeset approximate numbers in the input/data file
3) I do not read in a hydrogThetaFile, which should force the model to use the reference temperatures as a starting point
Sref is fixed at 35., and I've double checked a few times. Reference velocities are all zero.
I hadn't seen the timescale guidelines you mention for the sponge parameters and will try them later today to see if they help.
>> 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.
The viscosity and diffusion constants come from the input/data file in internal_wave
# ====================
# | Model parameters |
# ====================
#
# Continuous equation parameters
&PARM01
# Tref = 0.9684, 0.8665, 0.7645, 0.6626,
# 0.5607, 0.4587, 0.3568, 0.2548,
# 0.1529, 0.0510, -0.0510, -0.1529,
# -0.2548, -0.3568, -0.4587, -0.5607,
# -0.6626, -0.7645, -0.8665, -0.9684,
Treffile='Trefvar',
sRef= 20*35.,
viscAz=1.E-3,
viscAh=1.E-2,
no_slip_sides=.FALSE.,
no_slip_bottom=.FALSE.,
diffKhT=1.E-2
diffKzT=1.E-3,
f0=0.0,
beta=0.E-11,
eosType='LINEAR',
tAlpha=2.E-4,
sBeta =0.E-4,
gravity=9.81,
implicitFreeSurface=.TRUE.,
exactConserv=.TRUE.
nonHydrostatic=.FALSE.,
hFacMin=0.2,
implicSurfPress=0.5,
implicDiv2DFlow=0.5,
nonlinFreeSurf=3,
hFacInf=0.2,
hFacSup=1.8,
saltStepping=.FALSE.,
#- not safe to use globalFiles in multi-processors runs
#globalFiles=.TRUE.,
readBinaryPrec=64,
writeBinaryPrec=64,
writeStatePrec=64,
&
# Elliptic solver parameters
&PARM02
cg2dMaxIters=1000,
cg2dTargetResidual=1.E-13,
cg3dMaxIters=400,
cg3dTargetResidual=1.E-13,
&
# Time stepping parameters
&PARM03
nIter0=0,
nTimeSteps=15000,
deltaT=300.,
abEps=0.1,
pChkptFreq=0.,
chkptFreq=0.,
dumpFreq=3600.,
monitorFreq=2500.,
monitorSelect=2,
&
# Gridding parameters
&PARM04
usingCartesianGrid=.TRUE.,
delXfile='delXvar',
delYfile='delYvar',
delRfile='delZvar',
# delY=5.E3,
# delZ=20*100.,
&
# Input datasets
&PARM05
# hydrogThetaFile='T.init',
bathyFile='topog.flat',
&
*************
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. 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
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-support/attachments/20120208/563523d0/attachment-0001.htm>
More information about the MITgcm-support
mailing list