diff --git a/src/main.cpp b/src/main.cpp
index be703b8996425c34a09dcce7f9030495db8d3855..8a0dae3dec6bc199436379e5dad1b4750c5c7a64 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -153,7 +153,7 @@ int main(int argc, char *argv[])
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
- while((t<tmax) and (iteration<itermax)) {
+ while((t<tmax) and (iteration<1)) {
//double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
double dt = 0.8*finite_volumes_diffusion.diffusion_dt(rhoj, kj);
@@ -204,6 +204,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> uj = unknowns.uj();
double pi = 4.*std::atan(1.);
std::ofstream fout("comparaison u");
+ fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
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(-0.2) <<'\n'; // cas k non constant
@@ -215,6 +216,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const double*> Ej = unknowns.Ej();
double pi = 4.*std::atan(1.);
std::ofstream fout("comparaison E");
+ fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << Ej[j] << ' ' << (-(std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))+(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])))*0.5*(std::exp(-4.*pi*pi*0.2)-1.) + 2. <<'\n'; // cas k constant
}
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 4fb77fe85c07f23ef481ea855d965dd17da2c8cf..e9002afbe37e6dde7594a1c9d3da5f175de4c30b 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -50,7 +50,7 @@ public:
// pas constant
-
+ /*
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes())
@@ -60,11 +60,11 @@ public:
m_xr[r][0] = r*delta_x;
});
}
-
+*/
// pas non constant
- /*
+
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes())
@@ -85,7 +85,7 @@ public:
}
});
}
- */
+
~Mesh()
{
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 9eb4ce42fcd17ea2da2fe155a1637d73c5b20ac3..72e5c399659604e56e9b45c884c0075d3f560fbf 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -93,6 +93,7 @@ private:
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+ const Kokkos::View<const Rd*> xr = m_mesh.xr();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rdd sum = zero;
@@ -112,19 +113,32 @@ private:
});
// Conditions aux bords
- /*
- // Essai 1
+
+ /* // EN TRAVAUX
+ Rdd sum = zero;
+ double sum2 = 0.;
+ double sum3 = 0.;
+ for (int j=0; j<face_nb_cells(0); ++j) {
+ int cell_here = face_cells(l,j);
+ int local_face_number_in_cell = face_cell_local_face(l,j);
+ if (cell_here >= 0) {
+ sum -= tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell));
+ sum2 += kj(cell_here)*Vj(cell_here);
+ sum3 += Vj(cell_here);
+ */
+ // Essai 1 OK
int cell_here = face_cells(0,0);
int local_face_number_in_cell = face_cell_local_face(0,0);
- m_Fl(0) = -(kL(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
+ m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0);
- m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
-
- // Essai 2
- m_Fl(0) = -(kL(0) + kj(0))*0.5*(tensorProduct(uj(0), Cjr(0, 0)));
- m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(m_mesh.numberOfCells()-1))*0.5*(tensorProduct(uj(m_mesh.numberOfCells()-1), Cjr(m_mesh.numberOfCells()-1,1)));
- */
+ m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
+
+ double pi = 4.*std::atan(1.);;
+ std::ofstream fout("comparaison flux u");
+ for (int l=0; l<m_mesh.numberOfFaces(); ++l) {
+ fout << xr[l][0] << ' ' << m_Fl[l](0,0) << ' ' << 2.*pi*cos(pi*xr[l][0])*std::exp(-2*pi*pi*0.000015625) << std::endl;
+ }
return m_Fl ;
}
@@ -146,6 +160,8 @@ private:
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+ const Kokkos::View<const Rd*> xr = m_mesh.xr();
+
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rd sum = zero;
double sum2 = 0.;
@@ -163,6 +179,13 @@ private:
m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
+ // affichage
+ double pi = 4.*std::atan(1.);
+ std::ofstream fout("comparaison flux E");
+ for (int l=0; l<m_mesh.numberOfFaces(); ++l) {
+ fout << xr[l][0] << ' ' << m_Gl[l][0] << ' ' << sin(pi*xr[l][0])*2.*pi*cos(pi*xr[l][0])*std::exp(-4*pi*pi*0.000015625) << std::endl;
+ }
+
return m_Gl ;
}
@@ -288,9 +311,6 @@ public:
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
-
-
-
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR);