Snapshot port saving + DG LaxFriedrichs flux integrator completion (#39)

* SampleGenerator: save snapshot ports in sample generation.

* stokes_dg_mms: added direct solve

* ns_dg_mms: starting from stokes flow.

* with only domain integrator.

* checking LF scheme for steady ns.

* verifying TemamIntegrators. DGBdrTemamLFIntegrator needs debugging.

* temam integrators verified only for continuous galerkin. DGTemamFluxIntegrator needs debugging.

* us ns sketch initial loading

* initialize time integration

* navier solver step

* time integration loop and visualization

* sanity check at every time step.

* working version of DGLaxFriedrichsFluxIntegrator::AssembleFaceVector. other routines must be changed accordingly.

* steady ns with lax friedrichs flux verified.

* LaxFridriechsFlux ComputeFluxDotN.

* DGLaxFriedrichsFluxIntegrator::AssembleQuadVectorBase

* DGLaxFriedrichsFluxIntegrator::AssembleQuadGradBase

* DGLaxFriedrichsFluxIntegrator::AssembleQuadratureVector/Grad

* ROMNonlinearForm_fast.DGLaxFriedrichsFluxIntegrator

* DGLaxFriedrichsFlux::AssembleInterfaceVector/Grad

* DGLaxFriedrichsFlux::AssembleQuadratureVector/Grad for interface

* linalg_utils::AddwRtAP

* use AddwRtAP

* dg lf interface assemble vector/grad fast. this needs to be revisited.

include/dg_linear.hpp

@@ -80,10 +80,13 @@ class DGBdrLaxFriedrichsLFIntegrator : public LinearFormIntegrator
    Vector shape;
    VectorCoefficient &Q;
    DenseMatrix ELV;
+   Coefficient *Z = NULL;
+
+   double w;
 public:
    /// Constructs a boundary integrator with a given Coefficient QG
-   DGBdrLaxFriedrichsLFIntegrator(VectorCoefficient &QG)
-      : Q(QG) { }
+   DGBdrLaxFriedrichsLFIntegrator(VectorCoefficient &QG, Coefficient *ZG = NULL)
+      : Q(QG), Z(ZG) { }
 
    virtual bool SupportsDevice() { return false; }
 

include/interfaceinteg.hpp

@@ -70,6 +70,20 @@ class InterfaceNonlinearFormIntegrator : virtual public HyperReductionIntegrator
                                           const double &iw,
                                           const Vector &eltest1, const Vector &eltest2,
                                           Array2D<DenseMatrix*> &quadmats);
+
+   virtual void AddAssembleVector_Fast(const int s, const double qw,
+                                       FaceElementTransformations &Tr1,
+                                       FaceElementTransformations &Tr2,
+                                       const IntegrationPoint &ip,
+                                       const Vector &x1, const Vector &x2,
+                                       Vector &y1, Vector &y2);
+
+   virtual void AddAssembleGrad_Fast(const int s, const double qw,
+                                     FaceElementTransformations &Tr1,
+                                     FaceElementTransformations &Tr2,
+                                     const IntegrationPoint &ip,
+                                     const Vector &x1, const Vector &x2,
+                                     Array2D<SparseMatrix *> &jac);
 };
 
 class InterfaceDGDiffusionIntegrator : public InterfaceNonlinearFormIntegrator
@@ -260,25 +274,37 @@ class InterfaceDGElasticityIntegrator : public InterfaceNonlinearFormIntegrator
 
 /*
    Interior face, interface integrator
-   < [v], {uu \dot n} + \Lambda * [u] >
-   \Lambda = max( | u- \dot n |, | u+ \dot n | )
+   < [v], {uu \dot n} + \Lambda / 2 * [u] >
+   \Lambda = max( 2 * | u- \dot n |, 2 * | u+ \dot n | )
+
+   For boundary face,
+   (i) Neumann condition (UD == NULL)
+   < v-, (u-)(u-) \dot n >
+   (ii) Dirichlet condition (UD != NULL)
+   Same formulation, u+ = UD
 */
 class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator
 {
 private:
-   int dim, ndofs1, ndofs2, nvdofs;
-   double w, un1, un2, un;
+   int dim;
+   double un1, un2, un;
+   double u1mag, u2mag, normag;
+   Vector flux, u1, u2, nor;
+
+   int ndofs1, ndofs2, nvdofs;
+   double w;
    Coefficient *Q{};
+   VectorCoefficient *UD = NULL;
 
-   Vector nor, flux, shape1, shape2, u1, u2, tmp_vec;
+   Vector shape1, shape2;
    DenseMatrix udof1, udof2, elv1, elv2;
-   DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22, tmp;
+   DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22;
 
    // precomputed basis value at the sample point.
    Array<DenseMatrix *> shapes1, shapes2;
 public:
-   DGLaxFriedrichsFluxIntegrator(Coefficient &q, const IntegrationRule *ir = NULL)
-      : InterfaceNonlinearFormIntegrator(true, ir), Q(&q) {}
+   DGLaxFriedrichsFluxIntegrator(Coefficient &q, VectorCoefficient *ud = NULL, const IntegrationRule *ir = NULL)
+      : InterfaceNonlinearFormIntegrator(true, ir), Q(&q), UD(ud) {}
 
    void AssembleFaceVector(const FiniteElement &el1,
                            const FiniteElement &el2,
@@ -354,7 +380,43 @@ class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator
                                  const Vector &eltest1, const Vector &eltest2,
                                  Array2D<DenseMatrix*> &quadmats) override;
 
+   void AddAssembleVector_Fast(const int s, const double qw,
+                              FaceElementTransformations &Tr1,
+                              FaceElementTransformations &Tr2,
+                              const IntegrationPoint &ip,
+                              const Vector &x1, const Vector &x2,
+                              Vector &y1, Vector &y2) override;
+
+   void AddAssembleGrad_Fast(const int s, const double qw,
+                           FaceElementTransformations &Tr1,
+                           FaceElementTransformations &Tr2,
+                           const IntegrationPoint &ip,
+                           const Vector &x1, const Vector &x2,
+                           Array2D<SparseMatrix *> &jac) override;
+
 private:
+   void ComputeFluxDotN(const Vector &u1, const Vector &u2, const Vector &nor,
+                        const bool &eval2, Vector &flux);
+
+   void ComputeGradFluxDotN(const Vector &u1, const Vector &u2, const Vector &nor,
+                            const bool &eval2, const bool &ndofs2,
+                            DenseMatrix &gradu1, DenseMatrix &gradu2);
+
+   void AssembleQuadVectorBase(const FiniteElement &el1,
+                               const FiniteElement &el2,
+                               FaceElementTransformations *Tr1,
+                               FaceElementTransformations *Tr2,
+                               const IntegrationPoint &ip, const double &iw, const int &ndofs2,
+                               const DenseMatrix &elfun1, const DenseMatrix &elfun2,
+                               DenseMatrix &elvect1, DenseMatrix &elvect2);
+
+   void AssembleQuadGradBase(const FiniteElement &el1, const FiniteElement &el2,
+                             FaceElementTransformations *Tr1, FaceElementTransformations *Tr2,
+                             const IntegrationPoint &ip, const double &iw, const int &ndofs2,
+                             const DenseMatrix &elfun1, const DenseMatrix &elfun2,
+                             double &w, DenseMatrix &gradu1, DenseMatrix &gradu2,
+                             DenseMatrix &elmat11, DenseMatrix &elmat12, DenseMatrix &elmat21, DenseMatrix &elmat22);
+
    void AppendPrecomputeFaceCoeffs(const FiniteElementSpace *fes, 
                                     FaceElementTransformations *T,
                                     DenseMatrix &basis,

include/linalg_utils.hpp

@@ -136,6 +136,17 @@ SparseMatrix* SparseRtAP(DenseMatrix& R,
                         const Operator& A,
                         DenseMatrix& P);
 
+void AddwRtAP(DenseMatrix& R,
+            const Operator& A,
+            DenseMatrix& P,
+            DenseMatrix& RAP,
+            const double w=1.0);
+void AddwRtAP(DenseMatrix& R,
+            const Operator& A,
+            DenseMatrix& P,
+            SparseMatrix& RAP,
+            const double w=1.0);
+
 template<typename T>
 void PrintMatrix(const T &mat,
                 const std::string &filename);

include/multiblock_solver.hpp

@@ -128,6 +128,7 @@ friend class ParameterizedProblem;
    const bool IsNonlinear() const { return nonlinear_mode; }
    const bool UseRom() const { return use_rom; }
    ROMHandlerBase* GetROMHandler() const { return rom_handler; }
+   TopologyHandler* GetTopologyHandler() const { return topol_handler; }
    const TrainMode GetTrainMode() { return train_mode; }
    const bool IsVisualizationSaved() const { return save_visual; }
    const std::string GetSolutionFilePrefix() const { return sol_prefix; }

include/sample_generator.hpp

@@ -9,6 +9,7 @@
 #include "parameter.hpp"
 #include "linalg/BasisGenerator.h"
 #include "linalg_utils.hpp"
+#include "rom_handler.hpp"
 
 using namespace mfem;
 
@@ -19,6 +20,16 @@ enum SampleGeneratorType
    NUM_SAMPLE_GEN_TYPE
 };
 
+struct PortTag
+{
+   std::string Mesh1;
+   std::string Mesh2;
+   int Attr1;
+   int Attr2;
+};
+
+bool operator<(const PortTag &tag1, const PortTag &tag2);
+
 class SampleGenerator
 {
 protected:
@@ -53,6 +64,11 @@ class SampleGenerator
    std::vector<std::string> basis_tags;
    std::map<std::string, int> basis_tag2idx;
 
+   /* snapshot pairs per interface port, for nonlinear EQP */
+   std::vector<PortTag> port_tags;
+   std::map<PortTag, int> port_tag2idx;
+   Array<Array<int> *> port_colidxs;
+
 public:
    SampleGenerator(MPI_Comm comm);
 
@@ -89,9 +105,16 @@ class SampleGenerator
 
    const std::string GetSamplePath(const int& idx, const std::string &prefix = "");
 
-   void SaveSnapshot(BlockVector *U_snapshots, std::vector<std::string> &snapshot_basis_tags);
+   /*
+      Save each block of U_snapshots according to snapshot_basis_tags.
+      Number of blocks in U_snapshots must be equal to the size of snapshot_basis_tags.
+      The appended column indices of each basis tag are stored in col_idxs.
+   */
+   void SaveSnapshot(BlockVector *U_snapshots, std::vector<std::string> &snapshot_basis_tags, Array<int> &col_idxs);
+   void SaveSnapshotPorts(TopologyHandler *topol_handler, const TrainMode &train_mode, const Array<int> &col_idxs);
    void AddSnapshotGenerator(const int &fom_vdofs, const std::string &prefix, const std::string &basis_tag);
    void WriteSnapshots();
+   void WriteSnapshotPorts();
    const CAROM::Matrix* LookUpSnapshot(const std::string &basis_tag);
 
    void ReportStatus(const int &sample_idx);

sketches/CMakeLists.txt

@@ -27,6 +27,7 @@ add_executable(precond precond.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_mms ns_mms.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_dg_mms ns_dg_mms.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 add_executable(ns_rom ns_rom.cpp $<TARGET_OBJECTS:scaleupROMObj>)
+add_executable(usns usns.cpp $<TARGET_OBJECTS:scaleupROMObj>)
 
 file(COPY inputs/gen_interface.yml DESTINATION ${CMAKE_BINARY_DIR}/sketches/inputs)
 file(COPY meshes/2x2.mesh DESTINATION ${CMAKE_BINARY_DIR}/sketches/meshes)