diff --git a/src/main.cpp b/src/main.cpp
index c8b67fadecce48eed6e8569cce9aa754dcc7236b..afe1896fdce714b0ffcc25a61acf3e05e5f4887e 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();
@@ -262,7 +262,7 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
-
+ /*
// ETAPE 1 DU SPLITTING - EULER
double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
@@ -291,28 +291,28 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
+ */
- /*
// DIFFUSION PURE
- double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
+ double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
if (t+dt > tmax) {
dt = tmax-t;
}
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
- */
+
block_eos.updatePandCFromRhoE();
++iteration;
std::cout << "temps t : " << t << std::endl;
- /*
- // ECRITURE DANS UN FICHIER
-
- if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
+ // ECRITURE DANS UN FICHIER
+ /*
+ //if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
+
std::string ligne;
// rho
@@ -595,11 +595,11 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj();
double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
- std::ofstream fout("rho");
+ std::ofstream fout("rho1");
fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder
- //fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+ //fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder
+ fout << xj[j][0] << ' ' << rhoj[j] << '\n';
}
}
@@ -613,8 +613,9 @@ int main(int argc, char *argv[])
//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]*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';
}
}
@@ -630,13 +631,14 @@ int main(int argc, char *argv[])
//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] << ' ' << (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';
}
}
-
+ /*
{ // gnuplot output for entropy (gaz parfait en prenant cv = 1))
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj();
@@ -649,8 +651,8 @@ int main(int argc, char *argv[])
fout << xj[j][0] << ' ' << std::log(pj[j]*std::pow(rhoj[j],-gammaj[j])) << '\n';
}
}
-
-
+ */
+ /*
{ // gnuplot output for k
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> kj = unknowns.kj();
@@ -664,6 +666,8 @@ int main(int argc, char *argv[])
}
}
}
+ */
+
}
Kokkos::finalize();
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 470236c9e073d2670ebd4c2af54d74619bbc62e5..1f26bcac1ae9559bfe59af8a5b34be4904490654 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -80,8 +80,8 @@ public:
m_x0("x0", 1),
m_xmax("xmax", 1)
{
- double a = -1.;
- double b = 2.;
+ double a = 0.;
+ double b = 1.;
m_x0[0][0] = a;
m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells();
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 8232fbecd4a08740cfc347243190cbeda9de8065..e4ed0d0c837604c1de504584a80d5b617041ad33 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -182,8 +182,10 @@ private:
});
- m_ur[0]=zero;
- m_ur[m_mesh.numberOfNodes()-1]=zero;
+ //m_ur[0]=zero;
+ //m_ur[m_mesh.numberOfNodes()-1]=zero;
+ m_ur[0] = x0;
+ m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
// CL Kidder
/*
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 0147c3e7439eeb3d45fe3f2dd92576371fa7e38c..b2a714eccbce55886f990008273d949f283461c6 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -121,7 +121,7 @@ private:
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)));
-
+
// Kidder
@@ -236,7 +236,6 @@ 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));
-
// Kiddder
@@ -419,6 +418,8 @@ public:
// CL en diffusion pure
// TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
+ TL(0) = 1.+ t;
+ TR(0) = 3. + t;
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 84ce1bbb347bb69310d052b1b6f343ba09fd6405..2e137b352efedbea681a1fbdcf8c2b72c0cd5921 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -276,7 +276,7 @@ public:
}
-
+ /*
// --- Acoustic Solver ---
void initializeSod()
@@ -336,13 +336,13 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_kj[j] = xj[j][0];
- m_kj[j] = 0.2;
+ //m_kj[j] = 0.2;
// Sod
// k non regulier
- /*
+
if (xj[j][0]<0.7) {
m_kj[j]=0.;
} else {
@@ -373,13 +373,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];
- */
+
});
@@ -397,7 +397,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.;
+ m_nuj(j) = 0.5;
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -408,14 +408,14 @@ public:
m_TL[0] = m_ej[0];
m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
- m_nuL[0] = 0.;
- m_nuR[0] = 0.;
+ m_nuL[0] = 0.5;
+ m_nuR[0] = 0.5;
}
+ */
-
- /*
+
// DIFFUSION PURE
@@ -433,7 +433,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] = std::sin(pi*xj[j][0]);
+ m_uj[j] = xj[j][0];
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -445,7 +446,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] = 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();
@@ -458,16 +460,17 @@ public:
});
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] = 0.; // 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] = zero;
+ m_uR[0] = xj[0];
m_kL[0] = 0.;
- m_kR[0] = 1.;
+ m_kR[0] = 0.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
@@ -479,13 +482,13 @@ public:
// Conditions aux bords de Dirichlet sur T et nu
- m_TL[0] = 2.;
- m_TR[0] = 2.;
+ m_TL[0] = 1.;
+ m_TR[0] = 3.;
m_nuL[0] = 0.5;
m_nuR[0] = 1.;
}
- */
+
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
: m_mesh_data(mesh_data),