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

bonne version pour étudier la convergence pour test kidder avec transfert thermique

parent ddaafc81
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/main.cpp
+ 25
− 22
View file @ 1ee6b885
......@@ -142,7 +142,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=0.2;
const double tmax=1.5;
double t=0.;
int itermax=std::numeric_limits<int>::max();
......@@ -262,7 +262,7 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
/*
// ETAPE 1 DU SPLITTING - EULER
double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
......@@ -291,7 +291,7 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
*/
// AUTRE APPROCHE DU SPLITTING (PLUS LONG)
/*
......@@ -310,7 +310,7 @@ int main(int argc, char *argv[])
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
*/
/*
// DIFFUSION PURE
double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
......@@ -319,7 +319,7 @@ int main(int argc, char *argv[])
}
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
*/
block_eos.updatePandCFromRhoE();
......@@ -555,7 +555,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
/*
double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
......@@ -567,7 +567,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset
<< ": " << rang::fgB::green << error2 << rang::fg::reset << " \n";
*/
double error = 0.;
error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
......@@ -575,6 +575,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
<< ": " << rang::fgB::green << error << rang::fg::reset << " \n";
/*
double error4 = 0.;
error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);
......@@ -592,6 +593,7 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Erreur L infini E" << rang::style::reset
<< ": " << rang::fgB::green << error5 << rang::fg::reset << " \n";
*/
std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
<< ": " << rang::fgB::green << c << rang::fg::reset << " \n";
......@@ -606,15 +608,16 @@ int main(int argc, char *argv[])
//method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds();
{ // gnuplot output for density
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");
double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("rho1600");
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
fout << xj[j][0] << ' ' << rhoj[j] << '\n';
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';
}
}
......@@ -622,13 +625,13 @@ 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");
std::ofstream fout("u1600");
fout.precision(15);
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(-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] << ' ' << 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] << '\n';
}
......@@ -637,15 +640,15 @@ int main(int argc, char *argv[])
{ // gnuplot output for energy
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> Ej = unknowns.Ej();
double pi = 4.*std::atan(1.);
//double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("E");
//double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("E1600");
fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
//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] << ' ' << ((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] << '\n';
}
......
src/scheme/AcousticSolver.hpp
+ 4
− 3
View file @ 1ee6b885
......@@ -197,9 +197,10 @@ private:
m_ur[r]=invAr(r)*br(r);
});
/*
m_ur[0]=zero;
m_ur[m_mesh.numberOfNodes()-1]=zero;
*/
// Kidder
// Conditions aux limites dependant du temps
......@@ -207,11 +208,11 @@ private:
//m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
//R(t) = x*h(t) a la place de x(t)
/*
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;
}
......
src/scheme/FiniteVolumesDiffusion.hpp
+ 43
− 24
View file @ 1ee6b885
......@@ -115,7 +115,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)));
......@@ -124,7 +124,7 @@ 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) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
*/
// Kidder
......@@ -135,11 +135,11 @@ private:
// k = x
//m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
//m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
/*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0];
m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0];
*/
return m_Fl ;
}
......@@ -179,18 +179,19 @@ private:
});
// Conditions aux bords
m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
//m_Gl(0) = Fl(0)*uL(0);
//m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
// Kidder
// m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
//m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
/*
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 ;
......@@ -218,6 +219,9 @@ private:
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd*> x0 = m_mesh.x0();
const Kokkos::View<const Rd*> xmax = m_mesh.xmax();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rd sum = zero;
double sum2 = 0.;
......@@ -236,11 +240,20 @@ private:
// Conditions aux bords
// Diffusion pure
/*
int cell_here = face_cells(0,0);
m_Bl(0) = (nuL(0) + nuj(cell_here))*(1./(2*Vl(0)))*(Tj(cell_here) - TL(0));
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));
*/
// Kidder
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Bl(0) = ((1. + x0[0][0])*3.*x0[0][0])/(100.*h*h*h*h);
m_Bl(m_mesh.numberOfFaces()-1) = ((1. + xmax[0][0])*3.*xmax[0][0])/(100.*h*h*h*h);
return m_Bl ;
}
......@@ -394,12 +407,19 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const Rd*> x0 = m_mesh.x0();
const Kokkos::View<const Rd*> xmax = m_mesh.xmax();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
double pi = 4.*std::atan(1.);
TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
// la condition au bord a droite de T depend du temps
double h = std::sqrt(1. - (t*t)/(50./9.));
// Les CL de T dependent du temps
// Diffusion pure
//TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
......@@ -416,8 +436,8 @@ public:
// 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) {
//const int j = j0+1;
Kokkos::parallel_for(m_mesh.numberOfCells()-2, KOKKOS_LAMBDA(const int& j0) {
const int j = j0+1;
Rd momentum_fluxes = zero;
double energy_fluxes = 0.;
Rd trich = zero;
......@@ -433,16 +453,15 @@ public:
Ej[j] += (dt*inv_mj[j]) * energy_fluxes;
// test k non cst pour diff pure
uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi);
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);
//uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi);
//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)));
// 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))-(6*xj[j][0]+3.)/(100*(1-t*t/(50/9))*(1-t*t/(50/9))));
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))+(6*xj[j][0]+3.)/(100*(1-t*t/(50/9))*(1-t*t/(50/9))));
});
// Calcul de e par la formule e = E-0.5 u^2
......@@ -473,8 +492,8 @@ public:
double exact_u = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
//exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
exact_u = std::sin(pi*xj[j][0])*std::exp(-t); // solution exacte cas test k non constant
//exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
//exact_u = std::sin(pi*xj[j][0])*std::exp(-t); // solution exacte cas test k non constant
exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
}
err_u = std::sqrt(err_u);
......
src/scheme/FiniteVolumesEulerUnknowns.hpp
+ 29
− 11
View file @ 1ee6b885
......@@ -276,7 +276,7 @@ public:
}
/*
// --- Acoustic Solver ---
void initializeSod()
......@@ -284,25 +284,27 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
/*
if (xj[j][0]<0.5) {
m_rhoj[j]=1.;
} else {
m_rhoj[j]=0.125;
}
*/
//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){
/*
if (xj[j][0]<0.5) {
m_pj[j]=1;
} else {
m_pj[j]=0.1;
}
*/
//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.);
......@@ -311,8 +313,8 @@ public:
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_gammaj[j] = 1.4;
//m_gammaj[j] = 3.;
//m_gammaj[j] = 1.4;
m_gammaj[j] = 3.;
});
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
......@@ -333,9 +335,10 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_kj[j] = xj[j][0];
//m_kj[j] = 0.5;
//m_kj[j] = 0.5;
/*
// Sod
// k non regulier
......@@ -376,7 +379,7 @@ public:
int n = 10.;
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];
*/
});
......@@ -392,12 +395,26 @@ public:
m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));
m_S0(j) = m_entropy(j);
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_Tj[j] = m_ej[j];
});
// Conditions aux bords de Dirichlet sur T et nu
m_TL[0] = 1.;
m_TR[0] = 103./100.;
m_nuL[0] = 0.5;
m_nuR[0] = 1.;
}
*/
/*
// DIFFUSION PURE
......@@ -467,6 +484,7 @@ public:
m_nuR[0] = 1.;
}
*/
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
: m_mesh_data(mesh_data),
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment