From 80c28da039c12ab88140d49dc454570e90dc4c22 Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Mon, 4 Jun 2018 14:51:26 +0200
Subject: [PATCH] ajout autre approche du splitting
---
src/main.cpp | 23 +++++++++++++++++++++--
src/scheme/FiniteVolumesDiffusion.hpp | 22 ++++++++++++----------
src/scheme/FiniteVolumesEulerUnknowns.hpp | 8 ++++----
3 files changed, 37 insertions(+), 16 deletions(-)
diff --git a/src/main.cpp b/src/main.cpp
index 7065306c5..a25a9b1b8 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -180,7 +180,7 @@ int main(int argc, char *argv[])
// ETAPE 1 DU SPLITTING - EULER
- double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt_euler = 0.1*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt_euler > tmax) {
dt_euler = tmax-t;
@@ -192,7 +192,7 @@ int main(int argc, char *argv[])
double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
double t_diff = t-dt_euler;
-
+
if (dt_euler <= dt_diff) {
dt_diff = dt_euler;
finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
@@ -207,6 +207,25 @@ int main(int argc, char *argv[])
}
}
+
+ // AUTRE APPROCHE DU SPLITTING (PLUS LONG)
+ /*
+ double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
+ double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
+ double dt = 0.;
+ if (dt_euler < dt_diff) {
+ dt = dt_euler;
+ } else {
+ dt = dt_diff;
+ }
+ if (t+dt > tmax) {
+ dt = tmax-t;
+ }
+ acoustic_solver.computeNextStep(t,dt, unknowns);
+ finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
+ t += dt;
+ */
+
// DIFFUSION PURE
/*
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 0bfa990ca..ab50c2ea4 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -129,16 +129,16 @@ private:
// Kidder
// k = 0.5
- //m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
- //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5;
+ m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
+ m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5;
// 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];
+ //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 ;
}
@@ -346,11 +346,11 @@ public:
//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 = 0.5)
- //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+ 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)));
+ //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)));
});
@@ -362,7 +362,8 @@ public:
// Mise a jour de k
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- kj(j) = xj[j][0];
+ //kj(j) = xj[j][0];
+ kj(j) = 0.5;
/*
if (xj[j][0]<0.7) {
kj[j]=0.;
@@ -534,7 +535,8 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
entropy_new(j) = std::log(pj(j)*std::pow(rhoj(j),-gammaj(j)));
if (entropy_new(j) - entropy(j) < 0) {
- std::cout << "ATTENTION ENTROPIE NEGATIVE !";
+ std::cout << "ATTENTION ENTROPIE NEGATIVE !" << std::endl;
+ std::cout << "j = " << j << " difference = " << entropy_new(j) - entropy(j) << std::endl;
}
entropy(j) = std::log(pj(j)*std::pow(rhoj(j),-gammaj(j)));
});
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index a5455c799..f9bccfe27 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -261,8 +261,8 @@ 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] = xj[j][0];
+ m_kj[j] = 0.5;
/*
if (xj[j][0]<0.7) {
m_kj[j]=0.;
@@ -280,8 +280,8 @@ public:
m_uL[0] = zero;
m_uR[0] = zero;
- m_kL[0] = 0.;
- m_kR[0] = 1.;
+ m_kL[0] = 0.5;
+ m_kR[0] = 0.5;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));
--
GitLab