diff --git a/src/main.cpp b/src/main.cpp
index 2ad92a3662695c980356ef8c8a5a286e69466eff..c1a54433b47f96dbb07f386aa08c139fade277bb 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=1.5;
+ const double tmax=0.8;
double t=0.;
int itermax=std::numeric_limits<int>::max();
@@ -265,7 +265,7 @@ int main(int argc, char *argv[])
// 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;
@@ -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,nuj, cj);
+ double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj);
double t_diff = t-dt_euler;
if (dt_euler <= dt_diff) {
@@ -555,7 +555,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
+ /*
double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
@@ -575,7 +575,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
<< ": " << rang::fgB::green << error << rang::fg::reset << " \n";
- /*
+
double error4 = 0.;
error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);
@@ -612,28 +612,28 @@ int main(int argc, char *argv[])
{ // gnuplot output for density
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.));
+ //double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("rho");
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';
}
}
{ // 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(-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] << '\n';
+ fout << xj[j][0] << ' ' << uj[j][0] << '\n';
}
}
@@ -641,19 +641,20 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> Ej = unknowns.Ej();
//double pi = 4.*std::atan(1.);
- double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
+ //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] << ' ' << ((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] << '\n';
+ fout << xj[j][0] << ' ' << Ej[j] << '\n';
}
}
+ /*
{ // gnuplot output u^2*0.5 + T
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> Tj = unknowns.Tj();
@@ -665,7 +666,7 @@ int main(int argc, char *argv[])
fout << xj[j][0] << ' ' << Tj[j]+uj[j][0]*uj[j][0]*0.5 << ' ' << (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
}
}
-
+ */
/*
{ // gnuplot output for entropy (gaz parfait en prenant cv = 1))
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 1f26bcac1ae9559bfe59af8a5b34be4904490654..470236c9e073d2670ebd4c2af54d74619bbc62e5 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 = 0.;
- double b = 1.;
+ double a = -1.;
+ double b = 2.;
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 4370758c2544dc345f1ee1f51b9b0738bcdb5dab..309e5165404ce018961698bb6ed5fec367fbbecc 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -197,10 +197,10 @@ private:
m_ur[r]=invAr(r)*br(r);
});
- /*
+
m_ur[0]=zero;
m_ur[m_mesh.numberOfNodes()-1]=zero;
- */
+
// Kidder
// Conditions aux limites dependant du temps
@@ -208,11 +208,11 @@ private:
//m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
//R(t) = x*h(t) a la place de x(t)
-
+ /*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
-
+ */
return m_ur;
}
@@ -386,11 +386,11 @@ public:
});
// Mise a jour de k
-
+ /*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
kj(j) = xj[j][0];
- /*
+
//kj(j) = 0.5;
if (xj[j][0]<0.7) {
@@ -402,7 +402,7 @@ public:
kj[j]=0. ;
}
}
- */
+
});
@@ -411,7 +411,7 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
nuj(j) = 0.5*(1.+xj[j][0]);
});
-
+ */
}
};
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index d22213b63462a1773846327bc235bf96b5efbd08..7d7b9a825cb65fa5dd7b5f6d4dbadb2c8779fa55 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -115,7 +115,7 @@ private:
});
// Conditions aux bords
- /*
+
int cell_here = face_cells(0,0);
int local_face_number_in_cell = face_cell_local_face(0,0);
m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
@@ -124,7 +124,7 @@ private:
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)));
//m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
- */
+
// Kidder
@@ -135,11 +135,11 @@ private:
// k = x
//m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
//m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
-
+ /*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0];
m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
-
+ */
return m_Fl ;
}
@@ -179,19 +179,20 @@ private:
});
// Conditions aux bords
- //m_Gl(0) = Fl(0)*uL(0);
- //m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
+ m_Gl(0) = Fl(0)*uL(0);
+ m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
// Kidder
// m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
//m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
-
+ /*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0);
m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
-
+ */
+
return m_Gl ;
}
@@ -238,14 +239,14 @@ private:
});
// Conditions aux bords
- /*
+
int cell_here = face_cells(0,0);
m_Bl(0) = (nuL(0) + nuj(cell_here))*(1./(2*Vl(0)))*(Tj(cell_here) - TL(0));
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));
- */
+ /*
double h = std::sqrt(1. - (t*t)/(50./9.));
// nu = (1+x)*0.5
@@ -255,7 +256,8 @@ private:
// nu = 0.2
//m_Bl(0) = (0.2*3.*h*x0[0][0])/(50.*h*h*h*h);
//m_Bl(m_mesh.numberOfFaces()-1) = (0.2*3.*h*xmax[0][0])/(50.*h*h*h*h);
-
+ */
+
return m_Bl ;
}
@@ -359,14 +361,14 @@ public:
sum += kj(cell_nodes(j,m));
sum1 += nuj(cell_nodes(j,m));
}
-
+
if (sum > sum1) {
if (sum == 0.) {
dt_j[j] = std::numeric_limits<double>::max();
} else {
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl;
}
- } else {
+ } else {
if (sum1 == 0.) {
dt_j[j] = std::numeric_limits<double>::max();
} else {
@@ -415,25 +417,11 @@ public:
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
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.));
-
- // Les CL de T dependent du temps
+ // double pi = 4.*std::atan(1.);
+ // double h = std::sqrt(1. - (t*t)/(50./9.));
// Diffusion pure
//TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
-
- // Kidder
- /*
- TL(0) = (1./(100*h*h))*((3.*h*x0[0][0]*h*x0[0][0])/(h*h) + 100.);
- TR(0) = (1./(100*h*h))*((3.*h*xmax[0][0]*h*xmax[0][0])/(h*h) + 100.);
- nuL(0) = (h*x0[0][0]+1.)*0.5;
- nuR(0) = (h*xmax[0][0]+1.)*0.5;
- uL[0] = (-h*x0[0][0]*t)/((50./9.)-t*t);
- uR[0] = (-h*xmax[0][0]*t)/((50./9.)-t*t);
- kL[0] = h*x0[0][0];
- kR[0] = h*xmax[0][0] ;
- */
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
@@ -474,9 +462,9 @@ public:
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
// ajout second membre pour kidder (k = x)
- uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t));
+ //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))+(0.2*3.)/(50.*h*h*h*h));
- 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*h*h*h*h));
+ //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*h*h*h*h));
});
// Calcul de e par la formule e = E-0.5 u^2
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 0c8b9c01c22638530ae263a175285ea132794fe8..8d247083d757b56a0420dabea004b9a7fff4b9eb 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -284,27 +284,27 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- /*
+
if (xj[j][0]<0.5) {
m_rhoj[j]=1.;
} else {
m_rhoj[j]=0.125;
}
- */
+
//Kidder
- m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
+ //m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- /*
+
if (xj[j][0]<0.5) {
m_pj[j]=1;
} else {
m_pj[j]=0.1;
}
- */
+
//Kidder
- m_pj[j] = 2.*std::pow(m_rhoj[j],3);
+ //m_pj[j] = 2.*std::pow(m_rhoj[j],3);
});
double pi = 4.*std::atan(1.);
@@ -313,8 +313,8 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- //m_gammaj[j] = 1.4;
- m_gammaj[j] = 3.;
+ m_gammaj[j] = 1.4;
+ //m_gammaj[j] = 3.;
});
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
@@ -334,15 +334,15 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_kj[j] = xj[j][0];
+ //m_kj[j] = xj[j][0];
//m_kj[j] = 0.5;
- /*
+
// 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 = 10.;
+ 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];
- */
+
});
@@ -388,7 +388,7 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
m_kL[0] = 0.;
- m_kR[0] = 1.;
+ m_kR[0] = 0.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -396,8 +396,8 @@ 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.2;
+ //m_nuj(j) = 0.5*(1.+xj[j][0]);
+ m_nuj(j) = 0.5;
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -406,10 +406,10 @@ public:
// Conditions aux bords de Dirichlet sur T et nu
- m_TL[0] = 1.;
- m_TR[0] = 103./100.;
- 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] = m_nuj[0];
+ m_nuR[0] = m_nuj[m_mesh.numberOfCells()-1];
}