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

ajout pas non constant avec fonction x^2

parent dec901f6
No related branches found
No related tags found
No related merge requests found
......@@ -170,6 +170,7 @@ int main(int argc, char *argv[])
t += dt;
++iteration;
std::cout << "temps t : " << t << std::endl;
}
......@@ -197,7 +198,7 @@ int main(int argc, char *argv[])
double pi = 4.*std::atan(1.);
std::ofstream fout("comparaison u");
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.001) <<'\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
}
}
......
......@@ -50,15 +50,15 @@ public:
// pas constant
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes())
{
const double delta_x = 1./connectivity.numberOfCells();
Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
m_xr[r][0] = r*delta_x;
});
}
// Mesh(const Connectivity& connectivity)
//: m_connectivity(connectivity),
// m_xr("xr", connectivity.numberOfNodes())
//{
// const double delta_x = 1./connectivity.numberOfCells();
// Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
// m_xr[r][0] = r*delta_x;
// });
//}
// pas non constant
......@@ -77,6 +77,18 @@ public:
// });
//}
// pas non constant avec fonction x^2
Mesh(const Connectivity& connectivity)
: m_connectivity(connectivity),
m_xr("xr", connectivity.numberOfNodes())
{
const double delta_x = 1./connectivity.numberOfCells();
Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
m_xr[r][0] = (r*delta_x)*(r*delta_x);
});
}
~Mesh()
{
;
......
......@@ -214,10 +214,11 @@ public:
for (int ll=0; ll<cell_nb_faces(j); ++ll) {
minVl = std::min(minVl, Vl(cell_faces(j, ll)));
}
//k=2 => (kj(j+1) + 2*kj(j) + kj(j-1)) = 8
//dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./8.)*minVl;
// k=x
// k pas forcement constant
double sum = 0.;
for (int m = 0; m < cell_nb_nodes(j); ++m) {
sum += kj(cell_nodes(j,m));
......@@ -279,7 +280,7 @@ public:
momentum_fluxes += Fl(l)*Cjr(j,R);
energy_fluxes += (Gl(l), Cjr(j,R));
}
uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_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] += (dt*inv_mj[j]) * momentum_fluxes;
Ej[j] += (dt*inv_mj[j]) * energy_fluxes;
});
......@@ -313,13 +314,17 @@ public:
double erreur = 0.;
double exacte = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte
//exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.001); // solution exacte cas test k constant
exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
erreur += (exacte - uj[j][0])*(exacte - uj[j][0]);
}
erreur = std::sqrt(erreur);
return erreur;
}
// Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
// (quand la solution exacte est connue)
double error_Linf(UnknownsType& unknowns) {
Kokkos::View<Rd*> uj = unknowns.uj();
......@@ -327,10 +332,12 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
double pi = 4.*std::atan(1.);
double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2);
//double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.001); // k constant
double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant
double erreur = std::abs(exacte - uj[0][0]);
for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
//exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.001);
exacte = std::sin(pi*xj[j][0])*std::exp(-0.2);
if (std::abs(exacte - uj[j][0]) > erreur) {
......
......@@ -228,7 +228,7 @@ void initializeSod()
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_kj[j] = 2; // k constant
//m_kj[j] = 2.; // k constant
m_kj[j] = xj[j][0]; // k non constant, k = x
});
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment