README for the Streamfunction Vorticity Navier Stokes Solver These are the source files and examples for the streamfunction vorticity solution of the Navier Stokes equations. The directories are: -------------------- Src source files No Sflags.h file, all options are set in dimns.f90. To make the executable: make ns2dp Data read from file ns2d.dat Problem dependent file for initialization, examples are tgrvadj.f90 and theat.f90 Modify neumann.f90 to set the no-slip vorticity boundaries. Input contains some grid files and their partitions GrvAdj Gravitational adjustment problem OverFlow an overflow-like problem User-controlled flags in dimns.f90 ---------------------------------- "npts" refers to the number of vorticity points per element. There are two options for their distribution within an element, and is controlled by the logical flag "dgm_vort". If "dgm_vort=true" a Discontinuous Galerkin Method (DGM) is used for the vorticity, and the vorticity is collocated at the (interior) Gauss points. If "dgm_vort=false" the vorticity is discritized with the normal spectral element procedure: a Continuous Galerkin method (CGM), and the collocation points are the Gauss-Lobatto points. I found the latter to be more robust then the former. "nptp" is not really used anymore, I left it to avoid having to weed it out through the entire code, it is set equal to npts. "nptT" is the number of tracer (density) points within an element. These are collocated at the Gauss points. There are a bunch of logical flags "PARALLEL" If true the code runs in parallel "PERIOD_FLAG" Must be set to true if boundary conditions are periodic in x-direction. The wavelength of the periodicity must be also set in file period.f90 by modifying the parameter "wave_length". "SPHERE" Set to true if running on the spherical Earth. "ROTATE_POLE" Set to true if pole is within grid "TIMING_FLAG" Set to true if you want to perform timing on the code The timing module would then report the following statics 1. total wall clock time in the time-stepping routine 2. wall clock time per time step 3. various timing in the message passing portion of the code "LINEAR" If true the vorticity advection term (only) is omitted "dg_flag" Set to true if there is a tracer discretized with DGM "tracers_flag" Set to true if there is a tracer to be updated dg_flag and "tracers_flag" must be identical. "dgm_vort" If true vorticity is discretized with DGM, otherwise with CGM. see the explanation for "npts" for more information. No slip boundary conditions: ---------------------------- The assignment of the no-slip boundary condition is somewhat more involved. The input grid file holds Dirichlet boundary conditions on the vorticity and streamfunction. Dirichlet conditions on the vorticity are appropriate at inlets, free slip boundaries and symmetry planes, while streamfunction Dirichlet conditions are appropriate on no-normal flow boundaries. At no-slip walls the boundary vorticity is diagnosed from the definition; its calculation involves local boundary integrals and hence require extra information for its processing. This is all done in file "neumann.f90" in subroutine "Noslip_Vort_Boundaries". Actually, the file contains several example subroutines for the driven cavity, backward facing step, gravitational adjustment and overflow problems. These subroutines are called "Noslip_Vort_Boundaries_LABEL" where the label refers to the problem. The only relevant routine is the one without the _LABEL suffix. This subroutine checks all the boundary edges, and test their locations with respect to the no-slip boundaries. If they do coincide with one, the element and edge are added to the list of no-slip edges. The normal and tangential velocity components along the edge are also set; both components are zero for stationary boundaries. Input Files: ------------ The code reads the file "ns2d.dat" to set all its run-time parameters. There are sample files in the appropriate directories. The run-time parameters are: 1. The partition file necessary for the grid decomposition if "PARALLEL=true". The file must be omitted if the code is run in serial mode. 2. The vorticity file 3. the time step and the number of timesteps 4. the snapshot interval 5. flag to initialize variable if set to 1, variables are read from restart file 6. the initial time 7. the gravitational constant 8. the viscosity parameter 9. The diffusion parameter Output Files: ------------- fort.11: the vorticity at the specified snapshot interval fort.12: the streamfunction fort.13: the tracer field fort.21: the x-velocity component at the final time step fort.22: the y-velocity component at the final time step All field are output in elemental form, so there are "npts x npts x nelem" items per snapshot for the vorticity and streamfunction, and "nptT x nptT x nelem" items for the tracer. Post-processing: ---------------- The post-processing of the files for plotting purposes require slight modifications to the plotting and interpolation routines. Interpolation: The interpolation routine interp2d now requires an extra line to tell the code where the variables are collocated within an element, Gauss or Gauss-Lobatto points. if the "collocation points" flag below is set to 2, the points are the Gauss points, otherwise, they are the Gauss-Lobatto points. ================new interp2d.dat file=================================== F ! sphere logical T if on a sphere T ! elemental logical T if elemental /cluster/mohamed/Channel/channel_040x10.06 structured_grid 0 64.e3 201 ! xmin, xmax, nx -20 0 51 ! ymin, ymax, ny 2 5 ! collocation points (2 for Gauss), nptT seom_inputfile.dat structuredgrid_outputfile.dat 0 ================new interp2d.dat file=================================== Visualization: Drawing the streamfunction and vorticity fields directly from the seom data file is rather easy, any of the plotting scripts from the shallow water code will do. To plot the tracer field defined on the Gauss grid you need to change the iformat flag from -4 (for unformatted binary data on Gauss-Lobatto points in elemental form) to -11 (for unformatted binary data on Gauss points in elemental form, where the number of points is npts-1, npts is read from the grid file). Todo: 1. The main capability lacking in the code is the ability to use different viscosity and diffusion coefficients in the vertical and horizontal. Given the large discrepancies between these two scales it is an important capability. 2. Multiply connected domain cannot be handled. The streamfunction on internal "holes" must be determined as part of the solution, and is not currently implemented.