Select Git revision
DiscreteFunctionUtils.hpp
-
Stéphane Del Pino authoredStéphane Del Pino authored
FiniteVolumesDiffusion.hpp 15.93 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)));
// 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);
//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 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 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<Rdd*> m_Fl;
Kokkos::View<Rd*> m_Gl;
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())
{
;
}
// 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*>& 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.;
for (int m = 0; m < cell_nb_nodes(j); ++m) {
sum += kj(cell_nodes(j,m));
}
if (sum == 0.) {
dt_j[j] = Vj[j]/cj[j];
} else { // dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl;
dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*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<Rd*> Sj = unknowns.Sj();
Kokkos::View<double*> Ej = unknowns.Ej();
Kokkos::View<double*> ej = unknowns.ej();
Kokkos::View<double*> gammaj = unknowns.gammaj();
const Kokkos::View<const double*> kj = unknowns.kj();
Kokkos::View<Rd*> uL = unknowns.uL();
Kokkos::View<Rd*> uR = unknowns.uR();
Kokkos::View<double*> kL = unknowns.kL();
Kokkos::View<double*> kR = unknowns.kR();
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();
// Calcule les flux
computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, t);
const Kokkos::View<const Rdd*> Fl = m_Fl ;
const Kokkos::View<const Rd *> Gl = m_Gl ;
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) {
Rd momentum_fluxes = zero;
double energy_fluxes = 0.;
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);
energy_fluxes += (Gl(l), Cjr(j,R));
}
uj[j] += (dt*inv_mj[j]) * momentum_fluxes;
Ej[j] += (dt*inv_mj[j]) * energy_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
// 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]);
});
// met a jour les quantites (geometriques) associees au maillage
//m_mesh_data.updateAllData();
/*
// Calcul de rho avec la formule Mj = Vj rhoj
const Kokkos::View<const double*> mj = unknowns.mj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
rhoj[j] = mj[j]/Vj[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);
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(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;
}
// 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;
}
};
#endif // FINITE_VOLUMES_DIFFUSION_HPP