diff --git a/src/main.cpp b/src/main.cpp
index e16fdcc70f84a7b2903b4f19ad52f688fe2f3c60..e9db133da38122340acbc45578b80374093cf328 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -20,6 +20,7 @@
#include <AcousticSolver.hpp>
#include <FiniteVolumesDiffusion.hpp>
#include <NoSplitting.hpp>
+#include <Montee.hpp>
#include <TinyVector.hpp>
#include <TinyMatrix.hpp>
@@ -138,6 +139,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);
+ Montee<MeshDataType> montee(mesh_data, unknowns);
typedef TinyVector<MeshType::dimension> Rd;
@@ -301,7 +303,7 @@ int main(int argc, char *argv[])
// NAVIER-STOKES AVEC SPLITTING
/*
// Etape 1 du splitting - Euler
- double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt_euler = 0.5*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
}
@@ -337,7 +339,7 @@ int main(int argc, char *argv[])
*/
// NAVIER-STOKES SANS SPLITTING
-
+ /*
double dt = 0.1*acoustic_solver.acoustic_dt(Vj, cj); // pour le cas xi = 0
//double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
double dt_diff = 0.9*no_splitting.nosplitting_dt(Vj,cj,rhoj,kj);
@@ -349,6 +351,21 @@ int main(int argc, char *argv[])
}
no_splitting.computeNextStep(t,dt, unknowns);
t += dt;
+ */
+
+ // NAVIER-STOKES SANS SPLITTING (MONTEE EN ORDRE)
+
+ double dt = 0.1*acoustic_solver.acoustic_dt(Vj, cj); // pour le cas xi = 0
+ double dt_diff = 0.9*montee.montee_dt(Vj,cj,rhoj,kj);
+ if (dt_diff < dt) {
+ dt = dt_diff;
+ }
+ if (t+dt > tmax) {
+ dt = tmax-t;
+ }
+ montee.computeNextStep(t,dt, unknowns);
+ t += dt;
+
/*
// Fichier pas de temps
@@ -691,8 +708,8 @@ int main(int argc, char *argv[])
std::ofstream fout("rho");
fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder
- //fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+ //fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder
+ fout << xj[j][0] << ' ' << rhoj[j] << '\n';
}
}
@@ -736,10 +753,10 @@ int main(int argc, char *argv[])
//fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2) <<'\n'; //cas k constant
//fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-tmax) <<'\n'; // cas k non constant
- fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
+ //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
//fout << xj[j][0] << ' ' << uj[j][0] << ' ' << xj[j][0] << std::endl;
- //fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ fout << xj[j][0] << ' ' << uj[j][0] << '\n';
}
}
@@ -775,7 +792,7 @@ int main(int argc, char *argv[])
}
}
*/
- /*
+
{ // gnuplot output for k
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> kj = unknowns.kj();
@@ -789,7 +806,7 @@ int main(int argc, char *argv[])
}
}
}
- */
+
/*
{ // gnuplot output for vitesse
const Kokkos::View<const Rd*> xj = mesh_data.xj();
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index dd3b42900f4edde1c4d98c33a520cf103f782dbd..cee56e131148a1010e7813a8914c80b29607e1dd 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -533,9 +533,9 @@ public:
double err_u = 0.;
double exact_u = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
+ //exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
//exact_u = std::sin(pi*xj[j][0])*std::exp(-t); // solution exacte cas test k non constant
- //exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
+ exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
//exact_u = xj[j][0];
err_u += std::abs(exact_u - uj[j][0])*Vj(j);
}
@@ -639,7 +639,7 @@ public:
double h = std::sqrt(1. - (t*t)/(50./9.));
double exacte = std::sqrt((4.*((xj[0][0]*xj[0][0])/(h*h)) + 100.-(xj[0][0]*xj[0][0])/(h*h))/100.)/h;
double erreur = std::abs(exacte - rhoj[0]);
-
+
for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
exacte = std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h;
if (std::abs(exacte - rhoj[j]) > erreur) {
@@ -658,11 +658,13 @@ public:
double pi = 4.*std::atan(1.);
- double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*t);
+ //double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*t);
+ double exacte = -(xj[0][0]*t)/((50./9.)-t*t);
double erreur = std::abs(exacte - uj[0][0]);
for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
- exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*t);
+ //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*t);
+ exacte = -(xj[j][0]*t)/((50./9.)-t*t);
if (std::abs(exacte - uj[j][0]) > erreur) {
erreur = std::abs(exacte - uj[j][0]);
}
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index a3b67c8c4a29f317b54bf8f42d59c91bdc092348..b7ba56098c34816cd476696c88da7803a4475b4a 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -395,7 +395,7 @@ public:
/*
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];
+ m_kj[j] = 0.00014*m_kj[j];
*/
});
diff --git a/src/scheme/Montee.hpp b/src/scheme/Montee.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..80ad3cb3ef3be2859e03a40b345a5f3bb2b49c30
--- /dev/null
+++ b/src/scheme/Montee.hpp
@@ -0,0 +1,499 @@
+#ifndef MONTEE_HPP
+#define MONTEE_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 Montee
+
+template<typename MeshData>
+class Montee
+{
+ 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*>
+ computeRhoCj(const Kokkos::View<const double*>& kj,
+ const Kokkos::View<const double*>& Vj,
+ 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] + kj[j]/Vj[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_nb_cells
+ = m_connectivity.nodeNbCells();
+ const Kokkos::View<const unsigned short**>& node_cell_local_node
+ = m_connectivity.nodeCellLocalNode();
+ const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+ const Kokkos::View<const unsigned short*> cell_nb_nodes
+ = m_connectivity.cellNbNodes();
+
+ Kokkos::View<Rd*> xr = m_mesh.xr();
+ 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.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Rd& br = m_br(r);
+ br = zero;
+ for (int j=0; j<node_nb_cells(r); ++j) {
+
+ std::vector<double> stock_p(2);
+ std::vector<double> stock_u(2);
+
+ const int J = node_cells(r,j);
+ const int R = node_cell_local_node(r,j);
+
+ if ((r == 0) or (r == m_mesh.numberOfNodes()-1)) {
+
+ br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
+
+ } else {
+
+ for (int l=0; l<cell_nb_nodes(J); ++l) {
+ double sum_p = 0.;
+ double sum_u = 0.;
+ double sum = 0.;
+ int k = cell_nodes(J,l);
+
+ for (int i=0; i<node_nb_cells(k); ++i) {
+ int cell_here = node_cells(k,i);
+ sum_p += (1./Vj(cell_here))*pj(cell_here);
+ sum_u += (1./Vj(cell_here))*uj[cell_here][0];
+ sum += 1./Vj(cell_here);
+ }
+
+ stock_p[l] = sum_p/sum;
+ stock_u[l] = sum_u/sum;
+ }
+
+ br += Ajr(J,R)*(uj(J) + (stock_u[1]-stock_u[0])/Vj(J)*(xr[r]-xj[J])) + (pj(J) + (stock_p[1]-stock_p[0])/Vj(J)*(xr[r][0]-xj[J][0]))*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);
+ });
+
+ // --- CL ---
+
+ 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**>& node_cells = m_connectivity.nodeCells();
+ const Kokkos::View<const unsigned short*>& node_nb_cells
+ = m_connectivity.nodeNbCells();
+ const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+ const Kokkos::View<const unsigned short*> cell_nb_nodes
+ = m_connectivity.cellNbNodes();
+
+ Kokkos::View<Rd*> xr = m_mesh.xr();
+ 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) {
+
+ for (int r=0; r<cell_nb_nodes[j]; ++r) {
+ int R = cell_nodes(j,r);
+
+ if ((R == 0) or (R == m_mesh.numberOfNodes()-1)) {
+
+ m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(R))+pj(j)*Cjr(j,r);
+
+ } else {
+
+ std::vector<double> stock_p(2);
+ std::vector<double> stock_u(2);
+
+ for (int l=0; l<cell_nb_nodes(j); ++l) {
+ double sum_p = 0.;
+ double sum_u = 0.;
+ double sum = 0.;
+ int k = cell_nodes(j,l);
+
+ for (int i=0; i<node_nb_cells(k); ++i) {
+ int cell_here = node_cells(k,i);
+ sum_p += (1./Vj(cell_here))*pj(cell_here);
+ sum_u += (1./Vj(cell_here))*uj[cell_here][0];
+ sum += 1./Vj(cell_here);
+ }
+
+ stock_p[l] = sum_p/sum;
+ stock_u[l] = sum_u/sum;
+ }
+ m_Fjr(j,r) = Ajr(j,r)*((uj(j)+(stock_u[1]-stock_u[0])/Vj(j)*(xr[R]-xj[j])) - ur(R)) + (pj(j)+(stock_p[1]-stock_p[0])/Vj(j)*(xr[R][0]-xj[j][0]))*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*>& rhoj,
+ const Kokkos::View<const double*>& kj,
+ 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(kj, Vj, rhoj, cj);
+ const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
+
+ const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
+ const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj, 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<Rd*> m_ur0;
+ Kokkos::View<double*> m_rhocj;
+ Kokkos::View<double*> m_Vj_over_cj;
+
+public:
+ Montee(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_ur0("ur0", m_mesh.numberOfNodes()),
+ m_rhocj("rho_c", m_mesh.numberOfCells()),
+ m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+ {
+ ;
+ }
+
+ // Calcule une evaluation du pas de temps verifiant une CFL du type
+ // c*dt/dx<1 (c modifie)
+ KOKKOS_INLINE_FUNCTION
+ double montee_dt(const Kokkos::View<const double*>& Vj,
+ const Kokkos::View<const double*>& cj,
+ const Kokkos::View<const double*>& rhoj,
+ const Kokkos::View<const double*>& kj) const {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ m_Vj_over_cj[j] = Vj[j]/(cj[j]+kj[j]/(rhoj[j]*Vj[j]));
+ });
+
+ double dt = std::numeric_limits<double>::max();
+ Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(m_Vj_over_cj), dt);
+
+ return dt;
+ }
+
+
+ // 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();
+ Kokkos::View<double*> PTj = unknowns.PTj();
+ Kokkos::View<double*> kL = unknowns.kL();
+ Kokkos::View<double*> kR = unknowns.kR();
+ Kokkos::View<Rd*> uL = unknowns.uL();
+ Kokkos::View<Rd*> uR = unknowns.uR();
+
+ const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+ const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+ const Kokkos::View<const double*> Vl = m_mesh_data.Vl();
+ const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
+ Kokkos::View<Rd*> xr = m_mesh.xr();
+
+ 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 Rd**> Fjr = m_Fjr;
+ const Kokkos::View<const Rd*> ur = m_ur;
+ const Kokkos::View<const Rd*> ur0 = m_ur0;
+
+ // Calcul de PT (3eme essai, avec uR du solveur de Riemann)
+
+ // initialisation pour t = 0
+ if (t==0) {
+ for (int r = 0; r<m_mesh.numberOfNodes(); r++) {
+ m_ur[r][0] = 0.;
+ }
+ }
+
+ // iteration pour point fixe ur
+ for (int itconv=0; itconv<100; ++itconv){
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ double sum = 0;
+ for (int k=0; k<cell_nb_nodes(j); ++k) {
+ int node_here = cell_nodes(j,k);
+ sum += (ur(node_here), Cjr(j,k));
+ }
+ PTj(j) = pj(j) - kj(j)*sum/Vj(j);
+ });
+
+ // Copie
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ m_ur0[inode][0]=m_ur[inode][0];
+ }
+
+ // Calcule les flux
+ computeExplicitFluxes(xr, xj, rhoj, kj, uj, PTj, cj, Vj, Cjr, t);
+
+ // Relaxation de ur
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ m_ur[inode][0]=0.5*m_ur[inode][0]+0.5*m_ur0[inode][0];
+ }
+
+ // Calcul erreur L1
+ double sum=0.;
+ for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
+ sum+=std::abs(m_ur0[inode][0]-m_ur[inode][0]);
+ }
+ sum/=double(m_mesh.numberOfNodes());
+ std::cout << " it " << itconv << " sum " << sum << std::endl;
+ if(sum<1.e-6) break;
+ }
+
+
+ // 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;
+
+ // ajout second membre pour kidder (k cst)
+ //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+ // ajout second membre pour kidder (k = x)
+ //uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t));
+ //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+ });
+
+ // 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 // MONTEE_HPP
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index f25ddc0a631db2c9347c92698d27e118749b8122..d5533e1145c5c23c4451539e9bd0909c970ee7ca 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -81,6 +81,15 @@ private:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
//m_rhocj[j] = rhoj[j]*cj[j];
m_rhocj[j] = rhoj[j]*cj[j] + kj[j]/Vj[j];
+ /*
+ double tempo = kj[0];
+ for (int i=0; i<m_mesh.numberOfCells(); i++) {
+ if (kj[i]>tempo) {
+ tempo = kj[i];
+ }
+ }
+ m_rhocj[j] = rhoj[j]*cj[j] + tempo/Vj[j];
+ */
});
return m_rhocj;
}
@@ -163,10 +172,10 @@ private:
});
// --- CL ---
-
+
m_ur[0]=zero;
m_ur[m_mesh.numberOfNodes()-1]=zero;
-
+
//m_ur[0] = x0;
//m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
@@ -380,21 +389,19 @@ public:
});
-
- std::ofstream fout2("pj");
- fout2.precision(15);
- for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- fout2 << xj[j][0] << ' ' << pj[j] << '\n';
- }
- std::ofstream fout3("pTj");
- fout3.precision(15);
- for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- fout3 << xj[j][0] << ' ' << PTj[j] << '\n';
- }
+ // Calcule les flux
+ computeExplicitFluxes(xr, xj, rhoj, kj, uj, PTj, cj, Vj, Cjr, t);
*/
-
+
// Calcul de PT (3eme essai, avec uR du solveur de Riemann)
+ // initialisation pour t = 0
+ if (t==0) {
+ for (int r = 0; r<m_mesh.numberOfNodes(); r++) {
+ m_ur[r][0] = 0.;
+ }
+ }
+
// iteration pour point fixe ur
for (int itconv=0; itconv<100; ++itconv){
@@ -405,7 +412,6 @@ public:
sum += (ur(node_here), Cjr(j,k));
}
PTj(j) = pj(j) - kj(j)*sum/Vj(j);
-
});
// Copie
@@ -418,7 +424,7 @@ public:
// Relaxation de ur
for (int inode=0;inode<m_mesh.numberOfNodes();++inode){
- m_ur[inode][0]=0.5*m_ur[inode][0]+0.5*m_ur0[inode][0];
+ m_ur[inode][0]=0.1*m_ur[inode][0]+0.9*m_ur0[inode][0];
}
// Calcul erreur L1
@@ -430,6 +436,7 @@ public:
std::cout << " it " << itconv << " sum " << sum << std::endl;
if(sum<1.e-6) break;
}
+
// Mise a jour de la vitesse et de l'energie totale specifique
const Kokkos::View<const double*> inv_mj = unknowns.invMj();