Semi-BDF Time Scheme

examples/ns/CMakeLists.txt

@@ -112,7 +112,9 @@ SET( perigee_analysis_lib_src
   ${perigee_source}/Solver/PLinear_Solver_PETSc.cpp
   ${perigee_SOURCE_DIR}/src/PLocAssem_VMS_NS_GenAlpha.cpp
   ${perigee_SOURCE_DIR}/src/PLocAssem_VMS_NS_GenAlpha_WeakBC.cpp
+  ${perigee_SOURCE_DIR}/src/PLocAssem_VMS_NS_SemiBDF1.cpp
   ${perigee_SOURCE_DIR}/src/PGAssem_NS_FEM.cpp
+  ${perigee_SOURCE_DIR}/src/PGAssem_NS_SemiBDF1.cpp
   ${perigee_SOURCE_DIR}/src/PDNSolution_NS.cpp
   ${perigee_SOURCE_DIR}/src/PNonlinear_NS_Solver.cpp
   ${perigee_SOURCE_DIR}/src/PTime_NS_Solver.cpp
@@ -127,6 +129,7 @@ SET( perigee_postprocess_lib_src
   ${perigee_source}/Postproc_Tool/Interpolater.cpp
   ${perigee_SOURCE_DIR}/src/VisDataPrep_NS.cpp
   ${perigee_SOURCE_DIR}/src/VTK_Writer_NS.cpp
+  ${perigee_SOURCE_DIR}/src/Post_error_Cubic.cpp
   )
 
 INCLUDE_DIRECTORIES( ../../include )
@@ -155,21 +158,25 @@ INCLUDE(../../conf/message_compiler_setup.cmake)
 # ===================================================================
 ADD_EXECUTABLE( preprocess3d preprocess_ns.cpp)
 ADD_EXECUTABLE( ns3d ns_driver.cpp)
+ADD_EXECUTABLE( ns3d_semiBDF1 ns_semiBDF1_driver.cpp)
 ADD_EXECUTABLE( prepost3d prepost_ns.cpp)
 ADD_EXECUTABLE( vis_ns vis_ns.cpp)
 ADD_EXECUTABLE( vis_p1_wss vis_p1_wss.cpp)
 ADD_EXECUTABLE( vis_p2_wss vis_p2_wss.cpp)
 ADD_EXECUTABLE( vis_601_p1_wss vis_601_p1_wss.cpp)
 ADD_EXECUTABLE( vis_602_p2_wss vis_602_p2_wss.cpp)
+ADD_EXECUTABLE( post_comp post_compare_manu.cpp)
 
 TARGET_LINK_LIBRARIES( preprocess3d perigee_preprocess )
 TARGET_LINK_LIBRARIES( ns3d perigee_analysis )
+TARGET_LINK_LIBRARIES( ns3d_semiBDF1 perigee_analysis )
 TARGET_LINK_LIBRARIES( prepost3d perigee_preprocess )
 TARGET_LINK_LIBRARIES( vis_ns perigee_postprocess )
 TARGET_LINK_LIBRARIES( vis_p1_wss perigee_postprocess )
 TARGET_LINK_LIBRARIES( vis_p2_wss perigee_postprocess )
 TARGET_LINK_LIBRARIES( vis_601_p1_wss perigee_postprocess )
 TARGET_LINK_LIBRARIES( vis_602_p2_wss perigee_postprocess )
+TARGET_LINK_LIBRARIES( post_comp perigee_postprocess )
 
 if(OPENMP_CXX_FOUND)
   SET_TARGET_PROPERTIES( perigee_preprocess PROPERTIES COMPILE_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS} -DUSE_OPENMP" )

examples/ns/include/PDNSolution_NS.hpp

@@ -41,8 +41,14 @@ class PDNSolution_NS : public PDNSolution
     void Init_pipe_parabolic( const APart_Node * const &pNode_ptr,
         const FEANode * const &fNode_ptr );
 
-  private:
-    const bool is_print;
+    // case 3: generate flow parabolic for an arbitrary inlet face
+    // with v = v_max * (1 - r^3 / R^3), v_max = 5/3 * v_ave
+    void Init_flow_cubic( const APart_Node *const &pNode_ptr,
+        const FEANode *const &fNode_ptr,
+        const ALocal_InflowBC *const &infbc );
+
+private:
+  const bool is_print;
 };
 
 #endif

examples/ns/include/PGAssem_NS_SemiBDF1.hpp

@@ -0,0 +1,231 @@
+#ifndef PGASSEM_NS_SEMIBDF1_HPP
+#define PGASSEM_NS_SEMIBDF1_HPP
+// ==================================================================
+// PGAssem_NS_SemiBDF1.hpp
+//
+// Parallel golbal assembly based on PETSc, using AIJ matrix format.
+// The assembly routine is designed for classical C0 FEM method, which
+// means we do not need extraction operators and local mesh sizes.
+//
+// This assembly is for the NS Equations written in VMS formulation.
+// The time scheme is Semi-BDF. 
+// The input solution vectors contains
+//  [ pressure; velocity ].
+//
+// Author: Chi Ding 
+// Date  : Mar 11, 2024
+// ==================================================================
+#include "IPGAssem.hpp"
+#include "PETSc_Tools.hpp"
+#include "PDNSolution_NS.hpp"
+
+class PGAssem_NS_SemiBDF1 : public IPGAssem
+{
+  public:
+    // Constructor for NS equations
+    PGAssem_NS_SemiBDF1( 
+        IPLocAssem * const &locassem_ptr,
+        FEAElement * const &elements,
+        const IQuadPts * const &quads,
+        const IAGlobal_Mesh_Info * const &agmi_ptr,
+        const ALocal_Elem * const &alelem_ptr,
+        const ALocal_IEN * const &aien_ptr,
+        const APart_Node * const &pnode_ptr,
+        const ALocal_NBC * const &part_nbc,
+        const ALocal_EBC * const &part_ebc,
+        const int &in_nz_estimate=60 );
+
+    // Destructor
+    virtual ~PGAssem_NS_SemiBDF1();
+
+    // Nonzero pattern estimate for the NS equations
+    virtual void Assem_nonzero_estimate(
+        const ALocal_Elem * const &alelem_ptr,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &elements,
+        const IQuadPts * const &quad_s,
+        const ALocal_IEN * const &lien_ptr,
+        const APart_Node * const &node_ptr,
+        const ALocal_NBC * const &nbc_part,
+        const ALocal_EBC * const &ebc_part);
+
+    // Assem mass matrix and residual vector
+    virtual void Assem_mass_residual(
+        const PDNSolution * const &sol,
+        const ALocal_Elem * const &alelem_ptr,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &elementv,
+        FEAElement * const &elements,
+        const IQuadPts * const &quad_v,
+        const IQuadPts * const &quad_s,
+        const ALocal_IEN * const &lien_ptr,
+        const FEANode * const &fnode_ptr,
+        const ALocal_NBC * const &nbc_part,
+        const ALocal_EBC * const &ebc_part );
+
+    // Assembly the residual vector for the NS equations
+    virtual void Assem_residual(
+        const PDNSolution * const &sol_0,
+        const PDNSolution * const &sol,
+        const double &curr_time,
+        const double &dt,
+        const ALocal_Elem * const &alelem_ptr,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &elementv,
+        FEAElement * const &elements,
+        const IQuadPts * const &quad_v,
+        const IQuadPts * const &quad_s,
+        const ALocal_IEN * const &lien_ptr,
+        const FEANode * const &fnode_ptr,
+        const ALocal_NBC * const &nbc_part,
+        const ALocal_EBC * const &ebc_part );
+
+    // Assembly the residual vector and tangent matrix 
+    // for the NS equations
+    virtual void Assem_tangent_residual(
+        const PDNSolution * const &sol_0,
+        const PDNSolution * const &sol,
+        const double &curr_time,
+        const double &dt,
+        const ALocal_Elem * const &alelem_ptr,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &elementv,
+        FEAElement * const &elements,
+        const IQuadPts * const &quad_v,
+        const IQuadPts * const &quad_s,
+        const ALocal_IEN * const &lien_ptr,
+        const FEANode * const &fnode_ptr,
+        const ALocal_NBC * const &nbc_part,
+        const ALocal_EBC * const &ebc_part );
+
+    // Assembly routine for the surface integrals of flow rate and
+    // pressure
+    virtual double Assem_surface_flowrate(
+        const PDNSolution * const &sol,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &element_s,
+        const IQuadPts * const &quad_s,
+        const ALocal_EBC * const &ebc_part,
+        const int &ebc_id );
+
+    virtual double Assem_surface_flowrate(
+        const PDNSolution * const &sol,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &element_s,
+        const IQuadPts * const &quad_s,
+        const ALocal_InflowBC * const &infbc_part,
+        const int &nbc_id );
+
+    virtual double Assem_surface_ave_pressure(
+        const PDNSolution * const &sol,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &element_s,
+        const IQuadPts * const &quad_s,
+        const ALocal_EBC * const &ebc_part,
+        const int &ebc_id );
+
+    virtual double Assem_surface_ave_pressure(
+        const PDNSolution * const &sol,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &element_s,
+        const IQuadPts * const &quad_s,
+        const ALocal_InflowBC * const &infbc_part,
+        const int &nbc_id );
+
+  private:
+    // Private data
+    const int nLocBas, dof_sol, dof_mat, num_ebc, nlgn;
+
+    int snLocBas;
+
+    // Private function
+    // Essential boundary condition
+    void EssBC_KG( const ALocal_NBC * const &nbc_part, const int &field );
+
+    void EssBC_G( const ALocal_NBC * const &nbc_part, const int &field );
+
+    // Natural boundary condition
+    void NatBC_G( const double &curr_time, const double &dt,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &element_s,
+        const IQuadPts * const &quad_s,
+        const ALocal_NBC * const &nbc_part,
+        const ALocal_EBC * const &ebc_part );
+
+    // Backflow integral on outlet surfaces
+    void BackFlow_G(
+        const PDNSolution * const &sol_0,
+        IPLocAssem * const &lassem_ptr,
+        FEAElement * const &element_s,
+        const IQuadPts * const &quad_s,
+        const ALocal_NBC * const &nbc_part,
+        const ALocal_EBC * const &ebc_part );
+
+    // Resistance type boundary condition on outlet surfaces
+    // void NatBC_Resis_G( const double &curr_time, const double &dt,
+    //     const PDNSolution * const &dot_sol,
+    //     const PDNSolution * const &sol,
+    //     IPLocAssem * const &lassem_ptr,
+    //     FEAElement * const &element_s,
+    //     const IQuadPts * const &quad_s,
+    //     const ALocal_NBC * const &nbc_part,
+    //     const ALocal_EBC * const &ebc_part,
+    //     const IGenBC * const &gbc );
+
+    // void NatBC_Resis_KG( const double &curr_time, const double &dt,
+    //     const PDNSolution * const &dot_sol,
+    //     const PDNSolution * const &sol,
+    //     IPLocAssem * const &lassem_ptr,
+    //     FEAElement * const &element_s,
+    //     const IQuadPts * const &quad_s,
+    //     const ALocal_NBC * const &nbc_part,
+    //     const ALocal_EBC * const &ebc_part,
+    //     const IGenBC * const &gbc );
+
+    // Weak imposition of no-slip boundary condition on wall
+    // void Weak_EssBC_KG( const double &curr_time, const double &dt,
+    //     const PDNSolution * const &sol,
+    //     IPLocAssem * const &lassem_ptr,
+    //     FEAElement * const &element_vs,
+    //     const IQuadPts * const &quad_s,
+    //     const ALocal_IEN * const &lien_ptr,
+    //     const FEANode * const &fnode_ptr,
+    //     const ALocal_NBC * const &nbc_part,
+    //     const ALocal_WeakBC * const &wbc_part);
+
+    // void Weak_EssBC_G( const double &curr_time, const double &dt,
+    //     const PDNSolution * const &sol,
+    //     IPLocAssem * const &lassem_ptr,
+    //     FEAElement * const &element_vs,
+    //     const IQuadPts * const &quad_s,
+    //     const ALocal_IEN * const &lien_ptr,
+    //     const FEANode * const &fnode_ptr,
+    //     const ALocal_NBC * const &nbc_part,
+    //     const ALocal_WeakBC * const &wbc_part);
+
+    void GetLocal(const double * const &array, const int * const &IEN,
+        double * const &local_array) const
+    {
+      for(int ii=0; ii<nLocBas; ++ii)
+      {
+        const int offset1 = ii * dof_sol;
+        const int offset2 = IEN[ii] * dof_sol;
+        for(int jj=0; jj<dof_sol; ++jj)
+          local_array[offset1 + jj] = array[offset2 + jj];
+      }
+    }
+
+    void GetLocal( const double * const &array, const int * const &IEN,
+        const int &in_locbas, double * const &local_array) const
+    {
+      for(int ii=0; ii<in_locbas; ++ii)
+      {
+        const int offset1 = ii * dof_sol;
+        const int offset2 = IEN[ii] * dof_sol;
+        for(int jj=0; jj<dof_sol; ++jj)
+          local_array[offset1 + jj] = array[offset2 + jj];
+      }
+    }
+};
+
+#endif