NavierStokes3DJacobian.c

Go to the documentation of this file.

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

int NavierStokes3DJacobian(
                    double  *Jac, 
                    double  *u,   
                    void    *p,   
                    int     dir,  
                    int     nvars,
                    int     upw   
                   )
{
  NavierStokes3D *param = (NavierStokes3D*) 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 */
  _NavierStokes3DEigenvalues_      (u,_NavierStokes3D_stride_,D,param->gamma,dir);
  _NavierStokes3DLeftEigenvectors_ (u,_NavierStokes3D_stride_,L,param->gamma,dir);
  _NavierStokes3DRightEigenvectors_(u,_NavierStokes3D_stride_,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 = 6;  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 = 12; 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 = 18; 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 = 24; 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]) );

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

  return(0);
}

int NavierStokes3DStiffJacobian(
                          double  *Jac, 
                          double  *u,   
                          void    *p,   
                          int     dir,  
                          int     nvars,
                          int     upw   
                   )
{
  NavierStokes3D *param = (NavierStokes3D*) 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 */
  _NavierStokes3DEigenvalues_      (u,_NavierStokes3D_stride_,D,param->gamma,dir);
  _NavierStokes3DLeftEigenvectors_ (u,_NavierStokes3D_stride_,L,param->gamma,dir);
  _NavierStokes3DRightEigenvectors_(u,_NavierStokes3D_stride_,R,param->gamma,dir);

  int aupw = absolute(upw), k;
  if (dir == _XDIR_) {
    k = 0;  D[k] = 0.0;
    k = 6;  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 = 12; D[k] = 0.0;
    k = 18; D[k] = 0.0;
    k = 24; 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]) );
  } else if (dir == _YDIR_) {
    k = 0;  D[k] = 0.0;
    k = 6;  D[k] = 0.0;
    k = 12; 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 = 18; D[k] = 0.0;
    k = 24; 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]) );
  } else if (dir == _ZDIR_) {
    k = 0;  D[k] = 0.0;
    k = 6;  D[k] = 0.0;
    k = 12; D[k] = 0.0;
    k = 18; 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 = 24; 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]) );
  }

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

  return(0);
}