Skip to content

Commit

Permalink
Merge pull request #58 from CHLNDDEV/dev-testing
Browse files Browse the repository at this point in the history 
Merging in small improvements from the development branch.
  • Loading branch information
keith roberts authored Mar 4, 2020
2 parents d21ad48 + 360e909 commit ae52714
Show file tree
Hide file tree
Showing 10 changed files with 363 additions and 173 deletions.
152 changes: 70 additions & 82 deletions @edgefx/edgefx.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -466,92 +466,80 @@
if filtit == 1
tic
bs = NaN([obj.nx,obj.ny]);
div = 1e4; grav = 9.807;
nb = ceil([obj.nx,obj.ny]/div); n1s = 1;
for ii = 1:nb(1)
n1e = min(obj.nx,n1s + div - 1); n2s = 1;
for jj = 1:nb(2)
n2e = min(obj.ny,n2s + div - 1);
% Rossby radius of deformation filter
% See Shelton, D. B., et al. (1998): Geographical variability of the first-baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.
ygg = yg(n1s:n1e,n2s:n2e);
f = 2*7.29e-5*abs(sind(ygg));
if barot
% barotropic case
c = sqrt(grav*max(1,-tmpz(n1s:n1e,n2s:n2e)));
else
% baroclinic case (estimate Nm to be 2.5e-3)
Nm = 2.5e-3;
c = Nm*max(1,-tmpz(n1s:n1e,n2s:n2e))/pi;
end
rosb = c./f;
clear f;
% update for equatorial regions
I = abs(ygg) < 5; Re = 6371e3;
twobeta = 4*7.29e-5*cosd(ygg(I))/Re;
rosb(I) = sqrt(c(I)./twobeta);
clear twobeta
% limit rossby radius to 10,000 km for practical purposes
rosb(rosb > 1e7) = 1e7;
% Keep lengthscales rbfilt * barotropic
% radius of deformation
rosb = min(10,max(0,floor(log2(rosb/dy/rbfilt))));
edges = double(unique(rosb(:)));
bst = rosb*0;
for i = 1:length(edges)
if edges(i) > 0
mult = 2^edges(i);
xl = max(1,n1s-mult/2);
xu = min(obj.nx,n1e+mult/2);
if (max(xg(:)) > 179 && min(xg(:)) < -179) || ...
(max(xg(:)) > 359 && min(xg(:)) < 1)
% wraps around
if xu == obj.nx && xl == 1
xr = [obj.nx-mult/2+1:1:obj.nx xl:xu ...
1:mult/2+n1e-obj.nx];
elseif xl == 1
% go to otherside
xr = [obj.nx-mult/2+1:1:obj.nx xl:xu];
elseif xu == obj.nx
% go to otherside
xr = [xl:xu 1:mult/2+n1e-obj.nx];
else
xr = xl:xu;
end
else
xr = xl:xu;
end
yl = max(1,n2s-mult/2);
yu = min(obj.ny,n2e+mult/2);
if max(yg(:)) > 89 && yu == obj.ny
% create mirror around pole
yr = [yl:yu yu-1:-1:2*obj.ny-n2e-mult/2];
else
yr = yl:yu;
end
if mult == 2
tmpz_ft = filt2(tmpz(xr,yr),dy,...
dy*2.01,'lp');
else
tmpz_ft = filt2(tmpz(xr,yr),...
dy,dy*mult,'lp');
end
% delete the padded region
tmpz_ft(1:find(xr == n1s)-1,:) = [];
tmpz_ft(n1e-n1s+2:end,:) = [];
tmpz_ft(:,1:find(yr == n2s)-1) = [];
tmpz_ft(:,n2e-n2s+2:end) = [];
% break into 10 deg latitude chunks, or less if higher res
div = ceil(min(1e7/obj.nx,10*obj.ny/(max(yg(:))-min(yg(:)))));
grav = 9.807;
nb = ceil(obj.ny/div);
n2s = 1;
for jj = 1:nb
n2e = min(obj.ny,n2s + div - 1);
% Rossby radius of deformation filter
% See Shelton, D. B., et al. (1998): Geographical variability of the first-baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.
ygg = yg(:,n2s:n2e);
dxx = mean(dx(n2s:n2e));
f = 2*7.29e-5*abs(sind(ygg));
if barot
% barotropic case
c = sqrt(grav*max(1,-tmpz(:,n2s:n2e)));
else
% baroclinic case (estimate Nm to be 2.5e-3)
Nm = 2.5e-3;
c = Nm*max(1,-tmpz(:,n2s:n2e))/pi;
end
rosb = c./f;
clear f;
% update for equatorial regions
I = abs(ygg) < 5; Re = 6371e3;
twobeta = 4*7.29e-5*cosd(ygg(I))/Re;
rosb(I) = sqrt(c(I)./twobeta);
clear twobeta
% limit rossby radius to 10,000 km for practical purposes
rosb(rosb > 1e7) = 1e7;
% Keep lengthscales rbfilt * barotropic
% radius of deformation
rosb = min(10,max(0,floor(log2(rosb/dy/rbfilt))));
edges = double(unique(rosb(:)));
bst = rosb*0;
for i = 1:length(edges)
if edges(i) > 0
mult = 2^edges(i);
xl = 1; xu = obj.nx;
if (max(xg(:)) > 179 && min(xg(:)) < -179) || ...
(max(xg(:)) > 359 && min(xg(:)) < 1)
% wraps around
xr = [obj.nx-mult/2+1:1:obj.nx xl:xu 1:mult/2];
else
tmpz_ft = tmpz(n1s:n1e,n2s:n2e);
xr = xl:xu;
end
[by,bx] = EarthGradient(tmpz_ft,dy,dx(n2s:n2e)); % get slope in x and y directions
tempbs = sqrt(bx.^2 + by.^2); % get overall slope
bst(rosb == edges(i)) = tempbs(rosb == edges(i));
yl = max(1,n2s-mult/2);
yu = min(obj.ny,n2e+mult/2);
if max(yg(:)) > 89 && yu == obj.ny
% create mirror around pole
yr = [yl:yu yu-1:-1:2*obj.ny-n2e-mult/2];
else
yr = yl:yu;
end
if mult == 2
tmpz_ft = filt2(tmpz(xr,yr),...
[dxx dy],dy*2.01,'lp');
else
tmpz_ft = filt2(tmpz(xr,yr),...
[dxx dy],dy*mult,'lp');
end
% delete the padded region
tmpz_ft(1:find(xr == 1,1,'first')-1,:) = [];
tmpz_ft(obj.nx+1:end,:) = [];
tmpz_ft(:,1:find(yr == n2s,1,'first')-1) = [];
tmpz_ft(:,n2e-n2s+2:end) = [];
else
tmpz_ft = tmpz(:,n2s:n2e);
end
bs(n1s:n1e,n2s:n2e) = bst;
n2s = n2e + 1;
[by,bx] = EarthGradient(tmpz_ft,dy,dx(n2s:n2e)); % get slope in x and y directions
tempbs = sqrt(bx.^2 + by.^2); % get overall slope
bst(rosb == edges(i)) = tempbs(rosb == edges(i));
end
n1s = n1e + 1;
bs(:,n2s:n2e) = bst;
n2s = n2e + 1;
end
toc
% legacy filter
Expand Down
13 changes: 12 additions & 1 deletion @geodata/geodata.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,18 @@
% Class constructor to parse NetCDF DEM data, NaN-delimited vector,
% or shapefile that defines polygonal boundary of meshing
% domain.

% options
%addOptional(p,'bbox',defval);
%addOptional(p,'shp',defval);
%addOptional(p,'h0',defval);
%addOptional(p,'dem',defval);
%addOptional(p,'backupdem',defval);
%addOptional(p,'fp',defval);
%addOptional(p,'weirs',defval);
%addOptional(p,'pslg',defval);
%addOptional(p,'boubox',defval);
%addOptional(p,'window',defval);

% Check for m_map dir
M_MAP_EXISTS=0 ;
if exist('m_proj','file')==2
Expand Down
27 changes: 10 additions & 17 deletions @msh/msh.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1060,6 +1060,7 @@ function write(obj,fname,type)
if mq_l < opt.mqa && (opt.ds || LT ~= LTn)
% Need to clean it again
disp('Poor or overlapping elements, cleaning again')
disp(['(Min Qual = ' num2str(mq_l) ')'])
% repeat without projecting (already projected)
ii = find(strcmp(varargino,'proj'));
if ~isempty(ii)
Expand Down Expand Up @@ -1855,7 +1856,6 @@ function write(obj,fname,type)
error('Mesh 2 is invalid. Please execute msh.clean on object');
end


disp('Forming outer boundary for inset...')
try
[cell1] = extdom_polygon(extdom_edges2(t1,p1),p1,-1,0);
Expand Down Expand Up @@ -1924,28 +1924,20 @@ function write(obj,fname,type)
end

merge = msh() ; merge.p = pm; merge.t = tm ;
clear pm tm pmid p1 t1 p2 t2

% This step does some mesh smoothing around the intersection
% of the meshes.
l = cellfun(@length,cell1);
[~,s] = sort(l,'descend'); cell1 = cell1(s);
% just do for largest 10 polygons
in = false(size(merge.p,1),1);
for ii = 1:min(length(cell1),10)
k = boundary(cell1{ii}(1:end-1,1),cell1{ii}(1:end-1,2));
xout = cell1{ii}(k,1); yout = cell1{ii}(k,2);
[xe,ye] = enlargePoly(xout,yout,1.35);
int = inpoly(merge.p,[xe,ye]);
in(int) = true;
end

if ~isempty(cleanargin)
locked = [merge.p(~in,:); pfixx];
% lock anything vertices more than 2*dmax away
% from intersection
[~,dst] = ourKNNsearch(p1',p1',2);
dmax = max(dst(:,2));
[~,dst1] = ourKNNsearch(p1',merge.p',1);
[~,dst2] = ourKNNsearch(p2',merge.p',1);
locked = [merge.p(dst1 > 2*dmax | dst2 > 2*dmax,:); pfixx];
cleanargin{end} = locked;
% iteration is done in clean if smoothing creates neg quality
merge = clean(merge,cleanargin);
end
clear pm tm pmid p1 t1 p2 t2

merge.pfix = pfixx ;
% edges don't move but they are no longer valid after merger
Expand Down Expand Up @@ -2228,6 +2220,7 @@ function write(obj,fname,type)
end
bars = [obj.t(:,[1,2]); obj.t(:,[1,3]); obj.t(:,[2,3])]; % Interior bars duplicated
bars = unique(sort(bars,2),'rows'); % Bars as node pairs
if nargout == 1; return; end
if type == 0
% Compute based on Haversine spherical earth distances
long = zeros(length(bars)*2,1);
Expand Down
Loading

0 comments on commit ae52714

Please sign in to comment.