riemannSolverLinear.cpph

Go to the documentation of this file.

 
#include <algorithm>
 
#include "../../Utils/KernelUtils.h"
 
namespace kernels {
namespace aderdg {
namespace generic {
namespace c {
 
template <bool useFlux, bool useNCP, bool useMM, typename SolverType>
void riemannSolverLinear(SolverType& solver, double* FL, double* FR,
                         const double* const QL, const double* const QR,
                         const double t,
                         const double dt,
                         const tarch::la::Vector<Dimensions, double>& faceCentre,
                         const tarch::la::Vector<Dimensions, double>& dx,
                         const int direction) {
  /*
   * For the linear kernels, we need the material parameters in the
   * extrapolated predictor.
   * We compute the averages of the material parameters but
   * do not use them in the max eigenvalue calculation.
   * 
   * Consider numberOfParameters in:
   * QL,QR,QavL,QavR,Qdiff
   * 
   * Not necessary to consider numberOfParameters in
   * FL,FR,LL,LR
   */
 
  constexpr int numberOfVariables  = SolverType::NumberOfVariables;
  
  constexpr int numberOfParameters = SolverType::NumberOfParameters;
  constexpr int numberOfData       = numberOfVariables+numberOfParameters;
  
  constexpr int basisSize          = SolverType::Order+1;
  constexpr int basisSize2         = basisSize*basisSize;
  constexpr int order              = basisSize - 1;
 
  
  // Compute the average variables and parameters from the left and the right
  double QavL[numberOfData] = {0.0}; // size: numberOfVariables+numberOfParameters
  double QavR[numberOfData] = {0.0}; // size: numberOfVariables+numberOfParameters
  kernels::riemannsolvers::util::averageRiemannInputs<basisSize,numberOfData>(
      QL,exahype2::solvers::aderdg::ADERDGSolver::weights,QavL);
  kernels::riemannsolvers::util::averageRiemannInputs<basisSize,numberOfData>(
      QR,exahype2::solvers::aderdg::ADERDGSolver::weights,QavR);
 
  //get abs max eigenvalue
  const double maxHyperbolicEigenvalueL = solver.maxEigenvalue(QavL, faceCentre, dx, t, dt, direction);
  const double maxHyperbolicEigenvalueR = solver.maxEigenvalue(QavR, faceCentre, dx, t, dt, direction);
  const double smax = std::max(maxHyperbolicEigenvalueL, maxHyperbolicEigenvalueR);
 
  double Qdiff[numberOfData] = {0.0};
  double Qavg[numberOfData] = {0.0};
 
  idx2 idx_2d(Dimensions, numberOfVariables);
  
  double flux[numberOfVariables];
  
  double ncp[numberOfVariables] = {0.0};
  
  for(int l = 0; l < numberOfData; l++) {
    Qavg[l] = 0.5 *  (QavR[l] + QavL[l]);
  }
 
  //We copy the averaged material parameters to Qdiff as they are used in the flux term
  //These must not be overritten !
  for(int l = numberOfVariables; l < numberOfData; l++) {
    Qdiff[l] = Qavg[l];
  }
 
  {
    idx3 idx_FLR(basisSize,basisSize, numberOfVariables);
    idx3 idx_QLR(basisSize,basisSize, numberOfData);    
 
    std::fill_n (FL, basisSize2 * numberOfVariables, 0.0);
    std::fill_n (FR, basisSize2 * numberOfVariables, 0.0);
    for (int i = 0; i < basisSize; i++) {
      for (int j = 0; j < basisSize; j++) {
        for(int l = 0 ; l< numberOfVariables ; l++){ 
          Qdiff[l] = 0.5 * (QR[idx_QLR(i, j, l)] - QL[idx_QLR(i ,j, l)]);
        }
 
        if(useNCP) {
          
          solver.nonconservativeProduct(
            Qavg,
            Qdiff,
            faceCentre,
            dx,
            t,
            dt,
            direction,
            ncp
          );
/*        
          if(useMM){
            solver.multiplyMaterialParameterMatrix(Qavg, ncp);
          }
*/        
          for(int l=0; l < numberOfVariables; l++) {
            FL[idx_FLR(i, j, l)] += ncp[l]; 
            FR[idx_FLR(i, j, l)] += ncp[l]; 
          }
        }
           
            
        if(useFlux){
 
          solver.flux(
            Qdiff,
            faceCentre,
            dx,
            t,
            dt,
            direction,
            flux
          );
 
//          if(useMM){
//            solver.multiplyMaterialParameterMatrix(Qavg, flux);
//          }
        
          for(int l=0; l < numberOfVariables; l++) {
            FL[idx_FLR(i, j, l)] +=  flux[l];
            FR[idx_FLR(i, j, l)] +=  flux[l];
          }
        }
 
        for(int l=0; l < numberOfVariables; l++) {
          FL[idx_FLR(i, j, l)] -= smax*Qdiff[l];
          FR[idx_FLR(i, j, l)] += smax*Qdiff[l];
        }
      }
    }
  }
  
}
 
 
}  // namespace c
}  // namespace generic
}  // namespace aderdg
}  // namespace kernels