From 180a9c36715fb618f4880b4f6bf61f4d0627e90b Mon Sep 17 00:00:00 2001 From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr> Date: Tue, 29 May 2018 09:48:55 +0200 Subject: [PATCH] correction CFL et position R(t) a la place de x(t) dans Cl de computeUr --- src/main.cpp | 3 +-- src/mesh/Mesh.hpp | 35 ++++++++++++++++++++++++--- src/scheme/AcousticSolver.hpp | 33 ++++++++++++++++++++++--- src/scheme/FiniteVolumesDiffusion.hpp | 11 +++++---- 4 files changed, 68 insertions(+), 14 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 895ea2577..b2ca417a7 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -138,7 +138,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.; + const double tmax=0.7; double t=0; int itermax=std::numeric_limits<int>::max(); @@ -185,7 +185,6 @@ int main(int argc, char *argv[]) } t_diff += dt_diff; } - std::cout << "diff : " << t_diff << std::endl; } diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 462ef51e6..ef3176a03 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -15,6 +15,8 @@ public: private: const Connectivity& m_connectivity; Kokkos::View<Rd*> m_xr; + Kokkos::View<Rd*> m_x0; + Kokkos::View<Rd*> m_xmax; public: @@ -48,16 +50,43 @@ public: return m_xr; } + Kokkos::View<Rd*> x0() + { + return m_x0; + } + + const Kokkos::View<const Rd*> x0() const + { + return m_x0; + } + + Kokkos::View<Rd*> xmax() + { + return m_xmax; + } + + const Kokkos::View<const Rd*> xmax() const + { + return m_xmax; + } + + // pas constant Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), - m_xr("xr", connectivity.numberOfNodes()) + m_xr("xr", connectivity.numberOfNodes()), + m_x0("x0", 1), + m_xmax("xmax", 1) { - const double delta_x = (1.-0.1)/connectivity.numberOfCells(); + double a = 0.1; + double b = 1.; + m_x0[0][0] = a; + m_xmax[0][0] = b; + const double delta_x = (b-a)/connectivity.numberOfCells(); Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){ - m_xr[r][0] = 0.1+r*delta_x; + m_xr[r][0] = a + r*delta_x; }); } diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 78b6b6c62..a4b1250e1 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -129,10 +129,12 @@ private: computeBr(const Kokkos::View<const Rdd**>& Ajr, const Kokkos::View<const Rd**>& Cjr, const Kokkos::View<const Rd*>& uj, - const Kokkos::View<const double*>& pj) { + const Kokkos::View<const double*>& pj, + const double t) { const Kokkos::View<const unsigned int**>& node_cells = m_connectivity.nodeCells(); const Kokkos::View<const unsigned short**>& node_cell_local_node = m_connectivity.nodeCellLocalNode(); const Kokkos::View<const unsigned short*>& node_nb_cells = m_connectivity.nodeNbCells(); + Kokkos::View<Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { Rd& br = m_br(r); @@ -142,6 +144,22 @@ private: const int R = node_cell_local_node(r,j); br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R); } + // if (r==0){ + // for (int j=0; j<node_nb_cells(r); ++j) { + // const int J = node_cells(r,j); + // const int R = node_cell_local_node(r,j); + // double rho0=3./(ht*ht)*(xr(0),xr(0)); + // double p0=2*rho0*rho0*rho0; + // br-=p0*Cjr(J,R); + // } + // } + // if (r==(m_mesh.numberOfNodes()-1)){ + // for (int j=0; j<node_nb_cells(r); ++j) { + // const int J = node_cells(r,j); + // const int R = node_cell_local_node(r,j); + // br-=pmax*Cjr(J,R); + // } + // } }); return m_br; @@ -172,14 +190,21 @@ private: const Kokkos::View<const Rdd*> invAr = m_inv_Ar; Kokkos::View<Rd*> xr = m_mesh.xr(); + Kokkos::View<Rd*> x0 = m_mesh.x0(); + Kokkos::View<Rd*> xmax = m_mesh.xmax(); Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { m_ur[r]=invAr(r)*br(r); }); // Conditions aux limites dependant du temps - m_ur[0]=(-t/((50./9.)-t*t))*xr[0]; - m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1]; + //m_ur[0]=(-t/((50./9.)-t*t))*xr[0]; + //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; } @@ -236,7 +261,7 @@ private: const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr); const Kokkos::View<const Rdd*> Ar = computeAr(Ajr); - const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj); + const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj, t); Kokkos::View<Rd*> ur = m_ur; Kokkos::View<Rd**> Fjr = m_Fjr; diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 434af08ac..8d4c27674 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -124,13 +124,13 @@ private: //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))); */ - // 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]; - // k = 0.5 //m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5; //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5; + + // 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]; return m_Fl ; } @@ -264,7 +264,8 @@ public: sum += kj(cell_nodes(j,m)); } - dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl; + // dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl; + dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl; }); -- GitLab