Diagnostics don't add up!

[rduran]

I wrote a simple matlab code to check if momentum diagnostics add to zero, turns out that the u-momentum does and v-momentum do not. Here is the code I used so you can see it is not a bug (typo or so) from my part:

Code: Select all

for kk=1:tot
file=files{kk};
ncload(file)
vmom=makenan(v_accel-(v_cor+v_hadv+v_vadv+v_prsgrd+v_hvisc+v_vvisc));
umom=makenan(u_accel-(u_cor+u_hadv+u_vadv+u_prsgrd+u_hvisc+u_vvisc));
V(kk,:)=lims(vmom);
U(kk,:)=lims(umom);
end

the program lims gives me the absolute extrema of any array. the program makenan converts Fill_Value to NaN.

I ckecked the 2D momentum as well and its looks better but something still seems to be wrong, here are the corresponding code lines for scrutiny

Code: Select all

vmom=makenan(vbar_accel-(vbar_cor+vbar_hadv+vbar_bstr+vbar_prsgrd+vbar_hvisc+vbar_sstr));
umom=makenan(ubar_accel-(ubar_cor+ubar_hadv+ubar_bstr+ubar_prsgrd+ubar_hvisc+ubar_sstr));
V(kk,:)=lims(vmom);
U(kk,:)=lims(umom);

Finally I plot each with

Code: Select all

    
semilogy(time,V(:,1),'ro',time,V(:,2),'rx',time,U(:,1),'bo',time,U(:,2),'bx')
legend('MaxV','MinV','MaxU','MinU','Location','Best')

EDIT: Because I am interested in the order of magnitude and not the sign, and so I can plot in a Log plot, before plotting I do

Code: Select all

U=abs(U);
V=abs(V);

end edit===

the resulting figures are attached.

[pmaccc]

I just tried adding up the momentum diagnostics in one of my runs and couldn't find any obvious error, or any difference between u and v momentum budgets. Are you using "nudging to climatology" in your domain? I believe this would add an extra term to the momentum budgets that would not be in the assumed balance.

[rduran]

Parker,

Thank you very much for taking a look at this!

Are you using "nudging to climatology" in your domain?

I do not have M2CLM_NUDGING nor M3CLM_NUDGING defined (I found that with daily boundary conditions there was no need for nudging).

This was done with ROMS 3.6 revision 571, but I am getting similar problems with revision 656. I have the feeling a bug may have been introduced in the writing/accumulating of diagnostic terms, and that it was introduced some time not too long ago. Perhaps, if your ROMS version is different, I can try downloading the same revision you used (for the diagnostics you described above) to test my hypothesis.

My simulation is realistic with a domain from 41 to 48 North and from -129 to -123.4 West. The regular output look just fine, however the diagnostics don't.

I am attaching some plots where the order of magnitude for each term can be seen.
The first plot is a time series of the mean and standard deviations for the sum of the u and v momentum terms.
The second two plots are a time series of U and V momentum terms at a single point (away from boundaries) this one looks reasonable, and thus misleading.

The last two plots are a time series of the zonally-, meridionally- and depth-averaged momentum terms for u and for v.
Something is clearly wrong in the v-momentum, while u-momentum seems to be in a better shape.

Where could those huge vertical viscosity values be coming from?
I am including the code with which the plots were produced for scrutiny.


Here are my CPP options from the DIA files themselves.

Code: Select all

:CPP_options = "ORE2005, ANA_BSFLUX, ANA_BTFLUX, ASSUMED_SHAPE, AVERAGES, DIAGNOSTICS_UV, DJ_GRADPS, DOUBLE_PRECISION, MASKING, MIX_ISO_TS, MIX_GEO_UV, MPI, MY25_MIXING, NONLINEAR, NONLIN_EOS, POWER_LAW, PROFILE, K_GSCHEME, RADIATION_2D, !RST_SINGLE, SALINITY, SOLVE3D, SPLINES, TS_U3HADVECTION, TS_C4VADVECTION, TS_DIF2, TS_PSOURCE, UV_ADV, UV_COR, UV_U3HADVECTION, UV_C4VADVECTION, UV_QDRAG, UV_PSOURCE, UV_VIS2, VAR_RHO_2D" ;

Here is the code that computes and plots the values in the attached figures:

Code: Select all

zi=35; %z52(loi,lai,zi)=-128.4202 meters depth
loi=80;% lon(loi)=-126.3667 West
lai=190;%lat(lai)=45.5831 North
kk=1;
for kk=1:length(loopvec)
file=files{loopvec(kk)};
ncload(file)

v_accel=makenan(v_accel); % CONVERT FILL_VALUES TO NAN
v_cor=makenan(v_cor);
v_hadv=makenan(v_hadv);
v_vadv=makenan(v_vadv);
v_prsgrd=makenan(v_prsgrd);
v_hvisc=makenan(v_hvisc);
v_vvisc=makenan(v_vvisc);

u_accel=makenan(u_accel);
u_cor=makenan(u_cor);
u_hadv=makenan(u_hadv);
u_vadv=makenan(u_vadv);
u_prsgrd=makenan(u_prsgrd);
u_hvisc=makenan(u_hvisc);
u_vvisc=makenan(u_vvisc);

vmom=(v_accel-(v_cor+v_hadv+v_vadv+v_prsgrd+v_hvisc+v_vvisc)); %SUM OF MOMENTUM TERMS
umom=(u_accel-(u_cor+u_hadv+u_vadv+u_prsgrd+u_hvisc+u_vvisc));

vpaccel(kk)=v_accel(zi,lai,loi); %pointwise values for v-momentum
vpcor(kk)=v_cor(zi,lai,loi);
vphadv(kk)=v_hadv(zi,lai,loi);
vpvadv(kk)=v_vadv(zi,lai,loi);
vpprs(kk)=v_prsgrd(zi,lai,loi);
vphvisc(kk)=v_hvisc(zi,lai,loi);
vpvvisc(kk)=v_vvisc(zi,lai,loi);

upaccel(kk)=u_accel(zi,lai,loi); %pointwise values for u-momentum
upcor(kk)=u_cor(zi,lai,loi);
uphadv(kk)=u_hadv(zi,lai,loi);
upvadv(kk)=u_vadv(zi,lai,loi);
upprs(kk)=u_prsgrd(zi,lai,loi);
uphvisc(kk)=u_hvisc(zi,lai,loi);
upvvisc(kk)=u_vvisc(zi,lai,loi);


%V
v_cort(kk,:)=lims2(v_cor); %GET MAX, MIN, MEAN STD, MEDIAN SEE lims2 function below
v_prsgrdt(kk,:)=lims2(v_prsgrd);
v_accelt(kk,:)=lims2(v_accel);
v_hadvt(kk,:)=lims2(v_hadv);
v_vadvt(kk,:)=lims2(v_vadv);
v_hvisct(kk,:)=lims2(v_hvisc);
v_vvisct(kk,:)=lims2(v_vvisc);
%U
u_cort(kk,:)=lims2(u_cor);
u_prsgrdt(kk,:)=lims2(u_prsgrd);
u_accelt(kk,:)=lims2(u_accel);
u_hadvt(kk,:)=lims2(u_hadv);
u_vadvt(kk,:)=lims2(u_vadv);
u_hvisct(kk,:)=lims2(u_hvisc);
u_vvisct(kk,:)=lims2(u_vvisc);

%balance
V(kk,:)=lims2(vmom);
U(kk,:)=lims2(umom);

end %for kk

time=1:size(u_accelt,1); %generic time

%this is the plot for the mean of each u-momentum term
hh=plot(time,u_accelt(:,3),time,u_hadvt(:,3),time,u_prsgrdt(:,3),time,u_vvisct(:,3),time,u_cort(:,3),time,u_hvisct(:,3),time,u_vadvt(:,3)); 
set(hh,'MarkerSize',8,'Linewidth',1.5)
legend('MeanU accel','Mean U hadv','Mean U prsgrd','Mean U vvisct','Mean Ucor','Mean U hvisc','Mean U vadv','Location','Best')

%this is the plot for the mean of each v-momentum term
hh=plot(time,v_accelt(:,3),time,v_hadvt(:,3),time,v_prsgrdt(:,3),time,v_vvisct(:,3),time,v_cort(:,3),time,v_hvisct(:,3),time,v_vadvt(:,3)); 
set(hh,'MarkerSize',8,'Linewidth',1.5)
legend('MeanV accel','MeanV hadv','MeanV prsgrd','MeanV vvisct','MeanV cor','MeanV hvisc','MeanV vadv','Location','Best')

%this is plot for pointwise u-momentum
hh=plot(time,upaccel,time,uphadv,time,upprs,time,upvvisc,time,upcor,time,uphvisc,time,upvadv); 
set(hh,'MarkerSize',8,'Linewidth',1.5)
legend(' PointU accel','PointU hadv','PointU prsgrd','PointU vvisct','PointU cor','PointU hvisc','PointU vadv','Location','Best')

%this is plot for pointwise v-momentum
hh=plot(time,vpaccel,time,vphadv,time,vpprs,time,vpvvisc,time,vpcor,time,...
    vphvisc,time,vpvadv); 
set(hh,'MarkerSize',8,'Linewidth',1.5)
legend(' PointV accel','PointV hadv','PointV prsgrd','PointV vvisct','PointV cor','PointV hvisc','PointV vadv','Location','Best')


% this plots the mean and the standard deviation of the sum of v-momentum and u-momentum %terms
U=abs(U);
V=abs(V);
hh=semilogy(time,V(:,3),'ro',time,V(:,4),'rx',time,U(:,3),'bo',time,U(:,4),'bx');
set(hh,'MarkerSize',8,'Linewidth',1.5)
legend('MeanV','StdV','MeanU','StdU','Location','Best')

Here is lims2.m function used to get mean and standard deviation :

Code: Select all

function all=lims2(x)
% all=lims2(x)
% all=[absmax, absmin, mean, standard_deviation, median];
x=x(:);
absmax=max(x); %all of these are matlab built-in functions
absmin=min(x);
m=nanmean(x);
s=nanstd(x);
me=nanmedian(x);

if nargout>0
all=[absmax, absmin, m, s, me];
disp('    absMax   absMin   Mean     Std        Median')
disp(all)
lb %lb prints a line
else
disp('    absMax   absMin   Mean     Std        Median')
disp([absmax, absmin, m, s, me])
lb %lb prints a line
end
    
end

[arango]

I doubt that there is a problem with the accumulation (time integral) in the diagnostic terms. John Wilkin look extensible at the tracer budget last year. We fixed few things with respect to the biological tracers diagnostics in the Fennel model. Time-averaging of momentum and tracers terms are treated in the same way. If there is a problem in the accumulation of momentum terms, we also have it in the tracers terms.

Have you tried the diagnostic terms for (ubar, vbar)? Notice that we have surface and bottom stresses in them.

[rduran]

Hernan,

Thank you for your answer. In the first post of this thread you can see that the 2D diagnostics also look suspicious. You can also see I added the code for the sum of 2D momentum terms (I did include bottom and surface stress).

[rduran]

I verified the results in different matlabs, different computers and using matlab's official way to read netcdf (i.e. ncread in matlab R2012B) instead of using ncload and makenan. I did this for a run made with ROMS revision 571 and in run made with a revision 656 I am getting same results as above.

[apaloczy]

rduran, did you figure out why the momentum balance terms didn't add up? I am having the same problem.

I am computing the following sums for the 2D and 3D momentum balances, respectively (and analogously for the xi-direction balance):

residue_2D = vbar_hadv+vbar_cor+vbar_prsgrd+vbar_sstr+vbar_bstr+vbar_hvisc-vbar_accel
residue_3D = v_hadv+v_vadv+v_cor+v_prsgrd+v_hvisc+v_vvisc-v_accel

I am getting O(1e-12 m/s2) residues in my application, so I looked at the momentum balance in the upwelling test case (https://www.myroms.org/wiki/index.php/UPWELLING_CASE) to see if it closed, but still the same maximum residues of O(1e-12 m/s2) appear in all balances. Attached are plots for the 2D and 3D residues as functions of grid cell, for the last record in the diagnostic file (about 10 days run time). The dots are the residues, and the lines are the smallest non-zero terms in the budget. In the eta direction the smallest term is about the same order of magnitude as the residue, and in the xi direction it is smaller.

I also looked at the terms themselves, and they look generally good, I think. What seems to prevail is geostrophic balance in the across-isobath (eta) direction, and a more complicated balance in the along-isobath (xi) direction. The 2D momentum balance in the along-isobath (xi) direction seems to be frictional (sstr+bstr~0). So here looks like there is nothing physically unreasonable like in your case, but the residues are still too large compared to the smallest terms in the budget.

Originally I thought it could have something to do with the curvilinearity of my grid (even though it is nearly rectangular, since it is spherical and spans 4 degrees of longitude and 3 degrees of latitude, centered at 20.5 S). But then I found the same thing with the 100% rectangular grid of the upwelling test case.

Does anyone else see this in their applications? Did I make a mistake and forgot something?

[rduran]

Hi,

A few grid points at a river mouth were plain wrong. Checking for balanced momentum in an integrated fashion might be a good first step but I needed to look at the details to see where things went wrong. Your case seems to be different. Perhaps the question "what is zero?" is a relevant one. One way to answer that is to check for relative error, say residue / leading term.

If the leading terms are several orders of magnitude bigger than residue it may be that it's not a cause of concern. A safer approach would be to require the smaller terms to be considerably larger than the residue.

Perhaps users with more experience could chip in with the answer to what is zero, it could be your balance is not that bad.

p.s. where does that step-like structure come from? I would look into that, seems odd.

[apaloczy]

Hi rduran, thanks for your suggestion.

I looked at the ratio of the residue and the leading term, and the ratio of the residue and the smallest term.

The ratios between the residue and the leading term show that the residues are everywhere 9-7 orders of magnitude smaller than the leading terms. The ratios between the residue and the smallest term, however, show that the residue can be up to 8 orders of magnitude larger than the smallest term. This is really strange.

A few grid points at a river mouth were plain wrong.

I thought that all grid cells had to add up to zero. So once you identified the cells that were wrong, you removed them from your analysis?

The step-like structure seems to (I'm not sure) be a consequence of the symmetry of the domain, which is rectangular and has a linearly inclined topography (https://www.myroms.org/wiki/index.php/UPWELLING_CASE). Note that these step-like features exist only in the 2D momentum balance. I think it is each term repeating itself along all cells in the xi (along-isobath) direction, for each eta line.

[rduran]

A few grid points at a river mouth were plain wrong.

I thought that all grid cells had to add up to zero. So once you identified the cells that were wrong, you removed them from your analysis?

Yes, momentum terms adding up to zero (whatever zero means) should be pointwise. I guess it is possible to ignore a few bad cells if you can make sure it's not affecting the results you are looking into. In my case I needed to re-run the model anyway so I solved the problem with the river.

The ratios between the residue and the leading term show that the residues are everywhere 9-7 orders of magnitude smaller than the leading terms.

That doesn't seem too bad to me, if you were working with single precision (you are probably using double precision, but just to make a point) then that would be zero. It would still be good to know why its not more like 16 orders of magnitude (zero in double precision) which I do not know.

The ratios between the residue and the smallest term, however, show that the residue can be up to 8 orders of magnitude larger than the smallest term. This is really strange.

Not necessarily, if the smallest term is actually zero (say you are in steady state and u_accel=0) then for any residue ~=0 you get that residue/u_accel = inf

Similarly for a linear flow with u_adv=0 you would get residue / u_adv = inf ...

This shows that you have to consider each term separately. One thing is for sure: a residue similar in magnitude to the leading terms cannot be good.


On a separate note another problem with my plots showing that momentum balance does not add up is that plotting the maximum residue values (in a time series) can be misleading, see the attached plot for an example. It is the distribution of residue values for 2D momentum. Here q1 means first quartile, the blue points are the cumulative distribution, as you can see 50% of the data is around 10^-14 but an outlier (max value) would show a residue of 10^-8, it would be misleading if I just plotted max residues.

[apaloczy]

Thanks for the reply, rduran.

That doesn't seem too bad to me, if you were working with single precision (you are probably using double precision, but just to make a point) then that would be zero. It would still be good to know why its not more like 16 orders of magnitude (zero in double precision) which I do not know.

...if the smallest term is actually zero (say you are in steady state and u_accel=0) then for any residue ~=0 you get that residue/u_accel = inf

Similarly for a linear flow with u_adv=0 you would get residue / u_adv = inf ...

On a separate note another problem with my plots showing that momentum balance does not add up is that plotting the maximum residue values (in a time series) can be misleading

Thanks for bringing these points to my attention. I have been looking at the terms and residues pointwise and for each averaged record (instead of looking at their maximum values in a time series for the whole simulation, for example), and am using double precision. So I guess that for any analyses that rely on the first or second order balances I think this residue situation should be safe.

Thanks again,