FiniteVolumesEulerUnknowns.hpp

#ifndef FINITE_VOLUMES_EULER_UNKNOWNS_HPP
#define FINITE_VOLUMES_EULER_UNKNOWNS_HPP

template <typename TMeshData>
class FiniteVolumesEulerUnknowns
{
public:
  typedef TMeshData MeshDataType;
  typedef typename MeshDataType::MeshType MeshType;

  static constexpr size_t dimension = MeshType::dimension;
  typedef TinyVector<dimension> Rd;

private:
  const MeshDataType& m_mesh_data;
  const MeshType& m_mesh;

  Kokkos::View<double*> m_rhoj;
  Kokkos::View<Rd*> m_uj;
  Kokkos::View<double*> m_Ej;

  Kokkos::View<double*> m_ej;
  Kokkos::View<double*> m_pj;
  Kokkos::View<double*> m_gammaj;
  Kokkos::View<double*> m_cj;
  Kokkos::View<double*> m_mj;
  Kokkos::View<double*> m_inv_mj;
  Kokkos::View<double*> m_kj;

  Kokkos::View<Rd*> m_uL;
  Kokkos::View<Rd*> m_uR;
  Kokkos::View<double*> m_kL;
  Kokkos::View<double*> m_kR;

  Kokkos::View<double*> m_entropy;
  Kokkos::View<double*> m_entropy_new;
  Kokkos::View<double*> m_S0;

  Kokkos::View<double*> m_Tj;
  Kokkos::View<double*> m_nuj;
  Kokkos::View<double*> m_TL;
  Kokkos::View<double*> m_TR;
  Kokkos::View<double*> m_nuL;
  Kokkos::View<double*> m_nuR;

public:
  Kokkos::View<double*> rhoj()
  {
    return m_rhoj;
  }

  const Kokkos::View<const double*> rhoj() const
  {
    return m_rhoj;
  }

  Kokkos::View<Rd*> uj()
  {
    return m_uj;
  }

  const Kokkos::View<const Rd*> uj() const
  {
    return m_uj;
  }

  Kokkos::View<double*> Ej()
  {
    return m_Ej;
  }

  const Kokkos::View<const double*> Ej() const
  {
    return m_Ej;
  }

  Kokkos::View<double*> ej()
  {
    return m_ej;
  }

  const Kokkos::View<const double*> ej() const
  {
    return m_ej;
  }

  Kokkos::View<double*> pj()
  {
    return m_pj;
  }

  const Kokkos::View<const double*> pj() const
  {
    return m_pj;
  }

  Kokkos::View<double*> gammaj()
  {
    return m_gammaj;
  }

  const Kokkos::View<const double*> gammaj() const
  {
    return m_gammaj;
  }

  Kokkos::View<double*> cj()
  {
    return m_cj;
  }

  const Kokkos::View<const double*> cj() const
  {
    return m_cj;
  }

  Kokkos::View<double*> mj()
  {
    return m_mj;
  }

  const Kokkos::View<const double*> mj() const
  {
    return m_mj;
  }

  Kokkos::View<double*> invMj()
  {
    return m_inv_mj;
  }

  const Kokkos::View<const double*> invMj() const
  {
    return m_inv_mj;
  }

  Kokkos::View<double*> kj()
  {
    return m_kj;
  }

  const Kokkos::View<const double*> kj() const
  {
    return m_kj;
  }

  Kokkos::View<Rd*> uL()
  {
    return m_uL;
  }

  const Kokkos::View<const Rd*> uL() const
  {
    return m_uL;
  }

  Kokkos::View<Rd*> uR()
  {
    return m_uR;
  }

  const Kokkos::View<const Rd*> uR() const
  {
    return m_uR;
  }

  Kokkos::View<double*> kL()
  {
    return m_kL;
  }

  const Kokkos::View<const double*> kL() const
  {
    return m_kL;
  }

  Kokkos::View<double*> kR()
  {
    return m_kR;
  }

  const Kokkos::View<const double*> kR() const
  {
    return m_kR;
  }

  Kokkos::View<double*> nuj()
  {
    return m_nuj;
  }

  const Kokkos::View<const double*> nuj() const
  {
    return m_nuj;
  }

  Kokkos::View<double*> Tj()
  {
    return m_Tj;
  }

  const Kokkos::View<const double*> Tj() const
  {
    return m_Tj;
  }

  Kokkos::View<double*> TL()
  {
    return m_TL;
  }

  const Kokkos::View<const double*> TL() const
  {
    return m_TL;
  }

  Kokkos::View<double*> TR()
  {
    return m_TR;
  }

  const Kokkos::View<const double*> TR() const
  {
    return m_TR;
  }

  Kokkos::View<double*> nuL()
  {
    return m_nuL;
  }

  const Kokkos::View<const double*> nuL() const
  {
    return m_nuL;
  }

  Kokkos::View<double*> nuR()
  {
    return m_nuR;
  }

  const Kokkos::View<const double*> nuR() const
  {
    return m_nuR;
  }
  
  
  Kokkos::View<double*> entropy()
  {
    return m_entropy;
  }

  const Kokkos::View<const double*> entropy() const
  {
    return m_entropy;
  }

  Kokkos::View<double*> entropy_new()
  {
    return m_entropy_new;
  }

  const Kokkos::View<const double*> entropy_new() const
  {
    return m_entropy_new;
  }

  Kokkos::View<double*> S0()
  {
    return m_S0;
  }

  const Kokkos::View<const double*> S0() const
  {
    return m_S0;
  }


  /*
  // --- Acoustic Solver ---
  
  void initializeSod()
  {
    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	if (xj[j][0]<0.5) {
  	  m_rhoj[j]=1.;
  	} else {
  	  m_rhoj[j]=0.125;
	}
	
	//Kidder
	//m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	if (xj[j][0]<0.5) {
  	  m_pj[j]=1;
  	} else {
  	  m_pj[j]=0.1;
	}
	
	//Kidder
	//m_pj[j] = 2.*std::pow(m_rhoj[j],3);
      });
    
    double pi = 4.*std::atan(1.);
    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_uj[j] = zero;
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_gammaj[j] = 1.4;
	//m_gammaj[j] = 3.;
      });

    BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
    block_eos.updateEandCFromRhoP();

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
      });

    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_mj[j] = m_rhoj[j] * Vj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_inv_mj[j] = 1./m_mj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
  	m_kj[j] =  xj[j][0];
	//m_kj[j] = 0.5;

	
	// Sod

	// k non regulier
	
	if (xj[j][0]<0.7) {
  	  m_kj[j]=0.;
  	} else {
	  if (xj[j][0]<0.9){
	  
	    // Re = 28
	    // m_kj[j]=0.05;

	    // Re = 100
	    //m_kj[j] = 0.014;

	    // Re = 500
	    //m_kj[j] = 0.0028; 

	    // Re = 1000
	    //m_kj[j] = 0.0014;

	    // Re = 2000
	    //m_kj[j] = 0.0007;

	    // Re = 10 000
	    m_kj[j] = 0.00014;

	    // Re = 100 000
	    //m_kj[j] = 0.000014;

	  } else {
	    m_kj[j]=0. ;
	  }
	}
	

	// k regulier
	int n = 10.;
	m_kj[j] = std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.7)*(xj[j][0]<0.7+0.1/n) + std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.9-0.1/n)*(xj[j][0]<0.9) + (xj[j][0]>0.7+0.1/n)*(xj[j][0]<0.9-0.1/n);
	m_kj[j] = 0.014*m_kj[j];
	
	
      });

     // Conditions aux bords de Dirichlet sur u et k

    m_uL[0] = zero;
    m_uR[0] = zero;
    m_kL[0] = 0.;
    m_kR[0] = 1.;

    
    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){  	
	m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));
	m_S0(j) = m_entropy(j);
      });
    
  }
  */
  
  
   
  
  // DIFFUSION PURE

    void initializeSod()
    {
    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
    double pi = 4.*std::atan(1.);

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_rhoj[j]=1.;
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	  m_pj[j]=1.;
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_uj[j] = std::sin(pi*xj[j][0]);
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_gammaj[j] = 2.;
      });

    BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
    block_eos.updateEandCFromRhoP();

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	//m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
	m_Ej[j] = 2.;
      });

    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_mj[j] = m_rhoj[j] * Vj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_inv_mj[j] = 1./m_mj[j];
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	//m_kj[j] = 2.; // k constant
	m_kj[j] = xj[j][0]; // k non constant, k = x
      });

    // Conditions aux bords de Dirichlet sur u et k
    
    m_uL[0] = zero;
    m_uR[0] = zero;
    m_kL[0] = 0.;
    m_kR[0] = 1.;

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
      });

    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
	m_Tj[j] = m_ej[j];
      });

    // Conditions aux bords de Dirichlet sur T et nu
    
    m_TL[0] = 2.;
    m_TR[0] = 2.;
    m_nuL[0] = 0.5;
    m_nuR[0] = 1.;
  }


  FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
    : m_mesh_data(mesh_data),
      m_mesh(m_mesh_data.mesh()),
      m_rhoj("rhoj",m_mesh.numberOfCells()),
      m_uj("uj",m_mesh.numberOfCells()),
      m_Ej("Ej",m_mesh.numberOfCells()),
      m_nuj("nuj",m_mesh.numberOfCells()),
      m_Tj("Tj",m_mesh.numberOfCells()),
      m_ej("ej",m_mesh.numberOfCells()),
      m_pj("pj",m_mesh.numberOfCells()),
      m_gammaj("gammaj",m_mesh.numberOfCells()),
      m_cj("cj",m_mesh.numberOfCells()),
      m_mj("mj",m_mesh.numberOfCells()),
      m_inv_mj("inv_mj",m_mesh.numberOfCells()),
      m_kj("kj",m_mesh.numberOfCells()),
      m_uL("uL", 1),
      m_uR("uR", 1),
      m_kL("kL", 1),
      m_kR("kR", 1),
      m_TL("TL", 1),
      m_TR("TR", 1),
      m_nuL("nuL", 1),
      m_nuR("nuR", 1),
      m_entropy("entropy", m_mesh.numberOfCells()),
      m_entropy_new("entropy_new", m_mesh.numberOfCells()),
      m_S0("S0", m_mesh.numberOfCells())
  {
    ;
  }

};

#endif // FINITE_VOLUMES_EULER_UNKNOWNS_HPP