recommit with the right username

Author: jiajiaecho

src/PDE/CompFlow/DGCompFlow.hpp

@@ -438,6 +438,7 @@ class CompFlow {
     {
       const auto ndof = g_inputdeck.get< tag::ndof >();
       const auto rdof = g_inputdeck.get< tag::rdof >();
+      bool viscous=false;
 
       const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
 
@@ -468,8 +469,8 @@ class CompFlow {
 
       // compute internal surface flux integrals
       tk::surfInt( pref, 1, m_mat_blk, t, ndof, rdof, inpoel, solidx,
-                   coord, fd, geoFace, geoElem, m_riemann, velfn, U, P, ndofel,
-                   dt, R, riemannDeriv );
+                   coord, fd, geoFace, geoElem, m_riemann, visc_flux,
+                   velfn, U, P, ndofel, dt, R, riemannDeriv, viscous);
 
       // compute optional source term
       tk::srcInt( m_mat_blk, t, ndof, fd.Esuel().size()/4,
@@ -478,14 +479,14 @@ class CompFlow {
       if(ndof > 1)
         // compute volume integrals
         tk::volInt( 1, t, m_mat_blk, ndof, rdof,
-                    fd.Esuel().size()/4, inpoel, coord, geoElem, flux, velfn,
-                    U, P, ndofel, R );
+                    fd.Esuel().size()/4, inpoel, coord, geoElem, flux, 
+                    visc_flux, velfn, U, P, ndofel, R, viscous );
 
       // compute boundary surface flux integrals
       for (const auto& b : m_bc)
         tk::bndSurfInt( pref, 1, m_mat_blk, ndof, rdof, std::get<0>(b),
-                        fd, geoFace, geoElem, inpoel, coord, t, m_riemann,
-                        velfn, std::get<1>(b), U, P, ndofel, R, riemannDeriv );
+                        fd, geoFace, geoElem, inpoel, coord, t, m_riemann, visc_flux,
+                        velfn, std::get<1>(b), U, P, ndofel, R, riemannDeriv, viscous );
 
      // compute external (energy) sources
       const auto& ic = g_inputdeck.get< tag::ic >();
@@ -1024,6 +1025,17 @@ class CompFlow {
     {
       return {{ ul, Problem::initialize( ncomp, mat_blk, x, y, z, t ) }};
     }
+    
+static tk::FluxFn::result_type    
+   visc_flux( ncomp_t ncomp,
+          const std::vector< EOS >& mat_blk,
+          const std::vector< tk::real >& ugp,
+          const std::vector< std::array< tk::real, 3 > > & grad_all ) 
+    {
+    std::vector< std::array< tk::real, 3 > > fl( ugp.size(),
+                                             std::array<tk::real, 3 >{{0, 0, 0}}); 
+      return fl;
+    }
 
     //! \brief Boundary state function providing the left and right state of a
     //!   face at symmetry boundaries

src/PDE/Integrate/Basis.cpp

@@ -1052,3 +1052,57 @@ tk::massMatrixTaylorRefEl(std::size_t dof)
 
   return Mt;
 }
+
+std::vector< std::array< tk::real, 3 > > 
+tk::eval_state_gradient ( ncomp_t ncomp,
+                 const std::size_t ndof,
+                 const std::size_t ndof_el,
+                 const std::size_t e,
+                 const Fields& U,
+                 const std::array< std::vector<tk::real>, 3 >& dBdx )
+// *****************************************************************************
+//  Compute the state variables for the tetrahedron element
+//! \param[in] ncomp Number of scalar components in this PDE system
+//! \param[in] ndof Maximum number of degrees of freedom
+//! \param[in] ndof_el Number of degrees of freedom for the local element
+//! \param[in] e Index for the tetrahedron element
+//! \param[in] U Solution vector at recent time step
+//! \param[in] B Vector of basis functions
+//! \return Vector of state variable for tetrahedron element
+// *****************************************************************************
+{
+  // This is commented for now because that when p0/p1 adaptive with limiter
+  // applied, the size of basis will be 10. However, ndof_el will be 4 which
+  // leads to a size mismatch in limiter function.
+  //Assert( B.size() == ndof_el, "Size mismatch" );
+
+  if (U.empty()) return {};
+
+  
+  std::vector< std::array< tk::real, 3 > >  state_grad(ncomp, 
+                                            std::array< real, 3 >{{0, 0, 0}});
+
+  
+ for (std::size_t i=0; i<3; ++i) {
+   for (ncomp_t c=0; c<ncomp; ++c)
+    {
+      auto mark = c*ndof;
+      state_grad[c][i] += U( e, mark+1 ) * dBdx[i][1] 
+                       + U( e, mark+2 ) * dBdx[i][2]
+                       + U( e, mark+3 ) * dBdx[i][3];
+
+      if( ndof_el > 4 )
+    {
+      state_grad[c][i] += U( e, mark+4 ) * dBdx[i][4] 
+                       + U( e, mark+5 ) * dBdx[i][5]
+                       + U( e, mark+6 ) * dBdx[i][6]
+                       + U( e, mark+4 ) * dBdx[i][7] 
+                       + U( e, mark+5 ) * dBdx[i][8]
+                       + U( e, mark+6 ) * dBdx[i][9];
+    }
+  }
+}
+  
+ return state_grad;
+  
+}

src/PDE/Integrate/Basis.hpp

@@ -133,6 +133,14 @@ invMassMatTaylorRefEl( std::size_t dof );
 std::vector< std::vector< tk::real > >
 massMatrixTaylorRefEl(std::size_t dof);
 
+std::vector< std::array< tk::real, 3 > > 
+eval_state_gradient( ncomp_t ncomp,
+                 const std::size_t ndof,
+                 const std::size_t ndof_el,
+                 const std::size_t e,
+                 const Fields& U,
+                 const std::array< std::vector<tk::real>, 3 >& dBdx );
+
 } // tk::
 
 #endif // Basis_h

src/PDE/Integrate/Boundary.cpp

@@ -44,14 +44,16 @@ bndSurfInt( const bool pref,
             const UnsMesh::Coords& coord,
             real t,
             const RiemannFluxFn& flux,
+            const FluxFn& visc_flux,
             const VelFn& vel,
             const StateFn& state,
             const Fields& U,
             const Fields& P,
             const std::vector< std::size_t >& ndofel,
             Fields& R,
             std::vector< std::vector< tk::real > >& riemannDeriv,
-            int intsharp )
+            bool viscous,
+            int intsharp)
 // *****************************************************************************
 //! Compute boundary surface flux integrals for a given boundary type for DG
 //! \details This function computes contributions from surface integrals along
@@ -141,6 +143,8 @@ bndSurfInt( const bool pref,
 
         std::array< real, 3 >
           fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
+          
+
 
         // Gaussian quadrature
         for (std::size_t igp=0; igp<ng; ++igp)
@@ -189,6 +193,65 @@ bndSurfInt( const bool pref,
 
           // Compute the numerical flux
           auto fl = flux(mat_blk, fn, var, vel(ncomp, gp[0], gp[1], gp[2], t));
+          if (viscous){
+            std::vector< std::array< tk::real, 3 > > grad_all(2*ncomp, 
+                                            std::array< real, 3 >{{0, 0, 0}});
+            std::vector<real> fluxl(5, 0);
+            std::vector<real> fluxr(5, 0);
+ 
+            auto jacInv_l =
+          tk::inverseJacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
+           auto dBdx_l = tk::eval_dBdx_p1(dof_el, jacInv_l );
+           
+          if (dof_el > 4){
+              std::array< std::vector< real >, 3 > coordgp_3;
+              coordgp_3[0].resize( ng );
+              coordgp_3[1].resize( ng );
+              coordgp_3[2].resize( ng );
+              for(int i=0; i<3; i++)
+                coordgp_3[i][igp]=ref_gp_l[i];
+              eval_dBdx_p2( igp, coordgp_3, jacInv_l, dBdx_l );
+          }
+            
+            auto state_U_grad_l = eval_state_gradient (ncomp, ndof, dof_el, 
+                         el, U, dBdx_l );
+            auto state_P_grad_l = eval_state_gradient (nprim, ndof, dof_el, 
+                         el, P, dBdx_l );
+            for (ncomp_t c=0; c<ncomp; ++c){
+              grad_all.push_back(state_U_grad_l[c]);
+             }
+
+            for (ncomp_t c=0; c<ncomp; ++c){
+              grad_all.push_back(state_P_grad_l[c]);
+            } 
+            auto fl_vis_l = visc_flux(ncomp, mat_blk, var[0], grad_all); 
+
+            auto fl_vis_r = visc_flux(ncomp, mat_blk, var[1], grad_all) ; 
+      
+ 
+
+           // Flux functions
+           for (std::size_t i=0; i<3; ++i)
+            {
+             fluxl[1] = fluxl[1] + fl_vis_l[1][i]*fn[i];
+             fluxl[2] = fluxl[2] + fl_vis_l[2][i]*fn[i];
+             fluxl[3] = fluxl[3] + fl_vis_l[3][i]*fn[i];
+             fluxl[4] = fluxl[4] + fl_vis_l[4][i]*fn[i];
+             fluxr[1] = fluxr[1] + fl_vis_r[1][i]*fn[i];
+             fluxr[2] = fluxr[2] + fl_vis_r[2][i]*fn[i];
+             fluxr[3] = fluxr[3] + fl_vis_r[3][i]*fn[i];
+             fluxr[4] = fluxr[4] + fl_vis_r[4][i]*fn[i];
+            }
+    
+        
+       
+           // Numerical flux function
+           for(std::size_t c=1; c<5; ++c)
+            fl[c] = fl[c]+0.5 * (fluxl[c] + fluxr[c]);
+      
+      
+      
+           }
 
           // Code below commented until details about the form of these terms in
           // the \alpha_k g_k equations are sorted out.

src/PDE/Integrate/Boundary.hpp

@@ -42,14 +42,16 @@ bndSurfInt( const bool pref,
             const UnsMesh::Coords& coord,
             real t,
             const RiemannFluxFn& flux,
+            const FluxFn& visc_flux,
             const VelFn& vel,
             const StateFn& state,
             const Fields& U,
             const Fields& P,
             const std::vector< std::size_t >& ndofel,
             Fields& R,
             std::vector< std::vector< tk::real > >& riemannDeriv,
-            int intcompr=0 );
+            bool viscous,
+            int intcompr=0);
 
 //! Update the rhs by adding the boundary surface integration term
 void

src/PDE/Integrate/Surface.cpp

@@ -44,14 +44,17 @@ surfInt( const bool pref,
          const Fields& geoFace,
          const Fields& geoElem,
          const RiemannFluxFn& flux,
+         const FluxFn& visc_flux,
          const VelFn& vel,
          const Fields& U,
          const Fields& P,
          const std::vector< std::size_t >& ndofel,
          const tk::real /*dt*/,
          Fields& R,
          std::vector< std::vector< tk::real > >& riemannDeriv,
-         int intsharp )
+         bool viscous,
+         int intsharp
+         )
 // *****************************************************************************
 //  Compute internal surface flux integrals
 //! \param[in] pref Indicator for p-adaptive algorithm
@@ -226,6 +229,89 @@ surfInt( const bool pref,
 
       // compute flux
       auto fl = flux( mat_blk, fn, state, v );
+      
+      if (viscous){
+          std::vector< tk::real> fluxl(5, 0);
+          std::vector< tk::real> fluxr(5, 0);
+          //std::array< std::vector<real>, 3 > dBdx_l, dBdx_r;
+          std::vector< std::array< tk::real, 3 > > grad_all(2*ncomp, 
+                                     std::array< real, 3 >{{0, 0, 0}});
+          auto jacInv_l =
+          tk::inverseJacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
+          auto dBdx_l = tk::eval_dBdx_p1(dof_el, jacInv_l );
+          auto jacInv_r =
+          tk::inverseJacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
+          auto dBdx_r = tk::eval_dBdx_p1(dof_er, jacInv_r );
+          
+        if (dof_el > 4){          
+          std::array< std::vector< real >, 3 > coordgp_3;
+          coordgp_3[0].resize( ng );
+          coordgp_3[1].resize( ng );
+          coordgp_3[2].resize( ng );
+          for(int i=0; i<3; i++)
+             coordgp_3[i][igp]=ref_gp_l[i];
+          eval_dBdx_p2( igp, coordgp_3, jacInv_l, dBdx_l );
+        }  
+        if (dof_er > 4){          
+          std::array< std::vector< real >, 3 > coordgp_3;
+          coordgp_3[0].resize( ng );
+          coordgp_3[1].resize( ng );
+          coordgp_3[2].resize( ng );
+          for(int i=0; i<3; i++)
+             coordgp_3[i][igp]=ref_gp_r[i];
+          eval_dBdx_p2( igp, coordgp_3, jacInv_r, dBdx_r );
+        }  
+    
+      auto state_U_grad_l = eval_state_gradient (ncomp, ndof, dof_el, 
+                         el, U, dBdx_l );
+      auto state_P_grad_l = eval_state_gradient (nprim, ndof, dof_el, 
+                         el, P, dBdx_l );
+      for (ncomp_t c=0; c<ncomp; ++c){
+          grad_all.push_back(state_U_grad_l[c]);
+      }
+
+      for (ncomp_t c=0; c<ncomp; ++c){
+          grad_all.push_back(state_P_grad_l[c]);
+      } 
+      auto fl_vis_l = visc_flux(ncomp, mat_blk, state[0], grad_all); 
+      auto state_U_grad_r = eval_state_gradient (ncomp, ndof, dof_er, 
+                         er, U, dBdx_r );
+      auto state_P_grad_r = eval_state_gradient (ncomp, ndof, dof_er, 
+                         er, P, dBdx_r );
+      grad_all.clear();
+      for (ncomp_t c=0; c<ncomp; ++c){
+          grad_all.push_back(state_U_grad_r[c]);
+      }
+
+      for (ncomp_t c=0; c<ncomp; ++c){
+          grad_all.push_back(state_P_grad_r[c]);
+      } 
+      auto fl_vis_r = visc_flux(ncomp, mat_blk, state[1], grad_all); 
+      
+ 
+
+           // Flux functions
+    for (std::size_t i=0; i<3; ++i)
+    {
+    fluxl[1] = fluxl[1] + fl_vis_l[1][i]*fn[i];
+    fluxl[2] = fluxl[2] + fl_vis_l[2][i]*fn[i];
+    fluxl[3] = fluxl[3] + fl_vis_l[3][i]*fn[i];
+    fluxl[4] = fluxl[4] + fl_vis_l[4][i]*fn[i];
+    fluxr[1] = fluxr[1] + fl_vis_r[1][i]*fn[i];
+    fluxr[2] = fluxr[2] + fl_vis_r[2][i]*fn[i];
+    fluxr[3] = fluxr[3] + fl_vis_r[3][i]*fn[i];
+    fluxr[4] = fluxr[4] + fl_vis_r[4][i]*fn[i];
+    }
+    
+        
+       
+    // Numerical flux function
+    for(std::size_t c=1; c<5; ++c)
+      fl[c] = fl[c]+0.5 * (fluxl[c] + fluxr[c]);
+      
+      
+      
+      }
 
       // Add the surface integration term to the rhs
       update_rhs_fa( ncomp, nmat, ndof, ndofel[el], ndofel[er], wt, fn,

src/PDE/Integrate/Surface.hpp

@@ -42,14 +42,16 @@ surfInt( const bool pref,
          const Fields& geoFace,
          const Fields& geoElem,
          const RiemannFluxFn& flux,
+         const FluxFn& visc_flux,
          const VelFn& vel,
          const Fields& U,
          const Fields& P,
          const std::vector< std::size_t >& ndofel,
          const tk::real dt,
          Fields& R,
          std::vector< std::vector< tk::real > >& riemannDeriv,
-         int intcompr=0 );
+         bool viscous,
+         int intcompr=0);
 
 // Update the rhs by adding surface integration term
 void