From e0c466b31743a408b60f0821101ad8edba335670 Mon Sep 17 00:00:00 2001 From: William Pringle Date: Thu, 31 Jan 2019 12:34:45 -0500 Subject: [PATCH] Added bathy slope limiter method - msh class now has lim_bathy_slope function to limit bathymetric or topographic slope to user defined value (e.g. 0.25). Uses Mesh2D limgrad function --- @msh/msh.m | 57 +++++++++++++---- utilities/limgrad.m | 146 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 190 insertions(+), 13 deletions(-) create mode 100644 utilities/limgrad.m diff --git a/@msh/msh.m b/@msh/msh.m index e17092ce..186fd32b 100644 --- a/@msh/msh.m +++ b/@msh/msh.m @@ -297,9 +297,10 @@ function write(obj,fname,type) m_trimesh(obj.t,obj.p(:,1),obj.p(:,2),q); %m_trisurf(obj.t,obj.p(:,1),obj.p(:,2),q); else - %trimesh(obj.t,obj.p(:,1),obj.p(:,2),q); - trisurf(obj.t,obj.p(:,1),obj.p(:,2),q) - view(2); shading interp; + trimesh(obj.t,obj.p(:,1),obj.p(:,2),q,'facecolor','none'); + %trisurf(obj.t,obj.p(:,1),obj.p(:,2),q) + view(2); + %shading interp; end if exist('demcmap','file') demcmap(q) @@ -696,6 +697,33 @@ function write(obj,fname,type) end end + function obj = lim_bathy_slope(obj,dfdx,overland) + if nargin < 3 + overland = 0; + end + % Limit to topo or bathymetric slope to dfdx on the edges + imax = 100; + [edge,elen] = GetBarLengths(obj); + bt = obj.b; + if overland + I = bt > 0; + word = 'topographic'; + else + I = bt < 0; + word = 'bathymetric'; + end + bt(I) = 0; + [bnew,flag] = limgrad(edge,elen,bt,dfdx,imax); + if flag + obj.b(~I) = bnew(~I); + disp(['Successfully limited ' word ' slope to ' ... + num2str(dfdx) ' in limgrad function']) + else + warning(['Could not limit ' word ' slope to ' ... + num2str(dfdx) ' in limgrad function']) + end + end + % make nodestrings function obj = makens(obj,type,dir,cutlim,depthlim) if nargin < 2 @@ -1269,16 +1297,7 @@ function write(obj,fname,type) function [out1,barlen,bars] = CalcCFL(obj,dt) g = 9.81; % gravity - 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 - long = zeros(length(bars)*2,1); - lat = zeros(length(bars)*2,1); - long(1:2:end) = obj.p(bars(:,1),1); - long(2:2:end) = obj.p(bars(:,2),1); - lat(1:2:end) = obj.p(bars(:,1),2); - lat(2:2:end) = obj.p(bars(:,2),2); - % Get spherical earth distances for bars - barlen = m_lldist(long,lat); barlen = barlen(1:2:end)*1e3; % L = Bar lengths in meters + [bars,barlen] = GetBarLengths(obj); % sort bar lengths in ascending order [barlen,IA] = sort(barlen,'ascend'); bars = bars(IA,:); @@ -1305,6 +1324,18 @@ function write(obj,fname,type) end end + function [bars,barlen] = GetBarLengths(obj) + 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 + long = zeros(length(bars)*2,1); + lat = zeros(length(bars)*2,1); + long(1:2:end) = obj.p(bars(:,1),1); + long(2:2:end) = obj.p(bars(:,2),1); + lat(1:2:end) = obj.p(bars(:,1),2); + lat(2:2:end) = obj.p(bars(:,2),2); + % Get spherical earth distances for bars + barlen = m_lldist(long,lat); barlen = barlen(1:2:end)*1e3; % L = Bar lengths in meters + end function obj = CheckTimestep(obj,dt,varargin) %% Decimate mesh to achieve CFL condition for stability. diff --git a/utilities/limgrad.m b/utilities/limgrad.m new file mode 100644 index 00000000..1e48bf0c --- /dev/null +++ b/utilities/limgrad.m @@ -0,0 +1,146 @@ +function [ffun,flag] = limgrad(edge,elen,ffun,dfdx,imax) +%LIMGRAD impose "gradient-limits" on a function defined over +%an undirected graph. +% [FNEW] = LIMGRAD(EDGE,ELEN,FFUN,DFDX,ITER) computes a +% "gradient-limited" function FNEW on the undirected graph +% {EDGE,ELEN}, where EDGE is an NE-by-2 array of edge ind- +% ices, and ELEN is an NE-by-1 array of edge lengths. +% Gradients are limited over the graph edges, such that +% +% ABS(FNEW(N2)-FNEW(N1)) <= ELEN(II) * DFDX, +% +% where N1=EDGE(II,1) and N2=EDGE(II,2) are the two nodes +% in the II-TH edge. An iterative algorithm is used, swee- +% ping over an "active-set" of graph edges until converge- +% nce is achieved. A maximum of IMAX iterations are done. +% +% [FNEW,FLAG] = LIMGRAD(...) also returns a boolean FLAG, +% with FLAG=TRUE denoting convergence. +% +% See also LIMHFN2 + +% Darren Engwirda : 2017 -- +% Email : engwirda@mit.edu +% Last updated : 18/04/2017 + +%---------------------------------------------- basic checks + if ( ~isnumeric(edge) || ... + ~isnumeric(elen) || ... + ~isnumeric(ffun) || ... + ~isnumeric(dfdx) || ... + ~isnumeric(imax) ) + error('limgrad:incorrectInputClass' , ... + 'Incorrect input class.') ; + end +%---------------------------------------------- basic checks + if (ndims(edge) ~= +2 || ... + ndims(elen) > +2 || ... + ndims(ffun) > +2 || ... + numel(dfdx) ~= +1 || ... + numel(imax) ~= +1 ) + error('limgrad:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + if (size(edge,2) < +2 || ... + size(elen,2)~= +1 || ... + size(ffun,2)~= +1 || ... + size(edge,1)~= size(elen,1) ) + error('limgrad:incorrectDimensions' , ... + 'Incorrect input dimensions.'); + end + + nnod = size(ffun,1) ; + +%---------------------------------------------- basic checks + if (dfdx < +0. || imax < +0) + error('limgrad:invalidInputArgument', ... + 'Invalid input parameter.'); + end + if (min(min(edge(:,1:2))) < +1 || ... + max(max(edge(:,1:2))) > nnod ) + error('limgrad:invalidInputArgument', ... + 'Invalid EDGE input array.') ; + end + +%-- IVEC(NPTR(II,1):NPTR(II,2)) are edges adj. to II-TH node + nvec = [edge(:,1); edge(:,2)]; + ivec = [(1:size(edge,1))'; ... + (1:size(edge,1))'] ; + + [nvec,pidx] = sort (nvec) ; + ivec = ivec (pidx) ; + + mark = false(nnod,1) ; + mark(edge(:,1)) = true ; + mark(edge(:,2)) = true ; + + idxx = find(diff(nvec) > +0) ; + + nptr = zeros(nnod,2) ; + nptr(:,2) = -1 ; + nptr(mark,1) = [+1; idxx+1]; + nptr(mark,2) = [idxx; nnod]; + +%----------------------------- ASET=ITER if node is "active" + aset = zeros(size(ffun,1),1) ; + +%----------------------------- exhaustive 'til all satisfied + ftol = min(ffun) * sqrt(eps) ; + + for iter = +1 : imax + + %------------------------- find "active" nodes this pass + aidx = find(aset == iter - 1) ; + + if (isempty(aidx)), break; end + + %------------------------- reorder => better convergence + [~,idxx] = sort(ffun(aidx)) ; + + aidx = aidx(idxx); + + %------------------------- visit adj. edges and set DFDX + for ipos = 1 : length(aidx) + npos = aidx(ipos) ; + for jpos = nptr(npos,1) ... + : nptr(npos,2) + + epos = ivec(jpos,1) ; + + nod1 = edge(epos,1) ; + nod2 = edge(epos,2) ; + + %----------------- calc. limits about min.-value + if (ffun(nod1) > ffun(nod2)) + + fun1 = ffun(nod2) ... + + elen(epos) * dfdx ; + + if (ffun(nod1) > fun1+ftol) + ffun(nod1) = fun1; + aset(nod1) = iter; + end + + else + + fun2 = ffun(nod1) ... + + elen(epos) * dfdx ; + + if (ffun(nod2) > fun2+ftol) + ffun(nod2) = fun2; + aset(nod2) = iter; + end + + end + + end + end + + end + + flag = (iter < imax) ; + +end + + +