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

tests pas de temps

parent 593f6bf0
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/main.cpp
+ 19
− 18
View file @ cf5640e5
...@@ -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=1.5; 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();
...@@ -282,9 +282,9 @@ int main(int argc, char *argv[]) ...@@ -282,9 +282,9 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) { while((t<tmax) and (iteration<itermax)) {
// NAVIER-STOKES AVEC SPLITTING // NAVIER-STOKES AVEC SPLITTING
/*
// Etape 1 du splitting - Euler // Etape 1 du splitting - Euler
double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj); double dt_euler = 0.5*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) { if (t+dt_euler > tmax) {
dt_euler = tmax-t; dt_euler = tmax-t;
} }
...@@ -292,14 +292,15 @@ int main(int argc, char *argv[]) ...@@ -292,14 +292,15 @@ int main(int argc, char *argv[])
t += dt_euler; t += dt_euler;
// Etape 2 du splitting - Diffusion // Etape 2 du splitting - Diffusion
double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR); double param = 1.;
double dt_diff = param*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
double t_diff = t-dt_euler; double t_diff = t-dt_euler;
if (dt_euler <= dt_diff) { if (dt_euler <= dt_diff) {
dt_diff = dt_euler; dt_diff = dt_euler;
finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
} else { } else {
while (t > t_diff) { while (t > t_diff) {
dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj,nuL, nuR, kL,kR); dt_diff = param*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj,nuL, nuR, kL,kR);
if (t_diff+dt_diff > t) { if (t_diff+dt_diff > t) {
dt_diff = t-t_diff; dt_diff = t-t_diff;
} }
...@@ -307,7 +308,7 @@ int main(int argc, char *argv[]) ...@@ -307,7 +308,7 @@ int main(int argc, char *argv[])
t_diff += dt_diff; t_diff += dt_diff;
} }
} }
*/
// Diffusion pure // Diffusion pure
/* /*
double dt = 0.9*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);
...@@ -319,14 +320,14 @@ int main(int argc, char *argv[]) ...@@ -319,14 +320,14 @@ int main(int argc, char *argv[])
*/ */
// NAVIER-STOKES SANS SPLITTING // NAVIER-STOKES SANS SPLITTING
/*
double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR); double dt = 0.1*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
if (t+dt > tmax) { if (t+dt > tmax) {
dt = tmax-t; dt = tmax-t;
} }
no_splitting.computeNextStep(t,dt, unknowns); no_splitting.computeNextStep(t,dt, unknowns);
t += dt; t += dt;
*/
block_eos.updatePandCFromRhoE(); block_eos.updatePandCFromRhoE();
...@@ -595,7 +596,7 @@ int main(int argc, char *argv[]) ...@@ -595,7 +596,7 @@ int main(int argc, char *argv[])
// Erreurs sous differents normes // Erreurs sous differents normes
/*
// Density // Density
double error1 = 0.; double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax); error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
...@@ -612,7 +613,7 @@ int main(int argc, char *argv[]) ...@@ -612,7 +613,7 @@ int main(int argc, char *argv[])
error = finite_volumes_diffusion.error_L2_u(unknowns, tmax); error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
<< ": " << rang::fgB::green << error << rang::fg::reset << " \n"; << ": " << rang::fgB::green << error << rang::fg::reset << " \n";
/*
double error4 = 0.; double error4 = 0.;
error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax); error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);
std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
...@@ -661,8 +662,8 @@ int main(int argc, char *argv[]) ...@@ -661,8 +662,8 @@ int main(int argc, char *argv[])
std::ofstream fout("rho"); std::ofstream fout("rho");
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';
} }
} }
...@@ -696,10 +697,10 @@ int main(int argc, char *argv[]) ...@@ -696,10 +697,10 @@ int main(int argc, char *argv[])
//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';
} }
} }
...@@ -714,10 +715,10 @@ int main(int argc, char *argv[]) ...@@ -714,10 +715,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] << ' ' << (-(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] << ' ' << ((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] << ' ' << 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';
} }
} }
......
src/mesh/Mesh.hpp
+ 2
− 2
View file @ cf5640e5
...@@ -82,8 +82,8 @@ public: ...@@ -82,8 +82,8 @@ public:
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();
......
src/scheme/FiniteVolumesEulerUnknowns.hpp
+ 19
− 16
View file @ cf5640e5
...@@ -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);
...@@ -352,8 +352,8 @@ public: ...@@ -352,8 +352,8 @@ 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.5; //m_kj[j] = 0.014;
// TAC // TAC
...@@ -392,11 +392,11 @@ public: ...@@ -392,11 +392,11 @@ public:
*/ */
// k regulier // k regulier
/*
int n = 1; int n = 1;
m_kj[j] = std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.7)*(xj[j][0]<0.7+0.1/n) + std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.9-0.1/n)*(xj[j][0]<0.9) + (xj[j][0]>0.7+0.1/n)*(xj[j][0]<0.9-0.1/n); m_kj[j] = std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.7)*(xj[j][0]<0.7+0.1/n) + std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.9-0.1/n)*(xj[j][0]<0.9) + (xj[j][0]>0.7+0.1/n)*(xj[j][0]<0.9-0.1/n);
m_kj[j] = 0.014*m_kj[j]; m_kj[j] = 0.014*m_kj[j];
*/
}); });
...@@ -409,7 +409,7 @@ public: ...@@ -409,7 +409,7 @@ public:
m_uL[0] = zero; m_uL[0] = zero;
m_uR[0] = zero; m_uR[0] = zero;
m_kL[0] = 0.; m_kL[0] = 0.;
m_kR[0] = 1.; m_kR[0] = 0.;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
...@@ -422,13 +422,16 @@ public: ...@@ -422,13 +422,16 @@ public:
}); });
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_Tj[j] = m_ej[j]; //m_Tj[j] = m_ej[j];
m_Tj[j] = 0.;
}); });
// Conditions aux bords de Dirichlet sur T et nu // Conditions aux bords de Dirichlet sur T et nu
m_TL[0] = m_ej[0]; //m_TL[0] = m_ej[0];
m_TR[0] = m_ej[m_mesh.numberOfCells()-1]; //m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
m_TL[0] = 0.;
m_TR[0] = 0.;
m_nuL[0] = 0.; m_nuL[0] = 0.;
m_nuR[0] = 0.; m_nuR[0] = 0.;
......
src/scheme/NoSplitting.hpp
+ 12
− 12
View file @ cf5640e5
...@@ -190,19 +190,19 @@ private: ...@@ -190,19 +190,19 @@ private:
m_ur[r]=invAr(r)*br(r); m_ur[r]=invAr(r)*br(r);
}); });
/*
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;
} }
...@@ -350,8 +350,8 @@ public: ...@@ -350,8 +350,8 @@ public:
// 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)));
}); });
// Calcul de e par la formule e = E-0.5 u^2 // Calcul de e par la formule e = E-0.5 u^2
...@@ -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)*2.*sum/sum1; PTj(j) = pj(j) - kj(j)*2.*sum/sum1;
} }
...@@ -396,12 +396,12 @@ public: ...@@ -396,12 +396,12 @@ public:
}); });
// Mise a jour de k // Mise a jour de k
/*
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
kj(j) = xj[j][0]; kj(j) = xj[j][0];
}); });
*/
} }
}; };
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment