From 6f103b40b32b6fb408ba107b87d37def3cdebdca Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Tue, 17 Jul 2018 11:30:42 +0200
Subject: [PATCH] test kidder avec conduction thermique sans viscosite
cinematique
---
src/main.cpp | 22 +++++++-------
src/scheme/AcousticSolver.hpp | 13 +++++----
src/scheme/FiniteVolumesDiffusion.hpp | 24 ++++++++++------
src/scheme/FiniteVolumesEulerUnknowns.hpp | 35 +++++++++++------------
4 files changed, 50 insertions(+), 44 deletions(-)
diff --git a/src/main.cpp b/src/main.cpp
index afe1896fd..011cf8f72 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -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,8 +291,8 @@ int main(int argc, char *argv[])
t_diff += dt_diff;
}
}
- */
+ /*
// DIFFUSION PURE
double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj);
@@ -301,7 +301,7 @@ int main(int argc, char *argv[])
}
finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
t += dt;
-
+ */
block_eos.updatePandCFromRhoE();
@@ -595,11 +595,11 @@ int main(int argc, char *argv[])
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("rho1");
+ std::ofstream fout("rho");
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';
}
}
@@ -613,8 +613,8 @@ int main(int argc, char *argv[])
//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] << ' ' << xj[j][0] << std::endl;
+ fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
+ //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << xj[j][0] << std::endl;
//fout << xj[j][0] << ' ' << uj[j][0] << '\n';
}
@@ -631,8 +631,8 @@ int main(int argc, char *argv[])
//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]*xj[j][0]*0.5 + 2.*xj[j][0] + tmax + 1. << std::endl;
+ 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]*xj[j][0]*0.5 + 2.*xj[j][0] + tmax + 1. << std::endl;
//fout << xj[j][0] << ' ' << Ej[j] << '\n';
}
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index e4ed0d0c8..84c9168b1 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -184,15 +184,15 @@ private:
//m_ur[0]=zero;
//m_ur[m_mesh.numberOfNodes()-1]=zero;
- m_ur[0] = x0;
- m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
+ //m_ur[0] = x0;
+ //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
// CL Kidder
- /*
+
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;
}
@@ -371,12 +371,13 @@ public:
});
*/
+
// Mise a jour de nu
- /*
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
nuj(j) = 0.5*(1.+xj[j][0]);
});
- */
+
}
};
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index b2a714ecc..0c832568e 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -113,7 +113,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)));
@@ -121,7 +121,7 @@ private:
cell_here = face_cells(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)));
-
+ */
// Kidder
@@ -133,6 +133,10 @@ private:
//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];
+
+ // k = 0
+ m_Fl(0,0) = 0.;
+ m_Fl(m_mesh.numberOfFaces()-1,0) = 0.;
return m_Fl ;
@@ -230,19 +234,20 @@ private:
});
// Conditions aux bords
-
+ /*
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));
+ */
// Kiddder
// nu = (1+x)*0.5
- //double h = std::sqrt(1. - (t*t)/(50./9.));
- //m_Bl(0) = ((1.+h*x0[0][0])*3.*h*x0[0][0])/(100.*h*h*h*h);
- //m_Bl(m_mesh.numberOfFaces()-1) = ((1.+h*xmax[0][0])*3.*h*xmax[0][0])/(100.*h*h*h*h);
+ double h = std::sqrt(1. - (t*t)/(50./9.));
+ m_Bl(0) = ((1.+h*x0[0][0])*3.*h*x0[0][0])/(100.*h*h*h*h);
+ m_Bl(m_mesh.numberOfFaces()-1) = ((1.+h*xmax[0][0])*3.*h*xmax[0][0])/(100.*h*h*h*h);
// nu = 0
@@ -414,12 +419,12 @@ 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.);
- TL(0) = 1.+ t;
- TR(0) = 3. + t;
+ //TL(0) = 1.+ t;
+ //TR(0) = 3. + t;
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
@@ -464,6 +469,7 @@ public:
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); // nu = 0
//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))+(0.2*3.)/(50.*h*h*h*h)); // nu = 0.2
//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*h*h*h*h)); // nu = (1+x)*0.5
+ Ej[j] -= (dt*inv_mj[j])*Vj(j)*((6*xj[j][0]+3.)/(100*h*h*h*h));
});
// Calcul de e par la formule e = E-0.5 u^2
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 2e137b352..b0e7718bd 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -276,7 +276,7 @@ public:
}
- /*
+
// --- Acoustic Solver ---
void initializeSod()
@@ -284,27 +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.);
@@ -313,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);
@@ -336,9 +336,9 @@ public:
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
//m_kj[j] = xj[j][0];
- //m_kj[j] = 0.2;
+ m_kj[j] = 0.;
-
+ /*
// Sod
// k non regulier
@@ -379,7 +379,7 @@ public:
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 +396,8 @@ 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.5;
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -409,14 +409,13 @@ public:
m_TL[0] = m_ej[0];
m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
m_nuL[0] = 0.5;
- m_nuR[0] = 0.5;
+ m_nuR[0] = 1.;
}
- */
-
+ /*
// DIFFUSION PURE
void initializeSod()
@@ -487,7 +486,7 @@ public:
m_nuL[0] = 0.5;
m_nuR[0] = 1.;
}
-
+ */
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
--
GitLab