Skip to content
Snippets Groups Projects
Commit cbba1f4f authored by Fanny CHOPOT's avatar Fanny CHOPOT
Browse files

bonne version pour etude convergence kidder

parent 328825cd
Branches
Tags Kidder
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/main.cpp
+ 9
− 9
View file @ cbba1f4f
...@@ -141,7 +141,7 @@ int main(int argc, char *argv[]) ...@@ -141,7 +141,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const double*> Vj = mesh_data.Vj(); const Kokkos::View<const double*> Vj = mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.7; const double tmax=1.5;
double t=0.; double t=0.;
int itermax=std::numeric_limits<int>::max(); int itermax=std::numeric_limits<int>::max();
...@@ -280,7 +280,7 @@ int main(int argc, char *argv[]) ...@@ -280,7 +280,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
/*
double error1 = 0.; double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax); error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
...@@ -304,7 +304,7 @@ int main(int argc, char *argv[]) ...@@ -304,7 +304,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
<< ": " << rang::fgB::green << error4 << rang::fg::reset << " \n"; << ": " << rang::fgB::green << error4 << rang::fg::reset << " \n";
*/
/* /*
double error3 = 0.; double error3 = 0.;
error3 = finite_volumes_diffusion.error_L2_E(unknowns); error3 = finite_volumes_diffusion.error_L2_E(unknowns);
...@@ -329,12 +329,12 @@ int main(int argc, char *argv[]) ...@@ -329,12 +329,12 @@ int main(int argc, char *argv[])
{ // gnuplot output for density { // gnuplot output for density
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj(); 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"); std::ofstream fout("rho800");
fout.precision(15); fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) { 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] << ' ' << 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] << '\n';
} }
} }
...@@ -347,8 +347,8 @@ int main(int argc, char *argv[]) ...@@ -347,8 +347,8 @@ int main(int argc, char *argv[])
for (size_t j=0; j<mesh.numberOfCells(); ++j) { 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(-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] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\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';
} }
} }
......
src/mesh/Mesh.hpp
+ 2
− 2
View file @ cbba1f4f
...@@ -80,8 +80,8 @@ public: ...@@ -80,8 +80,8 @@ public:
m_x0("x0", 1), m_x0("x0", 1),
m_xmax("xmax", 1) m_xmax("xmax", 1)
{ {
double a = -1.; double a = 0.;
double b = 2.; double b = 1.;
m_x0[0][0] = a; m_x0[0][0] = a;
m_xmax[0][0] = b; m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells(); const double delta_x = (b-a)/connectivity.numberOfCells();
......
src/scheme/AcousticSolver.hpp
+ 7
− 7
View file @ cbba1f4f
...@@ -197,8 +197,8 @@ private: ...@@ -197,8 +197,8 @@ private:
m_ur[r]=invAr(r)*br(r); m_ur[r]=invAr(r)*br(r);
}); });
m_ur[0]=zero; //m_ur[0]=zero;
m_ur[m_mesh.numberOfNodes()-1]=zero; //m_ur[m_mesh.numberOfNodes()-1]=zero;
// Kidder // Kidder
...@@ -207,11 +207,11 @@ private: ...@@ -207,11 +207,11 @@ private:
//m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1]; //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) //R(t) = x*h(t) a la place de x(t)
/*
double h = std::sqrt(1. - (t*t)/(50./9.)); double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0]; 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]; m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
*/
return m_ur; return m_ur;
} }
...@@ -386,9 +386,9 @@ public: ...@@ -386,9 +386,9 @@ public:
// Mise a jour de k // Mise a jour de k
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
//kj(j) = xj[j][0]; kj(j) = xj[j][0];
//kj(j) = 0.5; //kj(j) = 0.5;
/*
if (xj[j][0]<0.7) { if (xj[j][0]<0.7) {
kj[j]=0.; kj[j]=0.;
} else { } else {
...@@ -398,7 +398,7 @@ public: ...@@ -398,7 +398,7 @@ public:
kj[j]=0. ; kj[j]=0. ;
} }
} }
*/
}); });
} }
......
src/scheme/FiniteVolumesDiffusion.hpp
+ 12
− 11
View file @ cbba1f4f
...@@ -115,7 +115,7 @@ private: ...@@ -115,7 +115,7 @@ private:
}); });
// Conditions aux bords // Conditions aux bords
/*
int cell_here = face_cells(0,0); int cell_here = face_cells(0,0);
int local_face_number_in_cell = face_cell_local_face(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))); 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: ...@@ -124,7 +124,7 @@ private:
local_face_number_in_cell = face_cell_local_face(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))); 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))); //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 // Kidder
...@@ -135,11 +135,11 @@ private: ...@@ -135,11 +135,11 @@ private:
// k = x // k = x
//m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0]; //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]; //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.)); double h = std::sqrt(1. - (t*t)/(50./9.));
m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0]; 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]; m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
*/
return m_Fl ; return m_Fl ;
} }
...@@ -178,18 +178,18 @@ private: ...@@ -178,18 +178,18 @@ private:
}); });
// Conditions aux bords // Conditions aux bords
m_Gl(0) = Fl(0)*uL(0); //m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0); //m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
// Kidder // 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); //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.)); 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(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); m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
*/
return m_Gl ; return m_Gl ;
...@@ -332,7 +332,8 @@ public: ...@@ -332,7 +332,8 @@ public:
// Mise a jour de la vitesse et de l'energie totale specifique // Mise a jour de la vitesse et de l'energie totale specifique
const Kokkos::View<const double*> inv_mj = unknowns.invMj(); const Kokkos::View<const double*> inv_mj = unknowns.invMj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { Kokkos::parallel_for(m_mesh.numberOfCells()-2, KOKKOS_LAMBDA(const int& j0) {
const int j = j0+1;
Rd momentum_fluxes = zero; Rd momentum_fluxes = zero;
double energy_fluxes = 0.; double energy_fluxes = 0.;
int l = 0; int l = 0;
...@@ -350,8 +351,8 @@ public: ...@@ -350,8 +351,8 @@ public:
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); //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) // 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))); Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
}); });
......
src/scheme/FiniteVolumesEulerUnknowns.hpp
+ 12
− 12
View file @ cbba1f4f
...@@ -223,27 +223,27 @@ public: ...@@ -223,27 +223,27 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
/*
if (xj[j][0]<0.5) { if (xj[j][0]<0.5) {
m_rhoj[j]=1.; m_rhoj[j]=1.;
} else { } else {
m_rhoj[j]=0.125; m_rhoj[j]=0.125;
} }
*/
//Kidder //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){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
/*
if (xj[j][0]<0.5) { if (xj[j][0]<0.5) {
m_pj[j]=1; m_pj[j]=1;
} else { } else {
m_pj[j]=0.1; m_pj[j]=0.1;
} }
*/
//Kidder //Kidder
//m_pj[j] = 2.*std::pow(m_rhoj[j],3); m_pj[j] = 2.*std::pow(m_rhoj[j],3);
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
...@@ -251,8 +251,8 @@ public: ...@@ -251,8 +251,8 @@ public:
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_gammaj[j] = 1.4; //m_gammaj[j] = 1.4;
//m_gammaj[j] = 3.; m_gammaj[j] = 3.;
}); });
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
...@@ -272,9 +272,9 @@ public: ...@@ -272,9 +272,9 @@ public:
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ 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; //m_kj[j] = 0.5;
/*
if (xj[j][0]<0.7) { if (xj[j][0]<0.7) {
m_kj[j]=0.; m_kj[j]=0.;
} else { } else {
...@@ -284,7 +284,7 @@ public: ...@@ -284,7 +284,7 @@ public:
m_kj[j]=0. ; m_kj[j]=0. ;
} }
} }
*/
}); });
// Conditions aux bords de Dirichlet sur u et k // Conditions aux bords de Dirichlet sur u et k
...@@ -292,7 +292,7 @@ public: ...@@ -292,7 +292,7 @@ public:
m_uL[0] = zero; m_uL[0] = zero;
m_uR[0] = zero; m_uR[0] = zero;
m_kL[0] = 0.; m_kL[0] = 0.;
m_kR[0] = 0.; m_kR[0] = 1.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j])); m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment