diff --git a/src/main.cpp b/src/main.cpp
index f75fcacee24820644f0a213d5c33c5c661fc0519..51f9e385bcbd32ce2bc98a89dc1c7e6c4326a399 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -160,6 +160,12 @@ int main(int argc, char *argv[])
double c = 0.;
c = finite_volumes_diffusion.conservatif(unknowns);
+
+ const Kokkos::View<const Rd*> xj = mesh_data.xj();
+ int size = 3000;
+ std::vector<std::vector<double>> x(size, std::vector<double>(mesh.numberOfCells()));
+ std::vector<double> tempo(size);
+ int i = 0;
/*
// Ecriture des valeurs initiales dans un fichier
@@ -260,12 +266,13 @@ int main(int argc, char *argv[])
diff.close();
*/
+
while((t<tmax) and (iteration<itermax)) {
// ETAPE 1 DU SPLITTING - EULER
- double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
@@ -307,11 +314,11 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
-
-
+
+ /*
// ECRITURE DANS UN FICHIER
- /*
- //if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
+
+ if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
std::string ligne;
@@ -500,7 +507,7 @@ int main(int argc, char *argv[])
tempo.precision(5);
tempo << std::fixed << t << '\n';
tempo.close();
-
+
// Fichier k indicateur
std::ifstream diffint("diffinter");
std::ofstream diffout("k", std::ios::trunc);
@@ -525,17 +532,54 @@ int main(int argc, char *argv[])
}
riffint.close();
riffout.close();
-
+
}
*/
// 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;
+
}
+ std::cout << "i = " << i << std::endl;
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+ // CREATION FICHIER x = f(t)
+ /*
+ std::ofstream fout("cara");
+ fout.precision(15);
+ for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ for (size_t k=0; k<i; ++k) {
+ fout << tempo[k] << ' ' << x[k][j] << '\n';
+ }
+ fout << ' ' << '\n';
+ fout << ' ' << '\n';
+ }
+ */
+ std::ofstream fout("cara1");
+ fout.precision(15);
+ for (int j=0; j<mesh.numberOfCells(); ++j) {
+ if (j%10 == 0) {
+ for (int k=0; k<i; ++k) {
+ if ( (k%10 == 0) or (k == i-1) ) {
+ fout << tempo[k] << ' ' << x[k][j] << '\n';
+ }
+ }
+ fout << ' ' << '\n';
+ fout << ' ' << '\n';
+ }
+ }
+
/*
double error1 = 0.;
@@ -668,6 +712,21 @@ 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();
+ const Kokkos::View<const double*> pj = unknowns.pj();
+ const Kokkos::View<const double*> gammaj = unknowns.gammaj();
+ const Kokkos::View<const Rd*> uj = unknowns.uj();
+ std::ofstream fout("vitesse");
+ fout.precision(15);
+ for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << uj[j][0] + std::sqrt((gammaj[j]*pj[j])/rhoj[j]) << '\n';
+ }
+ }
+
+
+
}
Kokkos::finalize();
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 86f7b09a4d0c64afcebbfd860e6038cd9dcf7d8b..daf8fd64e249fbd40664ed534bdde482fd16060d 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -419,7 +419,7 @@ public:
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
// double pi = 4.*std::atan(1.);
- //double h = std::sqrt(1. - (t*t)/(50./9.));
+ double h = std::sqrt(1. - (t*t)/(50./9.));
// CL en diffusion pure
// TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
@@ -502,6 +502,7 @@ public:
//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);
@@ -541,7 +542,8 @@ 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.*t)-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 += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j);
}
err_E = std::sqrt(err_E);
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index bdaa859b44c60ecf0f045d36f9d0674a15a2a019..88a9500e781f68c584b16f538f2082cc26fdc35d 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -274,8 +274,6 @@ public:
{
return m_S0;
}
-
-
// --- Acoustic Solver ---
@@ -336,7 +334,7 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_kj[j] = xj[j][0];
- m_kj[j] = 0.;
+ //m_kj[j] = 0.;
/*
// Sod
@@ -373,13 +371,13 @@ public:
m_kj[j]=0. ;
}
}
-
+ */
// k regulier
int n = 1;
m_kj[j] = std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.7)*(xj[j][0]<0.7+0.1/n) + std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.9-0.1/n)*(xj[j][0]<0.9) + (xj[j][0]>0.7+0.1/n)*(xj[j][0]<0.9-0.1/n);
m_kj[j] = 0.014*m_kj[j];
- */
+
});
@@ -396,8 +394,19 @@ public:
m_S0(j) = m_entropy(j);
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_nuj(j) = 0.5*(1.+xj[j][0]);
- //m_nuj(j) = 0.5;
+ //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 {
+ if (xj[j][0]<0.9){
+ m_nuj[j]=0.5;
+ } else {
+ m_nuj[j]=0. ;
+ }
+ }
+ */
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -408,8 +417,8 @@ public:
m_TL[0] = m_ej[0];
m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
- m_nuL[0] = 0.5;
- m_nuR[0] = 1.;
+ m_nuL[0] = 0.0005;
+ m_nuR[0] = 0.0005;
}
@@ -514,7 +523,7 @@ public:
m_nuR("nuR", 1),
m_entropy("entropy", m_mesh.numberOfCells()),
m_entropy_new("entropy_new", m_mesh.numberOfCells()),
- m_S0("S0", m_mesh.numberOfCells())
+ m_S0("S0", m_mesh.numberOfCells())
{
;
}