flowSolver.h

Go to the documentation of this file.

#ifndef FLOWSOLVER_H
#define FLOWSOLVER_H

// i and j are consistent with the physical domain
// i.e., same i means same x coordinate
//       same j means same y coordinate
// discretization stencil is:
//              (i,j+1)
//     (i-1,j)  (i, j )   (i+1,j)
//              (i,j-1)

//
// fields of primitive varibles in the original two-dimensional domain.
// thus, use matrix form.
//

#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include <cmath>
#include <vector>
#include <cstdlib> // for rand() and srand()
#include <ctime>   // for time()
#include "inputData.h"
#include "initCond.h"
#include "exactSol.h"
#include "outputData.h"
#include "initStage.h"
#include "NSSolver.h"
#include "setDt.h"
#include "constant.h"
using namespace std;

#ifdef SINGLE
#define real float
#else
#define real double
#endif

class flowSolver {

 public:

  flowSolver();

  ~flowSolver();


  void solve();

 private: 

  void setCoord();
  void setUnknowns();

  int * allocate1dIntArray(int size);
  real * allocate1dRealArray(int size);
  int ** allocate2dArray(int nMeshes, int dimension);
  real ** allocate2dRealArray(int ny, int nz);
  real *** allocate3dRealArray(int nMeshes, int dim, int **pointNumberMesh);
  real **** allocate4dRealArray(int nMeshes, int **pointNumberMesh);
  int **** allocate4dIntArray(int nMeshes, int **pointNumberMesh);
  real **** allocate4dRealArrayCage(int nMeshes, int **pointNumberMesh);
  int **** allocate4dIntArrayCage(int nMeshes, int **pointNumberMesh);

  void delete1dArray(int *x);
  void delete1dArray(real *x);
  void delete2dArray(int **x, int size);
  void delete2dRealArray(real **x, int ny);
  void delete3dArray(real ***x, int nMeshes, int dim);
  void delete4dArray(real ****x, int nMeshes, int **pointNumberMesh);
  void delete4dArray(int ****x, int nMeshes, int **pointNumberMesh);
  void delete4dArrayCage(real ****x, int nMeshes, int **pointNumberMesh);
  void delete4dArrayCage(int ****x, int nMeshes, int **pointNumberMesh);
  void releaseMemory();
  string iToString(int k);
  //data members


  int dim; // dimension of the problem: 2d or 3d
  real t;  //current time
  real dt; // time increment
  real Re; // Reynolds number based on the radius of a particle

  int Nc; // the order of spectral expansion to truncate

  int nPatch; // each patch carries a square fine mesh
  int NumberOfMeshes;   // including the coarse mesh
  // the coarse mesh is the zeroth mesh and fine meshes are from 1st to the last

  // cell number of coarse or fine meshes
  int **nCell; // cellNumberCoarse[NumberOfMeshes][dim]

  int nCoarseCellCompDomain;
  // number of coarse cells
  // the computational domain use (>=12):


// coarse cell number in the coarse mesh that a patch uses
  int cellNumberPerPatch; 

  int nCellPerParticle; // fine cell number used by a particle

  int **StartIndex;
  // in matrix form: StartIndex[NumberOfMeshes][dimension];
  // note: the first element with index 0
  // with respect to pressure point


  // i for index in 1d array or x direction in 2d cases
  // j for y direction in 2d cases
  // k for z direction in 2d cases
  int i, j, k;

  int tStep;

  // m for mesh iteration: meshes for patchs and coarse grid
  int m;


  int **pointNumberMeshP; // for pressure in x and y directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  int **pointNumberMeshU; // for velocity u  in x and y directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  int **pointNumberMeshV; // for velocity v  in x and y directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  int **pointNumberMeshW; // for velocity w  in x, y, and z directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  int **pointNumberMeshOmegaX; // for vorticity omega  in x and y directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  int **pointNumberMeshOmegaY; // for vorticity omega  in x and y directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  int **pointNumberMeshOmegaZ; // for vorticity omega  in x and y directions
  // in the form: pointNumberMesh[NumberOfMeshes][dimension]; 

  // nx and ny for wavenumber in x and y directions
  int nx, ny, nz;


  // in the form: h[NumberOfMeshes]; 
  real hCoarse, hFine; 
  real *h; 


  // pressure
  real ***xP; //coordinates for pressure
  real ****fp; //for for pressure poisson equation
  real ****exactP; // exact prrssure
  real ****p; // solution corresponding to pUnknowns
  int  ****pUnknowns;

  // velocity u
  real ***xU; //coordinates for u
  real ****exactU; // exact prrssure
  real ****u; // velocity u
  real ****uAux; // auxilliary velocity
  int  ****uUnknowns;

  //velocity v
  real ***xV; //coordinates for v
  real ****exactV; // exact prrssure
  real ****v; // velocity v
  real ****vAux; // velocity v
  int  ****vUnknowns;

  //velocity w
  real ***xW; //coordinates for v
  real ****exactW; // exact prrssure
  real ****w; // velocity v
  real ****wAux; // velocity v
  int  ****wUnknowns;

  //particles
  real *xParticleCenter;  // coordinates of particle center
  real *yParticleCenter; 
  real *zParticleCenter; 

  // only assign for particles
  int  ****pInsideParticle; // is p is inside the particle surface
  int  ****pCagePoint; 
  int  ****uCagePoint;  
  int  ****vCagePoint;
  int  ****wCagePoint;
  int  ****omegaXCagePoint;
  int  ****omegaYCagePoint;
  int  ****omegaZCagePoint;

  // vorticity omega
  real ***xOmegaX; //coordinates for vorticity
  real ****omegaX; // vorticity

  real ***xOmegaY; //coordinates for vorticity
  real ****omegaY; // vorticity

  real ***xOmegaZ; //coordinates for vorticity
  real ****omegaZ; // vorticity

  // if fixed particles, calculate U, S, V for one partilce
  // in the initial step and copy here
  // in late steps, don't need to calculate the expensive part of U, S, V 
  // for efficiency
  int fixedParticles;
  vector<real> UFixedParticles;
  vector<real> SFixedParticles;
  vector<real> VFixedParticles;

  // for rotating particles
  int rotatingParticles;
  real **angVel; //angular velocity
  real **angVelAcc;  //angular velocity acceleration
  real **torqueOld; // torque
  real **torqueNew; // torque


  string filename;
  ofstream ofParticle ;
  int restart;

};

#endif

---