diff --git a/src/main.cpp b/src/main.cpp index 011cf8f724ce2823cb01408456d51b42f39be322..f75fcacee24820644f0a213d5c33c5c661fc0519 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(); @@ -598,8 +598,8 @@ int main(int argc, char *argv[]) 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'; } } @@ -613,10 +613,10 @@ 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'; + fout << xj[j][0] << ' ' << uj[j][0] << '\n'; } } @@ -631,10 +631,10 @@ 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'; + fout << xj[j][0] << ' ' << Ej[j] << '\n'; } } 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 84c9168b13ea5ba1bf8f0da89aefcbe1ec9e8d78..1a02d8e0951dbe2728da0c26ff98888c10d6ab6d 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -182,17 +182,18 @@ 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 - + /* 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; } @@ -371,13 +372,13 @@ public: }); */ - + /* // Mise a jour de nu 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 0c832568e3ede72ae2d007188826f9c134ecaae2..86f7b09a4d0c64afcebbfd860e6038cd9dcf7d8b 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -113,7 +113,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))); @@ -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 @@ -135,8 +135,8 @@ private: //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0]; // k = 0 - m_Fl(0,0) = 0.; - m_Fl(m_mesh.numberOfFaces()-1,0) = 0.; + //m_Fl(0,0) = 0.; + //m_Fl(m_mesh.numberOfFaces()-1,0) = 0.; return m_Fl ; @@ -234,21 +234,21 @@ 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)); - */ + // Kiddder - + /* // nu = (1+x)*0.5 double h = std::sqrt(1. - (t*t)/(50./9.)); m_Bl(0) = ((1.+h*x0[0][0])*3.*h*x0[0][0])/(100.*h*h*h*h); m_Bl(m_mesh.numberOfFaces()-1) = ((1.+h*xmax[0][0])*3.*h*xmax[0][0])/(100.*h*h*h*h); - + */ // nu = 0 //m_Bl(0) = 0.; @@ -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.); @@ -469,7 +469,7 @@ public: //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); // nu = 0 //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)); // nu = 0.2 //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)); // nu = (1+x)*0.5 - Ej[j] -= (dt*inv_mj[j])*Vj(j)*((6*xj[j][0]+3.)/(100*h*h*h*h)); + //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((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 b0e7718bd29aaf8f115f6ab737c8781b18b077c1..bdaa859b44c60ecf0f045d36f9d0674a15a2a019 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);