THE CRUNCH GROUP
primary phone: +1 (401) 863 3694
secondary phone: +1 (401) 863 1594
web: http://www.cfm.brown.edu/crunch
Division of Applied Mathematics
182 George st, Box F
Brown University, Providence
Rhode Island, RI 02912, USA
home  ::  directions  ::  contact  ::  credits  ::  fluid mechanics  ::  applied math


We will start with the basic equations for the incompressible case. Assuming a divergence-free velocity field we have





The splitting scheme used is a 3-step one. The details of it can be found in the Spectral/hp element methods for CFD, page 247. This is a brief (and maybe incomplete) introduction: For the O(1) lowest accuracy, the following is the rough procedure used: Assume we are at the nth step. We treat the nonlinear term first

For this simple case this is explicit. We solve for u-hat, and we can move on to the second intermediate step. To recover the second step, we have to employ the incompressibility condition:

ie: we have to keep in mind that we want to solve for (3) below, so we take its divergence, and solve for the pressure first. The result is

and we obtain pn+1. Everything is ready for the actual second step, which is coupled with the divergence-free condition used:

Next solve for u-double-hat. After this is done, we move on to the viscous step

which is equivalent to the 3rd step:<

/p>




For higher accuracy we can employ different approaches (and this is what is done in the code): We can use a forward multi-step scheme of the Adams family:

Given this, the RHS of the nonlinear step can be re-written as

Keeping the following notation

the LHS of (4) and (5) are also modified as follows:

This completes the basic idea of the method employed. All of the above will be considered in drive.C as a guideline for the high-level routines that will appear.





It is pretty easy to get familiar with the above structure. The idea is the following:





Quick compilation quide

To start compiling the code, go the specific directory of your architecture (just do a uname) and type

gmake {opt/dbx}

in the prompt. For the Linux machines, the architecture has to be specified and enforced in the commmand line:

gmake ARCH=GCC {opt/dbx}

(no spaces on either side of the '=' sign).

Same for the Hlib. For the removal of the object files (*.o) the well-known useful command

gmake clean

will clean the directory and let you compile all the files again with either the dbx or the opt version (linking between optimized and debug object files is not recommended - works but would be slow).




Nektar2d / src

Among others the source code contains:



There are more files to be included only if you turn a specific flag on while you compile, but these are the basic ones we will consider for the moment.

drive.C

This file contains all the high-level functions that will be executed when the simulation starts. First it sets the domain defined in the header file

Omega = PreProcess(   |__|   ,   |__|  )

where the arguments are just the flags and arguments that you specify in the command line. This routine lives in prepost.C. It alocates memory for most variables, opens files (.chk, .his etc) - in other words does all the dirty I/O. All memories are set now.

The outline of the solver has as follows: 


 StartUp

  while (step < nsteps) 

  {
     set_order
     MakeF(Prep)
     MakeF(WaveProp)
  u: Integrate_SS(...) 
  v: Integrate_SS(...) 
     MakeF(Pressure)
     P_Solve(...)
     MakeF(Viscous)
  u: V_Solve(...)
  v: V_Solve(...)    
     MakeF(Post)
     Analyser
  }


We have to keep in mind that everything in Nektar can be done in either physical or modal (transformed) space. According to the expansion
we denote by 'u' the physical space and by 'û' the transformed space. The default state everywhere is the transformed space. The only time when the physical space is used is during the evaluation of the nonlinear term, due to the nature of this nonlinearity. The way to check in which space we are in is through 


 if(omega->U->fhead->state = 'p')
   omega->U->Trans(omega->U, Q_to_J)


(remember: *a.U is equivalent to a->U).

The function Trans(...) is a member function of Element_List and transforms the fields from Quadrature/Physical space (Q) to Jacobi/Transformed space (J).

The Element_List U is a member of Omega and contains all the elements in the domain Omega. Imagine a list of elements:

|... 11 ...|... 2 ...|... 10 ...|... 7 ...|... 8 ...|... 23 ...|... 19 ...|...


This creates an Element_List, and each cell is an Element. The first element of the Element_List U has a name: fhead. This is an element class and one of its members is 'state'. This can take value 'p' (physical) or 't' (transformed). All we need to do is to check one element in the field to find out its state. Usually this is the first one (fhead). Then the whole field is bound to be at the same state. Trans(...) takes arguments 


 J_to_Q : Jacobi to Quadrature 
 Q_to_J : Quadrature to Jacobi


Since the fields are by default in Transformed space, the field can be in physical space if, say, we assign a given initial condition at the beginning of the run (GIVEN parameter). Restart files are in Transformed space.




The splitting scheme