NavierStokes3DPreStep.c File Reference

Pre-step function for 3D Navier Stokes equations. More...

#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <matmult_native.h>
#include <physicalmodels/navierstokes3d.h>
#include <hypar.h>

Go to the source code of this file.

Functions
int	NavierStokes3DPreStep (double *u, void *s, void *m, double waqt)

Detailed Description

Pre-step function for 3D Navier Stokes equations.

Author

Debojyoti Ghosh

Definition in file NavierStokes3DPreStep.c.

Function Documentation

int NavierStokes3DPreStep	(	double *	u,
		void *	s,
		void *	m,
		double	waqt
	)

Pre-step function for the 3D Navier Stokes equations: This function is called at the beginning of each time step.

• The solution at the beginning of the step is copied into NavierStokes3D::solution for the linearized flux partitioning.
• The linearized "fast" Jacobian (representing the acoustic modes) NavierStokes3D::fast_jac is computed as:

  \begin{equation} A_f\left({\bf u}\right) = R\left({\bf u}\right)\Lambda_f\left({\bf u}\right)L\left({\bf u}\right) \end{equation}

  where \(R\) and \(L\) are the matrices of right and left eigenvectors, and,

  \begin{equation} \Lambda_f = diag\left[0,0,0,u+c,u-c \right] \end{equation}

  with \(c\) being the speed of sound.
  
  Reference:
• Ghosh, D., Constantinescu, E. M., Semi-Implicit Time Integration of Atmospheric Flows with Characteristic-Based Flux Partitioning, SIAM Journal on Scientific Computing, 38 (3), 2016, A1848-A1875, http://dx.doi.org/10.1137/15M1044369.

Parameters

u	Solution vector
s	Solver object of type HyPar
m	MPI object of type MPIVariables
waqt	Current simulation time

Definition at line 32 of file NavierStokes3DPreStep.c.

{
  HyPar             *solver = (HyPar*)   s;
  NavierStokes3D    *param  = (NavierStokes3D*) solver->physics;
  int               *dim    = solver->dim_local;
  int               ghosts  = solver->ghosts, dir, p;
  double            *A;

  static const int  ndims   = _MODEL_NDIMS_;
  static const int  JacSize = _MODEL_NVARS_*_MODEL_NVARS_;
  static int        index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
  static double     D[_MODEL_NVARS_*_MODEL_NVARS_],L[_MODEL_NVARS_*_MODEL_NVARS_],
                    R[_MODEL_NVARS_*_MODEL_NVARS_],DL[_MODEL_NVARS_*_MODEL_NVARS_];

  /* set bounds for array index to include ghost points */
  _ArrayAddCopy1D_(dim,(2*ghosts),bounds,ndims);
  /* set offset such that index is compatible with ghost point arrangement */
  _ArraySetValue_(offset,ndims,-ghosts);
  /* copy the solution to act as a reference for linearization */
  _ArrayCopy1D_(u,param->solution,(solver->npoints_local_wghosts*_MODEL_NVARS_));

  int done = 0; _ArraySetValue_(index,ndims,0);
  while (!done) {
    _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
    int q = _MODEL_NVARS_*p;

    dir = _XDIR_;
    A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
    /* get the eigenvalues, and left & right eigenvectors */
    _NavierStokes3DEigenvalues_      ((u+q),_NavierStokes3D_stride_,D,param->gamma,dir);
    _NavierStokes3DLeftEigenvectors_ ((u+q),_NavierStokes3D_stride_,L,param->gamma,dir);
    _NavierStokes3DRightEigenvectors_((u+q),_NavierStokes3D_stride_,R,param->gamma,dir);
    /* remove the entropy modes (corresponding to eigenvalues u) */
    D[0] = D[12] = D[18] = 0.0;
    /* assemble the Jacobian */
    MatMult5(_MODEL_NVARS_,DL,D,L );
    MatMult5(_MODEL_NVARS_,A ,R,DL);

    dir = _YDIR_;
    A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
    /* get the eigenvalues, and left & right eigenvectors */
    _NavierStokes3DEigenvalues_      ((u+q),_NavierStokes3D_stride_,D,param->gamma,dir);
    _NavierStokes3DLeftEigenvectors_ ((u+q),_NavierStokes3D_stride_,L,param->gamma,dir);
    _NavierStokes3DRightEigenvectors_((u+q),_NavierStokes3D_stride_,R,param->gamma,dir);
    /* remove the entropy modes (corresponding to eigenvalues v) */
    D[0] = D[6] = D[18] = 0.0;
    /* assemble the Jacobian */
    MatMult5(_MODEL_NVARS_,DL,D,L );
    MatMult5(_MODEL_NVARS_,A ,R,DL);

    dir = _ZDIR_;
    A = (param->fast_jac + _MODEL_NDIMS_*JacSize*p + dir*JacSize);
    /* get the eigenvalues, and left & right eigenvectors */
    _NavierStokes3DEigenvalues_      ((u+q),_NavierStokes3D_stride_,D,param->gamma,dir);
    _NavierStokes3DLeftEigenvectors_ ((u+q),_NavierStokes3D_stride_,L,param->gamma,dir);
    _NavierStokes3DRightEigenvectors_((u+q),_NavierStokes3D_stride_,R,param->gamma,dir);
    /* remove the entropy modes (corresponding to eigenvalues v) */
    D[0] = D[6] = D[12] = 0.0;
    /* assemble the Jacobian */
    MatMult5(_MODEL_NVARS_,DL,D,L );
    MatMult5(_MODEL_NVARS_,A ,R,DL);

    _ArrayIncrementIndex_(ndims,bounds,index,done);
  }
  return(0);
}