Select Git revision
Synchronizer.hpp
-
Stéphane Del Pino authored
In the case of subitems of greater dimension that items, synchronization was missing for SubItemValuePerItem and SubItemArrayPerItem. This has been missing for ages but it not that useful. Actually, it is more of a completeness functionality. Related tests have been added.
Stéphane Del Pino authoredIn the case of subitems of greater dimension that items, synchronization was missing for SubItemValuePerItem and SubItemArrayPerItem. This has been missing for ages but it not that useful. Actually, it is more of a completeness functionality. Related tests have been added.
FiniteVolumesDiffusion.hpp 20.68 KiB
#ifndef FINITE_VOLUMES_DIFFUSION_HPP
#define FINITE_VOLUMES_DIFFUSION_HPP
// --- INCLUSION fichiers headers ---
#include <Kokkos_Core.hpp>
#include <rang.hpp>
#include <BlockPerfectGas.hpp>
#include <TinyVector.hpp>
#include <TinyMatrix.hpp>
#include <Mesh.hpp>
#include <MeshData.hpp>
#include <FiniteVolumesEulerUnknowns.hpp>
// ---------------------------------
// Creation classe FiniteVolumesDiffusion
template<typename MeshData>
class FiniteVolumesDiffusion
{
typedef typename MeshData::MeshType MeshType;
typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType;
MeshData& m_mesh_data;
const MeshType& m_mesh;
const typename MeshType::Connectivity& m_connectivity;
constexpr static size_t dimension = MeshType::dimension;
typedef TinyVector<dimension> Rd;
typedef TinyMatrix<dimension> Rdd;
private:
// Sert a calculer les reductions (en gros calculer le min sur des
// vecteurs en parallele)
struct ReduceMin
{
private:
const Kokkos::View<const double*> x_;
public:
typedef Kokkos::View<const double*>::non_const_value_type value_type;
ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {}
typedef Kokkos::View<const double*>::size_type size_type;
KOKKOS_INLINE_FUNCTION void
operator() (const size_type i, value_type& update) const
{
if (x_(i) < update) {
update = x_(i);
}
}
KOKKOS_INLINE_FUNCTION void
join (volatile value_type& dst,
const volatile value_type& src) const
{
if (src < dst) {
dst = src;
}
}
KOKKOS_INLINE_FUNCTION void
init (value_type& dst) const
{ // The identity under max is -Inf.
dst= Kokkos::reduction_identity<value_type>::min();
}
};
// Calcule un Fl
Kokkos::View<Rdd*>
computeFl(const Kokkos::View<const Rd**>& Cjr,
const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& kj,
const Kokkos::View<const Rd*>& uL,
const Kokkos::View<const Rd*>& uR,
const Kokkos::View<const double*>& kL,
const Kokkos::View<const double*>& kR,
const double& t) {
const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
const Kokkos::View<const unsigned short*> face_nb_cells
= m_connectivity.faceNbCells();
const Kokkos::View<const unsigned short**> face_cell_local_face
= m_mesh.connectivity().faceCellLocalFace();
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd*> xr = m_mesh.xr();
const Kokkos::View<const Rd*> x0 = m_mesh.x0();
const Kokkos::View<const Rd*> xmax = m_mesh.xmax();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rdd sum = zero;
double sum2 = 0.;
double sum3 = 0.;
for (int j=0; j<face_nb_cells(l); ++j) {
int cell_here = face_cells(l,j);
int local_face_number_in_cell = face_cell_local_face(l,j);
sum -= tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell));
sum2 += (1./Vj(cell_here))*kj(cell_here);
sum3 += 1./Vj(cell_here);
}
m_Fl(l) = ((sum2/sum3)/Vl(l))*sum;
});
// 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)));
//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);
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) = -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)));
// 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;
// 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];
*/
return m_Fl ;
}
// Calcule un Gl
Kokkos::View<Rd*>
computeGl(const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const Rdd*>& Fl,
const Kokkos::View<const Rd*>& uL,
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 short*> face_nb_cells
= m_connectivity.faceNbCells();
const Kokkos::View<const unsigned short**> face_cell_local_face
= m_mesh.connectivity().faceCellLocalFace();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd*> xr = m_mesh.xr();
const Kokkos::View<const Rd*> x0 = m_mesh.x0();
const Kokkos::View<const Rd*> xmax = m_mesh.xmax();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rd sum = zero;
double sum2 = 0.;
for (int j=0; j<face_nb_cells(l); ++j) {
int cell_here = face_cells(l,j);
sum += (1./Vj(cell_here))*uj(cell_here);
sum2 += 1./Vj(cell_here);
}
m_Gl(l) = (1./sum2)*Fl(l)*sum;
});
// Conditions aux bords
m_Gl(0) = Fl(0)*uL(0);
m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
// Kidder
//m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
//m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
/*
double h = std::sqrt(1. - (t*t)/(50./9.));
m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0);
m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0);
*/
return m_Gl ;
}
// Calcule un Bl
Kokkos::View<double*>
computeBl(const Kokkos::View<const Rd**>& Cjr,
const Kokkos::View<const double*>& Tj,
const Kokkos::View<const double*>& nuj,
const Kokkos::View<const double*>& TL,
const Kokkos::View<const double*>& TR,
const Kokkos::View<const double*>& nuL,
const Kokkos::View<const double*>& nuR,
const double& t) {
const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
const Kokkos::View<const unsigned short*> face_nb_cells
= m_connectivity.faceNbCells();
const Kokkos::View<const unsigned short**> face_cell_local_face
= m_mesh.connectivity().faceCellLocalFace();
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
Rd sum = zero;
double sum2 = 0.;
double sum3 = 0.;
for (int j=0; j<face_nb_cells(l); ++j) {
int cell_here = face_cells(l,j);
int local_face_number_in_cell = face_cell_local_face(l,j);
sum -= Tj(cell_here)*Cjr(cell_here, local_face_number_in_cell);
sum2 += (1./Vj(cell_here))*nuj(cell_here);
sum3 += 1./Vj(cell_here);
}
m_Bl(l) = ((sum2/sum3)/Vl(l))*sum[0];
});
// 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));
return m_Bl ;
}
// Calcule la liste des inverse d'une liste de matrices
// (pour l'instant, juste 1x1)
void inverse(const Kokkos::View<const Rdd*>& A,
Kokkos::View<Rdd*>& inv_A) const {
Kokkos::parallel_for(A.size(), KOKKOS_LAMBDA(const int& r) {
inv_A(r) = Rdd{1./(A(r)(0,0))};
});
}
// Calcule la liste des inverses d'une liste de reels
void inverse(const Kokkos::View<const double*>& x,
Kokkos::View<double*>& inv_x) const {
Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
inv_x(r) = 1./x(r);
});
}
// Enchaine les operations pour calculer les flux (Fjr et Gjr) pour
// pouvoir derouler le schema
KOKKOS_INLINE_FUNCTION
void computeExplicitFluxes(const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const Rd**>& Cjr,
const Kokkos::View<const double*>& kj,
const Kokkos::View<const Rd*>& uL,
const Kokkos::View<const Rd*>& uR,
const Kokkos::View<const double*>& kL,
const Kokkos::View<const double*>& kR,
const Kokkos::View<const double*>& Tj,
const Kokkos::View<const double*>& nuj,
const Kokkos::View<const double*>& TL,
const Kokkos::View<const double*>& TR,
const Kokkos::View<const double*>& nuL,
const Kokkos::View<const double*>& nuR,
const double& t) {
Kokkos::View<Rdd*> Fl = m_Fl ;
Fl = computeFl (Cjr, uj, kj, uL, uR, kL, kR, t);
Kokkos::View<Rd*> Gl = m_Gl ;
Gl = computeGl (uj, Fl, uL, uR, t);
Kokkos::View<double*> Bl = m_Bl ;
Bl = computeBl (Cjr, Tj, nuj, TL, TR, nuL, nuR, t);
}
Kokkos::View<Rdd*> m_Fl;
Kokkos::View<Rd*> m_Gl;
Kokkos::View<double*> m_Bl;
public:
FiniteVolumesDiffusion(MeshData& mesh_data,
UnknownsType& unknowns)
: m_mesh_data(mesh_data),
m_mesh(mesh_data.mesh()),
m_connectivity(m_mesh.connectivity()),
m_Fl("Fl", m_mesh.numberOfFaces()),
m_Gl("Gl", m_mesh.numberOfFaces()),
m_Bl("Bl", m_mesh.numberOfFaces())
{
;
}
// Calcule une evaluation du pas de temps verifiant le CFL parabolique
// Utilise la reduction definie dans la structure ReduceMin.
KOKKOS_INLINE_FUNCTION
double diffusion_dt(const Kokkos::View<const double*>& rhoj,
const Kokkos::View<const double*>& kj,
const Kokkos::View<const double*>& nuj,
const Kokkos::View<const double*>& cj) const {
Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells());
const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
const Kokkos::View<const unsigned int**>& cell_faces = m_connectivity.cellFaces();
const Kokkos::View<const unsigned short*> cell_nb_faces
= m_connectivity.cellNbFaces();
const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
const Kokkos::View<const unsigned short*> cell_nb_nodes
= m_connectivity.cellNbNodes();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
double minVl = std::numeric_limits<double>::max();
for (int ll=0; ll<cell_nb_faces(j); ++ll) {
minVl = std::min(minVl, Vl(cell_faces(j, ll)));
}
double sum = 0.;
double sum1 = 0.;
for (int m = 0; m < cell_nb_nodes(j); ++m) {
sum += kj(cell_nodes(j,m));
sum1 += nuj(cell_nodes(j,m));
}
if (sum > sum1) {
if (sum == 0.) {
dt_j[j] = std::numeric_limits<double>::max();
} else {
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl;
}
} else {
if (sum1 == 0.) {
dt_j[j] = std::numeric_limits<double>::max();
} else {
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum1)*minVl;
}
}
});
double dt = std::numeric_limits<double>::max();
Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(dt_j), dt);
return dt;
}
// Avance la valeur des inconnues pendant un pas de temps dt
void computeNextStep(const double& t, const double& dt,
UnknownsType& unknowns)
{
Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<Rd*> uj = unknowns.uj();
Kokkos::View<double*> Ej = unknowns.Ej();
Kokkos::View<double*> Tj = unknowns.Tj();
Kokkos::View<double*> ej = unknowns.ej();
Kokkos::View<double*> gammaj = unknowns.gammaj();
Kokkos::View<double*> kj = unknowns.kj();
Kokkos::View<double*> nuj = unknowns.nuj();
Kokkos::View<Rd*> uL = unknowns.uL();
Kokkos::View<Rd*> uR = unknowns.uR();
Kokkos::View<double*> kL = unknowns.kL();
Kokkos::View<double*> kR = unknowns.kR();
Kokkos::View<double*> TL = unknowns.TL();
Kokkos::View<double*> TR = unknowns.TR();
Kokkos::View<double*> nuL = unknowns.nuL();
Kokkos::View<double*> nuR = unknowns.nuR();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
double pi = 4.*std::atan(1.);
TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.);
// la condition au bord a droite de T depend du temps
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
const Kokkos::View<const Rdd*> Fl = m_Fl ;
const Kokkos::View<const Rd *> Gl = m_Gl ;
const Kokkos::View<const double*> Bl = m_Bl ;
const Kokkos::View<const unsigned short*>& cell_nb_faces
= m_connectivity.cellNbFaces();
const Kokkos::View<const unsigned int*[2]>& cell_faces
= m_connectivity.cellFaces();
// Mise a jour de la vitesse et de l'energie totale specifique
const Kokkos::View<const double*> inv_mj = unknowns.invMj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
//const int j = j0+1;
Rd momentum_fluxes = zero;
double energy_fluxes = 0.;
Rd trich = zero;
int l = 0;
for (int R=0; R<cell_nb_faces(j); ++R) {
l = cell_faces(j,R);
momentum_fluxes += Fl(l)*Cjr(j,R);
trich = Bl(l)*Cjr(j,R);
energy_fluxes += (Gl(l), Cjr(j,R)) + trich[0];
}
uj[j] += (dt*inv_mj[j]) * momentum_fluxes;
Ej[j] += (dt*inv_mj[j]) * energy_fluxes;
// test k non cst pour diff pure
uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi);
Ej[j] -= ((pi*pi*0.5*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) + xj[j][0]*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]))*(std::exp(-2.*t)-1.) - pi*0.5*std::exp(-2.*t)*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0]) + (1.+xj[j][0])*((3.*pi*pi*pi*std::cos(pi*xj[j][0])*std::sin(pi*xj[j][0])-xj[j][0]*pi*pi*pi*pi*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])))*(std::exp(-2.*t)-1.)+0.5*pi*pi*std::exp(-2.*t)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))))*(dt*inv_mj[j])*Vj(j);
// 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)));
// 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)));
});
// Calcul de e par la formule e = E-0.5 u^2
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
});
// Calcul de T par la formule T = E-0.5 u^2
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
Tj[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
});
}
// Calcul erreur entre solution analytique et solution numerique en norme L2
// (quand la solution exacte est connue)
double error_L2_u(UnknownsType& unknowns, const double& t) {
Kokkos::View<Rd*> uj = unknowns.uj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
//double pi = 4.*std::atan(1.);
double err_u = 0.;
double exact_u = 0.;
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(-0.2); // solution exacte cas test k non constant
exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
}
err_u = std::sqrt(err_u);
return err_u;
}
double error_L2_rho(UnknownsType& unknowns, const double& t) {
Kokkos::View<double*> rhoj = unknowns.rhoj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
double h = std::sqrt(1. - (t*t)/(50./9.));
double err_rho = 0.;
double exact_rho = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
exact_rho = std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h;
err_rho += (exact_rho - rhoj[j])*(exact_rho - rhoj[j])*Vj(j);
}
err_rho = std::sqrt(err_rho);
return err_rho;
}
/*
double error_L2_E(UnknownsType& unknowns) {
Kokkos::View<double*> Ej = unknowns.Ej();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
double pi = 4.*std::atan(1.);
double err_E = 0.;
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.*0.2)-1.) + 2.;
err_E += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j);
}
err_E = std::sqrt(err_E);
return err_E;
}
*/
// Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
// (quand la solution exacte est connue)
double error_Linf_rho(UnknownsType& unknowns, const double& t) {
Kokkos::View<double*> rhoj = unknowns.rhoj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
//double pi = 4.*std::atan(1.);
double h = std::sqrt(1. - (t*t)/(50./9.));
double exacte = std::sqrt((4.*((xj[0][0]*xj[0][0])/(h*h)) + 100.-(xj[0][0]*xj[0][0])/(h*h))/100.)/h;
double erreur = std::abs(exacte - rhoj[0]);
for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
exacte = std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h;
if (std::abs(exacte - rhoj[j]) > erreur) {
erreur = std::abs(exacte - rhoj[j]);
}
}
return erreur;
}
// Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
// (quand la solution exacte est connue)
double error_Linf_u(UnknownsType& unknowns, const double& t) {
Kokkos::View<Rd*> uj = unknowns.uj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
double exacte = -(xj[0][0]*t)/((50./9.)-t*t);
double erreur = std::abs(exacte - uj[0][0]);
for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
exacte = -(xj[j][0]*t)/((50./9.)-t*t);
if (std::abs(exacte - uj[j][0]) > erreur) {
erreur = std::abs(exacte - uj[j][0]);
}
}
return erreur;
}
// Verifie la conservativite de E
double conservatif(UnknownsType& unknowns) {
Kokkos::View<double*> Ej = unknowns.Ej();
Kokkos::View<double*> rhoj = unknowns.rhoj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
double sum = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
sum += Ej[j]*rhoj[j]*Vj[j];
}
return sum;
}
// Verifie la conservativite de quantite de mvt
double conservatif_mvt(UnknownsType& unknowns) {
Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<Rd*> uj = unknowns.uj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
double sum = 0.;
for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
sum += Vj[j]*rhoj[j]*uj[j][0];
}
return sum;
}
// Teste la croissance de l entropie d un gaz parfait
void entropie(UnknownsType& unknowns) {
Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<double*> pj = unknowns.pj();
Kokkos::View<double*> gammaj = unknowns.gammaj();
Kokkos::View<double*> entropy = unknowns.entropy();
Kokkos::View<double*> entropy_new = unknowns.entropy_new();
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::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)));
});
}
};
#endif // FINITE_VOLUMES_DIFFUSION_HPP