Difference between revisions of "ADCIRC"

From GRASS-Wiki
Jump to: navigation, search
(Interfacing with the ADCIRC coastal ocean model: more pre-processing helper script)
Line 16: Line 16:
 
You can extract the tidal forcing factor amplitude and phase for boundary nodes needed in the fort.15 setup file from raster maps already loaded into GRASS. In the following example the FES2004 dataset is used. (FES2002 for the K2 constituent as K2 is reputed to be broken in the 2004 edition)
 
You can extract the tidal forcing factor amplitude and phase for boundary nodes needed in the fort.15 setup file from raster maps already loaded into GRASS. In the following example the FES2004 dataset is used. (FES2002 for the K2 constituent as K2 is reputed to be broken in the 2004 edition)
  
<pre>
+
<source lang="bash">
 
constits="M2
 
constits="M2
 
N2
 
N2
Line 56: Line 56:
 
# Manually replace with values from nearby raster cells using d.what.rast
 
# Manually replace with values from nearby raster cells using d.what.rast
 
grep '*' *.tidal.prn
 
grep '*' *.tidal.prn
</pre>
+
</source>
  
 
Then, in Matlab or Octave:
 
Then, in Matlab or Octave:
  
<pre>
+
<source lang="matlab">
 
M2 = load('fes2004.M2.tidal.prn');
 
M2 = load('fes2004.M2.tidal.prn');
 
N2 = load('fes2004.N2.tidal.prn');
 
N2 = load('fes2004.N2.tidal.prn');
Line 118: Line 118:
  
 
% then manually fix c2N2 -> 2N2
 
% then manually fix c2N2 -> 2N2
</pre>
+
</source>
  
 
=== Post processing output data ===
 
=== Post processing output data ===
Line 128: Line 128:
 
* <tt>v.in.adcirc_grid</tt> does not yet support loading in node type (boundary, etc.) as point attributes -- vector tables can be a bit slow for maps with potentially millions of features. In the mean time here is a little shell script to extract the boundary nodes on the ''first'' open boundary into a separate vector points map:
 
* <tt>v.in.adcirc_grid</tt> does not yet support loading in node type (boundary, etc.) as point attributes -- vector tables can be a bit slow for maps with potentially millions of features. In the mean time here is a little shell script to extract the boundary nodes on the ''first'' open boundary into a separate vector points map:
  
 +
<source lang="bash">
 
  infile=fort.14
 
  infile=fort.14
 
  inmap=mesh_grid_pts
 
  inmap=mesh_grid_pts
Line 137: Line 138:
 
     tail -n +2 | tr '\n' ',' | sed -e 's/ //g' -e 's/,$//'`
 
     tail -n +2 | tr '\n' ',' | sed -e 's/ //g' -e 's/,$//'`
 
  v.extract in=$inmap out=$bdymap list="$cat_list"
 
  v.extract in=$inmap out=$bdymap list="$cat_list"
 +
</source>
  
 
==== Importing elevation or velocity reporting stations ====
 
==== Importing elevation or velocity reporting stations ====
Line 147: Line 149:
 
If you wish to have a look at where the sub-domains of a parallel run are, you can use the following Bash script:
 
If you wish to have a look at where the sub-domains of a parallel run are, you can use the following Bash script:
  
<pre>
+
<source lang="bash">
 
mkdir fort14s
 
mkdir fort14s
 
for DIR in PE00?? ; do
 
for DIR in PE00?? ; do
Line 160: Line 162:
 
   v.in.adcirc_grid -p in="$file" out="submesh_${base}_pts" --quiet
 
   v.in.adcirc_grid -p in="$file" out="submesh_${base}_pts" --quiet
 
done
 
done
</pre>
+
</source>
  
 
Here is a loop for displaying 80 sub-meshes, with random dark colors:
 
Here is a loop for displaying 80 sub-meshes, with random dark colors:
<pre>
+
 
 +
<source lang="bash">
 
for PE in {00..79} ; do
 
for PE in {00..79} ; do
 
   echo -n "[$PE]"
 
   echo -n "[$PE]"
Line 176: Line 179:
  
 
d.vect station_map color=black fcolor=red size=8 icon=basic/circle
 
d.vect station_map color=black fcolor=red size=8 icon=basic/circle
</pre>
+
</source>
  
 
==== Importing unit 61, 62, 63, 64 results ====
 
==== Importing unit 61, 62, 63, 64 results ====

Revision as of 19:03, 25 July 2019

Interfacing with the ADCIRC coastal ocean model

Grid preparation

OceanMesh2D

OceanMesh2D is an open source (GPL3) generator of 2D triangular meshes suitable for use with ADCIRC. It likes to work in lat/lon coordinate system.

You can output your DEM as NetCDF format with r.out.gdal.

You can create a domain boundary polyline with r.contour. You can then output the domain boundary polyline either as Shapefiles with v.out.ogr, or as two column x,y ASCII data with the v.out.ascii.mat addon module for GRASS 6.

Nodal attribute preparation

You can extract the tidal forcing factor amplitude and phase for boundary nodes needed in the fort.15 setup file from raster maps already loaded into GRASS. In the following example the FES2004 dataset is used. (FES2002 for the K2 constituent as K2 is reputed to be broken in the 2004 edition)

constits="M2
N2
S2
K1
O1
P1
K2
2N2
Q1
M4"

for const in $constits ; do
  MAP=fes2004.$const
  if [ $const = K2 ] ; then
    MAP=fes2002.K2
  fi
  echo "[$MAP]"

  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.amplitude -pi --q | sort -n > $MAP.ampl.prn

  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.phase -p --q | sort -n >  $MAP.phase_raw.prn
  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.phase.u -pi --q | sort -n >  $MAP.phase.u.prn
  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.phase.v -pi --q | sort -n >  $MAP.phase.v.prn

  echo "% $const: grid open boundary nodes" > $MAP.tidal.prn
  echo "% node ampl phase phase.u phase.v" >> $MAP.tidal.prn
  paste -d'|' $MAP.ampl.prn $MAP.phase_raw.prn \
      $MAP.phase.u.prn $MAP.phase.v.prn | \
    cut -f1,2,4,6,8 -d'|' | \
    tr '|' '\t' >> $MAP.tidal.prn
done

# check for nulls in the source data.
# Manually replace with values from nearby raster cells using d.what.rast
grep '*' *.tidal.prn

Then, in Matlab or Octave:

M2 = load('fes2004.M2.tidal.prn');
N2 = load('fes2004.N2.tidal.prn');
S2 = load('fes2004.S2.tidal.prn');
K1 = load('fes2004.K1.tidal.prn');
O1 = load('fes2004.O1.tidal.prn');
K2 = load('fes2002.K2.tidal.prn');
P1 = load('fes2004.P1.tidal.prn');
Q1 = load('fes2004.Q1.tidal.prn');

c2N2 = load('fes2004.2N2.tidal.prn');
M4 = load('fes2004.M4.tidal.prn');


[M2 cart2pol(M2(:,4), M2(:,5)) * 180/pi]

delta = M2(:,3) - (cart2pol(M2(:,4), M2(:,5)) * 180/pi);
[M2(:,1) delta]
plot(delta), grid
xlabel('boundary node')
ylabel('degrees')

[M2(:,1) M2(:,3) (cart2pol(M2(:,4), M2(:,5)) * 180/pi) delta]


%%%
constits={
'M2'
'N2'
'S2'
'K1'
'O1'
'P1'
'K2'
'c2N2'
'Q1'
'M4'}

fid = fopen('grid_constits.prn', 'w')

for i = 1: length(constits)
  const = char(constits(i));
  eval(['data = ' const ';'])
  new_phase = cart2pol(data(:,4), data(:,5)) * 180/pi;
  for j = 1:length(new_phase)
    if (new_phase(j) < 0)
      new_phase(j) = new_phase(j) + 360;
    end
  end
  fprintf(fid, ' %s\n', const);
  for j = 1:length(new_phase)
    fprintf(fid, '   %8.6f  %9.5f\n', data(j,2), new_phase(j));
  end
end

fclose(fid)

% then manually fix c2N2 -> 2N2

Post processing output data

Importing the fort.14 mesh grid

The v.in.adcirc_grid addon module for GRASS 6 will import a 3D triangular mesh from the ADCIRC coastal ocean model into a GRASS vector map.

  • v.in.adcirc_grid does not yet support loading in node type (boundary, etc.) as point attributes -- vector tables can be a bit slow for maps with potentially millions of features. In the mean time here is a little shell script to extract the boundary nodes on the first open boundary into a separate vector points map:
 infile=fort.14
 inmap=mesh_grid_pts
 bdymap=mesh_bdy_pts
 
 v.in.adcirc_grid -p in=$infile out=$inmap
 num_pts=`grep ' = Number of nodes for open boundary' "$infile" | head -n 1 | awk '{print $1}'`
 cat_list=`grep -A"$num_pts" ' = Number of nodes for open boundary' "$infile" | \
    tail -n +2 | tr '\n' ',' | sed -e 's/ //g' -e 's/,$//'`
 v.extract in=$inmap out=$bdymap list="$cat_list"

Importing elevation or velocity reporting stations

With a little awk you can load in the elevation reporting stations for unit 61 with v.in.ascii.

example goes here

Visualizing sub-domains

If you wish to have a look at where the sub-domains of a parallel run are, you can use the following Bash script:

mkdir fort14s
for DIR in PE00?? ; do
   cp -a $DIR/fort.14 fort14s/${DIR}.14
done
 
cd fort14s
for file in *.14 ; do
   base=`echo "$file" | sed -e 's/PE00//' -e 's/\.14//'`
   echo -n "[$base]"
   v.in.adcirc_grid in="$file" out="submesh_$base" --quiet
   v.in.adcirc_grid -p in="$file" out="submesh_${base}_pts" --quiet
done

Here is a loop for displaying 80 sub-meshes, with random dark colors:

for PE in {00..79} ; do
  echo -n "[$PE]"
  R=`expr $RANDOM / 256`
  G=`expr $RANDOM / 256`
  B=`expr $RANDOM / 256`

  d.vect "submesh_$PE" color="$R:$G:$B"
done

d.vect US_medium_shoreline color=black width=2

d.vect station_map color=black fcolor=red size=8 icon=basic/circle

Importing unit 61, 62, 63, 64 results

- addme: bash code + r.in.xyz + r.surf.nnbathy

d.barb addon module for GRASS 6 for rendering u,v vector arrows