Skip to content

Commit

Permalink
Added bathy slope limiter method
Browse files Browse the repository at this point in the history 
- 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
  • Loading branch information
@WPringle
WPringle committed Jan 31, 2019
1 parent b05a3d9 commit e0c466b
Show file tree
Hide file tree
Showing 2 changed files with 190 additions and 13 deletions.
57 changes: 44 additions & 13 deletions @msh/msh.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,:);
Expand All @@ -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.
Expand Down
146 changes: 146 additions & 0 deletions utilities/limgrad.m
View file
Original file line number Diff line number Diff line change
@@ -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 : [email protected]
% 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



0 comments on commit e0c466b

Please sign in to comment.