From 9305c29d76793935792cf02cda5e7ebcfc234bf1 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Fri, 28 Feb 2020 15:07:00 -0600 Subject: [PATCH 1/6] Improving barotropic rossby radius filter, and locked points for mesh merging - Changing filt2 to have different dx and dy inputs - In rossby radius filter split up into only y-direction chunks to provide a mean dx for latitude bands into filt2 - in the mesh merger change method which lock points that are far away from the intersection zone. --- @edgefx/edgefx.m | 152 +++++++++++++++++++++------------------------ @geodata/geodata.m | 13 +++- @msh/msh.m | 25 +++----- utilities/filt2.m | 19 ++++-- 4 files changed, 105 insertions(+), 104 deletions(-) diff --git a/@edgefx/edgefx.m b/@edgefx/edgefx.m index 29a6b525..6a4f416c 100644 --- a/@edgefx/edgefx.m +++ b/@edgefx/edgefx.m @@ -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 diff --git a/@geodata/geodata.m b/@geodata/geodata.m index 26f41580..b65e5421 100644 --- a/@geodata/geodata.m +++ b/@geodata/geodata.m @@ -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 diff --git a/@msh/msh.m b/@msh/msh.m index 9189d242..b863de0c 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -1923,28 +1923,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 @@ -2227,6 +2219,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); diff --git a/utilities/filt2.m b/utilities/filt2.m index aa65e5a8..dbbc5ab4 100644 --- a/utilities/filt2.m +++ b/utilities/filt2.m @@ -82,9 +82,9 @@ narginchk(4,4) % assert(license('test','image_toolbox')==1,'Error: I''m sorry, the filt2 function requires the Image Processing Toolbox.') assert(ismatrix(Z)==1,'Input error: Z must be a 2d matrix.') -assert(isscalar(res)==1,'Input error: res must be a scalar value.') +%assert(isscalar(res)==1,'Input error: res must be a scalar value.') assert(ismember(lower(filtertype),{'lp','hp','bp','bs'}),'Input error: filtertype must be ''hp'', ''lp'', ''bp'', or ''bs''.') -if lambda<=(2*res) +if lambda<=(2*max(res)) warning('Nyquist says the wavelength should exceed two times the resolution of the dataset, which is an unmet condition based on these inputs. I''ll give you some numbers, but I would''t trust ''em if I were you.') end @@ -97,11 +97,11 @@ %% Design filter: % 2*pi*sigma is the wavelength at which the amplitude is multiplied by a factor of about 0.6 (more exactly, exp(-0.5)) -sigma = (lambda(1)/res) /(2*pi); +sigma = (lambda(1)./res) /(2*pi); %f = fspecial('gaussian',2*ceil(2.6*sigma)+1,sigma); % WJP implementation without image processing toolbox, subfunc at bottom of % this function -f = gaussian2D([2*ceil(2.6*sigma)+1 2*ceil(2.6*sigma)+1], sigma); +f = gaussian2D([2*ceil(2.6*sigma(1))+1 2*ceil(2.6*sigma(end))+1], sigma); %% Now filter the data @@ -619,8 +619,17 @@ siz = (siz-1)./2; [x,y] = meshgrid(-siz(2):siz(2),-siz(1):siz(1)); + + std = max(std); + sr = siz(2)/siz(1); + if sr > 1 + sy = sr^2; sx = 1; + else + sx = sr^2; sy = 1; + end + % analytic function - h = exp(-(x.*x + y.*y)/(2*std*std)); + h = exp(-(sx*x.*x + sy*y.*y)/(2*std*std)); % truncate very small values to zero h(h Date: Fri, 28 Feb 2020 15:16:32 -0600 Subject: [PATCH 2/6] updates from dev - GridData updates from dev for a local nan fill option and changing the slope calculation to and RMS instead of mean(abs) - introducing the Calc_tau0 for constant tau0 > 0 option (computes allowable tau0 based on your timestep) --- @msh/private/GridData.m | 149 +++++++++++++++++++++++++--------------- utilities/Calc_tau0.m | 33 +++++++-- 2 files changed, 122 insertions(+), 60 deletions(-) diff --git a/@msh/private/GridData.m b/@msh/private/GridData.m index e04b8485..0a673f05 100644 --- a/@msh/private/GridData.m +++ b/@msh/private/GridData.m @@ -1,5 +1,5 @@ function obj = GridData(geodata,obj,varargin) -% obj = GridData(geodata,obj,'K',K,'type',type); +% obj = GridData(geodata,obj,varargin); % GridData: Uses the cell-averaged approach to interpolate the nodes % on the unstructured grid to the DEM. % Input : geodata - either a geodata class or a filename of the netcdf DEM @@ -65,7 +65,8 @@ K = (1:nn)'; % K is all of the grid type = 'all'; interp = 'CA'; -NaNs = 'ignore'; +NaNs = 'ignore'; nanfill = false; +ignoreOL = 0 ; N = 1 ; mindepth = -inf ; maxdepth = +inf ; @@ -87,30 +88,27 @@ elseif ii == 5 N = varargin{ind*2} ; elseif ii ==6 - mindepth = varargin{ind*2} ; + mindepth = varargin{ind*2} ; elseif ii ==7 - maxdepth = varargin{ind*2} ; + maxdepth = varargin{ind*2} ; elseif ii ==8 - ignoreOL = varargin{ind*2} ; + ignoreOL = varargin{ind*2} ; end end end end -if ~exist('ignoreOL','var') - ignoreOL = 0 ; -end if ignoreOL - disp('NaNing overland data before interpolating') + disp('NaNing overland data before interpolating') end if N > 1 - disp(['Enlarging CA stencil by factor ',num2str(N)]) ; + disp(['Enlarging CA stencil by factor ',num2str(N)]) ; end if mindepth > -inf - disp(['Bounding minimum depth to ',num2str(mindepth), ' meters.']) ; + disp(['Bounding minimum depth to ',num2str(mindepth), ' meters.']) ; end if maxdepth < inf - disp(['Bounding maximum depth to ',num2str(maxdepth), ' meters.']) ; + disp(['Bounding maximum depth to ',num2str(maxdepth), ' meters.']) ; end if strcmp(type,'slope') @@ -118,44 +116,48 @@ 'calculating the slopes']) end +if strcmp(NaNs,'fill') || strcmp(NaNs,'fillinside') + disp('Fill in NaNs using nearest neighbour interpolation.') + nanfill = true; +end if strcmp(NaNs,'fill') - disp('Fill in NaNs using nearest neighbour interpolation.') warning(['Note that will try and put bathy everywhere on mesh even ' ... - 'outside of gdat/dem extents unless K logical is set.']) + 'outside of gdat/dem extents unless K logical is set. ' ... + 'Set to nan to "fillinside" to avoid this.']) end %% Let's read the LON LAT of DEM if not already geodata flipUD = 0; if ~isa(geodata,'geodata') - try - DEM_XA = double(ncread(geodata,'lon')); - DEM_YA = double(ncread(geodata,'lat')); - catch - DEM_XA = double(ncread(geodata,'x')); - DEM_YA = double(ncread(geodata,'y')); - end + [xvn, yvn, zvn] = getdemvarnames(geodata); + DEM_XA = double(ncread(geodata,xvn)); + DEM_YA = double(ncread(geodata,yvn)); DELTA_X = mean(diff(DEM_XA)); DELTA_Y = mean(diff(DEM_YA)); if DELTA_Y < 0 flipUD = 1; - DEM_YA = flipud(DEM_YA); - DELTA_Y = mean(diff(DEM_YA)); + DELTA_Y = -DELTA_Y; end else DEM_XA = geodata.Fb.GridVectors{1}; DEM_YA = geodata.Fb.GridVectors{2}; DEM_Z = -geodata.Fb.Values; + if ignoreOL + DEM_Z(DEM_Z <= 0) = NaN ; + end DELTA_X = mean(diff(DEM_XA)); DELTA_Y = mean(diff(DEM_YA)); [DEM_X,DEM_Y] = ndgrid(DEM_XA,DEM_YA); end if max(DEM_XA) > 180 - lon_change = obj.p(:,1) < 0; + lon_change = obj.p(:,1) < 0; lon_dir = 1; +elseif max(obj.p(:,1)) > 180 + lon_change = obj.p(:,1) > 180; lon_dir = -1; else - lon_change = false(length(obj.p),1); + lon_change = false(length(obj.p),1); lon_dir = 0; end -obj.p(lon_change,1) = obj.p(lon_change,1) + 360; +obj.p(lon_change,1) = obj.p(lon_change,1) + lon_dir*360; % kjr edit 20180320 if length(K) == length(obj.p) @@ -170,6 +172,9 @@ obj.p(:,2) < min(DEM_YA)-DELTA_Y | ... obj.p(:,2) > max(DEM_YA)+DELTA_Y; K(outside) = []; + if isempty(K) + warning('no mesh vertices contained within DEM bounds'); return; + end end end @@ -249,23 +254,22 @@ if exist('DEM_X','var') clear DEM_X DEM_Y DEM_Z end - [DEM_X,DEM_Y] = ndgrid(DEM_XA(I),DEM_YA(J)); - - finfo = ncinfo(geodata); - for ii = 1:length(finfo.Variables) - if length(finfo.Variables(ii).Size) == 2 - zvarname = finfo.Variables(ii).Name; - break - end - end - DEM_Z = single(ncread(geodata,zvarname,... + [DEM_X,DEM_Y] = ndgrid(DEM_XA(I),DEM_YA(J)); + DEM_Z = single(ncread(geodata,zvn,... [I(1) J(1)],[length(I) length(J)])); if flipUD + % handle DEMS from packed starting from the bottom left + DEM_YA = flipud(DEM_YA) ; + DEM_Y = fliplr(DEM_Y); DEM_Z = fliplr(DEM_Z); end % make into depths (ADCIRC compliant) DEM_Z = -DEM_Z; + + if ignoreOL + DEM_Z(DEM_Z <= 0) = NaN ; + end end % bound all depths below mindepth @@ -331,20 +335,17 @@ % Average for the depths if strcmp(type,'depth') || strcmp(type,'all') for ii = 1:length(K) - %%% MTC special option goes here - if(ignoreOL) - pts = reshape(DEM_Z(IDXL(ii):IDXR(ii),... - IDXB(ii):IDXT(ii)),[],1); - pts(pts < 0) = NaN ; - else - pts = reshape(DEM_Z(IDXL(ii):IDXR(ii),... - IDXB(ii):IDXT(ii)),[],1); - end - b(ii) = mean(pts,'omitnan'); + pts = reshape(DEM_Z(IDXL(ii):IDXR(ii),... + IDXB(ii):IDXT(ii)),[],1); + b(ii) = mean(pts,'omitnan'); + end + if sum(~isnan(b)) == 0 + warning('All depths were NaNs. Doing nothing and returning'); + return end % Try and fill in the NaNs - if strcmp(NaNs,'fill') - if ~isempty(find(isnan(b),1)) + if nanfill + if sum(isnan(b)) > 0 localcoord = obj.p(K,:); KI = knnsearch(localcoord(~isnan(b),:),localcoord(isnan(b),:)); bb = b(~isnan(b),:); @@ -353,15 +354,17 @@ end end - % Average for the slopes + % RMS for the slopes if strcmp(type,'slope') || strcmp(type,'all') for ii = 1:length(K) - bx(ii) = mean(reshape(DEM_ZX(IDXL(ii):IDXR(ii),... - IDXB(ii):IDXT(ii)),[],1),'omitnan'); - by(ii) = mean(reshape(DEM_ZY(IDXL(ii):IDXR(ii),... - IDXB(ii):IDXT(ii)),[],1),'omitnan'); + pts = reshape(DEM_ZX(IDXL(ii):IDXR(ii),... + IDXB(ii):IDXT(ii)),[],1); + bx(ii) = sqrt(mean(pts.^2,'omitnan')); + pts = reshape(DEM_ZY(IDXL(ii):IDXR(ii),... + IDXB(ii):IDXT(ii)),[],1); + by(ii) = sqrt(mean(pts.^2,'omitnan')); end - if strcmp(NaNs,'fill') + if nanfill % Try and fill in the NaNs if ~isempty(find(isnan(bx),1)) localcoord = obj.p(K,:); @@ -419,6 +422,42 @@ obj.bx(K_o) = sign(Hx(K_o)).*obj.bx(K_o); obj.by(K_o) = sign(Hy(K_o)).*obj.by(K_o); end -obj.p(lon_change,1) = obj.p(lon_change,1) - 360; +obj.p(lon_change,1) = obj.p(lon_change,1) - lon_dir*360; %EOF end + +function [xvn, yvn, zvn] = getdemvarnames(fname) + % Define well-known variables for longitude and latitude + % coordinates in Digital Elevation Model NetCDF file (CF + % compliant). + xvn = []; yvn = []; zvn = []; + wkv_x = {'x','Longitude','longitude','lon','lon_z'} ; + wkv_y = {'y','Latitude', 'latitude','lat','lat_z'} ; + finfo = ncinfo(fname); + for ii = 1:length(finfo.Variables) + if ~isempty(xvn) && ~isempty(yvn) && ~isempty(zvn); break; end + if length(finfo.Variables(ii).Size) == 1 + if isempty(xvn) && ... + any(strcmp(finfo.Variables(ii).Name,wkv_x)) + xvn = finfo.Variables(ii).Name; + end + if isempty(yvn) && ... + any(strcmp(finfo.Variables(ii).Name,wkv_y)) + yvn = finfo.Variables(ii).Name; + end + elseif length(finfo.Variables(ii).Size) == 2 + if isempty(zvn) + zvn = finfo.Variables(ii).Name; + end + end + end + if isempty(xvn) + error('Could not locate x coordinate in DEM') ; + end + if isempty(yvn) + error('Could not locate y coordinate in DEM') ; + end + if isempty(zvn) + error('Could not locate z coordinate in DEM') ; + end +end diff --git a/utilities/Calc_tau0.m b/utilities/Calc_tau0.m index f608e32a..725382d4 100644 --- a/utilities/Calc_tau0.m +++ b/utilities/Calc_tau0.m @@ -20,16 +20,16 @@ attrname = 'primitive_weighting_in_continuity_equation'; -if isempty(obj.b) - error('No bathymetry data in grid to calculate the tau0 coefficients') -end - %% Test optional arguments % default MinDepth = 10; % m Distance = 2; % km +opt = -3; % calculate tau0 using -3 option +mm = 2/3; % consistent mass matrix +kappa = 0.4; % GWCE weight +sf = 0.6; % suggested safety factor if ~isempty(varargin) - names = {'depth','distance'}; + names = {'depth','distance','opt','kappa','sf'}; for ii = 1:length(names) ind = find(~cellfun(@isempty,strfind(varargin(1:2:end),names{ii}))); if ~isempty(ind) @@ -37,11 +37,34 @@ MinDepth = varargin{ind*2}; elseif ii == 2 Distance = varargin{ind*2}; + elseif ii == 3 + opt = varargin{ind*2}; + elseif ii == 4 + kappa = varargin{ind*2}; + elseif ii == 5 + sf = varargin{ind*2}; end end end end +if opt > 0 + dt = opt; + obj.f15.tau0 = sf*4*(2-mm)*(3*kappa-1)/dt; + obj.f15.a00b00c00 = [kappa kappa 1-2*kappa]; + obj.f15.im = 511111; + disp(['Computed constant value of tau0 = ' num2str(obj.f15.tau0) ' based on...']) + disp(['mm = ' num2str(mm)]) + disp(['kappa = ' num2str(kappa)]) + disp(['sf = ' num2str(sf)]) + disp(['dt = ' num2str(dt)]) + return; +end + +if isempty(obj.b) || all(obj.b) == 0 + error('No bathymetry data in grid to calculate the tau0 coefficients') +end + %% Get all the unique bar edges and distances all_bars = [obj.t(:,[1,2]); obj.t(:,[1,3]); obj.t(:,[2,3])]; bars = unique(sort(all_bars,2),'rows'); From 798b304e48cf1e928e330bd04d824f5d1e681ddf Mon Sep 17 00:00:00 2001 From: William Pringle Date: Fri, 28 Feb 2020 15:20:14 -0600 Subject: [PATCH 3/6] more dev changes - Changing Calc_f13_inpoly to Calc_f13 and introducing the varargin type inputs. also adding more attributes - output for Calc_IT_Fric defval header of C_it is changed to be the same as the input C_it to the function. --- utilities/Calc_IT_Fric.m | 4 +- utilities/Calc_f13.m | 128 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 130 insertions(+), 2 deletions(-) create mode 100644 utilities/Calc_f13.m diff --git a/utilities/Calc_IT_Fric.m b/utilities/Calc_IT_Fric.m index 9f211594..4ac903ff 100644 --- a/utilities/Calc_IT_Fric.m +++ b/utilities/Calc_IT_Fric.m @@ -152,7 +152,7 @@ % Default Values obj.f13.defval.Atr(NA).AttrName = attrname; % We can just put in the options here -obj.f13.defval.Atr(NA).Unit = strcat(type,', C_it = ',num2str(C_it),... +obj.f13.defval.Atr(NA).Unit = strcat(type,', C_it = ',num2str(C_it*4*pi),... ', ',num2str(crit),', D',num2str(MinDepth)); if isempty(N_filename) % We want to just give it the gradients of bathy (and J for the tensor) @@ -225,4 +225,4 @@ end %EOF -end \ No newline at end of file +end diff --git a/utilities/Calc_f13.m b/utilities/Calc_f13.m new file mode 100644 index 00000000..5bc01a54 --- /dev/null +++ b/utilities/Calc_f13.m @@ -0,0 +1,128 @@ +function obj = Calc_f13(obj,attribute,varargin) +% obj = Calc_f13(obj,attribute,varargin) +% Input a msh class object and get back the specified +% attribute in a f13 structure of the obj +% +% Inputs: 1) A msh class obj with bathy on it +% 2) The attribute indicator +% Attributes currently supported: +% 'Cf' ('quadratic_friction_coefficient_at_sea_floor') +% 'Ev' ('average_horizontal_eddy_viscosity_in_sea_water_wrt_depth') +% 'Mn' ('mannings_n_at_sea_floor') +% 'Ss' ('surface_submergence_state') +% 'Re' ('initial_river_elevation') +% +% 3) then either: +% 'inpoly' followed by... +% - A cell-arry of polygons in which you would like to alter +% the attribute. +% - A set of attribute values that correspond 1-to-1 with the +% cell of polygons. +% - (optional) A set of 0 or 1's that correspond 1-to-1 with the +% cell of polygons as to whether the in or out polygon is +% selected +% +% or 'assign' followed by... +% - an array of values to assign with length the same as the +% number of vertices in the msh obj +% +% Outputs: 1) msh class obj with attribute values populating the f13 struct +% +% Author: Keith Roberts, WP to make it for general attribute +% Created: April 5 2018, July 5 2018, June 6 2019 (cleaning up) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +if strcmpi(attribute,'Cf') + attrname = 'quadratic_friction_coefficient_at_sea_floor'; + default_val = 0.0025; % Default Cf +elseif strcmpi(attribute,'Ev') + attrname = 'average_horizontal_eddy_viscosity_in_sea_water_wrt_depth'; + default_val = 0.05; % Default smagorinsky eddy viscosity +elseif strcmpi(attribute,'Mn') + attrname = 'mannings_n_at_sea_floor'; + default_val = 0.020; +elseif strcmpi(attribute,'Ss') + attrname = 'surface_submergence_state'; + default_val = 0; +elseif strcmpi(attribute,'Re') + attrname = 'initial_river_elevation'; + default_val = 0; +else + error(['Attribute ' attribute ' not currently supported']) +end + +if strcmpi(varargin{1},'inpoly') + polys = varargin{2}; + cfvals = varargin{3}; + if length(varargin) < 4 + inverse = 0*cfvals; + else + inverse = varargin{4}; + end +elseif strcmpi(varargin{1},'assign') + Cf_val_on_mesh = varargin{2}; +else + error(['First varargin entry (' varargin{1} ') unsupported']) +end + +%% Make into f13 struct +if isempty(obj.f13) + % Add this as first entry in f13 struct + obj.f13.AGRID = obj.title; + obj.f13.NumOfNodes = length(obj.p); + obj.f13.nAttr = 1; + NA = 1; +else + broken = 0; + for NA = 1:obj.f13.nAttr + if strcmp(attrname,obj.f13.defval.Atr(NA).AttrName) + broken = 1; + % overwrite existing tau0 + break + end + end + if ~broken + % add internal_tide to list + obj.f13.nAttr = obj.f13.nAttr + 1; + NA = obj.f13.nAttr; + end +end + +% Default Values +obj.f13.defval.Atr(NA).AttrName = attrname; +% We can just put in the options here +obj.f13.defval.Atr(NA).Unit = 'unitless'; +valpernode = 1; +obj.f13.defval.Atr(NA).ValuesPerNode = valpernode ; +obj.f13.defval.Atr(NA).Val = default_val ; + +% User Values +obj.f13.userval.Atr(NA).AttrName = attrname ; +cf = obj.p(:,1)*0 + default_val; +if strcmpi(varargin{1},'inpoly') + for i = 1 : length(polys) + in = inpoly([obj.p(:,1),obj.p(:,2)],polys{i}); + if inverse(i) + cf(~in) = cfvals(i); + else + cf(in) = cfvals(i); + end + end +else + cf = Cf_val_on_mesh; +end +numnodes = length(find(cf ~= default_val)); +obj.f13.userval.Atr(NA).usernumnodes = numnodes ; +% Print out list of nodes for each +K = find(cf ~= default_val); +obj.f13.userval.Atr(NA).Val = [K cf(K)]'; + +if ~isempty(obj.f15) + % Change attribute in obj.f15 + disp(['Adding on ' attrname ' into fort.15 struct']) + obj.f15.nwp = obj.f15.nwp + 1; + obj.f15.AttrName(obj.f15.nwp).name = attrname; +end + +%EOF +end From 6bbb9c48f9a5762bef00d2f8ce9187e3cdbf959f Mon Sep 17 00:00:00 2001 From: William Pringle Date: Fri, 28 Feb 2020 15:26:20 -0600 Subject: [PATCH 4/6] extract_boundary and subdomain small edits for dev - fixmesh functions handles the indexing for carrying over bathy and slope in ExtractSubDomain - Fix on the type2 input to extract_boundary to be 22 instead of 2 for river. --- utilities/ExtractSubDomain.m | 3 +-- utilities/extract_boundary.m | 8 ++++---- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/utilities/ExtractSubDomain.m b/utilities/ExtractSubDomain.m index 0fa7b6b2..f4ef2c64 100644 --- a/utilities/ExtractSubDomain.m +++ b/utilities/ExtractSubDomain.m @@ -25,8 +25,7 @@ else t(in,:) = []; end -[p1,t] = fixmesh(p,t); -[~,~,ind] = intersect(p1,p,'rows','stable'); +[p1,t,ind] = fixmesh(p,t); obj.p = p1; obj.t = t; if ~isempty(obj.b) obj.b = obj.b(ind); diff --git a/utilities/extract_boundary.m b/utilities/extract_boundary.m index 8f91ddcb..0b73441a 100644 --- a/utilities/extract_boundary.m +++ b/utilities/extract_boundary.m @@ -16,8 +16,6 @@ % counter-clockwise (0) or clockwise (1). % opendat: open boundary information from a pre-existing grid % boudat: land boundary information from a pre-exist -% type: flux (1) or elevation (2) nodestring -% type2: if flux, no-flux (20) or river (2) flux nodestring % OUTPUTS: % poly: the boundary of each enclosing polygon sorted in winding-order % poly is returned as a cell-array of length number of polys. @@ -25,6 +23,8 @@ % poly % opendat: open boundary information with appened open bou % boudat: land boundary information with appended land bou +% type: flux(1) or elevation (2) type bc +% type2: if flux, zero flux (20) or non-zero flux river (22) bc % % kjr,UND,CHL,2017 % @@ -187,7 +187,7 @@ if type==1 nbou = 0; nvel = 0; nvell = []; nbbv = [] ; if ~exist('type2','var') - type2 = input('What kind of flux boundary is it, 20(island),2(River)?'); + type2 = input('What kind of flux boundary is it, 20(island),22(River)?'); end for ii = 1 : length(poly) nbou = nbou + 1; @@ -215,4 +215,4 @@ area = area + det([xp(i), xp(i+1); yp(i), yp(i+1)]); end area = 1/2*area; -end \ No newline at end of file +end From e00943d5761358cd0d81b11c3d97c3fa1ff91786 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Fri, 28 Feb 2020 18:22:11 -0600 Subject: [PATCH 5/6] Update msh.m --- @msh/msh.m | 1 + 1 file changed, 1 insertion(+) diff --git a/@msh/msh.m b/@msh/msh.m index b863de0c..2d63d702 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -1059,6 +1059,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) From 360e909acb92b34334eb2f3936799eaa8cf293c7 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Fri, 28 Feb 2020 18:34:16 -0600 Subject: [PATCH 6/6] Update msh.m --- @msh/msh.m | 49 ++++++++++++++++++++++++------------------------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/@msh/msh.m b/@msh/msh.m index 56810b9b..167e5b85 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -1856,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); @@ -1907,31 +1906,31 @@ function write(obj,fname,type) nn = setdiff(I',[bdx; bdnodes]); DTbase.Points(unique(nn),:) = []; end - + pm = DTbase.Points; tm = DTbase.ConnectivityList; - + pmid = (pm(tm(:,1),:)+pm(tm(:,2),:)+pm(tm(:,3),:))/3; - + %in1 is inside the inset boundary polygon in1 = inpoly(pmid,poly_vec1,edges1); - + %in2 is inside the global boundary polygon - in2 = inpoly(pmid,poly_vec2,edges2); - - % remove triangles that aren't in the global mesh or - % aren't in the inset mesh - del = (~in1 & ~in2); - tm(del,:) = []; + in2 = inpoly(pmid,poly_vec2,edges2); + + % remove triangles that aren't in the global mesh or + % aren't in the inset mesh + del = (~in1 & ~in2); + tm(del,:) = []; end - + merge = msh() ; merge.p = pm; merge.t = tm ; - + if ~isempty(cleanargin) - % lock anything vertices more than 2*dmax away - % from intersection - [~,dst] = ourKNNsearch(p1',p1',2); + % 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); + [~,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; @@ -1939,29 +1938,29 @@ function write(obj,fname,type) 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 merge.egfix = []; - + % put the weirs from obj1 back into merge if ~isempty(obj1.bd) && any(obj1.bd.ibtype == 24) disp('Putting the weirs from obj1 back into merged mesh') wm = msh(); wm.p = pwin1; wm.t = twin1; - merge = cat(merge,wm); + merge = cat(merge,wm); end - - % put the weirs from obj2 back into merge + + % put the weirs from obj2 back into merge if ~isempty(obj2.bd) && any(obj2.bd.ibtype == 24) disp('Putting the weirs from obj2 back into merged mesh') wm = msh(); wm.p = pwin2; wm.t = twin2; - merge = cat(merge,wm); + merge = cat(merge,wm); end - + % convert back to lat-lon wgs84 [merge.p(:,1),merge.p(:,2)] = ... m_xy2ll(merge.p(:,1),merge.p(:,2)); - + if nfix > 0 [merge.pfix(:,1),merge.pfix(:,2)] = ... m_xy2ll(pfixx(:,1),pfixx(:,2));