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

modifs sur les CL, pas terrible

parent d5474e01
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/main.cpp
+ 10
− 23
View file @ d41c3067
...@@ -138,7 +138,7 @@ int main(int argc, char *argv[]) ...@@ -138,7 +138,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=1.;
double t=0; double t=0;
int itermax=std::numeric_limits<int>::max(); int itermax=std::numeric_limits<int>::max();
...@@ -159,12 +159,6 @@ int main(int argc, char *argv[]) ...@@ -159,12 +159,6 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
<< ": " << rang::fgB::green << c << rang::fg::reset << " \n"; << ": " << rang::fgB::green << c << rang::fg::reset << " \n";
double c1 = 0.;
c1 = finite_volumes_diffusion.conservatif_mvt(unknowns);
std::cout << "* " << rang::style::underline << "Resultat conservativite rho u temps = 0" << rang::style::reset
<< ": " << rang::fgB::green << c1 << rang::fg::reset << " \n";
while((t<tmax) and (iteration<itermax)) { while((t<tmax) and (iteration<itermax)) {
// ETAPE 1 DU SPLITTING - EULER // ETAPE 1 DU SPLITTING - EULER
...@@ -219,14 +213,14 @@ int main(int argc, char *argv[]) ...@@ -219,14 +213,14 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
/*
double error = 0.; double error = 0.;
error = finite_volumes_diffusion.error_L2_u(unknowns); error = finite_volumes_diffusion.error_L2_u(unknowns);
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 error2 = 0.; double error2 = 0.;
error2 = finite_volumes_diffusion.error_Linf(unknowns); error2 = finite_volumes_diffusion.error_Linf(unknowns);
...@@ -247,12 +241,6 @@ int main(int argc, char *argv[]) ...@@ -247,12 +241,6 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::style::underline << "Resultat conservativite E" << rang::style::reset std::cout << "* " << rang::style::underline << "Resultat conservativite E" << rang::style::reset
<< ": " << rang::fgB::green << cons << rang::fg::reset << " \n"; << ": " << rang::fgB::green << cons << rang::fg::reset << " \n";
double cons1 = 0.;
cons1 = finite_volumes_diffusion.conservatif_mvt(unknowns);
std::cout << "* " << rang::style::underline << "Resultat conservativite rho u" << rang::style::reset
<< ": " << rang::fgB::green << cons1 << rang::fg::reset << " \n";
//method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); //method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds(); method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds();
...@@ -260,8 +248,8 @@ int main(int argc, char *argv[]) ...@@ -260,8 +248,8 @@ int main(int argc, char *argv[])
{ // gnuplot output for density { // gnuplot output for density
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 h = std::sqrt(1. - (0.2*0.2)/(50./9.)); double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
std::ofstream fout("rho1"); 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] << '\n'; //fout << xj[j][0] << ' ' << rhoj[j] << '\n';
...@@ -273,28 +261,27 @@ int main(int argc, char *argv[]) ...@@ -273,28 +261,27 @@ 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("u1"); std::ofstream fout("u");
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(-0.2) <<'\n'; // cas k non constant //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\n'; // cas k non constant
//fout << xj[j][0] << ' ' << uj[j][0] << '\n'; //fout << xj[j][0] << ' ' << uj[j][0] << '\n';
fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*0.2)/((50./9.)-0.2*0.2) << '\n'; fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n';
} }
} }
{ // gnuplot output for energy { // gnuplot output for 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();
//double pi = 4.*std::atan(1.); //double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (0.2*0.2)/(50./9.)); double h = std::sqrt(1. - (0.2*0.2)/(50./9.));
std::ofstream fout("e"); std::ofstream fout("E");
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] << ' ' << (-(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.*0.2)-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.*0.2)-1.) + 2. <<'\n' ; // cas k non constant
//fout << xj[j][0] << ' ' << Ej[j] << '\n'; fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h)*(std::sqrt((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h) + (-(xj[j][0]*0.2)/((50./9.)-0.2*0.2))*(-(xj[j][0]*0.2)/((50./9.)-0.2*0.2))*0.5 << '\n';
fout << xj[j][0] << ' ' << ej[j] << ' ' << ( std::sqrt((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h)*( std::sqrt((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h) << '\n';
} }
} }
......
src/scheme/AcousticSolver.hpp
+ 2
− 2
View file @ d41c3067
src/scheme/FiniteVolumesDiffusion.hpp
+ 21
− 14
View file @ d41c3067
...@@ -116,9 +116,11 @@ private: ...@@ -116,9 +116,11 @@ private:
int cell_here = face_cells(0,0); int cell_here = face_cells(0,0);
int local_face_number_in_cell = face_cell_local_face(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))); 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)));
//m_Fl(0) = -xr[0][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
cell_here = face_cells(m_mesh.numberOfFaces()-1,0); cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0); 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) + 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)));
return m_Fl ; return m_Fl ;
} }
...@@ -128,7 +130,8 @@ private: ...@@ -128,7 +130,8 @@ private:
computeGl(const Kokkos::View<const Rd*>& uj, computeGl(const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const Rdd*>& Fl, const Kokkos::View<const Rdd*>& Fl,
const Kokkos::View<const Rd*>& uL, const Kokkos::View<const Rd*>& uL,
const Kokkos::View<const Rd*>& uR) { const Kokkos::View<const Rd*>& uR,
const double& t) {
const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells(); const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
...@@ -157,6 +160,9 @@ private: ...@@ -157,6 +160,9 @@ private:
// Conditions aux bords // Conditions aux bords
m_Gl(0) = Fl(0)*uL(0); m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0); m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
//const double delta_x = 1./100.;
//double brico = - ((xr[m_mesh.numberOfNodes()-1][0]+delta_x)*t)/((50./9.)-t*t);
//m_Gl[m_mesh.numberOfFaces()-1][0] = brico*Fl[m_mesh.numberOfFaces()-1](0,0);
return m_Gl ; return m_Gl ;
...@@ -189,11 +195,12 @@ private: ...@@ -189,11 +195,12 @@ private:
const Kokkos::View<const Rd*>& uL, const Kokkos::View<const Rd*>& uL,
const Kokkos::View<const Rd*>& uR, const Kokkos::View<const Rd*>& uR,
const Kokkos::View<const double*>& kL, const Kokkos::View<const double*>& kL,
const Kokkos::View<const double*>& kR) { const Kokkos::View<const double*>& kR,
const double& t) {
Kokkos::View<Rdd*> Fl = m_Fl ; Kokkos::View<Rdd*> Fl = m_Fl ;
Fl = computeFl (Cjr, uj, kj, uL, uR, kL, kR); Fl = computeFl (Cjr, uj, kj, uL, uR, kL, kR);
Kokkos::View<Rd*> Gl = m_Gl ; Kokkos::View<Rd*> Gl = m_Gl ;
Gl = computeGl (uj, Fl, uL, uR); Gl = computeGl (uj, Fl, uL, uR, t);
} }
Kokkos::View<Rdd*> m_Fl; Kokkos::View<Rdd*> m_Fl;
...@@ -278,20 +285,19 @@ public: ...@@ -278,20 +285,19 @@ public:
Kokkos::View<double*> kL = unknowns.kL(); Kokkos::View<double*> kL = unknowns.kL();
Kokkos::View<double*> kR = unknowns.kR(); Kokkos::View<double*> kR = unknowns.kR();
Kokkos::View<const Rd*> xr = m_mesh.xr();
//uR[0] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
// Premiere possibilite CL (le mieux a mon avis et le plus logique)
uR[0] = (-t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1]; uR[0] = (-t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1];
// Deuxieme possiblitie CL
Kokkos::View<const Rd*> xr = m_mesh.xr();
//uR[0] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
// Calcule les flux // Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR); computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, t);
const Kokkos::View<const Rdd*> Fl = m_Fl ; const Kokkos::View<const Rdd*> Fl = m_Fl ;
const Kokkos::View<const Rd *> Gl = m_Gl ; const Kokkos::View<const Rd *> Gl = m_Gl ;
...@@ -318,7 +324,7 @@ public: ...@@ -318,7 +324,7 @@ public:
Ej[j] += (dt*inv_mj[j]) * energy_fluxes; Ej[j] += (dt*inv_mj[j]) * energy_fluxes;
//uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant //uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant
// ajout second membre pour kidder (k non constant) // 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)));
...@@ -351,12 +357,13 @@ public: ...@@ -351,12 +357,13 @@ public:
const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
double pi = 4.*std::atan(1.); //double pi = 4.*std::atan(1.);
double err_u = 0.; double err_u = 0.;
double exact_u = 0.; double exact_u = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { 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(-2.*pi*pi*0.2); // solution exacte cas test k constant
exact_u = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant // exact_u = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
exact_u = -(xj[j][0]*0.2)/((50./9.)-0.2*0.2);
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j); err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
} }
err_u = std::sqrt(err_u); err_u = std::sqrt(err_u);
......
src/scheme/FiniteVolumesEulerUnknowns.hpp
+ 4
− 2
View file @ d41c3067
...@@ -241,10 +241,12 @@ public: ...@@ -241,10 +241,12 @@ public:
// Conditions aux bords de Dirichlet sur u et k // Conditions aux bords de Dirichlet sur u et k
m_uL[0] = zero; m_uL[0] = zero;
m_uR[0] = zero; m_uR[0] = zero;
m_kL[0] = 0.; m_kL[0] = xj[0][0];
m_kR[0] = 1.; m_kR[0] = xj[m_mesh.numberOfCells()-1][0];
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment