InternalWaveModelArbitraryStratification requires 'monotonic density'

[qyxiao]

Hi Jeffrey, I'm trying to use your matlab code 'InternalWaveModelArbitraryStratification' to generate IGW that may exist based on a specific GCM setup. What I did was calculating the the rho at each grid depth and interpolate the data as a function handle input for your function. However, though I think the function is monotonic, I kept getting error "Density must be monotonic to use this class. The function you provided is not monotonic." Can you give any suggestions about what I should do?

To be more specific, below is what I've done:

%%%% this part following your example on https://github.com/JeffreyEarly/GLOceanKit/tree/master/Matlab/InternalWaveModel except changing Lz

aspectRatio = 8;

Lx = 100e3; Ly = aspectRatio*100e3; Lz = 1000;

Nx = 64; Ny = aspectRatio*64; Nz = 65; % 2^n + 1 grid points, to match the Winters model, but 2^n ok too.

latitude = 31;

%%%% these are densities and corresponding depth ZZ = [-5.000000e-01, -1.570000e+00, -2.790000e+00, -4.185000e+00, -5.780000e+00, -7.595000e+00, -9.660000e+00, -1.201000e+01, -1.468000e+01, -1.770500e+01, -2.112500e+01, -2.499000e+01, -2.934500e+01, -3.424000e+01, -3.972500e+01, -4.585500e+01, -5.269000e+01, -6.028000e+01, -6.868500e+01, -7.796500e+01, -8.817500e+01, -9.937000e+01, -1.116000e+02, -1.249150e+02, -1.393650e+02, -1.549900e+02, -1.718250e+02, -1.899000e+02, -2.092350e+02, -2.298550e+02, -2.517700e+02, -2.749850e+02, -2.995050e+02, -3.253200e+02, -3.524200e+02, -3.807900e+02, -4.104100e+02, -4.412550e+02, -4.733050e+02, -5.065400e+02, -5.409350e+02, -5.764650e+02, -6.131100e+02, -6.508550e+02, -6.896850e+02, -7.295950e+02, -7.705850e+02, -8.126600e+02, -8.558350e+02, -9.001350e+02, -9.455950e+02, -9.922600e+02, -1.040180e+03, -1.089425e+03, -1.140080e+03, -1.192235e+03, -1.246005e+03, -1.301520e+03, -1.358920e+03, -1.418375e+03, -1.480075e+03, -1.544225e+03, -1.611060e+03, -1.680845e+03, -1.753875e+03, -1.830475e+03, -1.911015e+03, -1.995905e+03, -2.085595e+03, -2.180595e+03, -2.281470e+03, -2.388845e+03, -2.503415e+03, -2.625955e+03, -2.757325e+03, -2.912665e+03];

rho_full = [999.1796 , 999.1814 , 999.1829 , 999.18396, 999.18445, 999.1848 , 999.1855 , 999.1861 , 999.1865 , 999.18713, 999.1879 , 999.1893 , 999.1916 , 999.19617, 999.2013 , 999.2068 , 999.21313, 999.21954, 999.22705, 999.23474, 999.2434 , 999.2527 , 999.262 , 999.2729 , 999.28394, 999.2963 , 999.3087 , 999.32275, 999.3365 , 999.35144, 999.36597, 999.38196, 999.39746, 999.4143 , 999.431 , 999.44904, 999.4674 , 999.48395, 999.50146, 999.5169 , 999.5331 , 999.5506 , 999.56793, 999.5868 , 999.6062 , 999.62585, 999.64496, 999.6639 , 999.6809 , 999.69916, 999.71844, 999.7413 , 999.7682 , 999.7986 , 999.8273 , 999.8481 , 999.85864, 999.8604 , 999.8618 , 999.8639 , 999.87024, 999.87305, 999.8742 , 999.87494, 999.87537, 999.8768 , 999.8779 , 999.8806 , 999.88324, 999.88354, 999.88873, 999.8988 , 999.8996 , 999.9002 , 999.90063, 999.9044 ];

%%%%% the function calling is slightly different from your example; I guess you forgot to update the example after you change the parameter list of the function maxModes = 32; wavemodel = InternalWaveModelArbitraryStratification([Lx, Ly, Lz], [Nx, Ny, Nz], @(z)interp1(ZZ, rho_full, z), ZZ, latitude, 'nModes', maxModes);

At this point, the 'monotonic density' error appears, which traced back to InternalModesSpectral.m. I tried to play with part of the function InitializeWithFunction(self, rho, zMin, zMax), and it looks like the appearance of error somehow realted to values of zMin, zMax. Still I'm not sure about the detail about this function and parameters.

It'd be great if you can provide some suggestions about how to run the code in this situation. Thank you very much!

[JeffreyEarly]

I couldn't repeat the same error, although I did find some other issues.

There's no need to create a function handle with the gridded density data... just pass it the data directly. That said, the modes that come from that data aren't the smoothest. (see the last two lines of my code below). You'll have to see how well it works... otherwise maybe smooth it a bit more.

I made one small change to the initialization routine of the WaveModel, so pull that before trying the following.

aspectRatio = 8;

%%%% these are densities and corresponding depth ZZ = [-5.000000e-01, -1.570000e+00, -2.790000e+00, -4.185000e+00,-5.780000e+00, -7.595000e+00, -9.660000e+00, -1.201000e+01,-1.468000e+01, -1.770500e+01, -2.112500e+01, -2.499000e+01,-2.934500e+01, -3.424000e+01, -3.972500e+01, -4.585500e+01,-5.269000e+01, -6.028000e+01, -6.868500e+01, -7.796500e+01,-8.817500e+01, -9.937000e+01, -1.116000e+02, -1.249150e+02,-1.393650e+02, -1.549900e+02, -1.718250e+02, -1.899000e+02,-2.092350e+02, -2.298550e+02, -2.517700e+02, -2.749850e+02,-2.995050e+02, -3.253200e+02, -3.524200e+02, -3.807900e+02,-4.104100e+02, -4.412550e+02, -4.733050e+02, -5.065400e+02,-5.409350e+02, -5.764650e+02, -6.131100e+02, -6.508550e+02,-6.896850e+02, -7.295950e+02, -7.705850e+02, -8.126600e+02,-8.558350e+02, -9.001350e+02, -9.455950e+02, -9.922600e+02,-1.040180e+03, -1.089425e+03, -1.140080e+03, -1.192235e+03,-1.246005e+03, -1.301520e+03, -1.358920e+03, -1.418375e+03,-1.480075e+03, -1.544225e+03, -1.611060e+03, -1.680845e+03,-1.753875e+03, -1.830475e+03, -1.911015e+03, -1.995905e+03,-2.085595e+03, -2.180595e+03, -2.281470e+03, -2.388845e+03,-2.503415e+03, -2.625955e+03, -2.757325e+03, -2.912665e+03]; rho_full = [999.1796 , 999.1814 , 999.1829 , 999.18396, 999.18445, 999.1848,999.1855 , 999.1861 , 999.1865 , 999.18713, 999.1879 , 999.1893 ,999.1916 , 999.19617, 999.2013 , 999.2068 , 999.21313, 999.21954,999.22705, 999.23474, 999.2434 , 999.2527 , 999.262 , 999.2729 ,999.28394, 999.2963 , 999.3087 , 999.32275, 999.3365 , 999.35144,999.36597, 999.38196, 999.39746, 999.4143 , 999.431 , 999.44904,999.4674 , 999.48395, 999.50146, 999.5169 , 999.5331 , 999.5506 ,999.56793, 999.5868 , 999.6062 , 999.62585, 999.64496, 999.6639 ,999.6809 , 999.69916, 999.71844, 999.7413 , 999.7682 , 999.7986 ,999.8273 , 999.8481 , 999.85864, 999.8604 , 999.8618 , 999.8639 ,999.87024, 999.87305, 999.8742 , 999.87494, 999.87537, 999.8768 ,999.8779 , 999.8806 , 999.88324, 999.88354, 999.88873, 999.8988 ,999.8996 , 999.9002 , 999.90063, 999.9044 ];

Lx = 100e3; Ly = aspectRatio*100e3; Lz = max(ZZ)-min(ZZ);

Nx = 64; Ny = aspectRatio*64; Nz = length(ZZ); % 2^n + 1 grid points, to match the Winters model, but 2^n ok too.

latitude = 31;

maxModes = 32; wavemodel = InternalWaveModelArbitraryStratification([Lx, Ly, Lz], [Nx, Ny, Nz], rho_full, ZZ, latitude, 'nModes', maxModes);

im = wavemodel.internalModes; im.ShowLowestModesAtWavenumber(0)

[qyxiao]

Hi Jeffrey, I downloaded the new code and did as you suggested but got new error at the step of calling InternalWaveModelArbitraryStratification:

Undefined variable "InterpolatingSpline" or class "InterpolatingSpline.KnotPointsForPoints".

Error in InternalModesSpectral/InitializeWithGrid (line 418) z_knot = InterpolatingSpline.KnotPointsForPoints(zIn,K,1);

Error in InternalModesWKBSpectral/InitializeWithGrid (line 123) InitializeWithGrid@InternalModesSpectral(self,rho,zIn);

Error in InternalModesBase (line 186) self.InitializeWithGrid(rho,z_in);

Error in InternalModesSpectral (line 118) self@InternalModesBase(rho,zIn,zOut,latitude,varargin{:});

Error in InternalModesWKBSpectral (line 39) self@InternalModesSpectral(rho,z_in,z_out,latitude, varargin{:});

Error in InternalModesAdaptiveSpectral (line 49) self@InternalModesWKBSpectral(rho,z_in,z_out,latitude, varargin{:});

Error in InternalModes (line 175) self.internalModes = InternalModesAdaptiveSpectral(rho,zIn,zOut,latitude,extraargs{:});

Error in InternalWaveModelArbitraryStratification (line 87) im = InternalModes(rho,z,z,latitude, varargin{:});

Error in QuickStart (line 115) wavemodel = InternalWaveModelArbitraryStratification([Lx, Ly, Lz], [Nx, Ny, Nz], rho_full, ZZ, latitude, 'nModes', maxModes);

[JeffreyEarly]

Sounds like you need to add a path to GLNumericalModelingKit/Matlab...

[qyxiao]

You are right. The code is working now. Thank you very much!

[qyxiao]

Hi Jeffrey, While using the code above, I keep receiving warning like "You will miss 199.58% of the energy due to limited vertical modes." no matter how large the parameter 'nModes' is. Do you have any clue why this is happening? Thanks a lot.

[JeffreyEarly]

I'm not sure what's going on. I ran the exact code I had posted above, then did wavemodel.InitializeWithGMSpectrum(1.0);

and that gave me,

You will miss 3.11% of the energy due to limited vertical modes.
Solving the EVP for 495 unique wavenumbers.

After it computed that modes,

After distributing energy across frequency and mode, you still have 88.60% of reference GM energy.
Due to restricted domain size, the j=1,k=l=0 mode contains 15.21% the total energy.
The (gridded, external) wave field sums to (88.60%, 0.00%) GM given the scales, and the randomized field sums to (88.54%, 0.00%) GM

Are you doing anything different?

[qyxiao]

Hi Jeffrey, yes it looks like this is caused by my changing of latitude from 31 to -35...

[cqdbc]

Hi Jeffrey, I copied the codes above but it didn't work. I just changed ZZ and rho_full into 76*1 vectors. The error is as follows

Wrong use of round Too many input parameters.

Error InternalModesExponentialStratification.IsStratificationExponential (line 397) N0 = round(sqrt( exp(intercept) ),3,'significant');

Error InternalModes (line 160) [isStratificationExponential, rho_params] = InternalModesExponentialStratification.IsStratificationExponential(rho,zIn);

Error InternalWaveModelArbitraryStratification (line 126) im = InternalModes(rho,zIn,z,latitude, extraargs{:});

Error test (line 22) wavemodel = InternalWaveModelArbitraryStratification([Lx, Ly, Lz], [Nx, Ny,Nz], rho, ZZ, latitude, 'nModes', maxModes);

Can you give any suggestions about what I should do?