SPINS User Guide: Difference between revisions
No edit summary |
m (moved Creating a SPINS case to SPINS User Guide: Better title) |
(No difference)
|
Revision as of 11:25, 12 July 2012
Welcome to the SPINS case setup page.
Getting SPINS and compiling
(Chris: Is it a good idea to put it on GitHub?)
SPINS consists of a bunch of C++ source files and a bunch of case files, and it requires four libraries. UMFPack, AMD and Blitz++ are supplied with SPINS, and it uses the system-installed FFTW.
Directory structure:
- cpp/src - SPINS source files
- cpp/src/cases - A few dozen example case files
- cpp/AMD - AMD library
- cpp/UMFPack - UMFPack library
- cpp/UFconfig - Helper for compiling AMD/UMFPack
To compile, ToDo: write this! .
The basics
The SPINS model is a Navier-Stokes solver that gets parameters and initial/boundary conditions from calls to user-provided routines. The user-provided routines are encapsulated in class derived from BaseCase (see BaseCase.hpp).
Creating your own custom configuration involves supplying the user-provided routines in a derived class based on BaseCase. The case file cases/doc_minimal.cpp shows the structure of a case file. It usually makes sense to start with a similar case file and customise it.
Calling sequence
Chris: can you verify this?
To help understand where/when to do certain operations, it's good to know the order that the user routines are called. The order is as follows:
- The constructor
- The type_* routines
- The size_* and and length_* routines
- The tracer number routines (numActive, numPassive, numtracers)
- is_mapped and do_mapping
- init_time
- init_vels
- init_tracer (single tracer) or init_tracers (multiple tracers)
At each time step, the following routines are called
- forcing:
- If no active tracers, only passive_forcing is called
- If active tracers, first vel_forcing, then tracer_forcing
- analysis
Examples of common operations
You can find examples of how to do various operations by digging through the case files. Some of the common operations are reproduced here.
Generating the grid vectors
At the top of your file include the global tensor indices:
// Tensor variables for indexing
blitz::firstIndex ii;
blitz::secondIndex jj;
blitz::thirdIndex kk;
Define arrays for the grid vectors:
// Arrays to hold the x, y, and z points
Array<double,1> xx, yy, zz;
You can initialise them in the constructor. For Periodic/free-slip the grid spacing is even:
// The constructor
casename() : xx(split_range(Nx)), yy(Ny), zz(Nz) {
// Assumes periodic or free-slip (cosine)
xx = Lx*(ii+0.5)/Nx; // x-coordinate: periodic
yy = Ly*(ii+0.5)/Ny; // y-coordinate: periodic
zz = Lz*(ii+0.5)/Nz; // z-coordinate: free-slip
}
For no-slip (Chebyshev) the grid spacing is not even:
// The constructor
casename() : xx(split_range(Nx)), yy(Ny), zz(Nz) {
// Assumes no-clip (Chebychev) in all dimensions
xx = -Lx*cos(ii*M_PI/(Nx-1))/2; // x-coordinate: Chebyshev
yy = -Ly*cos(ii*M_PI/(Ny-1))/2; // y-coordinate: Chebyshev
zz = -Lz*cos(ii*M_PI/(Nz-1))/2; // z-coordinate: Chebyshev
}
Outputting the grid
Create array tmp, then generate each grid and write the array and the grid reader out:
DTArray *tmp = alloc_array(Nz,Ny,Nz);
tmp = xx(ii) + 0*jj + 0*kk;
write_array(tmp,"xgrid");
write_reader(tmp,"xgrid",false);
tmp = 0*ii + yy(jj) + 0*kk;
write_array(tmp,"ygrid");
write_reader(tmp,"ygrid",false);
tmp = 0*ii + 0*jj + zz(kk);
write_array(tmp,"zgrid");
write_reader(tmp,"zgrid",false);
Initialising velocities
To Do: write this!
Initialising tracers
To Do: write this!
Forcing
To Do: write this!
Mapped grid
To Do: write this!
Analysis: energy diagnostic
If you're using a periodic grid, use this for kinetic energy diagnostic
double ke = 0.5*rho_0*pssum(sum(u*u+v*v+w*w))/(Nx*Ny*Nz)*Lx*Ly*Lz; // KE
if (master()) {
FILE * en_record = fopen("energy_record.txt","a");
assert(en_record);
fprintf(en_record,"%.8f %.14g\n",time,ke);
fclose(en_record);
}
and if you're on a mapped chebyshev grid, you can use this for the KE computation
double ke = pssum(sum((*get_quad_x())(ii)*(*get_quad_y())(jj)*(*get_quad_z())(kk)*
(pow(u(ii,jj,kk),2)+pow(v(ii,jj,kk),2)+pow(w(ii,jj,kk),2))));
Analysis: saving snapshots
To Do: write this!
Analysis: compute vorticity
To Do: write this!