diff --git a/src/main.cpp b/src/main.cpp
index 8a5237884398a28b365d41d1d6b17720a1a5770a..41953c8c17247d9b0f8b8da6bd035ebd7995215f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -116,7 +116,7 @@ int main(int argc, char *argv[])
// }
- { // class for acoustic solver test
+ {
Kokkos::Timer timer;
timer.reset();
Connectivity1D connectivity(number);
@@ -161,6 +161,7 @@ int main(int argc, char *argv[])
dt=tmax-t;
}
+ std::cout << "t=" << t << " dt=" << dt << '\n';
//acoustic_solver.computeNextStep(t,dt, unknowns);
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 3e6b9eab76fccb565147662323a3b7a91a309f8f..8778e5a4ca4dae785615170a2353bc3e12c23cd5 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -53,7 +53,7 @@ private:
}
KOKKOS_INLINE_FUNCTION
- void _updateVolume()
+ void _updateVolume() // A reflechir pour multiD
{
const Kokkos::View<const unsigned int**>& cell_nodes
= m_mesh.connectivity().cellNodes();
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index f7c0517df6d652778b361c7d305213757cca8d60..9ae88e9fae927ace7b44b77032eacedc537d7d19 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -73,9 +73,9 @@ private:
}
};
- // Calcule un Fjr
- Kokkos::View<Rd**>
- computeFjr(const Kokkos::View<const Rd**>& Cjr,
+ // Calcule un Fl
+ Kokkos::View<Rdd*>
+ computeFl(const Kokkos::View<const Rd**>& Cjr,
const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& kj) {
const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
@@ -85,29 +85,29 @@ private:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
for (int r=0; r<cell_nb_nodes[j]; ++r) {
- m_Fjr(j,r) = ((kj(cell_nodes(j,r)) + kj(cell_nodes(j-1,r)))/(2*Vl(j)))
+ m_Fl(j,r) = ((kj(cell_nodes(j,r)) + kj(cell_nodes(j-1,r)))/(2*Vl(j)))
* ((uj(j,r),Cjr(j,r))+ (uj(j-1,r),Cjr(j-1,r))) ; //tensorProduct(uj(j,r),Cjr(j,r)) ?
}
});
- return m_Fjr;
+ return m_Fl ;
}
- // Calcule un Gjr
- Kokkos::View<double**>
- computeGjr(const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const Rd**>& Fjr) {
+ // Calcule un Gl
+ Kokkos::View<Rd*>
+ computeGl(const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const Rdd*>& Fl) {
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) {
for (int r=0; r<cell_nb_nodes[j]; ++r) {
- m_Gjr(j,r) = 0.5*((uj(cell_nodes(j,r)) + uj(cell_nodes(j-1,r))),Fjr(j,r));
+ m_Gl(j,r) = 0.5*((uj(cell_nodes(j,r)) + uj(cell_nodes(j-1,r))),Fjr(j,r));
}
});
- return m_Gjr;
+ return m_Gl ;
}
// Calcule la liste des inverse d'une liste de matrices
@@ -134,14 +134,14 @@ private:
void computeExplicitFluxes(const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const Rd**>& Cjr,
const Kokkos::View<const double*>& kj) {
- Kokkos::View<Rd**> Fjr = m_Fjr;
- Fjr = computeFjr(Cjr, uj, kj);
- Kokkos::View<double**> Gjr = m_Gjr;
- Gjr = computeGjr(uj, Fjr);
+ Kokkos::View<Rdd*> Fl = m_Fl ;
+ Fl = computeFl (Cjr, uj, kj);
+ Kokkos::View<double**> Gl = m_Gl ;
+ Gl = computeGl (uj, Fl );
}
- Kokkos::View<Rd**> m_Fjr;
- Kokkos::View<double**> m_Gjr;
+ Kokkos::View<Rdd*> m_Fl;
+ Kokkos::View<Rd*> m_Gl;
public:
FiniteVolumesDiffusion(MeshData& mesh_data,
@@ -149,8 +149,8 @@ public:
: m_mesh_data(mesh_data),
m_mesh(mesh_data.mesh()),
m_connectivity(m_mesh.connectivity()),
- m_Fjr("Fjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
- m_Gjr("Gjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell())
+ m_Fl("Fl", m_mesh.numberOfFaces()),
+ m_Gl("Gl", m_mesh.numberOfFaces())
{
;
}
@@ -160,28 +160,33 @@ public:
KOKKOS_INLINE_FUNCTION
double diffusion_dt(const Kokkos::View<const double*>& rhoj,
const Kokkos::View<const double*>& kj) const {
+
Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells());
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+ const Kokkos::View<const unsigned int**>& cell_faces = m_connectivity.cellFaces();
+ const Kokkos::View<const unsigned short*> cell_nb_faces
+ = m_connectivity.cellNbFaces();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./(kj(j+1) + 2*kj(j) + kj(j-1)))
- * std::min(Vl(j),Vl(j-1)); // ATTENTION : kj != 0 !
- // * std::min(xj(j+1)-xj(j), xj(j)-xj(j-1));
- // * std::min((Dxj(j+1),Cjr(j,1)) + (xj(j),Cjr(j,0)),
- //(xj(j),Cjr(j,1)) + (xj(j-1),Cjr(j,0)) );
- //dt_j[j] = 0.0001;
- std::cout << Vl(j) << std::endl;
- //std::cout << std::min(Dj(j),Dj(j-1)) << std::endl;
- // Le probleme vient aussi de Dj
+
+ double minVl = std::numeric_limits<double>::max();
+ for (int ll=0; ll<cell_nb_faces(j); ++ll) {
+ minVl = std::min(minVl, Vl(cell_faces(j, ll)));
+ }
+ // k=1 => (kj(j+1) + 2*kj(j) + kj(j-1)) = 4
+ dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./4)*minVl;
});
-
+ // for (int j=0; j<m_mesh.numberOfCells(); ++j) {
+ // std::cout << "dt_j[" << j << "]=" << dt_j[j] << '\n';
+ // }
double dt = std::numeric_limits<double>::max();
Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(dt_j), dt);
- //std::cout << dt << std::endl;
+ // std::cout << dt << std::endl;
+ // std::exit(0);
return dt;
}
@@ -205,8 +210,8 @@ public:
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj);
- const Kokkos::View<const Rd**> Fjr = m_Fjr;
- const Kokkos::View<const double**> Gjr = m_Gjr;
+ const Kokkos::View<const Rdd*> Fl = m_Fl ;
+ const Kokkos::View<const Rd *> Gl = m_Gl ;
const Kokkos::View<const unsigned short*> cell_nb_nodes
= m_connectivity.cellNbNodes();
@@ -216,8 +221,8 @@ public:
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
for (int R=0; R<cell_nb_nodes[j]; ++R) {
- momentum_fluxes += Fjr(j,R);
- energy_fluxes += Gjr(j,R);
+ momentum_fluxes += Fl (j,R);
+ energy_fluxes += Gl (j,R);
}
uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;