Euler1DJacobian.c

Go to the documentation of this file.

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

int Euler1DJacobian(double  *Jac, 
                    double  *u,   
                    void    *p,   
                    int     dir,  
                    int     nvars,
                    int     upw    )
{
  Euler1D       *param = (Euler1D*) 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 */
  _Euler1DEigenvalues_      (u,D,param,0);
  _Euler1DLeftEigenvectors_ (u,L,param,0);
  _Euler1DRightEigenvectors_(u,R,param,0);

  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 = 4; 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 = 8; 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]) );

  MatMult3(3,DL,D,L);
  MatMult3(3,Jac,R,DL);

  return(0);
}

int Euler1DStiffJacobian( double  *Jac, 
                          double  *u,   
                          void    *p,   
                          int     dir,  
                          int     nvars,
                          int     upw    )
{
  Euler1D       *param = (Euler1D*) 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 */
  _Euler1DEigenvalues_      (u,D,param,0);
  _Euler1DLeftEigenvectors_ (u,L,param,0);
  _Euler1DRightEigenvectors_(u,R,param,0);

  int aupw = absolute(upw), k;
  k = 0; D[k] = 0.0; /* remove the entropy eigenmode */
  k = 4; 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 = 8; 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]) );

  MatMult3(3,DL,D,L);
  MatMult3(3,Jac,R,DL);

  return(0);
}