Skip to content

Commit

Permalink
fix up python, matlab in tutorials/doc (#386)
Browse files Browse the repository at this point in the history
* fix up python, matlab in tutorials/doc

* @jahn fixup, add BC gendata.py, minor format edits

* minor edit to baroclinic tut README

* makeover of the python analysis code, not using xarray

* fix matlab/python codes; fix mnc_test_-text bits; describe stat-diags output

* replace plots (better colormaps, we had changed a parm and no longer matched up)

* add xarray version of python analysis code, edit other versions

* few minor tweaks, replace streamfn plot

* do rename Zmd00015->Z in .py file as per discussion

* more specific discussion/use of numpy and xarray broadcasting

* read time axis rather than hard-code; also pcolor discussion

* more cleanup, remove some hard-coded numbers

* move matlab colormaps to utils/matlab

* new version of python code with shorter lines

* remove inadventent tabs inserted in several lines

* replace hard-coded 62s with Nx or Ny

* minor edits

* edited and reformatted xarray-version python script

* code improvements, remove additional hard-coding

* minor edits to comments

* move comment about xgcm to more relevant place

* fix up matlab_plots code to match up with python script

* additional formatting cleanup

* minor comment edits, including note about figures sizing

* add some comments to BT gyre genmake scripts

* add comments to gendata files, make code clearer

* rewrite/make clearer section on MITgcm file formats (3.9)

* edit line spacing to improve readability

* fix python syntax

* make matlab, pthyon codes in tutorial easier to read

* add some spaces to increase .py and .m code readability

* No "degree Kelvin", just "K" for Kelvin

* minor edits following review

* move README to README.md

* Update README.md

* Update README.md

Co-authored-by: Jean-Michel Campin <jmc@ocean.mit.edu>
  • Loading branch information
jrscott and jm-c committed Apr 21, 2021
1 parent d814156 commit ce0d9af
Show file tree
Hide file tree
Showing 18 changed files with 1,457 additions and 339 deletions.
102 changes: 73 additions & 29 deletions doc/examples/baroclinic_gyre/baroclinic_gyre.rst

Large diffs are not rendered by default.

Binary file modified doc/examples/baroclinic_gyre/figs/baroc_psi.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified doc/examples/baroclinic_gyre/figs/baroc_thetaplots.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified doc/examples/baroclinic_gyre/figs/trelax_freesurface.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified doc/examples/baroclinic_gyre/figs/trelax_theta_timeseries.png
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
37 changes: 25 additions & 12 deletions doc/examples/barotropic_gyre/barotropic_gyre.rst
Expand Up @@ -517,7 +517,8 @@ value is set to 0 (i.e., “land”). As discussed in :numref:`sec_baro_num_co
(i.e., periodic in both :math:`x`- and :math:`y`-directions), so boundary walls
are necessary to set up our enclosed box domain.
The matlab program :filelink:`verification/tutorial_barotropic_gyre/input/gendata.m`
was used to generate this bathymetry file. By default, this file is assumed to
was used to generate this bathymetry file (alternatively, see python equivalent
:filelink:`gendata.py <verification/tutorial_barotropic_gyre/input/gendata.py>`). By default, this file is assumed to
contain 32-bit (single precision) binary numbers.
See :numref:`sec_mitgcm_inp_file_format` for additional information on MITgcm input data file format specifications.

Expand All @@ -532,9 +533,11 @@ The units are Nm\ :sup:`--2`. Although :math:`\tau_{x}` is only a function of :m
this file must still define a complete 2-D map in order
to be compatible with the standard code for loading forcing fields
in MITgcm. The matlab program :filelink:`verification/tutorial_barotropic_gyre/input/gendata.m`
was used to generate this wind stress file. To run the barotropic jet variation of this tutorial example (see :numref:`baro_jet_solution`),
you will in fact need to run this
matlab program to generate the file ``input/windx_siny.bin``.
was used to generate this wind stress file (alternatively, see python equivalent
:filelink:`gendata.py <verification/tutorial_barotropic_gyre/input/gendata.py>`).
To run the barotropic jet variation of this tutorial example (see :numref:`baro_jet_solution`),
you will in fact need to run one of these
programs to generate the file ``input/windx_siny.bin``.

.. _baro_gyre_build_run:

Expand Down Expand Up @@ -703,21 +706,31 @@ binary data in ``Eta.0000077760.001.001.data`` is as simple as:
::

addpath ../../../utils/matlab/
XC=rdmds('XC'); YC=rdmds('YC');
Eta=rdmds('Eta',77760);
contourf(XC/1000,YC/1000,Eta,[-.04:.01:.04]); colorbar;
colormap((flipud(hot))); set(gca,'XLim',[0 1200]); set(gca,'YLim',[0 1200])
XC=rdmds('XC');
YC=rdmds('YC');
Eta=rdmds('Eta', 77760);
contourf(XC/1000, YC/1000, Eta, [-.04:.01:.04])
colorbar
colormap((flipud(hot)))
set(gca, 'XLim', [0 1200])
set(gca, 'YLim', [0 1200])

or using python (you will need to install the MITgcmutils package, see :numref:`sec_python`):
or using python (you will need to install the MITgcmutils package, see :numref:`MITgcmutils`):

::

from MITgcmutils import mds
import numpy as np
import matplotlib.pyplot as plt
XC = mds.rdmds('XC'); YC = mds.rdmds('YC')
XC = mds.rdmds('XC')
YC = mds.rdmds('YC')
Eta = mds.rdmds('Eta', 77760)
plt.contourf(XC, YC, Eta, np.linspace(-0.02, 0.05,8), cmap='hot_r')
plt.colorbar(); plt.show()
plt.contourf(XC/1000, YC/1000, Eta, np.linspace(-0.02, 0.05,8), cmap='hot_r')
plt.colorbar()
plt.show()

(for a more involved example with detailed explanations how to read in model output,
perform calculations using these data, and make plots, see tutorial :ref:`Baroclinic Ocean Gyre <baroclinic_gyre_solution>`)

Let’s simplify the example by considering the linear problem where we neglect the advection of momentum terms.
In other words, replace :math:`\frac{Du}{Dt}` and :math:`\frac{Dv}{Dt}` with
Expand Down
159 changes: 87 additions & 72 deletions doc/getting_started/getting_started.rst
Expand Up @@ -173,7 +173,7 @@ about git (continue reading...)

3. **Fully embracing the power of git!**

Git offers many tools to help organize and track changes in your work.
Git offers many tools to help organize and track changes in your work.
For example, one might keep separate projects on different branches,
and update the code separately (using ``git pull``) on these separate branches.
You can even make changes to code in the MIT repo tree; when git then
Expand Down Expand Up @@ -366,7 +366,7 @@ This section describes details and capabilities of
:filelink:`tools` directory. :filelink:`genmake2 <tools/genmake2>` is a shell
script written to work in
`bash <https://en.wikipedia.org/wiki/Bash_(Unix_shell)>`_ (and with all
“sh”–compatible shells including
“sh”–compatible shells including
`Bourne <https://en.wikipedia.org/wiki/Bourne_shell>`_ shells).
Like many unix tools, there is a help option that is invoked
thru ``genmake2 -h``. :filelink:`genmake2 <tools/genmake2>` parses
Expand Down Expand Up @@ -771,7 +771,7 @@ optfiles must be written. To create a new optfile, it is generally
best to start with one of the defaults and modify it to suit your needs.
Like
:filelink:`genmake2 <tools/genmake2>`, the optfiles are all written in `bash <https://en.wikipedia.org/wiki/Bash_(Unix_shell)>`_
(or using a simple
(or using a simple
`sh–compatible <https://en.wikipedia.org/wiki/Bourne_shell>`_ syntax).
While nearly all
`environment variables <https://en.wikipedia.org/wiki/Environment_variable>`_
Expand Down Expand Up @@ -1200,7 +1200,7 @@ which is made of the following files:
pressure i.e., downward).

- ``T.00000nIter`` - potential temperature (ocean:
:math:`^{\circ}\mathrm{C}`, atmosphere: :math:`^{\circ}\mathrm{K}`).
:math:`^{\circ}\mathrm{C}`, atmosphere: :math:`\mathrm{K}`).

- ``S.00000nIter`` - ocean: salinity (g/kg), atmosphere: water vapor
(g/kg).
Expand Down Expand Up @@ -1243,7 +1243,7 @@ NetCDF output files
append `netCDF <http://www.unidata.ucar.edu/software/netcdf>`_ files.
Unlike the :filelink:`pkg/mdsio` output, the :filelink:`pkg/mnc`–generated
output is usually placed within a subdirectory with a name such as
``mnc_output_`` (by default,
``mnc_output_`` (by default,
`netCDF <http://www.unidata.ucar.edu/software/netcdf>`_ tries to append,
rather than overwrite, existing files,
so a unique output directory is helpful for each separate run).
Expand Down Expand Up @@ -3032,102 +3032,117 @@ exiting subroutines, etc.)
MITgcm Input Data File Format
=============================

MITgcm input files for grid-related data (e.g., :varlink:`delXFile`), forcing
fields (e.g., :varlink:`tauThetaClimRelax`), parameter fields
(e.g., :varlink:`viscAhZfile`), etc. are assumed to be in "flat" or "unblocked"
`binary format <https://en.wikipedia.org/wiki/Binary_file>`_ . For historical
reasons, MITgcm files use big-endian
`byte ordering <https://en.wikipedia.org/wiki/Endianness>`_, **NOT**
little-endian which is the more common default for today's computers. Thus,
some care is required to create MITgcm-readable input files.
MITgcm input files for grid-related data (e.g., :varlink:`delXFile`),
forcing fields (e.g., :varlink:`tauThetaClimRelax`),
parameter fields (e.g., :varlink:`viscAhZfile`), etc. are assumed to
be in "flat" or "unblocked" `binary format <https://en.wikipedia.org/wiki/Binary_file>`_.
For historical reasons, MITgcm files use big-endian
`byte ordering <https://en.wikipedia.org/wiki/Endianness>`_,
**NOT** little-endian which is the more common default for today's computers.
Thus, some care is required to create MITgcm-readable input files.


- Using `MATLAB <https://www.mathworks.com/products/matlab.html>`_:
When writing binary files, MATLAB's
`fopen <https://www.mathworks.com/help/matlab/ref/fopen.html>`_ command
includes a MACHINEFORMAT option \\'b\\' which instructs MATLAB to read or
write using big-endian byte ordering. 2-D arrays should be index-ordered in
MATLAB as (:math:`x`, :math:`y`) and 3-D arrays as (:math:`x`, :math:`y`,
:math:`z`); data is ordered from low to high in each index, with :math:`x`
varying most rapidly.
- Using `MATLAB <https://www.mathworks.com/products/matlab.html>`_:
When writing binary files, MATLAB's `fopen <https://www.mathworks.com/help/matlab/ref/fopen.html>`_
command includes a MACHINEFORMAT option ``'b'`` which instructs MATLAB
to read or write using big-endian byte ordering. 2-D arrays should be
index-ordered in MATLAB as (:math:`x`, :math:`y`) and 3-D arrays as
(:math:`x`, :math:`y`, :math:`z`); data is ordered from low to high in
each index, with :math:`x` varying most rapidly.

An example to create a bathymetry file (from tutorial :ref:`sec_eg_baro`,
a simple enclosed, flat-bottom domain) is as follows:
An example to create a bathymetry file of single-precision, floating
point values (from tutorial :ref:`sec_eg_baro`, a simple enclosed,
flat-bottom domain) is as follows:

::

ieee='b'; % big endian format
accuracy='real*4'; % this is single precision
ieee = 'b'; % big-endian format
accuracy = 'float32'; % this is single-precision (='real*4')

Ho=5000; % ocean depth in meters
nx=62; % number of gridpoints in x-direction
ny=62; % number of gridpoints in y-direction
nx=62; % number of gridpoints in x-direction
ny=62; % number of gridpoints in y-direction

% Flat bottom at z=-Ho
h=-Ho*ones(nx,ny);
% Flat bottom at z = -Ho
h = -Ho * ones(nx, ny);

% Walls (surrounding domain) - generate bathymetry file
h([1 end],:)=0;
h(:,[1 end])=0;
fid=fopen('bathy.bin','w',ieee); fwrite(fid,h,accuracy); fclose(fid);
% Walls (surrounding domain)
h([1 end], :) = 0; % set ocean depth to zero at east and west walls
h(:, [1 end]) = 0; % set ocean depth to zero at south and north walls

- Using `Python <https://www.python.org/>`_:
Any Python script used to generate MITgcm input files must manually swap the
byte ordering before writing. This can be accomplished with the command:
% save as single-precision (float32) with big-endian byte ordering
fid = fopen('bathy.bin', 'w', ieee);
fwrite(fid, h, accuracy);
fclose(fid);

To read this bathymetry file back into MATLAB, reshaped back to (nx, ny):

::

if sys.byteorder == 'little': data.byteswap(True)
fid = fopen('bathy.bin', 'r', ieee);
h = reshape(fread(fid, Inf, accuracy), nx, ny);
fclose(fid);

- Using `Python <https://www.python.org/>`_:

or, convert as follows while writing an array to a file:
A python version of the above script to create a bathymetry file is as follows:

::

data.astype('>f4').tofile('data.bin')
import numpy as np

Note that 2-D and 3-D arrays should be index-ordered as
(:math:`y`, :math:`x`) and (:math:`z`, :math:`y`, :math:`x`),
respectively, to be written in proper ordering for MITgcm.
Ho = 5000 # ocean depth in meters
nx = 62 # number of gridpoints in x-direction
ny = 62 # number of gridpoints in y-direction

The above MATLAB example translated to Python is as follows:
# Flat bottom at z = -Ho
h = -Ho * np.ones((ny, nx))

# Walls (surrounding domain)
h[:, [0,-1]] = 0 # set ocean depth to zero at east and west walls
h[[0,-1], :] = 0 # set ocean depth to zero at south and north walls

# save as single-precision (NumPy type float32) with big-endian byte ordering
h.astype('>f4').tofile('bathy.bin')

The dtype specification ``'>f4'`` above instructs Python to write the file with
big-endian byte ordering (specifically, due to the '>') as single-precision real
numbers (due to the 'f4' which is NumPy ``float32`` or equivalently,
Fortran ``real*4`` format).

To read this bathymetry file back into Python, reshaped back to (ny, nx):

::

import numpy as np
import sys
Ho=5000; # ocean depth in meters
nx=62; # number of gridpoints in x-direction
ny=62; # number of gridpoints in y-direction
h = np.fromfile('bathy.bin', '>f4').reshape(ny, nx)

# Flat bottom at z=-Ho
h=-Ho*np.ones((ny,nx));
where again the dtype spec instructs Python to read a big-endian
file of single-precision, floating point values.

# Walls (surrounding domain) - generate bathymetry file
h[:,(0,-1)]=0;
h[(0,-1),:]=0;
# save as single precision with big-endian byte-ordering
h.astype('>f4').tofile('bathy.bin')
Note that 2-D and 3-D arrays should be index-ordered as
(:math:`y`, :math:`x`) and (:math:`z`, :math:`y`, :math:`x`),
respectively, to be written in proper ordering for MITgcm.

A more complicated example of using Python to generate input date is provided
in :filelink:`verification/seaice_itd/input/gendata.py`.
A more complicated example of using Python to generate input date is provided in
:filelink:`verification/tutorial_baroclinic_gyre/input/gendata.py`.

- Using `Fortran <https://en.wikipedia.org/wiki/Fortran>`_:
To create flat binary files in Fortran, open with syntax
``OPEN(..., ACCESS='DIRECT', ...)`` (i.e., **NOT** ``ACCESS='SEQUENTIAL'``
which includes additional metadata). By default Fortran will use the local
computer system's native byte ordering for reading and writing binary files,
To create flat binary files in Fortran, open with
syntax ``OPEN(..., ACCESS='DIRECT', ...)`` (i.e., **NOT** ``ACCESS='SEQUENTIAL'``
which includes additional metadata). By default Fortran will use the
local computer system's native byte ordering for reading and writing binary files,
which for most systems will be little-endian. One therefore has two options:
after creating a binary file in Fortran, use MATLAB or Python (or some other
utility) to read in and swap the bytes in the process of writing a new file;
or, determine if your local Fortran has a compiler flag to control
byte-ordering of binary files. Similar to MATLAB, 2-D and 3-D arrays in
Fortran should be index-ordered as (:math:`x`, :math:`y`) and
(:math:`x`, :math:`y`, :math:`z`), respectively.
after creating a binary file in Fortran, use MATLAB or Python (or some
other utility) to read in and swap the bytes in the process of writing a new file;
or, determine if your local Fortran has
a compiler flag to control byte-ordering of binary files.
Similar to MATLAB, 2-D and 3-D arrays in Fortran should be index-ordered
as (:math:`x`, :math:`y`) and (:math:`x`, :math:`y`, :math:`z`), respectively.

Using `NetCDF <http://www.unidata.ucar.edu/software/netcdf>`_ format for input
files is only partially implemented at present in MITgcm, and use is thus
discouraged.
Using `NetCDF <http://www.unidata.ucar.edu/software/netcdf>`_ format for input files is only
partially implemented at present in MITgcm, and use is thus discouraged.

Input files are by default single-precision real numbers (32-bit, ``real*4``),
but can be switched to double precision by setting namelist parameter
:varlink:`readBinaryPrec` (``PARM01`` in file ``data``) to a value of 64.
but can be switched to double precision by setting
namelist parameter :varlink:`readBinaryPrec` (``PARM01`` in file ``data``)
to a value of 64.

0 comments on commit ce0d9af

Please sign in to comment.