From 2e57e9de4f7ade7fc7716e655438ebd19f19f8ba Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Mon, 7 May 2018 09:47:27 +0200
Subject: [PATCH] Ajout fonction calcul conservativite et solution exacte pour
E quand k constant
---
src/main.cpp | 25 ++++++++++--
src/scheme/FiniteVolumesDiffusion.hpp | 48 +++++++++++++++--------
src/scheme/FiniteVolumesEulerUnknowns.hpp | 21 ++++++++--
3 files changed, 69 insertions(+), 25 deletions(-)
diff --git a/src/main.cpp b/src/main.cpp
index babd25dcb..e5adc9e7b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -189,17 +189,34 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini" << rang::style::reset
<< ": " << rang::fgB::green << error2 << rang::fg::reset << " \n";
+ double cons = 0.;
+ cons = finite_volumes_diffusion.conservatif(unknowns);
+
+ std::cout << "* " << rang::style::underline << "Resultat conservativite" << rang::style::reset
+ << ": " << rang::fgB::green << cons << rang::fg::reset << " \n";
+
+
//method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds();
- { // gnuplot output for density
+ { // 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.);
- std::ofstream fout("comparaison_u_k_non_cst_maill_non_unif_320");
+ std::ofstream fout("comparaison u");
+ 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
+ }
+ }
+
+ { // 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.);
+ std::ofstream fout("comparaison E");
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
+ fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::cos(pi*xj[j][0]))*(std::cos(pi*xj[j][0]))*std::sin(pi*xj[j][0])*((pi*pi*pi)*0.5)*std::exp(-4.*pi*0.2) <<'\n';
}
}
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index c01f1883c..ecac5a30b 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -93,10 +93,10 @@ private:
Rdd sum = zero;
double sum2 = 0.;
for (int j=0; j<face_nb_cells(l); ++j) {
- int cell_here = face_cells(l,j);
- int local_face_number_in_cell = face_cell_local_face(l,j);
- sum -= tensorProduct(uj(cell_here, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell));
- sum2 += kj(cell_here);
+ int cell_here = face_cells(l,j);
+ int local_face_number_in_cell = face_cell_local_face(l,j);
+ sum -= tensorProduct(uj(cell_here, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell));
+ sum2 += kj(cell_here);
}
m_Fl(l) = ((sum2/face_nb_cells(l))/Vl(l))*sum;
@@ -127,7 +127,7 @@ private:
sum += uj(cell_here, local_face_number_in_cell);
}
- m_Gl(l) = 0.5*Fl(l)*sum;
+ m_Gl(l) = (1./face_nb_cells(l))*Fl(l)*sum;
});
@@ -209,10 +209,6 @@ public:
minVl = std::min(minVl, Vl(cell_faces(j, ll)));
}
- //k=2 => (kj(j+1) + 2*kj(j) + kj(j-1)) = 8
- //dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./8.)*minVl;
-
- // k pas forcement constant
double sum = 0.;
for (int m = 0; m < cell_nb_nodes(j); ++m) {
sum += kj(cell_nodes(j,m));
@@ -274,8 +270,8 @@ public:
momentum_fluxes += Fl(l)*Cjr(j,R);
energy_fluxes += (Gl(l), Cjr(j,R));
}
- uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant
- //uj[j] += (dt*inv_mj[j]) * momentum_fluxes;
+ //uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant
+ uj[j] += (dt*inv_mj[j]) * momentum_fluxes;
Ej[j] += (dt*inv_mj[j]) * energy_fluxes;
});
@@ -310,8 +306,8 @@ public:
double erreur = 0.;
double exacte = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
- //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
- exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
+ exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
+ //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
erreur += (exacte - uj[j][0])*(exacte - uj[j][0])*Vj(j);
}
erreur = std::sqrt(erreur);
@@ -328,13 +324,13 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
double pi = 4.*std::atan(1.);
- //double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant
- double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant
+ double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant
+ //double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant
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(-2.*pi*pi*0.2);
- exacte = std::sin(pi*xj[j][0])*std::exp(-0.2);
+ exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2);
+ //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2);
if (std::abs(exacte - uj[j][0]) > erreur) {
erreur = std::abs(exacte - uj[j][0]);
@@ -345,6 +341,24 @@ public:
return erreur;
}
+ // Verifie la conservativite de E
+
+ double conservatif(UnknownsType& unknowns) {
+
+ Kokkos::View<double*> Ej = unknowns.Ej();
+
+ const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+
+ double sum = 0.;
+
+ for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
+ sum += Ej[j]*Vj[j];
+ }
+
+ return sum;
+
+ }
+
};
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 5f8e2fcee..cf2d37158 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -204,7 +204,18 @@ void initializeSod()
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_uj[j] = std::sin(pi*xj[j][0]);
+ if (j == m_mesh.numberOfCells()) {
+ m_uj[j] = zero;
+ }
+ else {
+ if (j == 0) {
+ m_uj[j] = zero;
+ }
+ else {
+ m_uj[j] = std::sin(pi*xj[j][0]);
+ }
+ }
+
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -215,7 +226,9 @@ void initializeSod()
block_eos.updateEandCFromRhoP();
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] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
+ m_Ej[j] = 0.5*(m_uj[j],m_uj[j]);
+ //m_Ej[j] = 2.;
});
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
@@ -228,8 +241,8 @@ void initializeSod()
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- //m_kj[j] = 2.; // k constant
- m_kj[j] = xj[j][0]; // k non constant, k = x
+ m_kj[j] = 2.; // k constant
+ //m_kj[j] = xj[j][0]; // k non constant, k = x
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
--
GitLab