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

ajout code pour tracer x = f(t)

parent de7475f4
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/main.cpp
+ 66
− 7
View file @ e11acce1
......@@ -161,6 +161,12 @@ int main(int argc, char *argv[])
double c = 0.;
c = finite_volumes_diffusion.conservatif(unknowns);
const Kokkos::View<const Rd*> xj = mesh_data.xj();
int size = 3000;
std::vector<std::vector<double>> x(size, std::vector<double>(mesh.numberOfCells()));
std::vector<double> tempo(size);
int i = 0;
/*
// Ecriture des valeurs initiales dans un fichier
......@@ -260,12 +266,13 @@ int main(int argc, char *argv[])
diff.close();
*/
while((t<tmax) and (iteration<itermax)) {
// ETAPE 1 DU SPLITTING - EULER
double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
......@@ -308,10 +315,10 @@ int main(int argc, char *argv[])
++iteration;
std::cout << "temps t : " << t << std::endl;
// ECRITURE DANS UN FICHIER
/*
//if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
// ECRITURE DANS UN FICHIER
if ((std::fmod(t,0.01) < 0.0001) or (t == tmax)) {
std::string ligne;
......@@ -532,11 +539,48 @@ int main(int argc, char *argv[])
// ENTROPY TEST
//finite_volumes_diffusion.entropie(unknowns);
// STOCKAGE COORDONNES ET TEMPS
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
x[i][j] = xj[j][0];
}
tempo[i] = t;
i = i + 1;
}
std::cout << "i = " << i << std::endl;
std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
// CREATION FICHIER x = f(t)
/*
std::ofstream fout("cara");
fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
for (size_t k=0; k<i; ++k) {
fout << tempo[k] << ' ' << x[k][j] << '\n';
}
fout << ' ' << '\n';
fout << ' ' << '\n';
}
*/
std::ofstream fout("cara1");
fout.precision(15);
for (int j=0; j<mesh.numberOfCells(); ++j) {
if (j%10 == 0) {
for (int k=0; k<i; ++k) {
if ( (k%10 == 0) or (k == i-1) ) {
fout << tempo[k] << ' ' << x[k][j] << '\n';
}
}
fout << ' ' << '\n';
fout << ' ' << '\n';
}
}
/*
double error1 = 0.;
error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
......@@ -668,6 +712,21 @@ int main(int argc, char *argv[])
}
*/
{ // gnuplot output for vitesse
const Kokkos::View<const Rd*> xj = mesh_data.xj();
const Kokkos::View<const double*> rhoj = unknowns.rhoj();
const Kokkos::View<const double*> pj = unknowns.pj();
const Kokkos::View<const double*> gammaj = unknowns.gammaj();
const Kokkos::View<const Rd*> uj = unknowns.uj();
std::ofstream fout("vitesse");
fout.precision(15);
for (size_t j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << uj[j][0] + std::sqrt((gammaj[j]*pj[j])/rhoj[j]) << '\n';
}
}
}
Kokkos::finalize();
......
src/scheme/FiniteVolumesDiffusion.hpp
+ 4
− 2
View file @ e11acce1
......@@ -419,7 +419,7 @@ public:
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
// double pi = 4.*std::atan(1.);
//double h = std::sqrt(1. - (t*t)/(50./9.));
double h = std::sqrt(1. - (t*t)/(50./9.));
// CL en diffusion pure
// TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
......@@ -502,6 +502,7 @@ public:
//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 = xj[j][0];
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
}
err_u = std::sqrt(err_u);
......@@ -541,7 +542,8 @@ public:
double exact_E = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
//exact_E = (-(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.;
exact_E = ((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.*t)-1.) + 2.;
//exact_E = ((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.*t)-1.) + 2.;
exact_E = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1. + t;
err_E += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j);
}
err_E = std::sqrt(err_E);
......
src/scheme/FiniteVolumesEulerUnknowns.hpp
+ 19
− 10
View file @ e11acce1
......@@ -275,8 +275,6 @@ public:
return m_S0;
}
// --- Acoustic Solver ---
void initializeSod()
......@@ -336,7 +334,7 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_kj[j] = xj[j][0];
m_kj[j] = 0.;
//m_kj[j] = 0.;
/*
// Sod
......@@ -373,13 +371,13 @@ public:
m_kj[j]=0. ;
}
}
*/
// k regulier
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] = 0.014*m_kj[j];
*/
});
......@@ -396,8 +394,19 @@ public:
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]);
//m_nuj(j) = 0.5;
//m_nuj(j) = 0.5*(1.+xj[j][0]);
m_nuj(j) = 0.0005;
/*
if (xj[j][0]<0.7) {
m_nuj[j]=0.;
} else {
if (xj[j][0]<0.9){
m_nuj[j]=0.5;
} else {
m_nuj[j]=0. ;
}
}
*/
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
......@@ -408,8 +417,8 @@ public:
m_TL[0] = m_ej[0];
m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
m_nuL[0] = 0.5;
m_nuR[0] = 1.;
m_nuL[0] = 0.0005;
m_nuR[0] = 0.0005;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment