Commit d2663c1d authored by Pierre Cazenave's avatar Pierre Cazenave

Add support for a sponge radius based on the distance to the nearest element...

Add support for a sponge radius based on the distance to the nearest element to the specified boundary nodes (rather than a fixed radius away from the boundary). To enable, set the SpongeRadius value given to add_sponge_nodes_list.m to a negative value and then call calc_sponge_radius.m to assign sponge layer coefficients to the relevant nodes. Finally, call add_spong_nodes_list to add the nodes to the MATLAB mesh struct
parent 15fcf203
......@@ -24,18 +24,20 @@ function [Mobj] = add_sponge_nodes_list(Mobj,SpongeList,SpongeName,SpongeRadius
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Pierre Cazenave (Plymouth Marine Laboratory)
% Karen Thurston (National Oceanography Centre, Liverpool)
%
% Revision history
% Modifed from add_sponge_nodes to read in nodes from a supplied list.
% 2012-11-26 Add ability to turn off the figures.
% 2013-01-18 Added support for variable sponge radius
%
%==============================================================================
subname = 'add_sponge_nodes';
global ftbverbose
if(ftbverbose)
fprintf('\n')
fprintf(['begin : ' subname '\n'])
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end
% Do we want a figure showing how we're getting along?
......@@ -46,6 +48,7 @@ end
%------------------------------------------------------------------------------
% Plot the mesh
%------------------------------------------------------------------------------
if plotFig == 1
if(lower(Mobj.nativeCoords(1:3)) == 'car')
x = Mobj.x;
......@@ -63,7 +66,7 @@ if plotFig == 1
axis('equal','tight')
end
npts = size(SpongeList,2);
npts = length(SpongeList);
if(npts == 0)
fprintf('No points in given list')
......@@ -79,11 +82,15 @@ Mobj.nSponge = Mobj.nSponge + 1;
Mobj.nSpongeNodes(Mobj.nSponge) = npts;
Mobj.sponge_nodes(Mobj.nSponge,1:npts) = SpongeList;
Mobj.sponge_name{Mobj.nSponge} = SpongeName;
Mobj.sponge_rad(Mobj.nSponge) = SpongeRadius;
Mobj.sponge_fac(Mobj.nSponge) = SpongeCoeff;
if(ftbverbose)
fprintf(['end : ' subname '\n'])
if max(size(SpongeRadius))==1 % if you have a constant sponge radius
Mobj.sponge_rad(Mobj.nSponge) = SpongeRadius;
else % if you have a variable sponge radius
Mobj.sponge_rad(Mobj.nSponge,1:npts) = SpongeRadius;
end
if(ftbverbose)
fprintf(['end : ' subname '\n'])
end
function [spongeRadius] = calc_sponge_radius(Mobj,Nlist)
% Calculate a variable sponge radius based on distance to the boundary
% node's furthest neighbour.
% (Adapted from Phil Hall's 'produce_netcdf_input_data.py')
%
% spongeRadius = calc_sponge_radius(Mobj,Nlist)
%
% DESCRIPTION
% Calculates the sponge radius for each node on the open boundary, based
% on the minimum of either the distance to the node's furthest
% neighbour, or 100 km.
%
% INPUT
% Mobj = Matlab mesh object
% Nlist = List of nodes
%
% OUTPUT
% spongeRadius = List of variable sponge radii
%
% EXAMPLE USAGE
% spongeRadius = calc_sponge_radius(Mobj,Nlist)
%
% Author(s)
% Karen Thurston (National Oceanography Centre, Liverpool)
%
%==========================================================================
subname = 'calc_sponge_radius';
global ftbverbose
if(ftbverbose)
fprintf('\n')
fprintf(['begin : ' subname '\n'])
end
%--------------------------------------------------------------------------
% Get a unique list and make sure they are in the range of node numbers
%--------------------------------------------------------------------------
Nlist = unique(Nlist);
spongeRadius = 100000+zeros(size(Nlist));
% For each node on the open boundary
for i =1:length(Nlist)
% Find the neighbouring nodes
[r,c]=find(Mobj.tri==Nlist(i));
neighbours = unique(Mobj.tri(r,:));
% Remove the node of interest from the neighbours list
n = find(neighbours~=Nlist(i));
neighbours = neighbours(n);
% Calculate the arc length (in degrees) between the node and its
% neighbours
arclen = distance(Mobj.lat(Nlist(i)),Mobj.lon(Nlist(i)),...
Mobj.lat(neighbours),Mobj.lon(neighbours));
% Convert from degrees to whole metres
arclen = ceil(1000*deg2km(arclen));
% If the smallest distance is less than 100km, keep it
if min(arclen)<spongeRadius(i)
spongeRadius(i)=min(arclen);
end
end
if(ftbverbose)
fprintf(['end : ' subname '\n'])
end
......@@ -19,8 +19,10 @@ function write_FVCOM_sponge(Mobj,filename)
%
% Author(s):
% Geoff Cowles (University of Massachusetts Dartmouth)
% Karen Thurston (National Oceanography Centre, Liverpool)
%
% Revision history
% 2013-01-18 Added support for variable sponge radius
%
%==============================================================================
subname = 'write_FVCOM_sponge';
......@@ -48,9 +50,15 @@ else
Total_Sponge = sum(Mobj.nSpongeNodes(1:Mobj.nSponge));
fprintf(fid,'Sponge Node Number = %d\n',Total_Sponge);
for i=1:Mobj.nSponge
for j=1:Mobj.nSpongeNodes(i)
fprintf(fid,'%d %f %f \n',Mobj.sponge_nodes(i,j),Mobj.sponge_rad(i),Mobj.sponge_fac(i));
end;
if max(size(Mobj.sponge_rad))==1 % if you have a constant sponge radius
for j=1:Mobj.nSpongeNodes(i)
fprintf(fid,'%d %f %f \n',Mobj.sponge_nodes(i,j),Mobj.sponge_rad(i),Mobj.sponge_fac(i));
end;
else % if you have a variable sponge radius
for j=1:Mobj.nSpongeNodes(i)
fprintf(fid,'%d %f %f \n',Mobj.sponge_nodes(i,j),Mobj.sponge_rad(i,j),Mobj.sponge_fac(i));
end;
end
end;
end;
fclose(fid);
......
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