libROM / MFEM v4.6 updates + Advection-diffusion

CMakeLists.txt

@@ -209,6 +209,9 @@ set(scaleupROMObj_SOURCES
   include/steady_ns_solver.hpp
   src/steady_ns_solver.cpp
 
+  include/advdiff_solver.hpp
+  src/advdiff_solver.cpp
+
   include/input_parser.hpp
   src/input_parser.cpp
 

docker/Dockerfile

@@ -29,6 +29,8 @@ WORKDIR /usr/lib
 RUN arch=$(uname -m) && sudo ln -s /usr/lib/${arch}-linux-gnu/libscalapack-openmpi.so.2.1.0 ./
 
 # re-install mfem with mumps, with cmake
+WORKDIR $LIB_DIR/mfem
+RUN sudo git checkout v4.6
 WORKDIR $LIB_DIR/mfem/build
 RUN sudo -E cmake .. -DBUILD_SHARED_LIBS=YES -DMFEM_USE_MPI=YES -DMFEM_USE_GSLIB=${MG} -DMFEM_USE_METIS=YES -DMFEM_USE_METIS_5=YES -DMFEM_USE_MUMPS=YES -DGSLIB_DIR="$LIB_DIR/gslib/build" -DMUMPS_DIR="$LIB_DIR/mumps/build/local"
 RUN sudo -E make -j 16
@@ -54,7 +56,7 @@ WORKDIR ./libROM
 RUN sudo git pull && sudo git checkout mpi-io
 WORKDIR ./build
 # libROM is using the MFEM without MUMPS right now.
-RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DUSE_MFEM=${USE_MFEM} -DMFEM_USE_GSLIB=${MFEM_USE_GSLIB}
+RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=${BUILD_TYPE} -DUSE_MFEM=OFF -DMFEM_USE_GSLIB=${MFEM_USE_GSLIB}
 RUN sudo make -j 16
 
 ENV LIBROM_DIR=$LIB_DIR/libROM

include/advdiff_solver.hpp

@@ -0,0 +1,66 @@
+// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: MIT
+
+#ifndef SCALEUPROM_ADVDIFF_SOLVER_HPP
+#define SCALEUPROM_ADVDIFF_SOLVER_HPP
+
+#include "poisson_solver.hpp"
+#include "stokes_solver.hpp"
+
+// By convention we only use mfem namespace as default, not CAROM.
+using namespace mfem;
+
+class AdvDiffSolver : public PoissonSolver
+{
+
+friend class ParameterizedProblem;
+
+protected:
+   double Pe = 0.0;  // Peclet number
+
+   // coefficient for prescribed velocity field.
+   // can be analytic function or solution from Stokes/SteadyNS equation.
+   Array<VectorCoefficient *> flow_coeffs;
+
+   /*
+      flow solver to obtain the prescribed velocity field. both StokesSolver / SteadyNSSolver can be used.
+   */
+   StokesSolver *stokes_solver = NULL;
+   bool load_flow = false;
+   bool save_flow = false;
+   std::string flow_file = "";
+
+   /* grid functions for visualizaing flow field */
+   /* NOTE(kevin): this will be set up at SaveVisualization. */
+   Array<FiniteElementSpace *> flow_fes;
+   Array<GridFunction *> flow_visual;
+   FiniteElementSpace *global_flow_fes = NULL;
+   Array<GridFunction *> global_flow_visual;
+
+public:
+   AdvDiffSolver();
+
+   virtual ~AdvDiffSolver();
+
+   void BuildDomainOperators() override;
+
+   // Component-wise assembly
+   void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp) override;
+
+   bool Solve() override;
+
+   void SetFlowAtSubdomain(std::function<void(const Vector &, Vector &)> F, const int m=-1);
+
+   void SetParameterizedProblem(ParameterizedProblem *problem) override;
+
+   void SaveVisualization() override;
+
+protected:
+   void SetMUMPSSolver() override;
+
+private:
+   void GetFlowField(ParameterizedProblem *flow_problem);
+};
+
+#endif

include/multiblock_solver.hpp

@@ -137,6 +137,8 @@ friend class ParameterizedProblem;
    BlockVector* GetSolution() { return U; }
    BlockVector* GetSolutionCopy() { return new BlockVector(*U); }
 
+   void SetSolutionSaveMode(const bool save_sol_);
+
    void GetVariableVector(const int &var_idx, BlockVector &global, BlockVector &var);
    void SetVariableVector(const int &var_idx, BlockVector &var, BlockVector &global);
 

include/parameterized_problem.hpp

@@ -92,9 +92,24 @@ extern double rforce_f;
 void tip_force(const Vector &x, Vector &u);
 }
 
+namespace advdiff_problem
+{
+
+extern bool analytic_flow;
+
+namespace advdiff_flow_past_array
+{
+   extern double q0, dq, qoffset;
+   extern Vector qk;
+
+   double qbdr(const Vector &x);
+}  //  namespace advdiff_flow_past_array
+
+}  // namespace advdiff_problem
+
 }
 
-enum BoundaryType
+enum class BoundaryType
 { 
    ZERO,
    DIRICHLET,
@@ -169,8 +184,8 @@ class PoissonComponent : public PoissonProblem
 public:
    PoissonComponent();
    virtual ~PoissonComponent() {};
-   virtual void SetParams(const std::string &key, const double &value) override;
-   virtual void SetParams(const Array<int> &indexes, const Vector &values) override;
+   void SetParams(const std::string &key, const double &value) override;
+   void SetParams(const Array<int> &indexes, const Vector &values) override;
 
 private:
    void SetBattr();
@@ -206,15 +221,17 @@ class StokesComponent : public StokesProblem
 
 class StokesFlowPastArray : public StokesComponent
 {
+friend class AdvDiffFlowPastArray;
+
 public:
    StokesFlowPastArray();
 
-   virtual void SetParams(const std::string &key, const double &value);
-   virtual void SetParams(const Array<int> &indexes, const Vector &values);
+   virtual void SetParams(const std::string &key, const double &value) override;
+   virtual void SetParams(const Array<int> &indexes, const Vector &values) override;
 
-private:
+protected:
    Vector *u0;
-   void SetBattr();
+   virtual void SetBattr();
 };
 
 class LinElastProblem : public ParameterizedProblem
@@ -251,6 +268,43 @@ class LinElastForceCantilever : public LinElastProblem
    LinElastForceCantilever();
 };
 
+namespace function_factory
+{
+
+namespace advdiff_problem
+{
+/*
+   flow_problem will be passed down to StokesSolver/SteadyNSSolver for obtaining velocity field.
+   It must be set appropriately within each AdvDiffSolver problems.
+*/
+extern StokesProblem *flow_problem;
+
+}  // namespace advdiff_problem
+
+}  // namespace function_factory
+
+class AdvDiffFlowPastArray : public StokesFlowPastArray
+{
+protected:
+   /*
+      flow_problem shares the same pointers with this class.
+      Thus every parameter set by this class is reflected to StokesFlowPastArrayProblem as well.
+      flow_problem will be passed down to StokesSolver/SteadyNSSolver for obtaining velocity field.
+   */
+   StokesFlowPastArray *flow_problem = NULL;
+
+public:
+   AdvDiffFlowPastArray();
+   virtual ~AdvDiffFlowPastArray();
+
+protected:
+   void SetBattr() override
+   {
+      StokesFlowPastArray::SetBattr();
+      flow_problem->SetBattr();
+   }
+};
+
 ParameterizedProblem* InitParameterizedProblem();
 
 #endif

include/poisson_solver.hpp

@@ -99,6 +99,9 @@ friend class ParameterizedProblem;
    void SanityCheckOnCoeffs();
 
    virtual void SetParameterizedProblem(ParameterizedProblem *problem);
+
+protected:
+   virtual void SetMUMPSSolver();
 };
 
 #endif

include/steady_ns_solver.hpp

@@ -150,35 +150,34 @@ friend class SteadyNSOperator;
 
    using StokesSolver::GetVariableNames;
 
-   virtual void InitVariables();
+   void InitVariables() override;
 
-   virtual void BuildOperators() override;
-   virtual void BuildDomainOperators();
+   void BuildDomainOperators() override;
 
-   virtual void SetupRHSBCOperators() override;
-   virtual void SetupDomainBCOperators() override;
+   void SetupRHSBCOperators() override;
+   void SetupDomainBCOperators() override;
 
-   virtual void Assemble();
+   void AssembleOperator() override;
 
-   virtual void SaveROMOperator(const std::string input_prefix="");
-   virtual void LoadROMOperatorFromFile(const std::string input_prefix="");
+   void SaveROMOperator(const std::string input_prefix="") override;
+   void LoadROMOperatorFromFile(const std::string input_prefix="") override;
 
    // Component-wise assembly
-   virtual void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp);
+   void BuildCompROMElement(Array<FiniteElementSpace *> &fes_comp) override;
    // virtual void BuildBdrROMElement(Array<FiniteElementSpace *> &fes_comp);
    // virtual void BuildInterfaceROMElement(Array<FiniteElementSpace *> &fes_comp);
-   virtual void SaveCompBdrROMElement(hid_t &file_id) override;
-   virtual void LoadCompBdrROMElement(hid_t &file_id) override;
+   void SaveCompBdrROMElement(hid_t &file_id) override;
+   void LoadCompBdrROMElement(hid_t &file_id) override;
 
-   virtual bool Solve();
+   bool Solve() override;
 
-   virtual void ProjectOperatorOnReducedBasis();
+   void ProjectOperatorOnReducedBasis() override;
 
-   virtual void SolveROM() override;
+   void SolveROM() override;
 
-   virtual void TrainEQP(SampleGenerator *sample_generator) override;
-   virtual void SaveEQP() override;
-   virtual void LoadEQP() override;
+   void TrainEQP(SampleGenerator *sample_generator) override;
+   void SaveEQP() override;
+   void LoadEQP() override;
 
 private:
    DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace);

include/stokes_solver.hpp

@@ -179,7 +179,7 @@ friend class ParameterizedProblem;
    // For bilinear case.
    // system-specific.
    virtual void AssembleInterfaceMatrices();
-   virtual void SetupMUMPSSolver();
+   virtual void SetupMUMPSSolver(bool set_oper);
    virtual void SetupPressureMassMatrix();
 
    // Component-wise assembly

sketches/ns_dg_mms.cpp

@@ -413,7 +413,7 @@ int main(int argc, char *argv[])
    MUMPSSolver *J_mumps = NULL;
    if (direct_solve)
    {
-      J_mumps = new MUMPSSolver;
+      J_mumps = new MUMPSSolver(MPI_COMM_WORLD);
       J_mumps->SetPrintLevel(-1);
       J_mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
       J_solver = J_mumps;

sketches/ns_rom.cpp

@@ -242,7 +242,7 @@ class SteadyNavierStokes : public Operator
 
       if (direct_solve)
       {
-         J_mumps = new MUMPSSolver();
+         J_mumps = new MUMPSSolver(MPI_COMM_WORLD);
          J_mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
          J_solver = J_mumps;
       }
@@ -1768,7 +1768,7 @@ int main(int argc, char *argv[])
       J_gmres.SetMaxIter(maxIter);
       J_gmres.SetPrintLevel(-1);
 
-      MUMPSSolver J_mumps;
+      MUMPSSolver J_mumps(MPI_COMM_WORLD);
       J_mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
 
       Solver *J_solver;

sketches/stokes_mms4.cpp

@@ -322,7 +322,7 @@ int main(int argc, char *argv[])
       HYPRE_BigInt sys_row_starts[2] = {0, sys_mat->NumRows()};
       HypreParMatrix Aop(MPI_COMM_WORLD, sys_glob_size, sys_row_starts, sys_mat);
 
-      MUMPSSolver mumps;
+      MUMPSSolver mumps(MPI_COMM_WORLD);
       mumps.SetPrintLevel(1);
       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
       // mumps.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_POSITIVE_DEFINITE);