NavierStokes2DJacobian.c

Go to the documentation of this file.

#include <mathfunctions.h>
#include <matmult_native.h>
#include <physicalmodels/navierstokes2d.h>

int NavierStokes2DJacobian(
                    double  *Jac, 
                    double  *u,   
                    void    *p,   
                    int     dir,  
                    int     nvars,
                    int     upw   
                   )
{
  NavierStokes2D *param = (NavierStokes2D*) p;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];

  /* get the eigenvalues and left,right eigenvectors */
  _NavierStokes2DEigenvalues_      (u,D,param->gamma,dir);
  _NavierStokes2DLeftEigenvectors_ (u,L,param->gamma,dir);
  _NavierStokes2DRightEigenvectors_(u,R,param->gamma,dir);

  int aupw = absolute(upw), k;
  k = 0;  D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 5;  D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 10; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 15; D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );

  MatMult4(_MODEL_NVARS_,DL,D,L);
  MatMult4(_MODEL_NVARS_,Jac,R,DL);

  return(0);
}

int NavierStokes2DStiffJacobian(
                          double  *Jac, 
                          double  *u,   
                          void    *p,   
                          int     dir,  
                          int     nvars,
                          int     upw   
                   )
{
  NavierStokes2D *param = (NavierStokes2D*) p;
  static double R[_MODEL_NVARS_*_MODEL_NVARS_], D[_MODEL_NVARS_*_MODEL_NVARS_], 
                L[_MODEL_NVARS_*_MODEL_NVARS_], DL[_MODEL_NVARS_*_MODEL_NVARS_];

  /* get the eigenvalues and left,right eigenvectors */
  _NavierStokes2DEigenvalues_      (u,D,param->gamma,dir);
  _NavierStokes2DLeftEigenvectors_ (u,L,param->gamma,dir);
  _NavierStokes2DRightEigenvectors_(u,R,param->gamma,dir);

  int aupw = absolute(upw), k;
  k = 0;  D[k] = absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) );
  k = 5;  D[k] = ( dir == _YDIR_ ? 0.0 : absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) ) );
  k = 10; D[k] = ( dir == _XDIR_ ? 0.0 : absolute( (1-aupw)*D[k] + 0.5*aupw*(1+upw)*max(0,D[k]) + 0.5*aupw*(1-upw)*min(0,D[k]) ) );
  k = 15; D[k] = 0.0;

  MatMult4(_MODEL_NVARS_,DL,D,L);
  MatMult4(_MODEL_NVARS_,Jac,R,DL);

  return(0);
}