Skip to content
Snippets Groups Projects
Select Git revision
  • 0af5c409bdc53bf791f132fc41ba34b0103efc77
  • develop default protected
  • feature/gmsh-reader
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

Synchronizer.hpp

Blame
    • Stéphane Del Pino's avatar
      0af5c409
      Add missing synchronization capabilities · 0af5c409
      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.
      0af5c409
      History
      Add missing synchronization capabilities
      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.
    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