Merge branch 'master' into hpcftools/flowSolver

- added comments to Blf and Lf contributions
- split setup and step into vel, auxiliary and pressure part

miniapps/navier/CMakeLists.txt

@@ -11,9 +11,14 @@
 
 if (MFEM_USE_MPI)
   list(APPEND NAVIER_COMMON_SOURCES
-      navier_solver.cpp)
+      navier_solver.cpp
+      incompressible_navier_solver.cpp
+      stokes_solver.cpp)
+
   list(APPEND NAVIER_COMMON_HEADERS
-      navier_solver.hpp)
+      navier_solver.hpp
+      incompressible_navier_solver.hpp
+      stokes_solver.hpp)
 
   convert_filenames_to_full_paths(NAVIER_COMMON_SOURCES)
   convert_filenames_to_full_paths(NAVIER_COMMON_HEADERS)
@@ -52,6 +57,11 @@ if (MFEM_USE_MPI)
     ${NAVIER_COMMON_FILES}
     LIBRARIES mfem)
 
+  add_mfem_miniapp(incompNS_2Dtest
+    MAIN incompNS_2Dtest.cpp
+    ${NAVIER_COMMON_FILES}
+    LIBRARIES mfem)
+
   add_mfem_miniapp(navier_turbchan
     MAIN navier_turbchan.cpp
     ${NAVIER_COMMON_FILES}

miniapps/navier/incompNS_2Dtest.cpp

@@ -0,0 +1,134 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// 3D flow over a cylinder benchmark example
+
+#include "incompressible_navier_solver.hpp"
+#include <fstream>
+
+using namespace mfem;
+using namespace incompressible_navier;
+
+void vel(const Vector &x, real_t t, Vector &u)
+{
+   real_t xi = x(0);
+   real_t yi = x(1);
+
+   u = 0.0;
+}
+
+void vel_inlet(const Vector &x, real_t t, Vector &u)
+{
+    u = 0.0;
+    if (x(0) < 0.001) {
+ 
+      u(0) = -0.001 * (std::pow(x(1) - 0.5, 2.0) - 0.25);
+    }
+}
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   int serial_refinements = 1;
+   int vOrder = 2;
+   int pOrder = 1;
+   int tOrder = 1;
+   real_t kin_vis = 20.0;
+   real_t dt = 1e-2;
+   real_t t = 0.0;
+   real_t t_final = 1.0;
+   bool last_step = false;
+
+   //Mesh *mesh = new Mesh("box-cylinder.mesh");
+   Mesh mesh = Mesh::MakeCartesian2D(90, 30, mfem::Element::QUADRILATERAL, true, 3, 1);
+
+   for (int i = 0; i < serial_refinements; ++i)
+   {
+      mesh.UniformRefinement();
+   }
+
+   if (Mpi::Root())
+   {
+      std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+   }
+
+   auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
+
+   // Create the flow solver.
+   IncompressibleNavierSolver flowsolver(pmesh, vOrder, pOrder, tOrder, kin_vis);
+   flowsolver.EnablePA(false);
+
+   // // Set the initial condition.
+   // ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
+   // VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
+   // u_ic->ProjectCoefficient(u_excoeff);
+
+   // Add Dirichlet boundary conditions to velocity space restricted to
+   // selected attributes on the mesh.
+   Array<int> attr(pmesh->bdr_attributes.Max());           attr = 0;
+   Array<int> attr_inlet(pmesh->bdr_attributes.Max());     attr_inlet = 0;
+   // Inlet is attribute 1.
+   attr[0] = 1;
+   // Walls is attribute 3.
+   attr[2] = 1;
+   flowsolver.AddVelDirichletBC(vel, attr);
+
+   attr_inlet[3] = 1;
+   flowsolver.AddVelDirichletBC(vel_inlet, attr_inlet);
+
+   flowsolver.Setup(dt);
+
+   ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
+   ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
+   ParGridFunction *psi_gf = flowsolver.GetCurrentPsi();
+
+   ParaViewDataCollection pvdc("3dfoc", pmesh);
+   pvdc.SetDataFormat(VTKFormat::BINARY32);
+   //pvdc.SetHighOrderOutput(true);
+   pvdc.SetCycle(0);
+   pvdc.SetTime(t);
+   pvdc.RegisterField("velocity", u_gf);
+   pvdc.RegisterField("pressure", p_gf);
+   pvdc.RegisterField("psi", psi_gf);
+   pvdc.Save();
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+      }
+
+      flowsolver.Step(t, dt, step);
+
+      if (step % 1 == 0)
+      {
+         pvdc.SetCycle(step);
+         pvdc.SetTime(t);
+         pvdc.Save();
+      }
+
+      if (Mpi::Root())
+      {
+         printf("%11s %11s\n", "Time", "dt");
+         printf("%.5E %.5E\n", t, dt);
+         fflush(stdout);
+      }
+   }
+
+   // flowsolver.PrintTimingData();
+
+   delete pmesh;
+
+   return 0;
+}

miniapps/navier/incompressible_navier_solver.cpp

@@ -0,0 +1,515 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "incompressible_navier_solver.hpp"
+#include "../../general/forall.hpp"
+#include <fstream>
+#include <iomanip>
+
+using namespace mfem;
+using namespace incompressible_navier;
+
+IncompressibleNavierSolver::IncompressibleNavierSolver(ParMesh *mesh, int velorder, int porder, int torder_, real_t kin_vis)
+   : pmesh(mesh), velorder(velorder), porder(porder), torder(torder_), kin_vis(kin_vis),
+     gll_rules(0, Quadrature1D::GaussLobatto), velGF(torder_+1,nullptr), pGF(torder_+1,nullptr)
+{
+   vfec   = new H1_FECollection(velorder, pmesh->Dimension());
+   psifec = new H1_FECollection(porder);
+   pfec   = new H1_FECollection(porder);
+   vfes   = new ParFiniteElementSpace(pmesh, vfec, pmesh->Dimension());
+   psifes = new ParFiniteElementSpace(pmesh, pfec);
+   pfes   = new ParFiniteElementSpace(pmesh, pfec);
+
+   // Check if fully periodic mesh
+   if (!(pmesh->bdr_attributes.Size() == 0))
+   {
+      vel_ess_attr.SetSize(pmesh->bdr_attributes.Max());
+      vel_ess_attr = 0;
+
+      pres_ess_attr.SetSize(pmesh->bdr_attributes.Max());
+      pres_ess_attr = 0;
+   }
+
+   int vfes_truevsize = vfes->GetTrueVSize();
+   int pfes_truevsize = pfes->GetTrueVSize();
+
+   for( int i = 0; i<torder+1; i++)
+   {
+      velGF[i] = new ParGridFunction(vfes); *velGF[i] = 0.0;
+      pGF[i]   = new ParGridFunction(pfes); *pGF[i]   = 0.0;
+   }
+
+   psiGF.SetSpace(psifes);
+   DvGF.SetSpace(vfes);
+   divVelGF.SetSpace(pfes);
+   pRHS.SetSpace(pfes);
+}
+
+void IncompressibleNavierSolver::Setup(real_t dt)
+{
+   if (verbose && pmesh->GetMyRank() == 0)
+   {
+      mfem::out << "Setup" << std::endl;
+      if (partial_assembly)
+      {
+         mfem::out << "Using Partial Assembly" << std::endl;
+      }
+      else
+      {
+         mfem::out << "Using Full Assembly" << std::endl;
+      }
+   }
+
+   this->Setup_velocity( dt );
+
+   this->Setup_auxiliary( dt );
+
+   this->Setup_pressure( dt );
+}
+
+void IncompressibleNavierSolver::Setup_velocity(real_t dt)
+{
+   // GLL integration rule (Numerical Integration)
+   const IntegrationRule &ir_ni = gll_rules.Get(vfes->GetFE(0)->GetGeomType(),
+                                                2 * velorder - 1);
+
+   vfes->GetEssentialTrueDofs(vel_ess_attr, vel_ess_tdof);
+
+   //-------------------------------------------------------------------------
+
+   //Setup of coefficient for mass term of Eq(13)
+   dtCoeff      = new ConstantCoefficient(1.0/dt);   
+   auto *vmass_blfi = new VectorMassIntegrator(*dtCoeff);       
+
+   //Setup of coefficient for stiffness term of Eq(13)
+   kinvisCoeff  = new ConstantCoefficient(kin_vis);
+   auto *vdiff_blfi = new VectorDiffusionIntegrator(*kinvisCoeff);   
+
+   // setup of Bilinear form of Eq(13)
+   velBForm = new ParBilinearForm(vfes);
+   if (numerical_integ)
+   {
+       vmass_blfi->SetIntRule(&ir_ni); 
+       vdiff_blfi->SetIntRule(&ir_ni);
+   }
+   velBForm->AddDomainIntegrator(vmass_blfi);
+   velBForm->AddDomainIntegrator(vdiff_blfi);
+   if (partial_assembly)
+   {
+      velBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+
+   velBForm->Assemble();
+   velBForm->FormSystemMatrix(vel_ess_tdof, vOp);
+
+   //-------------------------------------------------------------------------
+   
+   //Setup of coefficient for Eq(18)
+   pUnitVectorCoeff = new UnitVectorGridFunctionCoeff(pmesh->Dimension());
+   auto *pvel_lfi = new VectorDomainLFGradIntegrator(*pUnitVectorCoeff);
+
+   //Setup of coefficient for Eq(20)
+   nonlinTermCoeff = new NonLinTermVectorGridFunctionCoeff(pmesh->Dimension());
+   auto *p_nonlintermlfi = new VectorDomainLFIntegrator(*nonlinTermCoeff);
+
+   //Setup of coefficient for Eq(21)
+   prevVelLoadCoeff = new PrevVelVectorGridFunctionCoeff(pmesh->Dimension());
+   auto *prevVelLoadLFi = new VectorDomainLFIntegrator(*prevVelLoadCoeff);
+
+   //Setup of linear form of Eq(13)
+   velLForm = new ParLinearForm(vfes);
+   if (numerical_integ)
+   {
+      prevVelLoadLFi->SetIntRule(&ir_ni);
+      pvel_lfi->SetIntRule(&ir_ni);
+      p_nonlintermlfi->SetIntRule(&ir_ni);
+   }
+   velLForm->AddDomainIntegrator(prevVelLoadLFi);
+   velLForm->AddDomainIntegrator(pvel_lfi);
+   velLForm->AddDomainIntegrator(p_nonlintermlfi);
+
+   //-------------------------------------------------------------------------
+
+   if (partial_assembly)
+   {
+      Vector diag_pa(vfes->GetTrueVSize());
+      velBForm->AssembleDiagonal(diag_pa);
+      velInvPC = new OperatorJacobiSmoother(diag_pa, vel_ess_tdof);
+   }
+   else
+   {
+      velInvPC = new HypreSmoother(*vOp.As<HypreParMatrix>());
+      dynamic_cast<HypreSmoother *>(velInvPC)->SetType(HypreSmoother::Jacobi, 1);
+   }
+
+   velInv = new CGSolver(vfes->GetComm());
+   velInv->iterative_mode = true;
+   velInv->SetOperator(*vOp);
+   velInv->SetPreconditioner(*velInvPC);
+   velInv->SetPrintLevel(pl_velsolve);
+   velInv->SetRelTol(rtol_velsolve);
+   velInv->SetMaxIter(1200);
+}
+
+void IncompressibleNavierSolver::Setup_auxiliary(real_t dt)
+{
+   // GLL integration rule (Numerical Integration)
+   const IntegrationRule &ir_ni = gll_rules.Get(vfes->GetFE(0)->GetGeomType(),
+                                                2 * velorder - 1);
+   Array<int> empty;
+
+   // setup of Bilinear form of Eq(14)
+   psiBForm = new ParBilinearForm(psifes);
+   auto *psidiff_blfi = new DiffusionIntegrator;
+
+   if (numerical_integ)
+   {
+       psidiff_blfi->SetIntRule(&ir_ni); 
+   }
+   psiBForm->AddDomainIntegrator(psidiff_blfi);
+   if (partial_assembly)
+   {
+      psiBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+
+   psiBForm->Assemble();
+   psiBForm->FormSystemMatrix(empty, psiOp);
+
+   //-------------------------------------------------------------------------
+
+   //Setup of coefficient for linear form in Eq(14)
+   DvelCoeff = new VectorGridFunctionCoefficient;
+   auto *Dvel_lfi = new DomainLFGradIntegrator(*DvelCoeff);
+      
+   //Setup of linear form of Eq(14)
+   psiLForm = new ParLinearForm(psifes);
+   
+   if (numerical_integ)
+   {
+      Dvel_lfi->SetIntRule(&ir_ni);
+   }
+   psiLForm->AddDomainIntegrator(Dvel_lfi);
+
+   //-------------------------------------------------------------------------
+   if (partial_assembly)
+   {
+      int psifes_truevsize = psifes->GetTrueVSize();
+      mfem::Vector psin(psifes_truevsize);      psin = 0.0;
+      mfem::Vector respsi(psifes_truevsize);    respsi = 0.0;
+
+      lor = new ParLORDiscretization(*psiBForm, empty);
+      psiInvPC = new HypreBoomerAMG(lor->GetAssembledMatrix());
+      psiInvPC->SetPrintLevel(0);
+      psiInvPC->Mult(respsi, psin);
+      SpInvOrthoPC = new OrthoSolver(psifes->GetComm());
+      SpInvOrthoPC->SetSolver(*psiInvPC);
+   }
+   else
+   {
+      psiInvPC = new HypreBoomerAMG(*psiOp.As<HypreParMatrix>());
+      psiInvPC->SetPrintLevel(0);
+      SpInvOrthoPC = new OrthoSolver(psifes->GetComm());
+      SpInvOrthoPC->SetSolver(*psiInvPC);
+   }
+
+   psiInv = new CGSolver(psifes->GetComm());
+   psiInv->iterative_mode = true;
+   psiInv->SetOperator(*psiOp);
+   psiInv->SetPreconditioner(*SpInvOrthoPC);
+   psiInv->SetPrintLevel(pl_psisolve);
+   psiInv->SetRelTol(rtol_psisolve);
+   psiInv->SetMaxIter(1000);
+}
+
+void IncompressibleNavierSolver::Setup_pressure(real_t dt)
+{
+   // GLL integration rule (Numerical Integration)
+   const IntegrationRule &ir_ni = gll_rules.Get(vfes->GetFE(0)->GetGeomType(),
+                                                2 * velorder - 1);
+   Array<int> empty;
+
+   //-------------------------------------------------------------------------
+
+   // setup of Bilinear form of Eq(15)
+   pBForm = new ParBilinearForm(pfes);
+   auto *pmass_blfi = new MassIntegrator;
+
+   if (numerical_integ)
+   {
+       pmass_blfi->SetIntRule(&ir_ni); 
+   }
+   pBForm->AddDomainIntegrator(pmass_blfi);
+   if (partial_assembly)
+   {
+      pBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+
+   pBForm->Assemble();
+   pBForm->FormSystemMatrix(empty, pOp);
+
+   //-------------------------------------------------------------------------
+
+   //Setup of divergence of velocity coefficient for linear form in Eq(15)
+   divVelCoeff = new DivergenceGridFunctionCoefficient(velGF[0]);
+
+   //Setup of coefficient for linear form in Eq(14)
+   pRHSCoeff = new GridFunctionCoefficient(&pRHS);
+   auto *p_lfi = new DomainLFIntegrator(*pRHSCoeff);
+
+   //Setup of linear form of Eq(15)
+   pLForm = new ParLinearForm(pfes); 
+   if (numerical_integ)
+   {
+      p_lfi->SetIntRule(&ir_ni);
+   }
+   pLForm->AddDomainIntegrator(p_lfi);
+
+   //-------------------------------------------------------------------------
+   if (partial_assembly)
+   {
+      Vector diag_pa(pfes->GetTrueVSize());
+      pBForm->AssembleDiagonal(diag_pa);
+      pInvPC = new OperatorJacobiSmoother(diag_pa, empty);
+   }
+   else
+   {
+      pInvPC = new HypreSmoother(*pOp.As<HypreParMatrix>());
+      dynamic_cast<HypreSmoother *>(pInvPC)->SetType(HypreSmoother::Jacobi, 1);
+   }
+
+   pInv = new CGSolver(pfes->GetComm());
+   pInv->iterative_mode = true;
+   pInv->SetOperator(*pOp);
+   pInv->SetPreconditioner(*pInvPC);
+   pInv->SetPrintLevel(pl_psolve);
+   pInv->SetRelTol(rtol_psolve);
+   pInv->SetMaxIter(1000);
+}
+
+void IncompressibleNavierSolver::UpdateTimestepHistory(real_t dt)
+{
+
+}
+
+void IncompressibleNavierSolver::Step(real_t &time, real_t dt, int current_step)
+{
+   this->Step_velocity(time, dt, current_step);
+
+   this->Step_auxiliary(time, dt, current_step);
+
+   this->Step_pressure(time, dt, current_step);
+
+   *velGF[1] = *velGF[0];
+   *pGF[1]   = *pGF[0];
+   
+   mfem::out << "It: " << iter << " | Iter_U: " << iter_vsolve << " | Iter_Psi: " << iter_psisolve << " | Iter_P: " << iter_psolve  << "\n";
+   mfem::out << "It: " << iter << " | Resid_U: " << res_vsolve << " | Resid_Psi: " << res_psisolve << " | Resid_P: " << res_psisolve  << "\n";
+
+   time += dt;
+   iter ++;
+}
+
+void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt, int current_step)
+{
+   for (auto &vel_dbc : vel_dbcs)
+   {
+      velGF[0]->ProjectBdrCoefficient(*vel_dbc.coeff, vel_dbc.attr);
+      velGF[1]->ProjectBdrCoefficient(*vel_dbc.coeff, vel_dbc.attr);
+   }
+
+   //Update state in coefficient for Eq(18)
+   pUnitVectorCoeff->SetGridFunction( pGF[1] );
+   
+   //Update state in coefficient for Eq(20)
+   nonlinTermCoeff->SetGridFunction( velGF[1] );
+
+   //Update state in coefficient for Eq(21)
+   prevVelLoadCoeff ->SetGridFunction( velGF[1], dt ); 
+
+   velLForm->Assemble();
+   velLForm->ParallelAssemble(velLF);
+
+   Vector X1, B1;
+
+   if (partial_assembly)
+   {
+      auto *vpC = vOp.As<ConstrainedOperator>();
+      EliminateRHS(*velBForm, *vpC, vel_ess_tdof, *velGF[0], velLF, X1, B1, 1);
+   }
+   else
+   {
+      velBForm->FormLinearSystem(vel_ess_tdof, *velGF[0], velLF, vOp  , X1, B1, 1);
+   }
+
+   velInv->Mult(B1, X1);
+   iter_vsolve = velInv->GetNumIterations();
+   res_vsolve = velInv->GetFinalNorm();
+   velBForm->RecoverFEMSolution(X1, velLF, *velGF[0]);
+}
+
+void IncompressibleNavierSolver::Step_auxiliary(real_t &time, real_t dt, int current_step)
+{
+   // Compute new increment GF for LF of Eq(14)  and update state in coefficient
+   subtract(1.0/dt, *velGF[0], *velGF[1], DvGF);
+   DvelCoeff->SetGridFunction( &DvGF );
+
+   psiLForm->Assemble();
+   psiLForm->ParallelAssemble(psiLF);
+
+   Vector X2, B2;   
+   Array<int> empty;
+   if (partial_assembly)
+   {
+      auto *psipC = psiOp.As<ConstrainedOperator>();
+      EliminateRHS(*psiBForm, *psipC, empty, psiGF, psiLF, X2, B2, 1);
+   }
+   else
+   {
+      psiBForm->FormLinearSystem(empty, psiGF, psiLF, psiOp, X2, B2, 1);
+   }
+
+   psiInv->Mult(B2, X2);
+   iter_psisolve = psiInv->GetNumIterations();
+   res_psisolve = psiInv->GetFinalNorm();
+   psiBForm->RecoverFEMSolution(X2, psiLF, psiGF);
+}
+
+void IncompressibleNavierSolver::Step_pressure(real_t &time, real_t dt, int current_step)
+{
+   Array<int> empty;
+
+   // Compute new GF for LF of Eq(15) and update state in coefficient
+   divVelCoeff->SetGridFunction( velGF[0]);
+   divVelGF.ProjectCoefficient( *divVelCoeff );
+
+   add( *pGF[1], psiGF, pRHS);
+   add( pRHS, -1.0*kin_vis, divVelGF, pRHS);
+   pRHSCoeff->SetGridFunction( &pRHS );
+
+   pLForm->Assemble();
+   pLForm->ParallelAssemble(pLF);
+
+   Vector X3, B3;
+
+   if (partial_assembly)
+   {
+      auto *ppC = pOp.As<ConstrainedOperator>();
+      EliminateRHS(*pBForm, *ppC, empty, *pGF[0], pLF, X3, B3, 1);
+   }
+   else
+   {
+      pBForm->FormLinearSystem(empty, *pGF[0] , pLF , pOp , X3, B3, 1);
+   }
+
+   pInv->Mult(B3, X3);
+   iter_psolve = pInv->GetNumIterations();
+   res_psisolve = pInv->GetFinalNorm();
+   pBForm->RecoverFEMSolution(X3, pLF, *pGF[0]);
+}
+
+
+
+void IncompressibleNavierSolver::EliminateRHS(Operator &A,
+                                ConstrainedOperator &constrainedA,
+                                const Array<int> &ess_tdof_list,
+                                Vector &x,
+                                Vector &b,
+                                Vector &X,
+                                Vector &B,
+                                int copy_interior)
+{
+   const Operator *Po = A.GetOutputProlongation();
+   const Operator *Pi = A.GetProlongation();
+   const Operator *Ri = A.GetRestriction();
+   A.InitTVectors(Po, Ri, Pi, x, b, X, B);
+   if (!copy_interior)
+   {
+      X.SetSubVectorComplement(ess_tdof_list, 0.0);
+   }
+   constrainedA.EliminateRHS(X, B);
+}
+
+real_t IncompressibleNavierSolver::ComputeCFL(ParGridFunction &u, real_t dt)
+{
+   
+   return 0;
+}
+
+void IncompressibleNavierSolver::AddVelDirichletBC(VectorCoefficient *coeff, Array<int> &attr)
+{
+   vel_dbcs.emplace_back(attr, coeff);
+
+   if (verbose && pmesh->GetMyRank() == 0)
+   {
+      mfem::out << "Adding Velocity Dirichlet BC to attributes ";
+      for (int i = 0; i < attr.Size(); ++i)
+      {
+         if (attr[i] == 1)
+         {
+            mfem::out << i << " ";
+         }
+      }
+      mfem::out << std::endl;
+   }
+
+   for (int i = 0; i < attr.Size(); ++i)
+   {
+      MFEM_ASSERT((vel_ess_attr[i] && attr[i]) == 0,
+                  "Duplicate boundary definition deteceted.");
+      if (attr[i] == 1)
+      {
+         vel_ess_attr[i] = 1;
+      }
+   }
+
+}
+
+void IncompressibleNavierSolver::AddVelDirichletBC(VecFuncT *f, Array<int> &attr)
+{
+   AddVelDirichletBC(new VectorFunctionCoefficient(pmesh->Dimension(), f), attr);
+}
+
+IncompressibleNavierSolver::~IncompressibleNavierSolver()
+{
+   delete velBForm;
+   delete psiBForm;
+   delete pBForm;
+
+   delete kinvisCoeff;
+   delete dtCoeff;
+
+   for( int i = 0; i<torder+1; i++)
+   {
+      delete velGF[i];
+      delete pGF[i];
+   }
+
+   delete DvelCoeff;
+   delete divVelCoeff;
+   delete pRHSCoeff;
+   delete pUnitVectorCoeff;
+
+   delete velInv;
+   delete velInvPC;
+   delete psiInv;
+   delete SpInvOrthoPC;
+   delete psiInvPC;
+   delete lor;
+   delete pInv;
+   delete pInvPC;
+
+   delete vfec;
+   delete psifec;
+   delete pfec;
+   delete vfes;
+   delete psifes;
+   delete pfes;
+}

miniapps/navier/incompressible_navier_solver.hpp

@@ -0,0 +1,362 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_INCOMP_NAVIER_SOLVER_HPP
+#define MFEM_INCOMP_NAVIER_SOLVER_HPP
+
+#define INCOMP_NAVIER_VERSION 0.1
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+namespace incompressible_navier
+{
+using VecFuncT = void(const Vector &x, real_t t, Vector &u);
+using ScalarFuncT = real_t(const Vector &x, real_t t);
+
+//Coefficient which computed contribution of Eq(18)
+class UnitVectorGridFunctionCoeff : public VectorCoefficient
+{
+public:
+   UnitVectorGridFunctionCoeff( int dim)
+   : VectorCoefficient(dim*dim)
+   { }
+
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      real_t coeffVal = gridfunc_->GetValue(T, ip);
+
+      V.SetSize(vdim); V = 0.0;                           // FIXME
+      V[0] = coeffVal;
+      V[3] = coeffVal;
+   }
+
+   void SetGridFunction( GridFunction * gridfunc )
+   {
+      gridfunc_ = gridfunc;
+   }
+
+   GridFunction *gridfunc_ = nullptr;
+};
+
+//Coefficient which computed contribution of Eq(21)
+class PrevVelVectorGridFunctionCoeff : public VectorCoefficient
+{
+public:
+   PrevVelVectorGridFunctionCoeff( int dim)
+   : VectorCoefficient(dim)
+   { }
+
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      V.SetSize(vdim);
+      gridFuncCoeff->Eval(V, T, ip);
+
+      V *= 1.0/dt_;
+   }
+
+   void SetGridFunction( GridFunction * gridfunc, real_t dt )
+   {
+      gridfunc_ = gridfunc;
+      dt_ = dt;
+      delete gridFuncCoeff;
+      gridFuncCoeff = new VectorGridFunctionCoefficient( gridfunc );
+   }
+
+   GridFunction *gridfunc_ = nullptr;
+   VectorGridFunctionCoefficient *gridFuncCoeff = nullptr;
+   real_t dt_;
+};
+
+//Coefficient which computed contribution of Eq(20)
+class NonLinTermVectorGridFunctionCoeff : public VectorCoefficient
+{
+public:
+   NonLinTermVectorGridFunctionCoeff( int dim)
+   : VectorCoefficient(dim)
+   { }
+
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      Vector val(vdim);
+      Vector resultVal(vdim);
+      DenseMatrix vecGrad;
+      V.SetSize(vdim);
+      gridFuncCoeff->Eval(val, T, ip);
+
+      gridfunc_->GetVectorGradient(T, vecGrad);
+
+      vecGrad.MultTranspose( val, V );
+
+      V *= -1.0;
+   }
+
+   void SetGridFunction( ParGridFunction * gridfunc )
+   {
+      delete gridFuncCoeff;
+      gridfunc_ = gridfunc;
+      gridFuncCoeff = new VectorGridFunctionCoefficient( gridfunc );
+   }
+
+   VectorGridFunctionCoefficient *gridFuncCoeff = nullptr;
+   ParGridFunction *gridfunc_ = nullptr;
+};
+
+/// Container for a Dirichlet boundary condition of the velocity field.
+class VelDirichletBC_T
+{
+public:
+   VelDirichletBC_T(Array<int> attr, VectorCoefficient *coeff)
+      : attr(attr), coeff(coeff)
+   {}
+
+   VelDirichletBC_T(VelDirichletBC_T &&obj)
+   {
+      // Deep copy the attribute array
+      this->attr = obj.attr;
+
+      // Move the coefficient pointer
+      this->coeff = obj.coeff;
+      obj.coeff = nullptr;
+   }
+
+   ~VelDirichletBC_T() { delete coeff; }
+
+   Array<int> attr;
+   VectorCoefficient *coeff;
+};
+
+/// Transient incompressible Navier Stokes solver in a split scheme formulation.
+/**
+ * This implementation of a transient incompressible Navier Stokes solver uses
+ * the non-dimensionalized formulation. The coupled momentum and
+ * incompressibility equations are decoupled using the split scheme described in
+ * [1]. This leads to three solving steps.
+ *
+ */
+class IncompressibleNavierSolver
+{
+public:
+   /// Initialize data structures, set FE space order and kinematic viscosity.
+   /**
+    * The ParMesh @a mesh can be a linear or curved parallel mesh. The @a order
+    * of the finite element spaces is
+    */
+   IncompressibleNavierSolver(ParMesh *mesh, int velorder, int porder, int tOrder, real_t kin_vis);
+
+   /// Initialize forms, solvers and preconditioners.
+   void Setup(real_t dt);
+
+   void Setup_velocity(real_t dt);
+
+   void Setup_auxiliary(real_t dt);
+
+   void Setup_pressure(real_t dt);
+
+   /// Compute solution at the next time step t+dt.
+   /**
+    * This method can 
+    */
+   void Step(real_t &time, real_t dt, int cur_step);
+
+   void Step_velocity(real_t &time, real_t dt, int cur_step);
+
+   void Step_auxiliary(real_t &time, real_t dt, int cur_step);
+
+   void Step_pressure(real_t &time, real_t dt, int cur_step);
+
+   /// Return a pointer to the provisional velocity ParGridFunction.
+   ParGridFunction *GetProvisionalVelocity() { return velGF[1]; }
+
+   /// Return a pointer to the current velocity ParGridFunction.
+   ParGridFunction *GetCurrentVelocity() { return velGF[0]; }
+
+   /// Return a pointer to the current pressure ParGridFunction.
+   ParGridFunction *GetCurrentPressure() { return pGF[0]; }
+
+   /// Return a pointer to the current pressure ParGridFunction.
+   ParGridFunction *GetCurrentPsi() { return &psiGF ; }
+
+
+   /// Add a Dirichlet boundary condition to the velocity field.
+   void AddVelDirichletBC(VectorCoefficient *coeff, Array<int> &attr);
+
+   void AddVelDirichletBC(VecFuncT *f, Array<int> &attr);
+
+   /// Add a Dirichlet boundary condition to the pressure field.
+   // void AddPresDirichletBC(Coefficient *coeff, Array<int> &attr);
+
+   // void AddPresDirichletBC(ScalarFuncT *f, Array<int> &attr);
+
+   /// Enable partial assembly for every operator.
+   void EnablePA(bool pa) { partial_assembly = pa; }
+
+   /// Enable numerical integration rules. This means collocated quadrature at
+   /// the nodal points.
+   void EnableNI(bool ni) { numerical_integ = ni; }
+
+   /// Print timing summary of the solving routine.
+   /**
+    * The summary shows the timing in seconds in the first row of
+    *
+    */
+   void PrintTimingData();
+
+   ~IncompressibleNavierSolver();
+
+   /// Rotate entries in the time step and solution history arrays.
+   void UpdateTimestepHistory(real_t dt);
+
+
+   /// Compute CFL
+   real_t ComputeCFL(ParGridFunction &u, real_t dt);
+
+protected:
+
+   /// Eliminate essential BCs in an Operator and apply to RHS.
+   void EliminateRHS(Operator &A,
+                     ConstrainedOperator &constrainedA,
+                     const Array<int> &ess_tdof_list,
+                     Vector &x,
+                     Vector &b,
+                     Vector &X,
+                     Vector &B,
+                     int copy_interior = 0);
+
+   /// Enable/disable debug output.
+   bool debug = false;
+
+   /// Enable/disable verbose output.
+   bool verbose = true;
+
+   /// Enable/disable partial assembly of forms.
+   bool partial_assembly = false;
+
+   /// Enable/disable numerical integration rules of forms.
+   bool numerical_integ = false;
+
+   /// The parallel mesh.
+   ParMesh *pmesh = nullptr;
+
+   /// The order of the velocity and pressure space.
+   int velorder;
+   int porder;
+   int torder;
+
+   /// Kinematic viscosity (dimensionless).
+   real_t kin_vis;
+   Coefficient * kinvisCoeff = nullptr;
+
+   Coefficient *dtCoeff = nullptr;
+
+   IntegrationRules gll_rules;
+
+   /// Velocity $H^1$ finite element collection.
+   FiniteElementCollection *vfec = nullptr;
+
+   /// Psi $H^1$ finite element collection.
+   FiniteElementCollection *psifec = nullptr;
+
+   /// Pressure $H^1$ finite element collection.
+   FiniteElementCollection *pfec = nullptr;
+
+   /// Velocity $(H^1)^d$ finite element space.
+   ParFiniteElementSpace *vfes = nullptr;
+
+   /// Psi $(H^1)^d$ finite element space.
+   ParFiniteElementSpace *psifes = nullptr;
+
+   /// Pressure $H^1$ finite element space.
+   ParFiniteElementSpace *pfes = nullptr;
+
+   ParBilinearForm *velBForm = nullptr;
+   ParBilinearForm *psiBForm = nullptr;
+   ParBilinearForm *pBForm = nullptr;
+
+   ParLinearForm *velLForm = nullptr;
+   ParLinearForm *psiLForm = nullptr;
+   ParLinearForm *pLForm = nullptr;
+
+   std::vector<ParGridFunction*> velGF;
+   std::vector<ParGridFunction*> pGF;
+   ParGridFunction psiGF;
+
+   ParGridFunction DvGF, divVelGF, pRHS;
+   VectorGridFunctionCoefficient * DvelCoeff = nullptr;
+   DivergenceGridFunctionCoefficient * divVelCoeff = nullptr;
+   GridFunctionCoefficient * pRHSCoeff = nullptr;
+   UnitVectorGridFunctionCoeff * pUnitVectorCoeff = nullptr;
+   NonLinTermVectorGridFunctionCoeff * nonlinTermCoeff = nullptr;
+   PrevVelVectorGridFunctionCoeff * prevVelLoadCoeff = nullptr;
+   
+   OperatorHandle vOp;
+   OperatorHandle psiOp;
+   OperatorHandle pOp;
+
+   Solver *velInvPC = nullptr;
+   CGSolver *velInv = nullptr;
+
+   ParLORDiscretization *lor = nullptr;
+   HypreBoomerAMG *psiInvPC = nullptr;
+   OrthoSolver *SpInvOrthoPC = nullptr;
+   CGSolver *psiInv = nullptr;
+
+   Solver *pInvPC = nullptr;
+   CGSolver *pInv = nullptr;
+
+   Vector velLF, psiLF, pLF;
+
+   // All essential attributes.
+   Array<int> vel_ess_attr;
+   Array<int> pres_ess_attr;
+
+   // All essential true dofs.
+   Array<int> vel_ess_tdof;
+   Array<int> pres_ess_tdof;
+
+   // Bookkeeping for velocity dirichlet bcs.
+   std::vector<VelDirichletBC_T> vel_dbcs;
+
+   // Print levels.
+   int pl_psolve = 0;
+   int pl_psisolve = 0;
+   int pl_velsolve = 0;
+   int pl_amg = 0;
+
+   #if defined(MFEM_USE_DOUBLE)
+   real_t rtol_psolve = 1e-10;
+   real_t rtol_psisolve = 1e-10;
+   real_t rtol_velsolve = 1e-12;
+#elif defined(MFEM_USE_SINGLE)
+   real_t rtol_psolve = 1e-9;
+   real_t rtol_psisolve = 1e-5;
+   real_t rtol_velsolve = 1e-7;
+#else
+#error "Only single and double precision are supported!"
+   real_t rtol_psolve = 1e-12;
+   real_t rtol_psisolve = 1e-6;
+   real_t rtol_velsolve = 1e-8;
+#endif
+
+   // Iteration counts.
+   int iter = 1, iter_vsolve = 0, iter_psolve = 0, iter_psisolve = 0;
+
+   // Residuals.
+   real_t res_vsolve = 0.0, res_psolve = 0.0, res_psisolve = 0.0;
+
+};
+
+} // namespace incompressible_navier
+
+} // namespace mfem
+
+#endif

miniapps/navier/stokes_solver.cpp

@@ -0,0 +1,119 @@
+#include "stokes_solver.hpp"
+
+namespace mfem {
+
+StokesOperator::StokesOperator(ParFiniteElementSpace &vel_fes,
+                               ParFiniteElementSpace &pres_fes):
+    Operator(vel_fes.GetTrueVSize()+pres_fes.GetTrueVSize()),
+    vfes(vel_fes),
+    pfes(pres_fes),
+    offsets({0, vel_fes.GetTrueVSize(), pres_fes.GetTrueVSize()}),
+    intrules(0, Quadrature1D::GaussLobatto),
+    zero_coeff(0.0)
+{
+    if (vel_fes.GetParMesh()->bdr_attributes.Size() > 0)
+    {
+        vel_ess_bdr.SetSize(vel_fes.GetParMesh()->bdr_attributes.Max());
+        vel_ess_bdr = 0.0;
+        pres_ess_bdr.SetSize(vel_fes.GetParMesh()->bdr_attributes.Max());
+        pres_ess_bdr = 0.0;
+    }
+
+    vfes.GetEssentialTrueDofs(vel_ess_bdr, vel_ess_tdofs);
+    pfes.GetEssentialTrueDofs(pres_ess_bdr, pres_ess_tdofs);
+
+    offsets.PartialSum();
+
+
+    vel_bc_gf.reset(new ParGridFunction(&vfes));
+    *vel_bc_gf = 0.0; //set the velocity grid function to zero
+
+    pres_bc_gf.reset(new ParGridFunction(&pfes));
+    *pres_bc_gf = 0.0; //set the pressure grid function to zero
+
+    // The nonlinear convective integrators use over-integration (dealiasing) as
+    // a stabilization mechanism.
+    ir_nl = intrules.Get(vfes.GetFE(0)->GetGeomType(),
+                         (int)(ceil(1.5 * 2*(vel_fes.GetOrder(0)+1) - 3)));
+
+    ir = intrules.Get(vfes.GetFE(0)->GetGeomType(),
+                      (int)(2*(vel_fes.GetOrder(0)+1) - 3));
+
+    ir_face = intrules.Get(vfes.GetFaceElement(0)->GetGeomType(),
+                           (int)(2*(vel_fes.GetOrder(0)+1) - 3));
+
+    b11_form=nullptr;
+    b22_form=nullptr;
+    b12_form=nullptr;
+    b21_form=nullptr;
+
+}
+
+
+void StokesOperator::SetVelBC(std::vector<VelDirichletBC>& vvbc)
+{
+    for(auto vbc=vvbc.begin();vbc!=vvbc.end();vbc++)
+    {
+        for (int i = 0; i < vbc->second->Size(); i++)
+        {
+            if (*(vbc->second)[i] == 1)
+            {
+                vel_ess_bdr[i] = 1;
+            }
+        }
+    }
+
+    vfes.GetEssentialTrueDofs(vel_ess_bdr, vel_ess_tdofs);
+}
+
+void StokesOperator::SetPressBC(std::vector<PresDirichletBC>& vpbc)
+{
+    for(auto pbc=vpbc.begin();pbc!=vpbc.end();pbc++)
+    {
+        for(int i=0;i<pbc->second->Size();i++){
+            if (*(pbc->second)[i] == 1)
+            {
+                vel_ess_bdr[i] = 1;
+            }
+        }
+    }
+
+    pfes.GetEssentialTrueDofs(pres_ess_bdr, pres_ess_tdofs);
+}
+
+void StokesOperator::Mult(const Vector &x, Vector &y) const
+{
+
+}
+
+void StokesOperator::Setup()
+{
+    BilinearFormIntegrator *integrator;
+
+    delete b11_form;
+    b11_form=new ParBilinearForm(&vfes);
+    integrator=new ElasticityIntegrator(zero_coeff,*viscosity);
+    integrator->SetIntRule(&ir);
+    b11_form->AddDomainIntegrator(integrator);
+
+    delete b12_form;
+    b12_form=new ParMixedBilinearForm(&pfes,&vfes);
+    integrator=new VectorDivergenceIntegrator();
+    integrator->SetIntRule(&ir);
+    b12_form->AddDomainIntegrator(integrator);
+
+    delete b21_form;
+    b21_form=new ParMixedBilinearForm(&vfes,&pfes);
+    integrator=new GradientIntegrator();
+    integrator->SetIntRule(&ir);
+    b21_form->AddDomainIntegrator(integrator);
+
+    if (matrix_free)
+    {
+        b11_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+        b12_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+        b21_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+    }
+}
+
+}

miniapps/navier/stokes_solver.hpp

@@ -0,0 +1,76 @@
+#ifndef STOKESSOLVER_H
+#define STOKESSOLVER_H
+
+#define STOKES_VERSION 0.1
+
+#include "mfem.hpp"
+
+
+namespace mfem {
+
+using VelDirichletBC = std::pair<VectorCoefficient *, Array<int> *>;
+using PresDirichletBC = std::pair<Coefficient *, Array<int> *>;
+
+class StokesOperator:public Operator
+{
+public:
+    StokesOperator(ParFiniteElementSpace &vel_fes,
+                   ParFiniteElementSpace &pres_fes);
+
+    void SetVelBC(std::vector<VelDirichletBC>& vvbc);
+    void SetPressBC(std::vector<PresDirichletBC>& vpbc);
+
+    virtual
+    void Mult(const Vector &x, Vector &y) const override;
+
+    const Array<int>& GetOffsets() const
+    {
+        return offsets;
+    }
+
+    void Setup();
+    void Assemble();
+
+
+private:
+    ParFiniteElementSpace &vfes;
+    ParFiniteElementSpace &pfes;
+
+// ParGridFunction &kinematic_viscosity;
+
+    std::unique_ptr<ParGridFunction> vel_bc_gf;
+    std::unique_ptr<ParGridFunction> pres_bc_gf;
+
+    Array<int> vel_ess_bdr;
+    Array<int> pres_ess_bdr;
+
+    Array<int> vel_ess_tdofs;
+    Array<int> pres_ess_tdofs;
+
+    bool matrix_free;
+
+    Array<int> offsets;
+
+
+    IntegrationRules intrules;
+    IntegrationRule ir; //general integraion rule
+    IntegrationRule ir_nl; //non-linear integration rule
+    IntegrationRule ir_face; //face integration rule
+
+    ConstantCoefficient zero_coeff;
+
+    std::unique_ptr<Coefficient> viscosity;
+
+    ParBilinearForm *b11_form; //velocity
+    ParBilinearForm *b22_form; //pressure
+    ParMixedBilinearForm *b12_form; //mixed (velocity,pressure)
+    ParMixedBilinearForm *b21_form; //mized (pressure,velocity)
+
+    BlockOperator* A;
+
+};
+
+
+}
+
+#endif // STOKESSOLVER_H