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_Sj;

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*> Sj()
  {
   return m_Sj;
  }

  // --- 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;
  //	}
  //  });

  //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;
  //	}
  //  });

  //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;
  //  });

  //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] = 1; 
  //  });
  //}

  // ------------------

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){
	if (j == m_mesh.numberOfCells()) {
	  m_uj[j] = zero;
	}
	else {
	  if (j == 0) {
	    m_uj[j] = zero;
	  }
	  else {
	    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] = (std::cos(pi*xj[j][0]))*(std::cos(pi*xj[j][0]))*std::sin(pi*xj[j][0])*((pi*pi*pi)*0.5);
	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
      });

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


  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_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_Sj("Sj",m_mesh.numberOfCells())
  {
    ;
  }

};

#endif // FINITE_VOLUMES_EULER_UNKNOWNS_HPP