diff --git a/src/main.cpp b/src/main.cpp
index fad893502dbbc95dff1c5a310ca43a9ef6e5dd45..4713c168ee33421c75cff992ae403c9ee38f607b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -292,24 +292,6 @@ int main(int argc, char *argv[])
}
}
-
- // AUTRE APPROCHE DU SPLITTING (PLUS LONG)
- /*
- double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
- double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
- double dt = 0.;
- if (dt_euler < dt_diff) {
- dt = dt_euler;
- } else {
- dt = dt_diff;
- }
- if (t+dt > tmax) {
- dt = tmax-t;
- }
- acoustic_solver.computeNextStep(t,dt, unknowns);
- finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
- t += dt;
- */
/*
// DIFFUSION PURE
@@ -654,20 +636,6 @@ int main(int argc, char *argv[])
}
}
- /*
- { // gnuplot output u^2*0.5 + T
- const Kokkos::View<const Rd*> xj = mesh_data.xj();
- const Kokkos::View<const double*> Tj = unknowns.Tj();
- const Kokkos::View<const Rd*> uj = unknowns.uj();
- double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
- std::ofstream fout("essai");
- fout.precision(15);
- for (size_t j=0; j<mesh.numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << Tj[j]+uj[j][0]*uj[j][0]*0.5 << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; // kidder
- }
- }
- */
-
/*
{ // gnuplot output for entropy (gaz parfait en prenant cv = 1))
const Kokkos::View<const Rd*> xj = mesh_data.xj();
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index eb13817957b5eacd4598d4b6b6afaf73595406ef..a1f4dd88f29a5717377b3f77a428333671b887b2 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -38,7 +38,7 @@ private:
// -------
// Sert a calculer les reductions (en gros calculer le min sur des
- // vecteurs en parallele) Ne pas regarder plus comment ca marche.
+ // vecteurs en parallele)
struct ReduceMin
{
private:
@@ -164,22 +164,6 @@ private:
return m_br;
}
-
- /*
- Kokkos::View<Rd*> // calcule u_r (vitesse au sommet du maillage et flux de vitesse)
- computeUr(const Kokkos::View<const Rdd*>& Ar,
- const Kokkos::View<const Rd*>& br) {
- inverse(Ar, m_inv_Ar);
- const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
- 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;
-
- return m_ur;
- }
- */
Kokkos::View<Rd*> // calcule u_r (vitesse au sommet du maillage et flux de vitesse)
computeUr(const Kokkos::View<const Rdd*>& Ar,
@@ -201,13 +185,7 @@ private:
//m_ur[0]=zero;
//m_ur[m_mesh.numberOfNodes()-1]=zero;
- // Kidder
-
- // Conditions aux limites dependant du temps
- //m_ur[0]=(-t/((50./9.)-t*t))*xr[0];
- //m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
-
- //R(t) = x*h(t) a la place de x(t)
+ // CL Kidder
double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
@@ -273,7 +251,6 @@ private:
Kokkos::View<Rd*> ur = m_ur;
Kokkos::View<Rd**> Fjr = m_Fjr;
- //ur = computeUr(Ar, br);
ur = computeUr(Ar, br, t);
Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
}
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 772038f4ef5a8010016d87c1cd77a9059f0bd417..b3720dd661bdce8c7d584fddbddcb0ac80da8749 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -35,8 +35,6 @@ class FiniteVolumesDiffusion
private:
- // Sert a calculer les reductions (en gros calculer le min sur des
- // vecteurs en parallele)
struct ReduceMin
{
private:
@@ -119,11 +117,11 @@ 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) = -xr[0][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) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
+
*/
// Kidder
@@ -133,9 +131,6 @@ private:
//m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5;
// k = x
- //m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
- //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
-
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0];
m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
@@ -179,15 +174,12 @@ private:
});
// Conditions aux bords
- // m_Gl(0) = Fl(0)*uL(0);
+
+ //m_Gl(0) = Fl(0)*uL(0);
//m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
// Kidder
-
- // m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
- //m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
-
-
+
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0);
m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
@@ -246,16 +238,19 @@ private:
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));
*/
-
- double h = std::sqrt(1. - (t*t)/(50./9.));
+
+ // Kiddder
// nu = (1+x)*0.5
+ //double h = std::sqrt(1. - (t*t)/(50./9.));
//m_Bl(0) = ((1.+h*x0[0][0])*3.*h*x0[0][0])/(100.*h*h*h*h);
//m_Bl(m_mesh.numberOfFaces()-1) = ((1.+h*xmax[0][0])*3.*h*xmax[0][0])/(100.*h*h*h*h);
+
+
+ // nu = 0
m_Bl(0) = 0.;
m_Bl(m_mesh.numberOfFaces()-1) = 0.;
-
// nu = 0.2
//m_Bl(0) = (0.2*3.*h*x0[0][0])/(50.*h*h*h*h);
//m_Bl(m_mesh.numberOfFaces()-1) = (0.2*3.*h*xmax[0][0])/(50.*h*h*h*h);
@@ -423,8 +418,8 @@ public:
// double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (t*t)/(50./9.));
- // Diffusion pure
- //TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
+ // CL en diffusion pure
+ // TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
@@ -466,9 +461,9 @@ public:
// 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)));
- //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(0.2*3.)/(50.*h*h*h*h));
- //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(6*xj[j][0]+3.)/(100*h*h*h*h));
+ Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); // nu = 0
+ //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(0.2*3.)/(50.*h*h*h*h)); // nu = 0.2
+ //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(6*xj[j][0]+3.)/(100*h*h*h*h)); // nu = (1+x)*0.5
});
// Calcul de e par la formule e = E-0.5 u^2