diff --git a/src/main.cpp b/src/main.cpp
index bc1424410fe2f7652da9ba489a860a6502fbf96b..90feecfff0bcef20d897cf777deabf334f644b61 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -19,6 +19,7 @@
#include <Mesh.hpp>
#include <AcousticSolver.hpp>
#include <FiniteVolumesDiffusion.hpp>
+#include <NoSplitting.hpp>
#include <TinyVector.hpp>
#include <TinyMatrix.hpp>
@@ -136,6 +137,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns);
FiniteVolumesDiffusion<MeshDataType> finite_volumes_diffusion(mesh_data, unknowns);
+ NoSplitting<MeshDataType> no_splitting(mesh_data, unknowns);
typedef TinyVector<MeshType::dimension> Rd;
@@ -166,6 +168,7 @@ int main(int argc, char *argv[])
double c = 0.;
c = finite_volumes_diffusion.conservatif(unknowns);
+ // Diagrammes de marche
/*
const Kokkos::View<const Rd*> xj = mesh_data.xj();
int size = 3000;
@@ -175,6 +178,8 @@ int main(int argc, char *argv[])
int i = 0;
*/
+ // Animations
+ /*
// Ecriture des valeurs initiales dans un fichier
const Kokkos::View<const Rd*> xj = mesh_data.xj();
@@ -271,27 +276,24 @@ int main(int argc, char *argv[])
}
}
diff.close();
-
+ */
while((t<tmax) and (iteration<itermax)) {
-
- // ETAPE 1 DU SPLITTING - EULER
-
- double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
-
+ // NAVIER-STOKES AVEC SPLITTING
+ /*
+ // Etape 1 du splitting - Euler
+ double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj)
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
}
acoustic_solver.computeNextStep(t,dt_euler, unknowns);
t += dt_euler;
- // ETAPE 2 DU SPLITTING - DIFFUSION
-
+ // Etape 2 du splitting - Diffusion
double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
double t_diff = t-dt_euler;
-
if (dt_euler <= dt_diff) {
dt_diff = dt_euler;
finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
@@ -305,11 +307,9 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
-
-
+ */
+ // Diffusion pure
/*
- // DIFFUSION PURE
-
double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj,nuL,nuR,kL,kR);
if (t+dt > tmax) {
dt = tmax-t;
@@ -317,6 +317,16 @@ int main(int argc, char *argv[])
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
*/
+
+ // NAVIER-STOKES SANS SPLITTING
+
+ double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+ if (t+dt > tmax) {
+ dt = tmax-t;
+ }
+ no_splitting.computeNextStep(t,dt, unknowns);
+ t += dt;
+
block_eos.updatePandCFromRhoE();
@@ -324,8 +334,8 @@ int main(int argc, char *argv[])
std::cout << "temps t : " << t << std::endl;
- // ECRITURE DANS UN FICHIER
-
+ // Animations
+ /*
if ((std::fmod(t,0.001) < 0.0001) or (t == tmax)) {
std::string ligne;
@@ -542,13 +552,13 @@ int main(int argc, char *argv[])
riffout.close();
}
-
+ */
- // ENTROPY TEST
+ // Entropy test
//finite_volumes_diffusion.entropie(unknowns);
+ // Diagrammes de marche
/*
- // STOCKAGE COORDONNES ET TEMPS
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
x[i][j] = xj[j][0];
rho_marche[i][j] = rhoj[j];
@@ -564,21 +574,11 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
- // CREATION FICHIER x = f(t)
+
+ // Diagrammes de marche, creation fichier x = f(t)
/*
std::ofstream fout("cara");
fout.precision(15);
- for (size_t j=0; j<mesh.numberOfCells(); ++j) {
- for (size_t k=0; k<i; ++k) {
- fout << tempo[k] << ' ' << x[k][j] << '\n';
- }
- fout << ' ' << '\n';
- fout << ' ' << '\n';
- }
- */
- /*
- std::ofstream fout("cararho");
- fout.precision(15);
for (int j=0; j<mesh.numberOfCells(); ++j) {
if (j%10 == 0) {
for (int k=0; k<i; ++k) {
@@ -592,56 +592,52 @@ int main(int argc, char *argv[])
}
*/
+
+ // Erreurs sous differents normes
/*
+
+ // Density
double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
-
std::cout << "* " << rang::style::underline << "Erreur L2 rho" << rang::style::reset
<< ": " << rang::fgB::green << error1 << rang::fg::reset << " \n";
double error2 = 0.;
error2 = finite_volumes_diffusion.error_Linf_rho(unknowns, tmax);
-
std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset
<< ": " << rang::fgB::green << error2 << rang::fg::reset << " \n";
-
-
- double err0 = 0.;
- err0 = finite_volumes_diffusion.error_L1_u(unknowns, tmax);
-
- std::cout << "* " << rang::style::underline << "Erreur L1 u" << rang::style::reset
- << ": " << rang::fgB::green << err0 << rang::fg::reset << " \n";
+ // Vitesse
double error = 0.;
error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
-
std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
<< ": " << rang::fgB::green << error << rang::fg::reset << " \n";
-
double error4 = 0.;
error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);
-
std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
<< ": " << rang::fgB::green << error4 << rang::fg::reset << " \n";
-
- double err1 = 0.;
- err1 = finite_volumes_diffusion.error_L1_E(unknowns, tmax);
-
- std::cout << "* " << rang::style::underline << "Erreur L1 E" << rang::style::reset
- << ": " << rang::fgB::green << err1 << rang::fg::reset << " \n";
-
+ double err0 = 0.;
+ err0 = finite_volumes_diffusion.error_L1_u(unknowns, tmax);
+ std::cout << "* " << rang::style::underline << "Erreur L1 u" << rang::style::reset
+ << ": " << rang::fgB::green << err0 << rang::fg::reset << " \n";
+
+ // Energy
double error3 = 0.;
error3 = finite_volumes_diffusion.error_L2_E(unknowns,tmax);
-
std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset
<< ": " << rang::fgB::green << error3 << rang::fg::reset << " \n";
double error5 = 0.;
error5 = finite_volumes_diffusion.error_Linf_E(unknowns, tmax);
-
std::cout << "* " << rang::style::underline << "Erreur L infini E" << rang::style::reset
<< ": " << rang::fgB::green << error5 << rang::fg::reset << " \n";
+
+ double err1 = 0.;
+ err1 = finite_volumes_diffusion.error_L1_E(unknowns, tmax);
+ std::cout << "* " << rang::style::underline << "Erreur L1 E" << rang::style::reset
+ << ": " << rang::fgB::green << err1 << rang::fg::reset << " \n";
+
*/
std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index adf4b739b6f8c5852f6affb41488335c72267441..4313b3f73e4d2ef0a37feb0312c041f67c859fd9 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -76,15 +76,14 @@ public:
// pas constant
-
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes()),
m_x0("x0", 1),
m_xmax("xmax", 1)
{
- double a = -1.;
- double b = 2.;
+ double a = 0.;
+ double b = 1.;
m_x0[0][0] = a;
m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells();
@@ -95,7 +94,6 @@ public:
// pas non constant
-
/*
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
@@ -127,7 +125,6 @@ public:
// pas aleatoire
-
/*
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
@@ -163,7 +160,6 @@ public:
*/
-
~Mesh()
{
;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 1a02d8e0951dbe2728da0c26ff98888c10d6ab6d..7351f938e0bb925acdec2838503ae00b71124885 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -182,18 +182,18 @@ private:
});
- m_ur[0]=zero;
- m_ur[m_mesh.numberOfNodes()-1]=zero;
+ //m_ur[0]=zero;
+ //m_ur[m_mesh.numberOfNodes()-1]=zero;
//m_ur[0] = x0;
//m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
// CL Kidder
- /*
+
double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
- */
+
return m_ur;
}
@@ -372,9 +372,9 @@ public:
});
*/
- /*
+
// Mise a jour de nu
-
+ /*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
nuj(j) = 0.5*(1.+xj[j][0]);
});
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 19079bb705b1794c29a83c404cfee47aba10478b..2b78943cdf57d9719d22932f528a58ea04d0e299 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -284,27 +284,29 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
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.);
+ */
+ // 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);
+ */
+ // Kidder
+ m_pj[j] = 2.*std::pow(m_rhoj[j],3);
});
double pi = 4.*std::atan(1.);
@@ -313,8 +315,10 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_gammaj[j] = 1.4;
- //m_gammaj[j] = 3.;
+ // TAC
+ // m_gammaj[j] = 1.4;
+ // Kidder
+ m_gammaj[j] = 3.;
});
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
@@ -334,20 +338,19 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- //m_kj[j] = xj[j][0];
-
- m_kj[j] = 0.;
+ // Differents k (xi)
+ m_kj[j] = xj[j][0];
+ //m_kj[j] = 0.;
- /*
- // Sod
+ // 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;
@@ -374,8 +377,9 @@ public:
}
}
*/
- /*
+
// 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];
@@ -388,7 +392,7 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
m_kL[0] = 0.;
- m_kR[0] = 0.;
+ m_kR[0] = 1.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -413,7 +417,6 @@ public:
}
-
/*
// DIFFUSION PURE
@@ -432,7 +435,9 @@ public:
});
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];
});
@@ -444,8 +449,9 @@ public:
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]);
+ // 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.;
});
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f046a3d8664ae65fd7d5c437c20537a6032db128
--- /dev/null
+++ b/src/scheme/NoSplitting.hpp
@@ -0,0 +1,376 @@
+#ifndef NO_SPLITTING_HPP
+#define NO_SPLITTING_HPP
+
+// --- INCLUSION fichiers headers ---
+
+#include <Kokkos_Core.hpp>
+
+#include <rang.hpp>
+#include <BlockPerfectGas.hpp>
+
+#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+#include <Mesh.hpp>
+#include <MeshData.hpp>
+#include <FiniteVolumesEulerUnknowns.hpp>
+
+// ---------------------------------
+
+// Creation classe NoSplitting
+
+template<typename MeshData>
+class NoSplitting
+{
+ typedef typename MeshData::MeshType MeshType;
+ typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType;
+
+ MeshData& m_mesh_data;
+ MeshType& m_mesh;
+ const typename MeshType::Connectivity& m_connectivity;
+
+ constexpr static size_t dimension = MeshType::dimension;
+
+ typedef TinyVector<dimension> Rd;
+ typedef TinyMatrix<dimension> Rdd;
+
+private:
+
+ struct ReduceMin
+ {
+ private:
+ const Kokkos::View<const double*> x_;
+
+ public:
+ typedef Kokkos::View<const double*>::non_const_value_type value_type;
+
+ ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {}
+
+ typedef Kokkos::View<const double*>::size_type size_type;
+
+ KOKKOS_INLINE_FUNCTION void
+ operator() (const size_type i, value_type& update) const
+ {
+ if (x_(i) < update) {
+ update = x_(i);
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION void
+ join (volatile value_type& dst,
+ const volatile value_type& src) const
+ {
+ if (src < dst) {
+ dst = src;
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION void
+ init (value_type& dst) const
+ { // The identity under max is -Inf.
+ dst= Kokkos::reduction_identity<value_type>::min();
+ }
+ };
+
+
+ // ----
+
+ KOKKOS_INLINE_FUNCTION
+ const Kokkos::View<const double*>
+ computePj(const Kokkos::View<const double*>& pj,
+ const Kokkos::View<const double*>& kj,
+ const Kokkos::View<const Rd*>& uj)
+ {
+ const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+
+ const Kokkos::View<const unsigned short*> cell_nb_nodes
+ = m_connectivity.cellNbNodes();
+
+ const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ double new_p = 0.;
+ double sum = 0;
+ for (int m=0; m<cell_nb_nodes(j); ++m) {
+ sum += uj[cell_nodes(j,m)][0];
+ }
+ new_p = pj(j) - kj(j)*sum/Vj(j);
+ pj(j) = new_p;
+ });
+ return pj;
+ }
+
+ // ----
+
+
+ KOKKOS_INLINE_FUNCTION
+ const Kokkos::View<const double*>
+ computeRhoCj(const Kokkos::View<const double*>& rhoj,
+ const Kokkos::View<const double*>& cj)
+ {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ m_rhocj[j] = rhoj[j]*cj[j];
+ });
+ return m_rhocj;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const Kokkos::View<const Rdd**>
+ computeAjr(const Kokkos::View<const double*>& rhocj,
+ const Kokkos::View<const Rd**>& Cjr) {
+ const Kokkos::View<const unsigned short*> cell_nb_nodes
+ = m_connectivity.cellNbNodes();
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ for (int r=0; r<cell_nb_nodes[j]; ++r) {
+ m_Ajr(j,r) = tensorProduct(rhocj(j)*Cjr(j,r), Cjr(j,r));
+ }
+ });
+
+ return m_Ajr;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const Kokkos::View<const Rdd*>
+ computeAr(const Kokkos::View<const Rdd**>& Ajr) {
+ const Kokkos::View<const unsigned int**> node_cells = m_connectivity.nodeCells();
+ const Kokkos::View<const unsigned short**> node_cell_local_node = m_connectivity.nodeCellLocalNode();
+ const Kokkos::View<const unsigned short*> node_nb_cells = m_connectivity.nodeNbCells();
+
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Rdd sum = zero;
+ for (int j=0; j<node_nb_cells(r); ++j) {
+ const int J = node_cells(r,j);
+ const int R = node_cell_local_node(r,j);
+ sum += Ajr(J,R);
+ }
+ m_Ar(r) = sum;
+ });
+
+ return m_Ar;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const Kokkos::View<const Rd*>
+ computeBr(const Kokkos::View<const Rdd**>& Ajr,
+ const Kokkos::View<const Rd**>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const double*>& pj,
+ const double t) {
+ const Kokkos::View<const unsigned int**>& node_cells = m_connectivity.nodeCells();
+ const Kokkos::View<const unsigned short**>& node_cell_local_node = m_connectivity.nodeCellLocalNode();
+ const Kokkos::View<const unsigned short*>& node_nb_cells = m_connectivity.nodeNbCells();
+ Kokkos::View<Rd*> xr = m_mesh.xr();
+
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Rd& br = m_br(r);
+ br = zero;
+ for (int j=0; j<node_nb_cells(r); ++j) {
+ const int J = node_cells(r,j);
+ const int R = node_cell_local_node(r,j);
+ br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
+ }
+ });
+
+ return m_br;
+ }
+
+ Kokkos::View<Rd*>
+ computeUr(const Kokkos::View<const Rdd*>& Ar,
+ const Kokkos::View<const Rd*>& br,
+ const double& t) {
+
+ inverse(Ar, m_inv_Ar);
+ const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
+
+ Kokkos::View<Rd*> xr = m_mesh.xr();
+ Kokkos::View<Rd*> x0 = m_mesh.x0();
+ Kokkos::View<Rd*> xmax = m_mesh.xmax();
+
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ m_ur[r]=invAr(r)*br(r);
+ });
+
+
+ //m_ur[0]=zero;
+ //m_ur[m_mesh.numberOfNodes()-1]=zero;
+
+ //m_ur[0] = x0;
+ //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
+
+ // CL Kidder
+
+ double h = std::sqrt(1. - (t*t)/(50./9.));
+ m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
+ m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
+
+
+ return m_ur;
+ }
+
+ Kokkos::View<Rd**>
+ computeFjr(const Kokkos::View<const Rdd**>& Ajr,
+ const Kokkos::View<const Rd*>& ur,
+ const Kokkos::View<const Rd**>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const double*>& pj) {
+ const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+ const Kokkos::View<const unsigned short*> cell_nb_nodes
+ = m_connectivity.cellNbNodes();
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ for (int r=0; r<cell_nb_nodes[j]; ++r) {
+ m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
+ }
+ });
+
+ return m_Fjr;
+ }
+
+ // Calcul la liste des inverses d'une liste de matrices (pour
+ // l'instant seulement $R^{1\times 1}$)
+ void inverse(const Kokkos::View<const Rdd*>& A,
+ Kokkos::View<Rdd*>& inv_A) const {
+ Kokkos::parallel_for(A.size(), KOKKOS_LAMBDA(const int& r) {
+ inv_A(r) = Rdd{1./(A(r)(0,0))};
+ });
+ }
+
+ // Calcul la liste des inverses d'une liste de reels
+ void inverse(const Kokkos::View<const double*>& x,
+ Kokkos::View<double*>& inv_x) const {
+ Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
+ inv_x(r) = 1./x(r);
+ });
+ }
+
+ // Enchaine les operations pour calculer les flux (Fjr et ur) pour
+ // pouvoir derouler le schema
+ KOKKOS_INLINE_FUNCTION
+ void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
+ const Kokkos::View<const Rd*>& xj,
+ const Kokkos::View<const double*>& kj,
+ const Kokkos::View<const double*>& rhoj,
+ const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const double*>& pj,
+ const Kokkos::View<const double*>& cj,
+ const Kokkos::View<const double*>& Vj,
+ const Kokkos::View<const Rd**>& Cjr,
+ const double& t) {
+ const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
+ const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
+
+ const Kokkos::View<const double*> PTj = computePj(pj, kj, uj);
+
+ const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
+ const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, PTj, t);
+
+ Kokkos::View<Rd*> ur = m_ur;
+ Kokkos::View<Rd**> Fjr = m_Fjr;
+ ur = computeUr(Ar, br, t);
+ Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
+ }
+
+ Kokkos::View<Rd*> m_br;
+ Kokkos::View<Rdd**> m_Ajr;
+ Kokkos::View<Rdd*> m_Ar;
+ Kokkos::View<Rdd*> m_inv_Ar;
+ Kokkos::View<Rd**> m_Fjr;
+ Kokkos::View<Rd*> m_ur;
+ Kokkos::View<double*> m_rhocj;
+ Kokkos::View<double*> m_PTj;
+ Kokkos::View<double*> m_Vj_over_cj;
+
+public:
+ NoSplitting(MeshData& mesh_data,
+ UnknownsType& unknowns)
+ : m_mesh_data(mesh_data),
+ m_mesh(mesh_data.mesh()),
+ m_connectivity(m_mesh.connectivity()),
+ m_br("br", m_mesh.numberOfNodes()),
+ m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
+ m_Ar("Ar", m_mesh.numberOfNodes()),
+ m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
+ m_Fjr("Fjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
+ m_ur("ur", m_mesh.numberOfNodes()),
+ m_rhocj("rho_c", m_mesh.numberOfCells()),
+ m_PTj("PTj", m_mesh.numberOfCells()),
+ m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+ {
+ ;
+ }
+
+ // Avance la valeur des inconnues pendant un pas de temps dt
+ void computeNextStep(const double& t, const double& dt,
+ UnknownsType& unknowns)
+ {
+ Kokkos::View<double*> rhoj = unknowns.rhoj();
+ Kokkos::View<Rd*> uj = unknowns.uj();
+ Kokkos::View<double*> Ej = unknowns.Ej();
+
+ Kokkos::View<double*> ej = unknowns.ej();
+ Kokkos::View<double*> pj = unknowns.pj();
+ Kokkos::View<double*> gammaj = unknowns.gammaj();
+ Kokkos::View<double*> cj = unknowns.cj();
+ Kokkos::View<double*> kj = unknowns.kj();
+ Kokkos::View<double*> nuj = unknowns.nuj();
+
+ const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+ const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+ const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
+ Kokkos::View<Rd*> xr = m_mesh.xr();
+
+ // Calcule les flux
+ computeExplicitFluxes(xr, xj, kj, rhoj, uj, pj, cj, Vj, Cjr, t);
+
+ const Kokkos::View<const Rd**> Fjr = m_Fjr;
+ const Kokkos::View<const Rd*> ur = m_ur;
+ const Kokkos::View<const unsigned short*> cell_nb_nodes
+ = m_connectivity.cellNbNodes();
+ const Kokkos::View<const unsigned int**>& cell_nodes
+ = m_connectivity.cellNodes();
+
+ // Mise a jour de la vitesse et de l'energie totale specifique
+ const Kokkos::View<const double*> inv_mj = unknowns.invMj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (int R=0; R<cell_nb_nodes[j]; ++R) {
+ const int r=cell_nodes(j,R);
+ momentum_fluxes += Fjr(j,R);
+ energy_fluxes += (Fjr(j,R), ur[r]);
+ }
+ uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
+ Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
+ });
+
+ // Calcul de e par la formule e = E-0.5 u^2
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
+ });
+
+ // deplace le maillage (ses sommets) en utilisant la vitesse
+ // donnee par le schema
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
+ xr[r] += dt*ur[r];
+ });
+
+ // met a jour les quantites (geometriques) associees au maillage
+ m_mesh_data.updateAllData();
+
+ // Calcul de rho avec la formule Mj = Vj rhoj
+ const Kokkos::View<const double*> mj = unknowns.mj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ rhoj[j] = mj[j]/Vj[j];
+ });
+
+ // Mise a jour de k
+ /*
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ kj(j) = xj[j][0];
+
+ });
+ */
+ }
+};
+
+#endif // NO_SPLITTING_HPP