diff --git a/src/main.cpp b/src/main.cpp
index b3e40857bb4859da61409dfb2a0fe364fde07a05..8b357fe5660acb63320be84a7afd891259f7195b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -275,7 +275,7 @@ int main(int argc, char *argv[])
// ETAPE 2 DU SPLITTING - DIFFUSION
- double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
+ double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj);
double t_diff = t-dt_euler;
if (dt_euler <= dt_diff) {
@@ -283,7 +283,7 @@ int main(int argc, char *argv[])
finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
} else {
while (t > t_diff) {
- dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
+ dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj);
if (t_diff+dt_diff > t) {
dt_diff = t-t_diff;
}
@@ -310,7 +310,7 @@ int main(int argc, char *argv[])
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
*/
-
+
// DIFFUSION PURE
double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
@@ -325,7 +325,7 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
-
+
/*
// ECRITURE DANS UN FICHIER
@@ -567,6 +567,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset
<< ": " << rang::fgB::green << error2 << rang::fg::reset << " \n";
+ */
double error = 0.;
error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
@@ -579,15 +580,19 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
<< ": " << rang::fgB::green << error4 << rang::fg::reset << " \n";
- */
- /*
+
double error3 = 0.;
- error3 = finite_volumes_diffusion.error_L2_E(unknowns);
+ 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";
+
std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
<< ": " << rang::fgB::green << c << rang::fg::reset << " \n";
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index d53830e188debc24a161f09707bd95c2b1f21536..c074e8e6ef8bd2f73fe6c26fb531bc68d7d4a306 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -441,7 +441,7 @@ 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))-(6*xj[j][0]+3.)/(100*(1-t*t/(50/9))*(1-t*t/(50/9))));
});
@@ -468,13 +468,13 @@ public:
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
- //double pi = 4.*std::atan(1.);
+ double pi = 4.*std::atan(1.);
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(-0.2); // solution exacte cas test k non constant
- exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
+ 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
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
}
err_u = std::sqrt(err_u);
@@ -500,8 +500,8 @@ public:
return err_rho;
}
- /*
- double error_L2_E(UnknownsType& unknowns) {
+
+ double error_L2_E(UnknownsType& unknowns, const double& t) {
Kokkos::View<double*> Ej = unknowns.Ej();
@@ -514,14 +514,14 @@ public:
double exact_E = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
//exact_E = (-(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.;
- exact_E = ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2.;
+ exact_E = ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*t)-1.) + 2.;
err_E += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j);
}
err_E = std::sqrt(err_E);
return err_E;
}
- */
-
+
+
// Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
// (quand la solution exacte est connue)
@@ -555,11 +555,13 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
- double exacte = -(xj[0][0]*t)/((50./9.)-t*t);
+ double pi = 4.*std::atan(1.);
+
+ double exacte = std::sin(pi*xj[0][0])*std::exp(-t);
double erreur = std::abs(exacte - uj[0][0]);
for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
- exacte = -(xj[j][0]*t)/((50./9.)-t*t);
+ exacte = std::sin(pi*xj[j][0])*std::exp(-t);
if (std::abs(exacte - uj[j][0]) > erreur) {
erreur = std::abs(exacte - uj[j][0]);
}
@@ -568,6 +570,29 @@ public:
return erreur;
}
+ // Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
+ // (quand la solution exacte est connue)
+
+ double error_Linf_E(UnknownsType& unknowns, const double& t) {
+
+ Kokkos::View<double*> Ej = unknowns.Ej();
+
+ const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+
+ double pi = 4.*std::atan(1.);
+
+ double exacte = ((xj[0][0]*pi*pi*0.5)*(std::sin(pi*xj[0][0])*std::sin(pi*xj[0][0]) - std::cos(xj[0][0]*pi)*std::cos(pi*xj[0][0])) - pi*0.5*std::sin(pi*xj[0][0])*std::cos(pi*xj[0][0]))*(std::exp(-2.*t)-1.) + 2.;
+ double erreur = std::abs(exacte - Ej[0]);
+
+ for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
+ exacte = ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*t)-1.) + 2.;
+ if (std::abs(exacte - Ej[j]) > erreur) {
+ erreur = std::abs(exacte - Ej[j]);
+ }
+ }
+
+ return erreur;
+ }
// Verifie la conservativite de E
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index bcf1443d12e955dce68cd964a5e195d0cf78b5c6..fdf9c24e74513fc0cc44dcd15b469fe16ef3c163 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -456,7 +456,7 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_Tj[j] = 2. - 0.5*std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]); // k = x
+ m_Tj[j] = m_ej[j];
});
// Conditions aux bords de Dirichlet sur T et nu