advect_1d_FD.f90
Solves the 1D linear advection equation
in a finite domain subject to periodic boundary conditions with first order finite differencing schemes.
User can choose to use either an upwind method, or FTCS (forward time centred space). One of those choices is a bad idea!
The main use of this code is to
- Demonstrate that FTCS is unconditionally unstable regardless of CFL condition (i.e. choice on timestep)
- Demonstrate that upwinding is stable for C <= 1
- Demonstrate that C=1 for upwinding produces a perfect translation of the initial condition
- Demonstrate that C<1 introduces numerical diffusion of the profile during the evolution
The code is more "complicated" than it needs to be for such a simple task, which could other wise be coded as a rather simple script. It is broken into various modules and subroutines in anticipation of the requirements of a more sophisticated hydrodynamics code, where such a structure would be useful.
All modules and subroutines are self-contained within the single text file since it is so short.
It is intended the user change user_control and (if you like) initial_conditions
It is fairly self explanatory:
- INT nx : number of computational cells to use
- REAL x_min and x_max : set the size of the 1D domain
- REAL u : set characteristic speed
- REAL cfl : set a courant number. For upwinding, 1 is perfect, <1 shows dissipation, >1 is unstable
- REAL t_end: time to end the simulation.
- INT nsteps : if >=0 then the code will output at this step. Set to -1 to run to t_end instead
- LOGICAL ftcs and upwind : set one to true and the other false to choose the scheme
- LOGICAL verbose : set to true to print some basic diagnostic stuff during evolution
Since its all in one file something like this should do:
gfortran advect_1d_FD.f90 -o advect_1d_FD
./advect_1d_FD
However, there is also a simple makefile.
The code saves the coordinates xb and the function a(x) as .dat files at the step it is terminating (either according to nsteps or t_end) to the current working directory. It also calls a simple python script to automatically render a simple plot of xb vs a(x). This is called "output.png"
You should be able to quickly clean this stuff up with "make clea" and "make cleanall"
This particular code's history is as an independent Fortran implementation of the algorithms discussed in Chapter 4 of this excellent book http://bender.astro.sunysb.edu/hydro_by_example/CompHydroTutorial.pdf (he also provides python scripts to do the same sort of solves if you are interested to see the different approaches).