Some comments to make it clear what some newer parts are doing.

miniapps/navier/makefile

@@ -24,7 +24,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 NAVIER_COMMON_SRC = navier_solver.cpp ortho_solver.cpp
 NAVIER_COMMON_OBJ = $(NAVIER_COMMON_SRC:.cpp=.o)
 
-PAR_MINIAPPS = navier_mms navier_kovasznay navier_tgv navier_shear navier_3dfoc
+PAR_MINIAPPS = navier_mms navier_kovasznay navier_tgv navier_shear navier_3dfoc navier_ldc
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS =

miniapps/navier/navier_ldc.cpp

@@ -0,0 +1,342 @@
+// Copyright (c) 2010-2020, 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.
+//
+// Navier MMS example
+//
+// A manufactured solution is defined as
+//
+// u = [pi * sin(t) * sin(pi * x)^2 * sin(2 * pi * y),
+//      -(pi * sin(t) * sin(2 * pi * x)) * sin(pi * y)^2].
+//
+// p = cos(pi * x) * sin(t) * sin(pi * y)
+//
+// The solution is used to compute the symbolic forcing term (right hand side),
+// of the equation. Then the numerical solution is computed and compared to the
+// exact manufactured solution to determine the error.
+
+#ifdef WINDOWS
+#include <direct.h>
+#define GetCurrentDir _getcwd
+#else
+#include <unistd.h>
+#define GetCurrentDir getcwd
+#endif
+
+#include "navier_solver.hpp"
+#include <fstream>
+#include <iostream>
+#include <string>
+
+std::string get_current_dir() {
+   char buff[FILENAME_MAX]; //create string buffer to hold path
+   GetCurrentDir( buff, FILENAME_MAX );
+   std::string current_working_dir(buff);
+   return current_working_dir;
+}
+
+
+
+using namespace mfem;
+using namespace navier;
+
+
+struct s_NavierContext
+{
+   int ser_ref_levels = 1;
+   int order = 5;
+   double kinvis = .0001;
+   double t_final = .5;
+   double dt = 0.25e-3;
+   bool pa = true;
+   bool ni = false;
+   bool visualization = false;
+   bool checkres = false;
+   double lid = .1;
+   double top = 0.0;
+   double rel_tol = 1e-6;
+} ctx;
+
+void vel(const Vector &x, double t, Vector &u)
+{
+   double xi = x(0);
+   double yi = x(1);
+
+   if (yi==ctx.top)
+   {
+      u(0) = ctx.lid;
+   }
+   else
+   {
+      u(0) = 0.0;
+   }
+   
+   u(1) = 0.0;
+}
+
+// double p(const Vector &x, double t)
+// {
+//    double xi = x(0);
+//    double yi = x(1);
+
+//    return cos(M_PI * xi) * sin(t) * sin(M_PI * yi);
+// }
+
+// void accel(const Vector &x, double t, Vector &u)
+// {
+//    double xi = x(0);
+//    double yi = x(1);
+
+//    u(0) = M_PI * sin(t) * sin(M_PI * xi) * sin(M_PI * yi)
+//           * (-1.0
+//              + 2.0 * pow(M_PI, 2.0) * sin(t) * sin(M_PI * xi)
+//              * sin(2.0 * M_PI * xi) * sin(M_PI * yi))
+//           + M_PI
+//           * (2.0 * ctx.kinvis * pow(M_PI, 2.0)
+//              * (1.0 - 2.0 * cos(2.0 * M_PI * xi)) * sin(t)
+//              + cos(t) * pow(sin(M_PI * xi), 2.0))
+//           * sin(2.0 * M_PI * yi);
+
+//    u(1) = M_PI * cos(M_PI * yi) * sin(t)
+//           * (cos(M_PI * xi)
+//              + 2.0 * ctx.kinvis * pow(M_PI, 2.0) * cos(M_PI * yi)
+//              * sin(2.0 * M_PI * xi))
+//           - M_PI * (cos(t) + 6.0 * ctx.kinvis * pow(M_PI, 2.0) * sin(t))
+//           * sin(2.0 * M_PI * xi) * pow(sin(M_PI * yi), 2.0)
+//           + 4.0 * pow(M_PI, 3.0) * cos(M_PI * yi) * pow(sin(t), 2.0)
+//           * pow(sin(M_PI * xi), 2.0) * pow(sin(M_PI * yi), 3.0);
+// }
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   OptionsParser args(argc, argv);
+   args.AddOption(&ctx.ser_ref_levels,
+                  "-rs",
+                  "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&ctx.order,
+                  "-o",
+                  "--order",
+                  "Order (degree) of the finite elements.");
+   args.AddOption(&ctx.kinvis, "-k", "--kinematic-viscosity", 
+                  "Kinematic viscosity of the simulated fluid.");
+   args.AddOption(&ctx.lid, "-lid", "--lidspeed", 
+                  "Lid speed of the cavity.");
+   args.AddOption(&ctx.dt, "-dt", "--time-step", "Time step.");
+   args.AddOption(&ctx.t_final, "-tf", "--final-time", "Final time.");
+   args.AddOption(&ctx.rel_tol, "-rt", "--relative-tolerance", "Relative tolerance for stopping criterion.");
+   args.AddOption(&ctx.pa,
+                  "-pa",
+                  "--enable-pa",
+                  "-no-pa",
+                  "--disable-pa",
+                  "Enable partial assembly.");
+   args.AddOption(&ctx.ni,
+                  "-ni",
+                  "--enable-ni",
+                  "-no-ni",
+                  "--disable-ni",
+                  "Enable numerical integration rules.");
+   args.AddOption(&ctx.visualization,
+                  "-vis",
+                  "--visualization",
+                  "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(
+      &ctx.checkres,
+      "-cr",
+      "--checkresult",
+      "-no-cr",
+      "--no-checkresult",
+      "Enable or disable checking of the result. Returns -1 on failure.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(mfem::out);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(mfem::out);
+   }
+
+   std::string meshName(mesh_file);
+
+   //Determining where the lid is located.
+   Mesh *mesh = new Mesh(mesh_file);
+   for (int m = 0; m < mesh->GetNV(); m++){
+      double& vert = *(mesh->GetVertex(m)+1);
+      if (vert > ctx.top) {
+         ctx.top = vert;
+      }
+   }
+
+   // mesh->EnsureNodes();
+   // GridFunction *nodes = mesh->GetNodes();
+   // *nodes *= 2.0;
+   // *nodes -= 1.0;
+
+   for (int i = 0; i < ctx.ser_ref_levels; ++i)
+   {
+      mesh->UniformRefinement();
+   }
+
+   if (mpi.Root())
+   {
+      std::cout << "Number of elements: " << mesh->GetNE() << std::endl;
+   }
+
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Create the flow solver.
+   NavierSolver naviersolver(pmesh, ctx.order, ctx.kinvis);
+   naviersolver.EnablePA(ctx.pa);
+   naviersolver.EnableNI(ctx.ni);
+
+   // Set the initial condition.
+   ParGridFunction *u_ic = naviersolver.GetCurrentVelocity();
+   VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
+   u_ic->ProjectCoefficient(u_excoeff);
+
+   // FunctionCoefficient p_excoeff(p);
+
+   // Add Dirichlet boundary conditions to velocity space restricted to
+   // selected attributes on the mesh.
+   Array<int> attr(pmesh->bdr_attributes.Max());
+   attr = 1;
+   naviersolver.AddVelDirichletBC(vel, attr);
+
+   // Array<int> domain_attr(pmesh->attributes.Max());
+   // domain_attr = 1;
+   // naviersolver.AddAccelTerm(accel, domain_attr);
+
+   double t = 0.0;
+   double dt = ctx.dt;
+   double t_final = ctx.t_final;
+   bool last_step = false;
+
+   naviersolver.Setup(dt);
+
+   // double err_u = 0.0;
+   // double err_p = 0.0;
+   ParGridFunction *u_gf = nullptr;
+   ParGridFunction *p_gf = nullptr;
+   u_gf = naviersolver.GetCurrentVelocity();
+   p_gf = naviersolver.GetCurrentPressure();
+
+   std::string outputPath = get_current_dir();
+
+   ParaViewDataCollection pvdc("ldc_output", pmesh);
+   //pvdc.SetDataFormat(VTKFormat::BINARY32);
+   pvdc.SetDataFormat(VTKFormat::ASCII);
+   pvdc.SetHighOrderOutput(true);
+   pvdc.SetLevelsOfDetail(ctx.order);
+   pvdc.SetCycle(0);
+   pvdc.SetTime(t);
+   pvdc.RegisterField("velocity", u_gf);
+   pvdc.RegisterField("pressure", p_gf);
+   //pvdc.RegisterField("vorticity", &w_gf);
+   pvdc.Save();
+
+   ParGridFunction prev_u;
+   ParGridFunction prev_p;
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (t + dt >= t_final - dt / 2)
+      {
+         last_step = true;
+         std::ofstream myFile;
+         myFile.open("NoConverge.txt");
+         myFile <<  ctx.t_final << "  #This is final time allowed. \n";
+         myFile << ctx.rel_tol << "  #This is the relative tolerance between steps. \n";
+         myFile << ctx.kinvis << "  #This is the kinematic viscosity. \n";
+         myFile.close();
+      }
+
+      naviersolver.Step(t, dt, step);
+      //Take out this if statement if you want a pure transient solver.
+      if (step > 0)
+      {
+         double err_u = u_gf->DistanceTo(prev_u);
+         double err_p = p_gf->DistanceTo(prev_p);
+         if (err_u < ctx.rel_tol && err_p < ctx.rel_tol) //Relative tolerances between steps
+         {
+            last_step = true;
+         }
+      }
+
+      prev_u = *u_gf;
+      prev_p = *p_gf;
+
+      //Modify this if statement if you want more data saved.
+      //if (step % 50 == 0)
+      if (last_step)
+      {
+         pvdc.SetCycle(step);
+         pvdc.SetTime(t);
+         pvdc.Save();    
+
+      }
+      // Compare against exact solution of velocity and pressure.
+      // u_excoeff.SetTime(t);
+      // p_excoeff.SetTime(t);
+      // err_u = u_gf->ComputeL2Error(u_excoeff);
+      // err_p = p_gf->ComputeL2Error(p_excoeff);
+
+      // if (mpi.Root())
+      // {
+      //    printf("%11s %11s %11s %11s\n", "Time", "dt", "err_u", "err_p");
+      //    printf("%.5E %.5E %.5E %.5E err\n", t, dt, err_u, err_p);
+      //    fflush(stdout);
+      // }
+   }
+
+   if (ctx.visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank()
+               << "\n";
+      sol_sock << "solution\n" << *pmesh << *u_ic << std::flush;
+   }
+
+   naviersolver.PrintTimingData();
+
+   // Test if the result for the test run is as expected.
+   // if (ctx.checkres)
+   // {
+   //    double tol = 1e-3;
+   //    if (err_u > tol || err_p > tol)
+   //    {
+   //       if (mpi.Root())
+   //       {
+   //          mfem::out << "Result has a larger error than expected."
+   //                    << std::endl;
+   //       }
+   //       return -1;
+   //    }
+   // }
+
+   delete pmesh;
+
+   return 0;
+}