SPINS User Guide

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 the class BaseCase (see BaseCase.hpp).

Creating your own custom configuration involves building a derived class based on BaseCase. The case file cases/doc_minimal.cpp shows the structure of a case file. It makes sense to start with another similar case file and customise it.

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
}

but 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 and generate the grid, then 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 on 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!