diff --git a/src/main.cpp b/src/main.cpp
index 8b357fe5660acb63320be84a7afd891259f7195b..e34d3e45f93f7bdb9262c6922d8759bd85d2fb30 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.2;
+ const double tmax=1.5;
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,7 +291,7 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
- */
+
// AUTRE APPROCHE DU SPLITTING (PLUS LONG)
/*
@@ -310,7 +310,7 @@ int main(int argc, char *argv[])
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
*/
-
+ /*
// DIFFUSION PURE
double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
@@ -319,7 +319,7 @@ int main(int argc, char *argv[])
}
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
-
+ */
block_eos.updatePandCFromRhoE();
@@ -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);
@@ -567,14 +567,15 @@ 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 error = 0.;
error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
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);
@@ -592,7 +593,8 @@ 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";
@@ -606,15 +608,16 @@ int main(int argc, char *argv[])
//method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds();
+
{ // 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.));
- std::ofstream fout("rho");
+ double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
+ std::ofstream fout("rho1600");
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';
}
}
@@ -622,13 +625,13 @@ int main(int argc, char *argv[])
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("u");
+ std::ofstream fout("u1600");
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] << ' ' << 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] << '\n';
}
@@ -637,15 +640,15 @@ int main(int argc, char *argv[])
{ // 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 h = std::sqrt(1. - (tmax*tmax)/(50./9.));
- std::ofstream fout("E");
+ //double pi = 4.*std::atan(1.);
+ double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
+ std::ofstream fout("E1600");
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] << ' ' << ((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] << '\n';
}
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index f9cb525e09d5b8559605cb76385c7185bd73b95e..4370758c2544dc345f1ee1f51b9b0738bcdb5dab 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -197,9 +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
@@ -207,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;
}
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index c074e8e6ef8bd2f73fe6c26fb531bc68d7d4a306..49d24648542a2c2d26f6e9ee11ab566976063f92 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,18 +179,19 @@ 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(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 ;
@@ -217,7 +218,10 @@ private:
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
-
+
+ const Kokkos::View<const Rd*> x0 = m_mesh.x0();
+ const Kokkos::View<const Rd*> xmax = m_mesh.xmax();
+
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rd sum = zero;
double sum2 = 0.;
@@ -236,12 +240,21 @@ private:
// Conditions aux bords
+ // Diffusion pure
+ /*
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));
-
+ */
+
+ // Kidder
+
+ double h = std::sqrt(1. - (t*t)/(50./9.));
+ m_Bl(0) = ((1. + x0[0][0])*3.*x0[0][0])/(100.*h*h*h*h);
+ m_Bl(m_mesh.numberOfFaces()-1) = ((1. + xmax[0][0])*3.*xmax[0][0])/(100.*h*h*h*h);
+
return m_Bl ;
}
@@ -393,13 +406,20 @@ public:
Kokkos::View<double*> nuR = unknowns.nuR();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
-
+
+ const Kokkos::View<const Rd*> x0 = m_mesh.x0();
+ const Kokkos::View<const Rd*> xmax = m_mesh.xmax();
+
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.);
- TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
- // la condition au bord a droite de T depend du temps
+ double h = std::sqrt(1. - (t*t)/(50./9.));
+
+ // Les CL de T dependent du temps
+
+ // Diffusion pure
+ //TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
@@ -416,8 +436,8 @@ public:
// Mise a jour de la vitesse et de l'energie totale specifique
const Kokkos::View<const double*> inv_mj = unknowns.invMj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- //const int j = j0+1;
+ Kokkos::parallel_for(m_mesh.numberOfCells()-2, KOKKOS_LAMBDA(const int& j0) {
+ const int j = j0+1;
Rd momentum_fluxes = zero;
double energy_fluxes = 0.;
Rd trich = zero;
@@ -426,23 +446,22 @@ public:
l = cell_faces(j,R);
momentum_fluxes += Fl(l)*Cjr(j,R);
trich = Bl(l)*Cjr(j,R);
- energy_fluxes += (Gl(l), Cjr(j,R)) + trich[0];
+ energy_fluxes += (Gl(l), Cjr(j,R)) + trich[0];
}
uj[j] += (dt*inv_mj[j]) * momentum_fluxes;
Ej[j] += (dt*inv_mj[j]) * energy_fluxes;
// test k non cst pour diff pure
- uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi);
- Ej[j] -= ((pi*pi*0.5*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) + xj[j][0]*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]))*(std::exp(-2.*t)-1.) - pi*0.5*std::exp(-2.*t)*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]) + (1.+xj[j][0])*((3.*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0])-xj[j][0]*pi*pi*pi*pi*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])))*(std::exp(-2.*t)-1.)+0.5*pi*pi*std::exp(-2.*t)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))))*(dt*inv_mj[j])*Vj(j);
+ //uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi);
+ //Ej[j] -= ((pi*pi*0.5*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) + xj[j][0]*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]))*(std::exp(-2.*t)-1.) - pi*0.5*std::exp(-2.*t)*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]) + (1.+xj[j][0])*((3.*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0])-xj[j][0]*pi*pi*pi*pi*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])))*(std::exp(-2.*t)-1.)+0.5*pi*pi*std::exp(-2.*t)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))))*(dt*inv_mj[j])*Vj(j);
// ajout second membre pour kidder (k = 0.5)
//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));
- //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*(1-t*t/(50/9))*(1-t*t/(50/9))));
-
+ 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))+(6*xj[j][0]+3.)/(100*(1-t*t/(50/9))*(1-t*t/(50/9))));
});
// Calcul de e par la formule e = E-0.5 u^2
@@ -473,8 +492,8 @@ public:
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 = 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
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
}
err_u = std::sqrt(err_u);
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index fdf9c24e74513fc0cc44dcd15b469fe16ef3c163..c9dcf1625b8d31bdfc61cee2acac742ac3735203 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -276,7 +276,7 @@ public:
}
- /*
+
// --- Acoustic Solver ---
void initializeSod()
@@ -284,25 +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.);
@@ -311,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);
@@ -333,9 +335,10 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_kj[j] = xj[j][0];
+
//m_kj[j] = 0.5;
-
+ /*
// Sod
// k non regulier
@@ -376,7 +379,7 @@ public:
int n = 10.;
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];
-
+ */
});
@@ -392,12 +395,26 @@ public:
m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));
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]); // k = x
+ });
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ m_Tj[j] = m_ej[j];
+ });
+
+ // 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.;
+
}
- */
-
+
+ /*
// DIFFUSION PURE
@@ -467,6 +484,7 @@ public:
m_nuR[0] = 1.;
}
+ */
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
: m_mesh_data(mesh_data),