diff --git a/src/main.cpp b/src/main.cpp
index 51f9e385bcbd32ce2bc98a89dc1c7e6c4326a399..7c5afdf2917f902781151604f7370878bc963dc2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -142,7 +142,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const double*> Vj = mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
- const double tmax=0.8;
+ const double tmax=0.2;
double t=0.;
int itermax=std::numeric_limits<int>::max();
@@ -269,7 +269,7 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
-
+ /*
// ETAPE 1 DU SPLITTING - EULER
double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
@@ -298,8 +298,9 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
-
- /*
+ */
+
+
// DIFFUSION PURE
double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
@@ -308,7 +309,7 @@ int main(int argc, char *argv[])
}
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
- */
+
block_eos.updatePandCFromRhoE();
@@ -593,8 +594,14 @@ 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 err0 = 0.;
+ err0 = finite_volumes_diffusion.error_L1_u(unknowns, tmax);
+
+ std::cout << "* " << rang::style::underline << "Erreur L1 u" << rang::style::reset
+ << ": " << rang::fgB::green << err0 << rang::fg::reset << " \n";
+
double error = 0.;
error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
@@ -607,7 +614,13 @@ 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 err1 = 0.;
+ err1 = finite_volumes_diffusion.error_L1_E(unknowns, tmax);
+
+ std::cout << "* " << rang::style::underline << "Erreur L1 E" << rang::style::reset
+ << ": " << rang::fgB::green << err1 << rang::fg::reset << " \n";
+
double error3 = 0.;
error3 = finite_volumes_diffusion.error_L2_E(unknowns,tmax);
@@ -619,7 +632,7 @@ int main(int argc, char *argv[])
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";
@@ -650,35 +663,35 @@ int main(int argc, char *argv[])
{ // gnuplot output for vitesse
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const Rd*> uj = unknowns.uj();
- //double pi = 4.*std::atan(1.);
+ double pi = 4.*std::atan(1.);
std::ofstream fout("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(-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] << std::endl;
- fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ //fout << xj[j][0] << ' ' << uj[j][0] << '\n';
}
}
{ // gnuplot output for energy
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> Ej = unknowns.Ej();
- //double pi = 4.*std::atan(1.);
+ double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("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
+ 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
//fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((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.*tmax)-1.) + 2. <<'\n' ; // cas k non constant
//fout << xj[j][0] << ' ' << Ej[j] << ' ' << (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
//fout << xj[j][0] << ' ' << Ej[j] << ' ' << xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + tmax + 1. << std::endl;
- fout << xj[j][0] << ' ' << Ej[j] << '\n';
+ //fout << xj[j][0] << ' ' << Ej[j] << '\n';
}
}
@@ -711,7 +724,7 @@ int main(int argc, char *argv[])
}
}
*/
-
+ /*
{ // gnuplot output for vitesse
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj();
@@ -724,7 +737,7 @@ int main(int argc, char *argv[])
fout << xj[j][0] << ' ' << uj[j][0] + std::sqrt((gammaj[j]*pj[j])/rhoj[j]) << '\n';
}
}
-
+ */
}
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 470236c9e073d2670ebd4c2af54d74619bbc62e5..cacfa8cc76bdaa8164056bcf448e50e12767a767 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -3,6 +3,8 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
+#include <random>
+#include <iostream>
template <typename ConnectivityType>
class Mesh
@@ -73,7 +75,7 @@ public:
// pas constant
-
+ /*
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes()),
@@ -89,7 +91,7 @@ public:
m_xr[r][0] = a + r*delta_x;
});
}
-
+ */
// pas non constant
@@ -101,7 +103,7 @@ public:
m_xmax("xmax", 1)
{
double a = 0.;
- double b = 2.;
+ double b = 1.;
m_x0[0][0] = a;
m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells();
@@ -122,6 +124,36 @@ public:
}
*/
+
+ // pas aleatoire
+
+
+ Mesh(const Connectivity& connectivity)
+ : m_connectivity(connectivity),
+ m_xr("xr", connectivity.numberOfNodes()),
+ m_x0("x0", 1),
+ m_xmax("xmax", 1)
+ {
+ double a = 0.;
+ double b = 1.;
+ m_x0[0][0] = a;
+ m_xmax[0][0] = b;
+ const double h = (b-a)/connectivity.numberOfCells();
+ m_xr[0][0] = a;
+ m_xr[connectivity.numberOfNodes()-1] = b;
+ std::random_device rd;
+ std::mt19937 mt(rd());
+ std::uniform_real_distribution<double> dist(-h/2.1,h/2.1);
+ 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;
+ }
+ }
+
+
+
+
~Mesh()
{
;
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index daf8fd64e249fbd40664ed534bdde482fd16060d..c83a21bcc67a53c78c366c44e4ff22ebf6839b7e 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -484,10 +484,11 @@ public:
}
- // Calcul erreur entre solution analytique et solution numerique en norme L2
- // (quand la solution exacte est connue)
+ // Calcul erreur entre solution analytique et solution numerique pour differentes normes
- double error_L2_u(UnknownsType& unknowns, const double& t) {
+ // Norme L1
+
+ double error_L1_u(UnknownsType& unknowns, const double& t) {
Kokkos::View<Rd*> uj = unknowns.uj();
@@ -499,16 +500,38 @@ 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 += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
+ err_u += std::abs(exact_u - uj[j][0])*Vj(j);
}
- err_u = std::sqrt(err_u);
return err_u;
}
+ double error_L1_E(UnknownsType& unknowns, const double& t) {
+
+ Kokkos::View<double*> Ej = unknowns.Ej();
+
+ const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+
+ const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+
+ double pi = 4.*std::atan(1.);
+ double err_E = 0.;
+ 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.*t)-1.) + 2.;
+ //exact_E = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1. + t;
+ err_E += std::abs(exact_E - Ej[j])*Vj(j);
+ }
+ return err_E;
+ }
+
+
+ // Norme L2
+
double error_L2_rho(UnknownsType& unknowns, const double& t) {
Kokkos::View<double*> rhoj = unknowns.rhoj();
@@ -528,7 +551,28 @@ public:
return err_rho;
}
-
+ double error_L2_u(UnknownsType& unknowns, const double& t) {
+
+ Kokkos::View<Rd*> uj = unknowns.uj();
+
+ const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+
+ const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+
+ 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(-t); // solution exacte cas test k non constant
+ //exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
+ //exact_u = xj[j][0];
+ err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
+ }
+ err_u = std::sqrt(err_u);
+ return err_u;
+ }
+
double error_L2_E(UnknownsType& unknowns, const double& t) {
Kokkos::View<double*> Ej = unknowns.Ej();
@@ -541,18 +585,17 @@ public:
double err_E = 0.;
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 = (-(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.*t)-1.) + 2.;
- exact_E = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1. + t;
+ //exact_E = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1. + t;
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)
+
+ // Norme Linf
double error_Linf_rho(UnknownsType& unknowns, const double& t) {
@@ -575,9 +618,6 @@ 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_u(UnknownsType& unknowns, const double& t) {
Kokkos::View<Rd*> uj = unknowns.uj();
@@ -586,11 +626,11 @@ public:
double pi = 4.*std::atan(1.);
- double exacte = std::sin(pi*xj[0][0])*std::exp(-t);
+ double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*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(-t);
+ exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*t);
if (std::abs(exacte - uj[j][0]) > erreur) {
erreur = std::abs(exacte - uj[j][0]);
}
@@ -599,9 +639,6 @@ 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();
@@ -610,11 +647,11 @@ public:
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 exacte = (-(std::cos(pi*xj[0][0])*std::cos(pi*xj[0][0]))+(std::sin(pi*xj[0][0])*std::sin(pi*xj[0][0])))*0.5*(std::exp(-4.*pi*pi*0.2)-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.;
+ exacte = (-(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.;
if (std::abs(exacte - Ej[j]) > erreur) {
erreur = std::abs(exacte - Ej[j]);
}
@@ -643,6 +680,7 @@ public:
return sum;
}
+
// Verifie la conservativite de quantite de mvt
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 88a9500e781f68c584b16f538f2082cc26fdc35d..0e3719b8a4ab773074e395ccc6ca4a5eeb8d0c45 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -275,6 +275,8 @@ public:
return m_S0;
}
+
+ /*
// --- Acoustic Solver ---
void initializeSod()
@@ -336,7 +338,7 @@ public:
//m_kj[j] = 0.;
- /*
+
// Sod
// k non regulier
@@ -371,7 +373,7 @@ public:
m_kj[j]=0. ;
}
}
- */
+
// k regulier
int n = 1;
@@ -396,7 +398,7 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_nuj(j) = 0.5*(1.+xj[j][0]);
m_nuj(j) = 0.0005;
- /*
+
if (xj[j][0]<0.7) {
m_nuj[j]=0.;
} else {
@@ -406,7 +408,7 @@ public:
m_nuj[j]=0. ;
}
}
- */
+
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -421,10 +423,10 @@ public:
m_nuR[0] = 0.0005;
}
-
+ */
- /*
+
// DIFFUSION PURE
void initializeSod()
@@ -441,8 +443,8 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- //m_uj[j] = std::sin(pi*xj[j][0]);
- m_uj[j] = xj[j][0];
+ m_uj[j] = std::sin(pi*xj[j][0]);
+ //m_uj[j] = xj[j][0];
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -454,8 +456,8 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
- //m_Ej[j] = 2.;
- m_Ej[j] = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1.;
+ m_Ej[j] = 2.;
+ //m_Ej[j] = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1.;
});
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
@@ -468,20 +470,21 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_kj[j] = 0.; // k constant, k = 2
+ m_kj[j] = 2.; // k constant, k = 2
//m_kj[j] = xj[j][0]; // k non constant, k = x
});
// Conditions aux bords de Dirichlet sur u et k
m_uL[0] = zero;
- //m_uR[0] = zero;
- m_uR[0] = xj[0];
- m_kL[0] = 0.;
- m_kR[0] = 0.;
+ m_uR[0] = zero;
+ //m_uR[0] = xj[0];
+ m_kL[0] = 2.;
+ m_kR[0] = 2.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
+ m_nuj(j) = 0.;
+ //m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -489,13 +492,14 @@ public:
});
// Conditions aux bords de Dirichlet sur T et nu
-
- m_TL[0] = 1.;
- m_TR[0] = 3.;
- m_nuL[0] = 0.5;
- m_nuR[0] = 1.;
+
+ m_TL[0] = m_ej[0];
+ m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
+ m_nuL[0] = 0.;
+ m_nuR[0] = 0.;
+
}
- */
+
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)