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

test TAC

parent 992a21b5
No related branches found
No related tags found
No related merge requests found
...@@ -144,7 +144,7 @@ int main(int argc, char *argv[]) ...@@ -144,7 +144,7 @@ int main(int argc, char *argv[])
const Kokkos::View<const double*> Vj = mesh_data.Vj(); const Kokkos::View<const double*> Vj = mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.2; const double tmax=0.8;
double t=0.; double t=0.;
int itermax=std::numeric_limits<int>::max(); int itermax=std::numeric_limits<int>::max();
...@@ -658,18 +658,18 @@ int main(int argc, char *argv[]) ...@@ -658,18 +658,18 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj(); const Kokkos::View<const double*> rhoj = unknowns.rhoj();
double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("rho"); std::ofstream fout("rho_no_split");
fout.precision(15); fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) { 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] << ' ' << 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] << '\n';
} }
} }
{ // gnuplot output for pression { // gnuplot output for pression
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> pj = unknowns.pj(); const Kokkos::View<const double*> pj = unknowns.pj();
std::ofstream fout("p"); std::ofstream fout("p_no_split");
fout.precision(15); fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) { for (size_t j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << pj[j] << '\n'; fout << xj[j][0] << ' ' << pj[j] << '\n';
...@@ -679,7 +679,7 @@ int main(int argc, char *argv[]) ...@@ -679,7 +679,7 @@ int main(int argc, char *argv[])
{ // gnuplot output for internal energy { // gnuplot output for internal energy
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> ej = unknowns.ej(); const Kokkos::View<const double*> ej = unknowns.ej();
std::ofstream fout("e"); std::ofstream fout("e_no_split");
fout.precision(15); fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) { for (size_t j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << ej[j] << '\n'; fout << xj[j][0] << ' ' << ej[j] << '\n';
...@@ -690,16 +690,16 @@ int main(int argc, char *argv[]) ...@@ -690,16 +690,16 @@ int main(int argc, char *argv[])
const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const Rd*> uj = unknowns.uj(); const Kokkos::View<const Rd*> uj = unknowns.uj();
double pi = 4.*std::atan(1.); double pi = 4.*std::atan(1.);
std::ofstream fout("u"); std::ofstream fout("u_no_split");
fout.precision(15); fout.precision(15);
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'; //cas k constant //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] << ' ' << 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] << ' ' << xj[j][0] << std::endl;
//fout << xj[j][0] << ' ' << uj[j][0] << '\n'; fout << xj[j][0] << ' ' << uj[j][0] << '\n';
} }
} }
......
...@@ -75,15 +75,15 @@ public: ...@@ -75,15 +75,15 @@ public:
// pas constant // pas constant
/*
Mesh(const Connectivity& connectivity) Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity), : m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes()), m_xr("xr", connectivity.numberOfNodes()),
m_x0("x0", 1), m_x0("x0", 1),
m_xmax("xmax", 1) m_xmax("xmax", 1)
{ {
double a = 0.; double a = -1.;
double b = 1.; double b = 2.;
m_x0[0][0] = a; m_x0[0][0] = a;
m_xmax[0][0] = b; m_xmax[0][0] = b;
const double delta_x = (b-a)/connectivity.numberOfCells(); const double delta_x = (b-a)/connectivity.numberOfCells();
...@@ -91,7 +91,7 @@ public: ...@@ -91,7 +91,7 @@ public:
m_xr[r][0] = a + r*delta_x; m_xr[r][0] = a + r*delta_x;
}); });
} }
*/
// pas non constant // pas non constant
/* /*
...@@ -125,15 +125,15 @@ public: ...@@ -125,15 +125,15 @@ public:
// pas aleatoire // pas aleatoire
/*
Mesh(const Connectivity& connectivity) Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity), : m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes()), m_xr("xr", connectivity.numberOfNodes()),
m_x0("x0", 1), m_x0("x0", 1),
m_xmax("xmax", 1) m_xmax("xmax", 1)
{ {
double a = 0.; double a = -1.;
double b = 1.; double b = 2.;
m_x0[0][0] = a; m_x0[0][0] = a;
m_xmax[0][0] = b; m_xmax[0][0] = b;
const double h = (b-a)/connectivity.numberOfCells(); const double h = (b-a)/connectivity.numberOfCells();
...@@ -145,7 +145,7 @@ public: ...@@ -145,7 +145,7 @@ public:
for (int r=1; r<connectivity.numberOfNodes()-1; ++r){ for (int r=1; r<connectivity.numberOfNodes()-1; ++r){
const double delta_xr = dist(mt); const double delta_xr = dist(mt);
m_xr[r][0] = r*h + delta_xr; m_xr[r][0] = a + r*h + delta_xr;
} }
// creation fichier pour tracer h en fonction de x // creation fichier pour tracer h en fonction de x
...@@ -160,7 +160,7 @@ public: ...@@ -160,7 +160,7 @@ public:
//std::exit(0); //std::exit(0);
} }
*/
~Mesh() ~Mesh()
......
...@@ -182,18 +182,18 @@ private: ...@@ -182,18 +182,18 @@ private:
}); });
//m_ur[0]=zero; m_ur[0]=zero;
//m_ur[m_mesh.numberOfNodes()-1]=zero; m_ur[m_mesh.numberOfNodes()-1]=zero;
//m_ur[0] = x0; //m_ur[0] = x0;
//m_ur[m_mesh.numberOfNodes()-1] = xmax[0]; //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
// CL Kidder // CL Kidder
/*
double h = std::sqrt(1. - (t*t)/(50./9.)); double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0]; 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]; m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
*/
return m_ur; return m_ur;
} }
......
...@@ -299,28 +299,28 @@ public: ...@@ -299,28 +299,28 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// TAC // TAC
/*
if (xj[j][0]<0.5) { if (xj[j][0]<0.5) {
m_rhoj[j]=1.; m_rhoj[j]=1.;
} else { } else {
m_rhoj[j]=0.125; m_rhoj[j]=0.125;
} }
*/
// Kidder // 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){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// TAC // TAC
/*
if (xj[j][0]<0.5) { if (xj[j][0]<0.5) {
m_pj[j]=1; m_pj[j]=1.;
} else { } else {
m_pj[j]=0.1; m_pj[j]=0.1;
} }
*/
// Kidder // 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.); double pi = 4.*std::atan(1.);
...@@ -330,9 +330,9 @@ public: ...@@ -330,9 +330,9 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// TAC // TAC
// m_gammaj[j] = 1.4; m_gammaj[j] = 1.4;
// Kidder // Kidder
m_gammaj[j] = 3.; //m_gammaj[j] = 3.;
}); });
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
...@@ -353,7 +353,7 @@ public: ...@@ -353,7 +353,7 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
// Differents k (xi) // Differents k (xi)
//m_kj[j] = xj[j][0]; //m_kj[j] = xj[j][0];
m_kj[j] = 0.005; m_kj[j] = 0.014;
// TAC // TAC
...@@ -408,8 +408,8 @@ public: ...@@ -408,8 +408,8 @@ public:
m_uL[0] = zero; m_uL[0] = zero;
m_uR[0] = zero; m_uR[0] = zero;
m_kL[0] = 0.005; m_kL[0] = 0.014;
m_kR[0] = 0.005; m_kR[0] = 0.014;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
......
...@@ -191,18 +191,18 @@ private: ...@@ -191,18 +191,18 @@ private:
}); });
//m_ur[0]=zero; m_ur[0]=zero;
//m_ur[m_mesh.numberOfNodes()-1]=zero; m_ur[m_mesh.numberOfNodes()-1]=zero;
//m_ur[0] = x0; //m_ur[0] = x0;
//m_ur[m_mesh.numberOfNodes()-1] = xmax[0]; //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
// CL Kidder // CL Kidder
/*
double h = std::sqrt(1. - (t*t)/(50./9.)); double h = std::sqrt(1. - (t*t)/(50./9.));
m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0]; 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]; m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
*/
return m_ur; return m_ur;
} }
...@@ -348,7 +348,7 @@ public: ...@@ -348,7 +348,7 @@ public:
Ej[j] -= (dt*inv_mj[j]) * energy_fluxes; Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
// ajout second membre pour kidder (k cst) // 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))); //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) // ajout second membre pour kidder (k = x)
//uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t)); //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))); //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
...@@ -384,11 +384,11 @@ public: ...@@ -384,11 +384,11 @@ public:
} }
if (j == 0) { if (j == 0) {
//PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0); PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t)); //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
} else if (j == m_mesh.numberOfCells()-1) { } else if (j == m_mesh.numberOfCells()-1) {
PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t)); //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); PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
} else { } else {
//PTj(j) = pj(j) - kj(j)*sum/Vl(j); //PTj(j) = pj(j) - kj(j)*sum/Vl(j);
PTj(j) = pj(j) - kj(j)*2.*sum/sum1; PTj(j) = pj(j) - kj(j)*2.*sum/sum1;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment