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

ajout terme source Sj, cas k non constant, non OK

parent d6141b90
No related branches found
No related tags found
No related merge requests found
...@@ -184,7 +184,8 @@ int main(int argc, char *argv[]) ...@@ -184,7 +184,8 @@ int main(int argc, char *argv[])
double pi = 4.*std::atan(1.); double pi = 4.*std::atan(1.);
std::ofstream fout("comparaison u"); std::ofstream fout("comparaison u");
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'; // 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
} }
} }
......
...@@ -91,18 +91,21 @@ private: ...@@ -91,18 +91,21 @@ private:
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) { Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rdd sum = zero; Rdd sum = zero;
double sum2 = 0.;
for (int j=0; j<face_nb_cells(l); ++j) { for (int j=0; j<face_nb_cells(l); ++j) {
int cell_here = face_cells(l,j); int cell_here = face_cells(l,j);
int local_face_number_in_cell = face_cell_local_face(l,j); int local_face_number_in_cell = face_cell_local_face(l,j);
sum += tensorProduct(uj(cell_here, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell)); sum += tensorProduct(uj(cell_here, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell));
sum2 += kj(cell_here);
} }
//Fl(l) = ((kj(cell_nodes(j,r)) + kj(cell_nodes(j-1,r)))/(2*Vl(j))) // k = x
//* ((uj(j,r),Cjr(j,r))+ (uj(j-1,r),Cjr(j-1,r))) ; //tensorProduct(uj(j,r),Cjr(j,r)) ?
m_Fl(l) = ((sum2*0.5)/Vl(l))*sum;
// k = 2 // k = 2
m_Fl(l)= (2./Vl(l))*sum; // m_Fl(l)= (2./Vl(l))*sum;
}); });
...@@ -201,13 +204,27 @@ public: ...@@ -201,13 +204,27 @@ public:
const Kokkos::View<const unsigned short*> cell_nb_faces const Kokkos::View<const unsigned short*> cell_nb_faces
= m_connectivity.cellNbFaces(); = m_connectivity.cellNbFaces();
const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
const Kokkos::View<const unsigned short*> cell_nb_nodes
= m_connectivity.cellNbNodes();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
double minVl = std::numeric_limits<double>::max(); double minVl = std::numeric_limits<double>::max();
for (int ll=0; ll<cell_nb_faces(j); ++ll) { for (int ll=0; ll<cell_nb_faces(j); ++ll) {
minVl = std::min(minVl, Vl(cell_faces(j, ll))); minVl = std::min(minVl, Vl(cell_faces(j, ll)));
} }
// k=2 => (kj(j+1) + 2*kj(j) + kj(j-1)) = 8 // k=2 => (kj(j+1) + 2*kj(j) + kj(j-1)) = 8
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./8.)*minVl; // dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./8.)*minVl;
// k=x
double sum = 0.;
for (int m = 0; m < cell_nb_nodes(j); ++m) {
sum += kj(cell_nodes(j,m));
}
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl;
}); });
// for (int j=0; j<m_mesh.numberOfCells(); ++j) { // for (int j=0; j<m_mesh.numberOfCells(); ++j) {
...@@ -227,6 +244,7 @@ public: ...@@ -227,6 +244,7 @@ public:
{ {
Kokkos::View<double*> rhoj = unknowns.rhoj(); Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<Rd*> uj = unknowns.uj(); Kokkos::View<Rd*> uj = unknowns.uj();
Kokkos::View<Rd*> Sj = unknowns.Sj();
Kokkos::View<double*> Ej = unknowns.Ej(); Kokkos::View<double*> Ej = unknowns.Ej();
Kokkos::View<double*> ej = unknowns.ej(); Kokkos::View<double*> ej = unknowns.ej();
...@@ -261,7 +279,7 @@ public: ...@@ -261,7 +279,7 @@ public:
momentum_fluxes += Fl(l)*Cjr(j,R); momentum_fluxes += Fl(l)*Cjr(j,R);
energy_fluxes += (Gl(l), Cjr(j,R)); energy_fluxes += (Gl(l), Cjr(j,R));
} }
uj[j] -= (dt*inv_mj[j]) * momentum_fluxes; uj[j] -= std::exp(-t)*dt*Vj(j)*Sj(j)+(dt*inv_mj[j]) * momentum_fluxes;
Ej[j] -= (dt*inv_mj[j]) * energy_fluxes; Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
}); });
......
...@@ -26,6 +26,7 @@ private: ...@@ -26,6 +26,7 @@ private:
Kokkos::View<double*> m_mj; Kokkos::View<double*> m_mj;
Kokkos::View<double*> m_inv_mj; Kokkos::View<double*> m_inv_mj;
Kokkos::View<double*> m_kj; Kokkos::View<double*> m_kj;
Kokkos::View<Rd*> m_Sj;
public: public:
Kokkos::View<double*> rhoj() Kokkos::View<double*> rhoj()
...@@ -128,6 +129,11 @@ public: ...@@ -128,6 +129,11 @@ public:
return m_kj; return m_kj;
} }
Kokkos::View<Rd*> Sj()
{
return m_Sj;
}
// --- Acoustic Solver --- // --- Acoustic Solver ---
// void initializeSod() // void initializeSod()
...@@ -222,7 +228,12 @@ void initializeSod() ...@@ -222,7 +228,12 @@ void initializeSod()
}); });
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] = 2; //m_kj[j] = 2; // k constant
m_kj[j] = xj[j][0]; // k non constant, k = x
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_Sj[j] = std::sin(pi*xj[j][0])*(xj[j][0]-1.) - std::cos(xj[j][0]*pi);
}); });
} }
...@@ -239,7 +250,8 @@ void initializeSod() ...@@ -239,7 +250,8 @@ void initializeSod()
m_cj("cj",m_mesh.numberOfCells()), m_cj("cj",m_mesh.numberOfCells()),
m_mj("mj",m_mesh.numberOfCells()), m_mj("mj",m_mesh.numberOfCells()),
m_inv_mj("inv_mj",m_mesh.numberOfCells()), m_inv_mj("inv_mj",m_mesh.numberOfCells()),
m_kj("kj",m_mesh.numberOfCells()) m_kj("kj",m_mesh.numberOfCells()),
m_Sj("Sj",m_mesh.numberOfCells())
{ {
; ;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment