MSC321 Programming Projects
- Particle Tracking
Particle tracking using downloadable winds velocities from atmospheric prediction web-sites.
This can be useful to simulate pollution tracking.
The project includes potentially developing routines to:
- Phase 1 of the project can be downloaded.
- Read wind data from weather prediction websites.
Two files and some
libraries are now available to read NCEP forecast winds.
The code to read NCEP data and do the interpolation are now bundled
in a directory on metofis:~mohamed/PartSphere. The two source
files that you need are ncio.f90 and gfsread.f90.
The NCEP data is residing
on NCEP servers and you will read it over the web; this is what they call
cloud computing. However, this entails using quite a few libraries.
Fortunately these are available on metofis, but prior to using them you
need to add the following lines to your .bashrc file:
LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/netcdf4.1.2/lib
export LD_LIBRARY_PATH
- Interpolate the wind data in space and time.
Module lagint.f90
contains fortran code to do polynomial Lagrange interpolation. The velocity
at the particle location may not be available as an explicit expression. Instead
the velocity is known only on a discrete grid. One must then use interpolation
to estimate the velocity at the particle location.
- The example code provided in directory ~mohamed/PartSphere
includes a number of files. Here is a summary:
- ncio.f90: handles netcdf files sitting on NCEP servers
- gfsread.f90: calls ncio.f90 routines to get needed data
- timing.f90: computes timing of select code segments
- fileio.f90: does fortran type I/O in ASCII or binary form
- lagint.f90: does space-time interpolation
- project.f90: is the main program that coordinates the above tasks
- Makefile: automates the construction of the executables.
It is important that you use since the compilation requires extra
options you have not used before. To create the executable just
say:
make project
and an executable called project will be created.
- update the strength of the advected field if needed (radioactive decay say).
- Include stochastic perturbations to the wind field to model the effect of turbulence on
particle trajectories (optional)
- Plotting the evolution of the particles as they treck across the globe
This matlab plotting script
helps you visualize your output data on a sphere. This script uses a
toolbox called M_map to perform visualization on the sphere. The plotting
command are the same as if you are plotting in cartesian coordinates but
with a prefix of "m_". Look at the script for examples. The toolbox is
freely downloadable and a copy of it is available from "~mohamed/Matlab"
on metofis. The addpath command in the script should be modified to
reflect the location of the toolbox on the system. You don't have to do
anything if you running matlab from metofis.
- Parallelizing the particle evolution so as to handle large numbers of them.
The problem is also embarrassingly parallel with each processor tracking a set number of particles.
- Control the release of particle to model continuous release, sudden one-time burst, etc..
- Phase 2:
Your task is now to integrate the NCEP data, and the
interpolation routine into your code so as to track particles in a
real wind field on the sphere. The project consists of tracking
86400 particles released either from Fukichima Japan, or Miami
for an entire week as they carried by the winds. These particles
should be released one every second for the first 24 hours at the
level of 800mb. The non-dimensional pollutant they carry should be
normalized to 1 and decay like a Gaussian blob with a standard
deviation of 1/8 degree. Plot the contours of pollution every 12 hours.
The following write-upshould help with the
last phase of the project.
- Vortex Methods
Vortex methods are a simple extension to particle method. The vorticity
is represented by particles whose position and circulation are
known, these vorticity blobs are advected by the flow and can be
tracked. For 2D problems the vorticity carried by the particle is
conserved and remains constant, the only update needed is for the position
of the particle. In addition, a simple formula need to be programmed to
compute the flow field induced by the vortex blobs. One can then model the
evolution of an entire flow field with particle tracking only. The model
should be easily parallelizable with OpenMP directives. One interesting
application of a vortex model is the simulation of a vortex merger,
whereby a small vortex located close to a larger vortex would eventually
merge. The programming tasks will include some of the items listed below
- Phase 1 of the vortex method project.
(It is the same as for the particle tracking code).
- Initialize a set of vortex blobs
- Calculate the velocity induced by each one of these vortex blobs (this involves a number of geometric information)
- Update the blobs position using a time-stepping scheme like RK4
- Repeat the process for the next time-level.
- The post-processing would include plotting the position of the particle and the evolution of the vorticity/streamfunction fields.
- For large number of particles, a paralell version of the velocity computations can be developed to reduce turn around time.
Once the validation of the code components is done, you can move
on to simulate the axisymmetrization of an elliptical vortex. The
settings of the problem are described in the following
document
- Passive Tracer Advection
The purpose of this project is to write a parallel 2D tracer advection
code using a finite volume method. The tracer is passive (in
that it does not affect the flow field) and can represent temperature,
humidity, or a pollutant concentration. This is the Eulerian counterpart
to the Lagrangian particle tracking scheme. The programming tasks
would include:
- A discretization of the computational domain into rectangular finite volume cells
- An initialization routine to set the initial conditions of the tracer fields and
the wind field.
- A time-stepping procedure, to update the cell-averaged tracer in each cell e.g. RK4.
The time stepping procedure needs to do the following at each time-step:
- A Finite Volume reconstruction procedure to compute the tracers at the edges of the cell
from the cell averages.
- Calculate the advective flux of tracers at the edges of the cells
- Calculate the net fluxes (divergence) entering or leaving the cells
- The message passing paradigm necessary to communicate between different processors.
- The 1D -component of the finite volume project.
- The 2D -component of the finite volume project.
- Directory ~mohamed/FV2D on metofis contains a number of files
for the final projects they are
- fileio.f90 for reading/writing binary fortran files
- testio.f90 for testing and illustrating usage of module fileio.
- flow.f90 for setting the flow field
- initialize.f90 for setting the initial condition on tracer
- project.f90 skeleton (empty) main program for 2D Finite
Volume advection project.
- fv2d.m A matlab script to read the 2D FV grid, the
tracer field and performing contouring for visualization.
- Iterative poisson solver
1D/2D parallel iterative solver for the Poisson equation using finite difference methods.
The code components would be:
- The project outline
- A discretization of the computational domain using finite differences
- Calculations of the right hand side forcing term and initialization of the guess
- Iterate the guess until a convergence criterion is reached.
- The iterations can be quite simple and based on a weighed Jacobi or red-black
Gauss-Seidel iteration. The changes in the guess need to be recorded and in the residual.
As the number of points increases, the number of iterations increases substantially.
One way to counter the cost increase is to solve the problem in parallel.
- Implement a parallel decomposition of the solver and establish the communication
between processors to solve the problem in faster time.
- Boundary Element Method
Pure vortex methods permit the divergent component of the velocity field to be represented. Their simplicity however, comes at a price in that one cannot take into account the effects of boundaries (object)
on the flow. For that, one needs to compute the irrotational part of the flow, and that cannot be done with particles. For 2D flows, a boundary element method, which requires only the description of the boundaries
through line segments can be used. The Boundary Element code is already written, the list of tasks to
perform are:
- Compute the velocity imparted at the boundary by the vortex blobs.
- Solve the associated irrotation problem that would yield the opposite velocity,
the solver code would be provided.
- Add the divergent and irrotational velocity components.