SteadyNSSolver: LF integrator should be added for neumann condition a…

…s well.

include/parameterized_problem.hpp

@@ -233,6 +233,10 @@ friend class RandomSampleGenerator;
    Array<int> battr;
    Array<BoundaryType> bdr_type; // abstract boundary type
 
+   /* initial condition for time-dependent problem */
+   /* size with number of variables */
+   Array<function_factory::GeneralVectorFunction *> ic_ptr;
+
    Array<function_factory::GeneralScalarFunction *> general_scalar_ptr;
    Array<function_factory::GeneralVectorFunction *> general_vector_ptr;
 

include/unsteady_ns_solver.hpp

@@ -76,6 +76,8 @@ friend class SteadyNSOperator;
    using MultiBlockSolver::SaveVisualization;
    void SaveVisualization(const int step, const double time) override;
 
+   void SetParameterizedProblem(ParameterizedProblem *problem) override;
+
    void ProjectOperatorOnReducedBasis() override
    { mfem_error("UnsteadyNSSolver::ProjectOperatorOnReducedBasis is not implemented yet!\n"); }
 

src/parameterized_problem.cpp

@@ -170,6 +170,22 @@ void ubdr(const Vector &x, Vector &y)
    y(0) = - u0 * 4.0 / ygap / ygap * (y0 - x(1)) * (y1 - x(1));
 }
 
+void uic(const Vector &x, Vector &y)
+{
+   const int dim = x.Size();
+   y.SetSize(dim);
+   y = 0.0;
+
+   double ygap = (y1 - y0);
+   y(0) = (x(1) < y0) ? 0.0 : - u0 * 4.0 / ygap / ygap * (y0 - x(1)) * (y1 - x(1));
+}
+
+void pic(const Vector &x, Vector &y)
+{
+   y.SetSize(1);
+   y = 0.0;
+}
+
 }  // namespace backward_facing_step
 
 } // namespace flow_problem
@@ -795,6 +811,11 @@ BackwardFacingStep::BackwardFacingStep()
    /* technically only vector_bdr_ptr[0] will be used. */
    vector_bdr_ptr = &(function_factory::flow_problem::backward_facing_step::ubdr);
 
+   /* pointer to initial condition */
+   ic_ptr.SetSize(2);
+   ic_ptr[0] = &(function_factory::flow_problem::backward_facing_step::uic);
+   ic_ptr[1] = &(function_factory::flow_problem::backward_facing_step::pic);
+
    param_num = 4;
 
    // Default values.

src/steady_ns_solver.cpp

@@ -385,12 +385,15 @@ void SteadyNSSolver::SetupDomainBCOperators()
       {
          int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]);
          if (idx < 0) continue;
-         if (!BCExistsOnBdr(b)) continue;
+
          // TODO: Non-homogeneous Neumann stress bc
          if (bdr_type[b] == BoundaryType::NEUMANN)
-            continue;
-
-         hs[m]->AddBdrFaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta, ud_coeffs[b]), *bdr_markers[b]);
+            hs[m]->AddBdrFaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta), *bdr_markers[b]);
+         else
+         {
+            assert(BCExistsOnBdr(b));
+            hs[m]->AddBdrFaceIntegrator(new DGLaxFriedrichsFluxIntegrator(*minus_zeta, ud_coeffs[b]), *bdr_markers[b]);
+         }
       }
    }
 

src/unsteady_ns_solver.cpp

@@ -184,6 +184,21 @@ void UnsteadyNSSolver::SaveVisualization(const int step, const double time)
    MultiBlockSolver::SaveVisualization(step, time);
 }
 
+void UnsteadyNSSolver::SetParameterizedProblem(ParameterizedProblem *problem)
+{
+   SteadyNSSolver::SetParameterizedProblem(problem);
+
+   /* set up initial condition */
+   VectorFunctionCoefficient u_ic(vdim[0], problem->ic_ptr[0]);
+   VectorFunctionCoefficient p_ic(1, problem->ic_ptr[1]);
+
+   for (int m = 0; m < numSub; m++)
+   {
+      vels[m]->ProjectCoefficient(u_ic);
+      ps[m]->ProjectCoefficient(p_ic);
+   }
+}
+
 double UnsteadyNSSolver::ComputeCFL(const double dt_)
 {
    Vector ux, uy, uz;

utils/gmsh/bfs.geo

@@ -0,0 +1,54 @@
+// -----------------------------------------------------------------------------
+//
+// Backward facing step mesh
+//
+// -----------------------------------------------------------------------------
+
+// As usual, we start by defining some variables:
+
+Lc1 = 0.125;
+x1 = 2.0; xL = 6.0;
+yb = -1.0; yt = 1.0;
+
+Point(1) = {0, 0, 0, Lc1};
+Point(2) = {x1, 0, 0, Lc1};
+Point(3) = {x1, yb, 0, Lc1};
+Point(4) = {xL, yb, 0, Lc1};
+Point(5) = {xL, yt, 0, Lc1};
+Point(6) = {0, yt, 0, Lc1};
+
+Line(1)  = {1, 2};
+Line(2)  = {2, 3};
+Line(3)  = {3, 4};
+Line(4)  = {4, 5};
+Line(5)  = {5, 6};
+Line(6)  = {6, 1};
+
+// The third elementary entity is the surface. In order to define a simple
+// rectangular surface from the four curves defined above, a curve loop has
+// first to be defined. A curve loop is also identified by a tag (unique amongst
+// curve loops) and defined by an ordered list of connected curves, a sign being
+// associated with each curve (depending on the orientation of the curve to form
+// a loop):
+
+Curve Loop(1) = {1, 2, 3, 4, 5, 6};
+
+// We can then define the surface as a list of curve loops (only one here,
+// representing the external contour, since there are no holes--see `t4.geo' for
+// an example of a surface with a hole):
+
+Plane Surface(1) = {1};
+
+Physical Curve(1) = {6};
+Physical Curve(2) = {1,2,3,5};
+Physical Curve(3) = {4};
+Physical Surface(1) = {1};
+
+Mesh.Algorithm = 1;
+Mesh.AnisoMax = 0;
+Mesh.ElementOrder = 3;
+Mesh.MshFileVersion = 2.2;
+// Mesh.LcIntegrationPrecision = 1.0e-15;
+Mesh 2;
+Save "bfs.msh";
+// Exit;