#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<double*> m_PTj; 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<double*> PTj() { return m_PTj; } const Kokkos::View<const double*> PTj() const { return m_PTj; } 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(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ // TAC 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){ // TAC 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){ // TAC m_gammaj[j] = 1.4; // Kidder //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]); }); 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){ // Differents k (xi) //m_kj[j] = xj[j][0]; //m_kj[j] = 0.014; // TAC // 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 = 1; 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]; }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ m_PTj(j) = m_pj(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] = 0.; 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); }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ //m_nuj(j) = 0.5*(1.+xj[j][0]); m_nuj(j) = 0.; }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ //m_Tj[j] = m_ej[j]; m_Tj[j] = 0.; }); // Conditions aux bords de Dirichlet sur T et nu //m_TL[0] = m_ej[0]; //m_TR[0] = m_ej[m_mesh.numberOfCells()-1]; m_TL[0] = 0.; m_TR[0] = 0.; m_nuL[0] = 0.; m_nuR[0] = 0.; } /* // 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){ // Sans conduction thermique m_uj[j] = std::sin(pi*xj[j][0]); // Avec conduction thermique //m_uj[j] = 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){ // Sans conduction thermique m_Ej[j] = 2.; // Avec conduction thermique //m_Ej[j] = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1.; }); 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, k = 2 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_uR[0] = xj[0]; m_kL[0] = 0.; m_kR[0] = 1.; Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ //m_nuj(j) = 0.; 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] = m_ej[0]; //m_TR[0] = m_ej[m_mesh.numberOfCells()-1]; 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()), m_PTj("PTj",m_mesh.numberOfCells()) { ; } }; #endif // FINITE_VOLUMES_EULER_UNKNOWNS_HPP