Commit 8696e74d authored by Ricardo Torres's avatar Ricardo Torres 💬

added nesting and hybrid_coordinate functions as well as funtions to...

 added nesting and hybrid_coordinate functions as well as funtions to calculate distance to BC and coast that uses the matlab tool HJB_Solver_Package which can be downloaded from matlab central -- or we can add it to the repository but have not done it yet
parent eb348a3b
function [Nested]=find_nesting_region(conf,Mobj)
% Creates a nesting structure array for direct/indirect or weighted nesting
%
% function [Nested]=find_nesting_region(conf,M)
%
% DESCRIPTION:
% Uses the Mesh object M and the conf variable to search for the nodes
% and elements originating from the open boundaries into de domain
% interior for a specified number of levels.
%
% Optionally specify nesting type:
% 1/2: DIRECT/INDIRECT nesting:
% - Full variables/no surface elevation respectively.
% 3: RELAXATION nesting:
% - Nesting with a relaxation method.
%
% INPUT:
% conf = struct whose field names are the model configuration options
% as set by the user. These are generally assigned at the top of a
% make_input.m type of script. Minimum fields are
% Nesting type (1, 2 == direct nesting, 3 == weighted)
% conf.type = 3;
% If we're doing type 3, we can specify the number of levels of nested
% boundaries to use. The minimum valid value is 1. For Indirect or direct
% nesting use 1
% conf.levels = 5;
% conf.power = determines drop of weights from 1 [0 is linear, 1-4 is 1/conf.levels.^conf.power]
% M = Mesh object
%
% OUTPUT:
% Nested Mesh object with added nesting variables .
%
% EXAMPLE USAGE:
% conf.type = type of nesting [1, 2 == direct nesting, 3 == weighted]
% conf.levels = number of boundary bands to use for weighted option
% conf.power = determines drop of weights from 1 [0 is linear, 1-4 is 1/conf.levels.^conf.power]
% Mobj.tri = Triangulation table as pupulated by read_sms_grid
% Mobj.x = Nodes x coordinate [we should make it possible to use lon instead]
% Mobj.y = Nodes y coordinate [we should make it possible to use lat instead]
% Mobj.nObcNodes = number of nodes as defined by SMS and read through read_sms_grid
% Mobj.read_obc_nodes = nodes indices at boundaries as read in read_sms_grid
% Mobj.obc_type = Type of OBC as defined in mod_obcs.F [ values of 1-10]
% and parsed to Mobj from conf by add_obc_nodes_list
% Mobj.xc = Nodes x coordinate [we should make it possible to use lonc instead]
% Mobj.yc = Elements y coordinate [we should make it possible to use latc instead]
% Mobj.obc_nodes = matrix with node indices. Each row is a boundary level.
% Mobj.nObs = total number of open boundary levels (I think this is set
% in setup_metrics.m
%
% the global variable ftbverbose shows information at run time as well as
% generating a figure with the location of the nesting nodes and elements
%
% [Nested]=find_nesting_region(conf,Mobj)
%
% Nested.nObcNodes = number of nodes in new nesting zone
% Nested.read_obc_nodes = nodes indices in nesting zone
% Nested.obc_type = Type of OBC as defined in mod_obcs.F [ values of 1-10]
% Nested.obc_nodes = matrix with node indices. Each row is a boundary level.
% Nested.nObs = total number of open boundary levels
% Nested.weight_node = weights for nodes if using weighted type [0-1] see
% FVCOM manual for further info (Chapter 6 section 4 in version 3.2)
% Nested.weight_cell = weights for elements if using weighted type [0-1]
% Author(s):
% Ricardo Torres (Plymouth Marine Laboratory)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Darren Price (CH2MHill)
% Hakeem Johnson (CH2MHill)
%
% Revision history:
% 2015-11-01 First version based on Hakeem and Darren code as provided to
% Torres by Pierre
% script.
% 2016-01-19 Updated to a stand alone function and general tidy up
% general tidy up.
%
%==========================================================================
dump = dbstack;
subname = dump.name;
clear dump
global ftbverbose
if ftbverbose
fprintf('\nbegin : %s\n', subname)
end
M=Mobj;
TR = triangulation(M.tri, [M.x, M.y]);
et = cell(conf.levels, 1);
M.weight_cell = cell(1,length(M.nObcNodes)*conf.levels);
for oo=1:length(M.read_obc_nodes)
M.weight_node{oo} = ones(1, M.nObcNodes(oo));
end
nn=0;
if ftbverbose
figure,cla
patch('Vertices', [M.x, M.y], 'Faces', M.tri, ...
'Cdata', -M.h, 'edgecolor', 'b', 'facecolor', 'interp');
end
%% The initial OB are first in the order. The nested OB
% follow
etdx=0;
for oo=1:length(Mobj.read_obc_nodes)
for n = 1:conf.levels
nn=nn+1;etdx=etdx+1;
if (oo>1 && n == 1)
M.read_obc_nodes{nn+1} = Mobj.read_obc_nodes{oo};
M.obc_type(nn+1)=Mobj.obc_type(oo);
M.nObcNodes(nn+1) = length(Mobj.read_obc_nodes{oo});
M.obc_nodes(nn+1, 1:Mobj.nObcNodes(oo)) = Mobj.read_obc_nodes{oo};
ti = vertexAttachments(TR, double(M.read_obc_nodes{nn+1})');
et{etdx} = setdiff(unique([ti{:}]),[et{1:end}]);
nn=nn+1;
M.read_obc_nodes{nn+1} = int32(setdiff(unique(M.tri(et{etdx}, :)), ...
[M.read_obc_nodes{1:nn}]))';
M.obc_type(nn+1)=M.obc_type(nn);
M.nObs = M.nObs+1;
M.nObcNodes(nn+1) = length(M.read_obc_nodes{nn+1});
M.obc_nodes(nn+1, 1:M.nObcNodes(nn+1)) = M.read_obc_nodes{nn+1};
else
ti = vertexAttachments(TR, double(M.read_obc_nodes{nn})');
et{etdx} = setdiff(unique([ti{:}]),[et{1:end}]);
M.read_obc_nodes{nn+1} = int32(setdiff(unique(M.tri(et{etdx}, :)), ...
[M.read_obc_nodes{1:nn}]))';
M.obc_type(nn+1)=M.obc_type(nn);
M.nObs = M.nObs+1;
M.nObcNodes(nn+1) = length(M.read_obc_nodes{nn+1});
M.obc_nodes(nn+1, 1:M.nObcNodes(nn+1)) = M.read_obc_nodes{nn+1};
end
% Plot nesting region.
if ftbverbose
if n == 1 | nn==1
axis('equal', 'tight')
colormap('gray')
hold on
plot(M.x(M.read_obc_nodes{nn}), M.y(M.read_obc_nodes{nn}), 'wo')
plot(M.x(M.read_obc_nodes{nn+1}), M.y(M.read_obc_nodes{nn+1}), 'ro')
plot(M.xc(et{etdx}), M.yc(et{etdx}), 'wx')
else
plot(M.x(M.read_obc_nodes{nn+1}), M.y(M.read_obc_nodes{nn+1}), 'ro')
plot(M.xc(et{etdx}), M.yc(et{etdx}), 'wx')
end
end
end
end
if conf.type == 3
% probably best to have them as a linear decrease?
if conf.power ==0
% one should be the outermost boundary level which is why the vectors
% are flipped
weights_nodes = fliplr((1:conf.levels)./conf.levels);
weights_nodes(end+1) =0;
weights_elems = fliplr((1:conf.levels-1)./conf.levels-1);
weights_elems(end+1) =0;
else
weights_nodes = 1:conf.levels+1;
weights_nodes = 1./weights_nodes.^conf.power;
weights_elems = 1:conf.levels;
weights_elems = 1./weights_elems.^conf.power;
end
nn=0;etdx=0;
for oo=1:length(Mobj.read_obc_nodes)
nn=nn+1;
M.weight_node{nn} = repmat(weights_nodes(1), ...
1, M.nObcNodes(nn));
for n = 1:conf.levels
etdx=etdx+1;
M.weight_node{nn+1} = repmat(weights_nodes(n+1), ...
1, M.nObcNodes(nn+1));
M.weight_cell{etdx} = repmat(weights_elems(n), ...
1, length(et{etdx}), 1);
nn=nn+1;
end
end
end
Nested=M;
Nested.read_obc_elems=et;
return
% calculate and write hybrid coordinate layer
function Mobj=hybrid_coordinate(conf,Mobj)
% defaults
% conf.sigma_file='coord_hybrid.sig'
% conf.nsigma=41;DU=25;DL=25;Hmin=200;KU=5;KL=5;ZKU=[.5 .5 .5 .5 .5];ZKL=[.5 .5 .5 .5 .5 ];
% conf.H0=100;conf.nlev=20;conf.DU=25;conf.DL=25;conf.KU=5;conf.KL=5;
myset=optimset('MaxFunEvals',5000,'MaxIter',5000,'TolFun',10e-5,'TolX',1e-7);
% solve for z1-z2 to find Hmin parameter
H0=conf.H0;
nlev=conf.nlev;DU=conf.DU;DL=conf.DL;
KU=conf.KU;KL=conf.KL;
ZKU=repmat(DU./KU,1,KU);ZKL=repmat(DL./KL,1,KL);
myfunparams=@(H)hybrid_coordinate_hmin(H,nlev,DU,DL,KU,KL,ZKU,ZKL);
% myfunparams(200)
[Hmin,fval] = fminsearch(myfunparams,H0,myset);
% print hybrid_coord.sigma file
fout=fopen(conf.sigma_file,'wt');
fprintf(fout,'%s\n',['NUMBER OF SIGMA LEVELS = ',num2str(nlev)]);
fprintf(fout,'%s\n',['SIGMA COORDINATE TYPE = GENERALIZED']);
fprintf(fout,'DU = %4.1f\n',(DU));
fprintf(fout,'DL = %4.1f\n',(DL));
fprintf(fout,'MIN CONSTANT DEPTH = %10.1f\n',round(Hmin));
fprintf(fout,'%s\n',['KU = ',num2str(KU)]);
fprintf(fout,'%s\n',['KL = ',num2str(KL)]);
fprintf(fout,'ZKU = ');
for ii=1:length(ZKU)
fprintf(fout,'%4.1f',ZKU(ii));
end
fprintf(fout,'\n');
fprintf(fout,'ZKL = ');
for ii=1:length(ZKU)
fprintf(fout,'%4.1f',ZKL(ii));
end
fprintf(fout,'\n');
fclose(fout);
for xx=1:length(Mobj.h)
% loop through all nodes to create sig_coor
Mobj.siglev(xx,:) = sigma_gen(nlev,DL,DU,KL,KU,ZKL,ZKU,Mobj.h(xx),Hmin);
end
Mobj.siglay = Mobj.siglev(:,1:end-1) + (diff(Mobj.siglev,1,2)/2);
for zz=1:nlev-1
Mobj.siglevc(:,zz) = nodes2elems(Mobj.siglev(:,zz), Mobj);
% Create a siglay variable (i.e. midpoint in the sigma levels).
Mobj.siglayc(:,zz) = nodes2elems(Mobj.siglay(:,zz), Mobj);
end
Mobj.siglevc(:,zz+1) = nodes2elems(Mobj.siglev(:,zz+1), Mobj);
% Create a depth array for the element centres.
if ~isfield(Mobj,'hc')
Mobj.hc = nodes2elems(Mobj.h, Mobj);
end
Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* Mobj.siglev;
Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* Mobj.siglay;
Mobj.siglevzc = repmat(Mobj.hc, 1, nlev) .* Mobj.siglevc;
Mobj.siglayzc = repmat(Mobj.hc, 1, nlev-1) .* Mobj.siglayc;
% % Create a siglay variable (i.e. midpoint in the sigma levels).
% zlay = z(1:end-1) + (diff(z)/2);
%
% % Create a depth array for the element centres.
% hc = nodes2elems(Mobj.h, Mobj);
%
% Mobj.siglevz = repmat(Mobj.h, 1, nlev) .* repmat(z, Mobj.nVerts, 1);
% Mobj.siglayz = repmat(Mobj.h, 1, nlev-1) .* repmat(zlay, Mobj.nVerts, 1);
% Mobj.siglevzc = repmat(hc, 1, nlev) .* repmat(z, Mobj.nElems, 1);
% Mobj.siglayzc = repmat(hc, 1, nlev-1) .* repmat(zlay, Mobj.nElems, 1);
%
% % Add the sigma levels and layers to the Mobj.
% Mobj.siglev = z;
% Mobj.siglay = zlay;
return
%
% % test with made up bathymetr
% Hmax=200;
% y=0:100;B=100;
% H = Hmax .*exp( -((y./B-0.15).^2 ./(0.5).^2) );
% H=120
% Z0=zeros(length(H),nsigma);Z2=zeros(length(H),nsigma);
% for hh=1:length(H)
% Z0(hh,1)=0;DL2=0.001;DU2=0.001;KBM1=nsigma-1;
% for nn=1:nsigma-1
% X1=DL2+DU2;
% X1=X1*(KBM1-nn)/KBM1;
% X1=X1-DL2;
% X1=tanh(X1);
% X2=tanh(DL2);
% X3=X2+tanh(DU2);
%
% Z0(hh,nn+1)=((X1+X2)/X3)-1.0;
% end
%
% % s-coordinates
% X1=(H(hh)-DU-DL);
% X2=X1./H(hh);
% DR=X2./(nsigma-KU-KL-1);
%
% Z2(hh,1)=0.0;
%
% for K=2:KU+1
% Z2(hh,K)=Z2(hh,K-1)-(ZKU(K-1)./H(hh));
% end
%
% for K=KU+2:nsigma-KL
% Z2(hh,K)=Z2(hh,K-1)-DR;
% end
%
% KK=0;
% for K=nsigma-KL+1:nsigma
% KK=KK+1;
% Z2(hh,K)=Z2(hh,K-1)-(ZKL(KK)./H(hh));
% end
% end
function Mobj=distance_along_BC(Mobj,BCnodes)
CD=pwd;
% setup HPJ solver to calculate the distance function for the SMS mesh
[~,~,~,bnd] = connectivity([Mobj.x,Mobj.y],Mobj.tri);
% remove nodestring from coast.
% BCnodes=[Mobj.read_obc_nodes{:}];
coast_ind=(BCnodes);
% % calculate distance function
myParam.Max_Tri_per_Star = 20;
myParam.NumGaussSeidel = 40;
myParam.INF_VAL = 1000000;
myParam.TOL = 1e-12;
% in this case, we will assume the standard Euclidean metric
myMetric = [];
myTM.Vtx=[Mobj.x(:),Mobj.y(:)];
myTM.DoFmap=[Mobj.tri];
myTM.NegMask=false(size(bnd));
myBdy.Nodes=coast_ind(:);
myBdy.Data=zeros(size(myBdy.Nodes));
%
cd /users/modellers/rito/matlab/HJB_Solver_Package/HJB_Solver_Package
%
%
SEmex = SolveEikonalmex(myTM,myBdy,myParam,myMetric);
tic
Mobj.distOB = SEmex.Compute_Soln;
cd(CD)
end
......@@ -2,11 +2,12 @@ function Mobj=distance_to_coast(Mobj)
CD=pwd;
% setup HPJ solver to calculate the distance function for the SMS mesh
[~,~,~,bnd] = connectivity([Mobj.x,Mobj.y],Mobj.tri);
pv = [Mobj.x(find(bnd)) Mobj.y(find(bnd));Mobj.x(find(bnd(1))) Mobj.y(find(bnd(1)))];
% remove nodestring from coast.
BCnodes=Mobj.read_obc_nodes{1};
BCnodes=[Mobj.read_obc_nodes{:}];
bnd(BCnodes)=0;
coast_ind=find(bnd);
coast_ind(1:max(BCnodes))=[];
% % calculate distance function
myParam.Max_Tri_per_Star = 20;
......@@ -26,7 +27,7 @@ cd /users/modellers/rito/matlab/HJB_Solver_Package/HJB_Solver_Package
%
SEmex = SolveEikonalmex(myTM,myBdy,myParam,myMetric);
tic
Mobj.dist = SEmex.Compute_Soln;
Mobj.dist =SEmex.Compute_Soln;
cd(CD)
end
......@@ -134,7 +134,7 @@ switch ND
end
end
for jj=1:nLayers
[Plots(plotOPTS.figure).handles{jj}]=quiverm(double(plotOPTS.mesh.latc(igood)),double(plotOPTS.mesh.lonc(igood)),u,v,'k',plotOPTS.vel_sca,'filled'),
[Plots(plotOPTS.figure).handles{jj}]=quiver(double(plotOPTS.mesh.lonc(igood)),double(plotOPTS.mesh.latc(igood)),u,v,plotOPTS.vel_sca,'k','filled'),
% [Plots(plotOPTS.figure).handles{jj}]=quiverm(double(FVCOM.latc(igood)),double(FVCOM.lonc(igood)),u./100,v./100,'k',0,'filled'),
% [Plots(plotOPTS.figure).handles(jj),~]=m_vec(plotOPTS.vel_sca,x(igood),y(igood),...
% u(:,jj),v(:,jj),colourSpec(jj,:),...
......
......@@ -60,11 +60,12 @@ for i = 1:nb
[C2,~] = ismember(Mobj.tri(:,2), nodeIDs(:),'rows');
[C3,~] = ismember(Mobj.tri(:,3), nodeIDs(:),'rows');
obc_elems{i}= unique([find(C1);find(C2);find(C3)]);
if iscolumn( obc_elems{i});obc_elems{i}=obc_elems{i}';end
nObcElements(i) = numel(obc_elems{i}(:));
end
Mobj.relaxBC_nodes={reshape([Mobj.read_obc_nodes{:}],[],1)};
Mobj.relaxBC_elems={reshape([obc_elems{:}],[],1)};
Mobj.relaxBC_nodes={[Mobj.read_obc_nodes{:}]};
Mobj.relaxBC_elems={[obc_elems{:}]};
for bb=2:bc_width
nodeIDs = Mobj.tri(Mobj.relaxBC_elems{bb-1},:);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment