diff --git a/src/main.cpp b/src/main.cpp
index 7c5afdf2917f902781151604f7370878bc963dc2..2478cef554b7b62b35d731c26069aa731a296afe 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -158,6 +158,11 @@ int main(int argc, char *argv[])
Kokkos::View<Rd*> uj = unknowns.uj();
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+ Kokkos::View<double*> nuL = unknowns.nuL();
+ Kokkos::View<double*> nuR = unknowns.nuR();
+ Kokkos::View<double*> kL = unknowns.kL();
+ Kokkos::View<double*> kR = unknowns.kR();
+
double c = 0.;
c = finite_volumes_diffusion.conservatif(unknowns);
@@ -303,7 +308,7 @@ int main(int argc, char *argv[])
// DIFFUSION PURE
- double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
+ double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj,nuL,nuR,kL,kR);
if (t+dt > tmax) {
dt = tmax-t;
}
@@ -540,14 +545,14 @@ int main(int argc, char *argv[])
// ENTROPY TEST
//finite_volumes_diffusion.entropie(unknowns);
-
+ /*
// STOCKAGE COORDONNES ET TEMPS
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
x[i][j] = xj[j][0];
}
tempo[i] = t;
i = i + 1;
-
+ */
}
@@ -567,6 +572,7 @@ int main(int argc, char *argv[])
fout << ' ' << '\n';
}
*/
+ /*
std::ofstream fout("cara1");
fout.precision(15);
for (int j=0; j<mesh.numberOfCells(); ++j) {
@@ -580,7 +586,7 @@ int main(int argc, char *argv[])
fout << ' ' << '\n';
}
}
-
+ */
/*
double error1 = 0.;
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index cacfa8cc76bdaa8164056bcf448e50e12767a767..49bc72950e56f50305afb270f7e918c98316413f 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -140,15 +140,16 @@ public:
m_xmax[0][0] = b;
const double h = (b-a)/connectivity.numberOfCells();
m_xr[0][0] = a;
- m_xr[connectivity.numberOfNodes()-1] = b;
+ m_xr[connectivity.numberOfNodes()-1][0] = b;
std::random_device rd;
std::mt19937 mt(rd());
- std::uniform_real_distribution<double> dist(-h/2.1,h/2.1);
+ std::uniform_real_distribution<double> dist(h/1.01,h);
for (int r=1; r<connectivity.numberOfNodes()-1; ++r){
const double delta_xr = dist(mt);
m_xr[r][0] = m_xr[r-1][0] + std::abs(delta_xr);
std::cout << m_xr[r][0] << std::endl;
- }
+ }
+ std::exit(0);
}
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index c83a21bcc67a53c78c366c44e4ff22ebf6839b7e..074d0f9572446c866192f578e1d4d724bd2517f9 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -117,11 +117,12 @@ private:
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))*(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)));
+ //m_Fl(0) = -((kL(0)/Vj(0) + kj(cell_here)/Vj(cell_here))/(1./Vj(0) + 1./Vj(cell_here)))*(1./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))*(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)));
-
+ //m_Fl(m_mesh.numberOfFaces()-1) = -((kR(0)/Vj(m_mesh.numberOfCells()-1) + kj(cell_here)/Vj(cell_here))/(1./Vj(cell_here) + 1./Vj(m_mesh.numberOfCells()-1)))*(1./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)));
// Kidder
@@ -237,10 +238,12 @@ private:
int cell_here = face_cells(0,0);
m_Bl(0) = (nuL(0) + nuj(cell_here))*(1./(2*Vl(0)))*(Tj(cell_here) - TL(0));
+ //m_Bl(0) = (nuL(0)/Vj(0) + nuj(cell_here)/Vj(cell_here))/(1./Vj(0) + 1./Vj(cell_here))*((Tj(cell_here)-TL(0))/Vl(0));
+
cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
m_Bl(m_mesh.numberOfFaces()-1) = -(nuR(0) + nuj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(Tj(cell_here) - TR(0));
-
+ //m_Bl(m_mesh.numberOfFaces()-1) = (nuR(0)/Vj(m_mesh.numberOfCells()-1) + nuj(cell_here)/Vj(cell_here))/(1./Vj(m_mesh.numberOfCells()-1) + 1./Vj(cell_here))*((Tj(cell_here)-TR(0))/Vl(m_mesh.numberOfFaces()-1));
// Kiddder
/*
@@ -331,7 +334,12 @@ public:
double diffusion_dt(const Kokkos::View<const double*>& rhoj,
const Kokkos::View<const double*>& kj,
const Kokkos::View<const double*>& nuj,
- const Kokkos::View<const double*>& cj) const {
+ const Kokkos::View<const double*>& cj,
+ const Kokkos::View<const double*>& nuL,
+ const Kokkos::View<const double*>& nuR,
+ const Kokkos::View<const double*>& kL,
+ const Kokkos::View<const double*>& kR
+ ) const {
Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells());
@@ -350,6 +358,7 @@ public:
= m_connectivity.cellNbNodes();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+
double minVl = std::numeric_limits<double>::max();
for (int ll=0; ll<cell_nb_faces(j); ++ll) {
minVl = std::min(minVl, Vl(cell_faces(j, ll)));
@@ -357,12 +366,34 @@ public:
double sum = 0.;
double sum1 = 0.;
+
+ double vol = 0.;
for (int m = 0; m < cell_nb_nodes(j); ++m) {
- sum += kj(cell_nodes(j,m));
- sum1 += nuj(cell_nodes(j,m));
+ if (cell_nodes(j,m) == 0) {
+ sum += kL(0)/Vj(0);
+ sum1 += nuL(0)/Vj(0);
+ vol += 1./Vj(0);
+ } else if (cell_nodes(j,m) == m_mesh.numberOfCells()) {
+ sum += kR(0)/Vj(m_mesh.numberOfCells()-1);
+ sum1 = nuR(0)/Vj(m_mesh.numberOfCells()-1);
+ vol += 1./Vj(m_mesh.numberOfCells()-1);
+ } else {
+ sum += kj(cell_nodes(j,m))/Vj(cell_nodes(j,m));
+ sum1 += nuj(cell_nodes(j,m))/Vj(cell_nodes(j,m));
+ vol += 1./Vj(cell_nodes(j,m));
+ }
+ }
+ sum = sum/vol;
+ sum1 = sum1/vol;
+
+ /*
+ for (int m = 0; m < cell_nb_nodes(j); ++m) {
+ sum += kj(cell_nodes(j,m));
+ sum1 += nuj(cell_nodes(j,m));
}
-
+ */
+
if (sum > sum1) {
if (sum == 0.) {
dt_j[j] = std::numeric_limits<double>::max();