diff --git a/src/main.cpp b/src/main.cpp
index 2ab1f8d09538c6c6b87086930b462f9cfa57666b..f31a11b9a5c4319a70a571ba6072f231dde0bf8f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -144,7 +144,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();
@@ -167,6 +167,24 @@ int main(int argc, char *argv[])
double c = 0.;
c = finite_volumes_diffusion.conservatif(unknowns);
+
+ /*
+ // xj
+ const Kokkos::View<const Rd*> xj = mesh_data.xj();
+ std::ofstream fout("xj");
+ fout.precision(15);
+ for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << std::endl;
+ }
+
+ // xr
+ Kokkos::View<Rd*> xr = mesh.xr();
+ std::ofstream fout1("xr");
+ fout1.precision(15);
+ for (size_t j=0; j<mesh.numberOfNodes(); ++j) {
+ fout1 << xr[j][0] << ' ' << std::endl;
+ }
+ */
// Diagrammes de marche
/*
@@ -276,15 +294,14 @@ int main(int argc, char *argv[])
}
}
diff.close();
- */
-
+ */
while((t<tmax) and (iteration<itermax)) {
// NAVIER-STOKES AVEC SPLITTING
/*
// Etape 1 du splitting - Euler
- double dt_euler = 0.5*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt_euler = 0.1*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
}
@@ -321,14 +338,20 @@ int main(int argc, char *argv[])
// NAVIER-STOKES SANS SPLITTING
- double dt = 0.5*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+ double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
if (t+dt > tmax) {
dt = tmax-t;
}
no_splitting.computeNextStep(t,dt, unknowns);
t += dt;
-
-
+ /*
+ // Fichier pas de temps
+ std::ofstream past("pas_no_split", std::ios::app);
+ past.precision(5);
+ //past << std::fixed << dt_euler << std::endl; // split
+ past << std::fixed << dt << std::endl; // no split
+ past.close();
+ */
block_eos.updatePandCFromRhoE();
++iteration;
@@ -659,7 +682,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj();
double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
- std::ofstream fout("rho_3");
+ 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
@@ -676,7 +699,17 @@ int main(int argc, char *argv[])
fout << xj[j][0] << ' ' << pj[j] << '\n';
}
}
-
+
+ { // gnuplot output for pression total (no split)
+ const Kokkos::View<const Rd*> xj = mesh_data.xj();
+ const Kokkos::View<const double*> PTj = unknowns.PTj();
+ std::ofstream fout("pT");
+ fout.precision(15);
+ for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << PTj[j] << '\n';
+ }
+ }
+
{ // gnuplot output for internal energy
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> ej = unknowns.ej();
@@ -691,7 +724,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const Rd*> uj = unknowns.uj();
double pi = 4.*std::atan(1.);
- std::ofstream fout("u_3");
+ std::ofstream fout("u");
fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
@@ -715,10 +748,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 7440cdc15cd916932eb23549675990fa14b3a292..d5a608a3a2133e47a2cfed91e743d9d228b3b06c 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -145,9 +145,13 @@ public:
for (int r=1; r<connectivity.numberOfNodes()-1; ++r){
const double delta_xr = dist(mt);
- m_xr[r][0] = a + r*h + delta_xr;
+ //if (r< connectivity.numberOfCells()/2) {
+ m_xr[r][0] = a + r*h + delta_xr;
+ //} else {
+ //m_xr[r][0] = a + r*h;
+ //}
}
-
+
// creation fichier pour tracer h en fonction de x
//std::ofstream fout("alea");
//for (int r=1; r<connectivity.numberOfNodes()-1; ++r){
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 40298dc5c01217cc5002c3f14f78e12559406ffb..03cf62de92d8da737bdd4a45f5fec418a09cc891 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -181,19 +181,19 @@ private:
m_ur[r]=invAr(r)*br(r);
});
-
+ /*
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;
}
@@ -379,6 +379,13 @@ public:
nuj(j) = 0.5*(1.+xj[j][0]);
});
*/
+
+ // gnuplot output for vitesse riemann
+ std::ofstream fout("ur_split");
+ fout.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
+ fout << xr[j][0] << ' ' << ur[j][0] << '\n';
+ }
}
};
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index fe547951659e34fcc5561c13be2c3cc4f52058ae..281000de92e68e564f57c1ce8b7604c78bb4d406 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)));
@@ -123,13 +123,13 @@ private:
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)/Vj(m_mesh.numberOfCells()-1) + kj(cell_here)/Vj(cell_here))/(1./Vj(cell_here) + 1./Vj(m_mesh.numberOfCells()-1)))*(1./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
// 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;
+ 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
double h = std::sqrt(1. - (t*t)/(50./9.));
@@ -179,16 +179,16 @@ private:
});
// Conditions aux bords
-
+ /*
m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
-
+ */
// Kidder
- /*
+
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(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
- */
+
return m_Gl ;
@@ -236,7 +236,7 @@ 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));
//m_Bl(0) = (nuL(0)/Vj(0) + nuj(cell_here)/Vj(cell_here))/(1./Vj(0) + 1./Vj(cell_here))*((Tj(cell_here)-TL(0))/Vl(0));
@@ -245,7 +245,7 @@ private:
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));
//m_Bl(m_mesh.numberOfFaces()-1) = (nuR(0)/Vj(m_mesh.numberOfCells()-1) + nuj(cell_here)/Vj(cell_here))/(1./Vj(m_mesh.numberOfCells()-1) + 1./Vj(cell_here))*((Tj(cell_here)-TR(0))/Vl(m_mesh.numberOfFaces()-1));
-
+ */
// Kiddder
/*
@@ -256,8 +256,8 @@ private:
*/
// nu = 0
- //m_Bl(0) = 0.;
- //m_Bl(m_mesh.numberOfFaces()-1) = 0.;
+ m_Bl(0) = 0.;
+ m_Bl(m_mesh.numberOfFaces()-1) = 0.;
// nu = 0.2
//m_Bl(0) = (0.2*3.*h*x0[0][0])/(50.*h*h*h*h);
@@ -495,7 +495,7 @@ public:
//Ej[j] -= ((pi*pi*0.5*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) + xj[j][0]*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]))*(std::exp(-2.*t)-1.) - pi*0.5*std::exp(-2.*t)*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]) + (1.+xj[j][0])*((3.*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0])-xj[j][0]*pi*pi*pi*pi*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])))*(std::exp(-2.*t)-1.)+0.5*pi*pi*std::exp(-2.*t)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))))*(dt*inv_mj[j])*Vj(j);
// ajout second membre pour kidder (k = 0.5)
- //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)
//uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t));
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index eade1009bcd074db87372da6056c5985e06b7ef3..646eb3ea5e7ad2c9784f6ea134d701a8ddb8369b 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -71,37 +71,6 @@ private:
}
};
- /*
- // ----
-
- KOKKOS_INLINE_FUNCTION
- const Kokkos::View<const double*>
- computePj(const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& kj,
- const Kokkos::View<const Rd*>& uj)
- {
- const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
-
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_connectivity.cellNbNodes();
-
- const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
-
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- double new_p = 0.;
- double sum = 0;
- for (int m=0; m<cell_nb_nodes(j); ++m) {
- sum += uj[cell_nodes(j,m)][0];
- }
- new_p = pj(j) - kj(j)*sum/Vj(j);
- pj(j) = new_p;
- });
- return pj;
- }
-
- // ----
- */
-
KOKKOS_INLINE_FUNCTION
const Kokkos::View<const double*>
computeRhoCj(const Kokkos::View<const double*>& rhoj,
@@ -154,7 +123,7 @@ 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*>& PTj,
+ 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();
@@ -167,7 +136,7 @@ private:
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 += Ajr(J,R)*uj(J) + PTj(J)*Cjr(J,R);
+ br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
}
});
@@ -190,10 +159,10 @@ private:
m_ur[r]=invAr(r)*br(r);
});
- /*
- 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];
@@ -212,14 +181,14 @@ private:
const Kokkos::View<const Rd*>& ur,
const Kokkos::View<const Rd**>& Cjr,
const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& PTj) {
+ const Kokkos::View<const double*>& pj) {
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_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+PTj(j)*Cjr(j,r);
+ m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
}
});
@@ -248,10 +217,9 @@ private:
KOKKOS_INLINE_FUNCTION
void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
const Kokkos::View<const Rd*>& xj,
- const Kokkos::View<const double*>& kj,
const Kokkos::View<const double*>& rhoj,
const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& PTj,
+ const Kokkos::View<const double*>& pj,
const Kokkos::View<const double*>& cj,
const Kokkos::View<const double*>& Vj,
const Kokkos::View<const Rd**>& Cjr,
@@ -260,12 +228,12 @@ 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, PTj, t);
+ const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj, t);
Kokkos::View<Rd*> ur = m_ur;
Kokkos::View<Rd**> Fjr = m_Fjr;
ur = computeUr(Ar, br, t);
- Fjr = computeFjr(Ajr, ur, Cjr, uj, PTj);
+ Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
}
Kokkos::View<Rd*> m_br;
@@ -275,7 +243,6 @@ private:
Kokkos::View<Rd**> m_Fjr;
Kokkos::View<Rd*> m_ur;
Kokkos::View<double*> m_rhocj;
- Kokkos::View<double*> m_PTj;
Kokkos::View<double*> m_Vj_over_cj;
public:
@@ -291,7 +258,6 @@ public:
m_Fjr("Fjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
m_ur("ur", m_mesh.numberOfNodes()),
m_rhocj("rho_c", m_mesh.numberOfCells()),
- m_PTj("PTj", m_mesh.numberOfCells()),
m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
{
;
@@ -328,52 +294,9 @@ public:
const Kokkos::View<const unsigned short*> cell_nb_nodes
= m_connectivity.cellNbNodes();
- // Calcule les flux
- computeExplicitFluxes(xr, xj, kj, rhoj, uj, PTj, cj, Vj, Cjr, t);
-
const Kokkos::View<const Rd**> Fjr = m_Fjr;
const Kokkos::View<const Rd*> ur = m_ur;
- // Mise a jour de la vitesse et de l'energie totale specifique
- const Kokkos::View<const double*> inv_mj = unknowns.invMj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- Rd momentum_fluxes = zero;
- double energy_fluxes = 0;
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- const int r=cell_nodes(j,R);
- momentum_fluxes += Fjr(j,R);
- energy_fluxes += (Fjr(j,R), ur[r]);
- }
- uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
- Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
-
- // ajout second membre pour kidder (k cst)
- Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
- // ajout second membre pour kidder (k = x)
- //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)));
- });
-
- // Calcul de e par la formule e = E-0.5 u^2
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
- });
-
- // deplace le maillage (ses sommets) en utilisant la vitesse
- // donnee par le schema
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
- });
-
- // met a jour les quantites (geometriques) associees au maillage
- m_mesh_data.updateAllData();
-
- // Calcul de rho avec la formule Mj = Vj rhoj
- const Kokkos::View<const double*> mj = unknowns.mj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- rhoj[j] = mj[j]/Vj[j];
- });
-
/*
// Calcul de PT (1er essai)
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
@@ -401,6 +324,22 @@ public:
const Kokkos::View<const unsigned int**>& cell_faces = m_connectivity.cellFaces();
const Kokkos::View<const unsigned short*> cell_nb_faces
= m_connectivity.cellNbFaces();
+
+ // gnuplot output for vitesse
+ std::ofstream fout("uj");
+ fout.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ }
+
+ // gnuplot output for vitesse riemann
+ std::ofstream fout1("ur");
+ fout1.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
+ fout1 << xr[j][0] << ' ' << ur[j][0] << '\n';
+ }
+
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
std::vector<double> stock(2);
@@ -424,9 +363,25 @@ public:
} else {
PTj(j) = pj(j) - kj(j)*(stock[1]-stock[0])/Vj(j);
}
- */
+ */
// Calcul de PT (3eme essai, avec uR du solveur de Riemann)
+ computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, t);
+
+ // gnuplot output for vitesse
+ std::ofstream fout("uj");
+ fout.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
+ fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+ }
+
+ // gnuplot output for vitesse riemann
+ std::ofstream fout1("ur");
+ fout1.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
+ fout1 << xr[j][0] << ' ' << ur[j][0] << '\n';
+ }
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
if (j == 0) {
@@ -436,17 +391,77 @@ public:
PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
//PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
} else {
- PTj(j) = pj(j) - kj(j)*(ur[cell_nodes(j,1)][0]-ur[cell_nodes(j,0)][0])/Vj(j);
- }
+
+ double sum = 0;
+ for (int k=0; k<cell_nb_nodes(j); ++k) {
+ int node_here = cell_nodes(j,k);
+ sum += (ur(node_here), Cjr(j,k));
+ }
+ PTj(j) = pj(j) - kj(j)*sum/Vj(j);
+ //std::cout << PTj(j) << std::endl;
+ //PTj(j) = pj(j) - kj(j)*(ur[cell_nodes(j,0)][0]-ur[cell_nodes(j,1)][0])/Vj(j);
+ }
});
-
+
+
+ // Calcule les flux
+ computeExplicitFluxes(xr, xj, rhoj, uj, PTj, cj, Vj, Cjr, t);
+
+ // Mise a jour de la vitesse et de l'energie totale specifique
+ const Kokkos::View<const double*> inv_mj = unknowns.invMj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (int R=0; R<cell_nb_nodes[j]; ++R) {
+ const int r=cell_nodes(j,R);
+ momentum_fluxes += Fjr(j,R);
+ energy_fluxes += (Fjr(j,R), ur[r]);
+ }
+ uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
+ Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
+
+ // ajout second membre pour kidder (k cst)
+ Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+ // ajout second membre pour kidder (k = x)
+ //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)));
+ });
+
+ // Calcul de e par la formule e = E-0.5 u^2
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
+ });
+
+ // deplace le maillage (ses sommets) en utilisant la vitesse
+ // donnee par le schema
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
+ xr[r] += dt*ur[r];
+ });
+
+ // met a jour les quantites (geometriques) associees au maillage
+ m_mesh_data.updateAllData();
+
+ // Calcul de rho avec la formule Mj = Vj rhoj
+ const Kokkos::View<const double*> mj = unknowns.mj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ rhoj[j] = mj[j]/Vj[j];
+ });
+
// Mise a jour de k
/*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
kj(j) = xj[j][0];
});
*/
+ /*
+ // gnuplot output for vitesse riemann
+ std::ofstream fout("ur_no_split");
+ fout.precision(15);
+ for (size_t j=0; j<m_mesh.numberOfNodes(); ++j) {
+ fout << xr[j][0] << ' ' << ur[j][0] << '\n';
+ }
+ */
}
};