jeffery_hamel.cc

// An FEM formulation for solving the Jeffery-Hamel ODE:
// f''' + 2 * alpha * Re * f * f' + 4 * alpha^2 * f' = 0
// using C1 Hermite cubic and quartic finite elements.

// Libmesh includes
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/fe.h"
#include "libmesh/getpot.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/function_base.h"
#include "libmesh/exact_solution.h"

#include <ctype.h> // isspace

// Bring in everything from the libMesh namespace
using namespace libMesh;

// --------------------------------------------------------------------------------
// Helper classes

/**
 * This class is responsible for reading CSV data files with reference
 * solution data generated by the bvp_solver code, and linearly
 * interpolating it to provide function and gradient values for error
 * estimation routines.
 */
class JefferyHamelData
{
public:
  /**
   * The constructor reads the CSV file and stores everything in
   * vectors internally.
   */
  JefferyHamelData(const std::string & filename);

  /**
   * Returns true if there was no data file available.
   */
  bool is_dummy() { return _is_dummy; }

  // Returns the linearly interpolated value of f at the given point
  // x.  This assumes the bvp_solver data is sorted in increasing
  // order.
  Real f(Real x) const { return interpolate_value(x, _f); }

  // Returns the linearly interpolated value of f at the given point
  // x.  This assumes the bvp_solver data is sorted in increasing
  // order.
  Real fprime(Real x) const { return interpolate_value(x, _fprime); }

private:
  // Linearly interpolates a value from 'data_vector' at the point x.
  Real interpolate_value(Real x, const std::vector<Real> & data_vector) const;

  bool _is_dummy;
  std::vector<Real> _x;
  std::vector<Real> _f;
  std::vector<Real> _fprime;
};



/**
 * FunctionBase-derived object that can be attached to an
 * ExactSolution object and used for estimating error.  Has a
 * reference to a JefferyHamelData object that allows you to linearly
 * interpolate between tabulated values.
 */
class JefferyHamelFunction : public FunctionBase<Number>
{
public:
  JefferyHamelFunction(const JefferyHamelData & data) :
    FunctionBase<Number>(),
    _data(data)
  {}

  /**
   * Returns a new deep copy of the function.
   */
  virtual UniquePtr<FunctionBase<Number> > clone () const libmesh_override
  {
    return UniquePtr<FunctionBase<Number> >(new JefferyHamelFunction(_data));
  }

  /**
   * Returns the function's value at requested point and time.
   * Uses the underlying _data object to do interpolation.
   */
  virtual Number operator() (const Point & p,
                             const Real /*time*/) libmesh_override
  {
    return _data.f(p(0));
  }

  /**
   * Interface for vector-valued functions, fills up the passed in
   * 'output' vector with function values.  This is the function that
   * actually gets called by the ExactSolution object.
   */
  virtual void operator() (const Point & p,
                           const Real /*time*/,
                           DenseVector<Number> & output) libmesh_override
  {
    output.resize(1);
    output(0) = _data.f(p(0));
  }

private:
  const JefferyHamelData & _data;
};



/**
 * FunctionBase-derived object that can be attached to an
 * ExactSolution object and used for estimating error.  Has a
 * reference to a JefferyHamelData object that allows you to linearly
 * interpolate between tabulated values.
 */
class JefferyHamelGradient : public FunctionBase<Gradient>
{
public:
  JefferyHamelGradient(const JefferyHamelData & data) :
    FunctionBase<Gradient>(),
    _data(data)
  {}

  /**
   * Returns a new deep copy of the function.
   */
  virtual UniquePtr<FunctionBase<Gradient> > clone () const libmesh_override
  {
    return UniquePtr<FunctionBase<Gradient> >(new JefferyHamelGradient(_data));
  }

  /**
   * Returns the function's value at requested point and time.
   * Uses the underlying _data object to do interpolation.
   */
  virtual Gradient operator() (const Point & p,
                               const Real /*time*/) libmesh_override
  {
    return Gradient(_data.fprime(p(0)), 0., 0.);
  }

  /**
   * Interface for vector-valued functions, fills up the passed in
   * 'output' vector with function values.  This is the function that
   * actually gets called by the ExactSolution object.
   */
  virtual void operator() (const Point & p,
                           const Real /*time*/,
                           DenseVector<Gradient> & output) libmesh_override
  {
    output.resize(1);
    output(0) = Gradient(_data.fprime(p(0)), 0., 0.);
  }

private:
  const JefferyHamelData & _data;
};



// --------------------------------------------------------------------------------
// Misc. function prototypes

// Assembly function prototype.
void assemble_1D(EquationSystems & es, const std::string & system_name);

// Function used for setting the initial guess.
Number initial_value (const Point & p, const Parameters &, const std::string &, const std::string &);

// Function used for setting the initial guess gradient.
Gradient initial_gradient (const Point & p, const Parameters &, const std::string &, const std::string &);

// Write a Matplotlib script which generates a PDF line plot of the solution.
void write_lineplot(System & system);



// --------------------------------------------------------------------------------
// Main Program
int main(int argc, char ** argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

  // Create a mesh, with dimension to be overridden later, on the
  // default MPI communicator.
  Mesh mesh(init.comm());

  // Object for parsing command line arguments.
  GetPot command_line (argc, argv);

  // Number of elements, can be set on the command line.
  int n = 10;
  if (command_line.search(1, "-n"))
    n = command_line.next(n);

  // Set up constants used in the problem.
  // Positive Re -> diverging flow (away from the origin)
  // Negative Re -> converging flow (toward the origin)

  // Case I
  const Real Re = 30;
  const Real alpha_degrees = 15;

  // Case II - the .csv file for this case is not checked in, you will
  // have to regenerate it with the bvp_solver code.
  // const Real Re = 110;
  // const Real alpha_degrees = 3;

  // Case III
  // const Real Re = -80;
  // const Real alpha_degrees = 5;

  // Build 1D mesh with the requested number of elements.
  MeshTools::Generation::build_line(mesh, n, 0., 1., EDGE2);

  // Build up filename to be read in based on Re and alpha.  This
  // currently assumes that both Re and alpha that we want to test can
  // be expressed as integers, e.g. "bvp_solver_data_Re_30_alpha_15.csv"
  std::stringstream data_file;
  data_file << "bvp_solver_data_Re_"
            << static_cast<int>(Re) << "_alpha_"
            << static_cast<int>(alpha_degrees) << ".csv";

  // Create a data reader object that will be shared by both
  JefferyHamelData jhdata(data_file.str());

  // Construct a FunctionBase-derived object that can be used in
  // conjunction with the ExactSolution class.
  JefferyHamelFunction jhf(jhdata);
  JefferyHamelGradient jhg(jhdata);

  // Define the equation systems object and the system we are going
  // to solve.
  EquationSystems equation_systems(mesh);
  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem>("Jeffery-Hamel");

  // Add a variable "u" to the system.
  system.add_variable("u", THIRD, HERMITE);

  // Give the system a pointer to the matrix assembly function. This
  // will be called when needed by the library.
  system.attach_assemble_function(assemble_1D);

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Set the linear solver tolerance used by the underlying linear
  // solver.  For nonlinear problems, this should really adapt to the
  // size of the current nonlinear residual.
  equation_systems.parameters.set<Real> ("linear solver tolerance") = 1.e-3;

  // Set problem parameters on the equation_systems object
  equation_systems.parameters.set<Real> ("Re") = Re;
  equation_systems.parameters.set<Real> ("alpha_degrees") = alpha_degrees;

  // We'll set this flag when convergence is (hopefully) achieved.
  bool converged = false;

  // The maximum number of nonlinear steps.
  const unsigned int n_nonlinear_steps = 10;

  // Project the initial guess into the current_local_solution vector.
  // This seems to give convergence in slightly fewer Newton iterations
  // than an initial guess of zero.
  system.project_solution(initial_value, initial_gradient, equation_systems.parameters);

  // Extra solution vector used to determine convergence.
  UniquePtr<NumericVector<Number> > last_nonlinear_soln (system.solution->clone());

  // The finer the grid, the less tightly we seem to be able to
  // converge the residual...
  const Real nonlinear_residual_absolute_tol = 1.e-9;

  // Begin the Newton iterations.
  for (unsigned int l=0; l<n_nonlinear_steps; ++l)
    {
      // Always use an initial guess of zero for the solution to the
      // nonlinear update.
      system.solution->zero();

      // Solve the linear system of equations.
      system.solve();

      // How many iterations were required to solve the linear system?
      const unsigned int n_linear_iterations = system.n_linear_iterations();

      // Doing 0 linear iterations typically implies the solver failed for some reason...
      if (n_linear_iterations == 0)
        {
          libMesh::out << "Previous solve likely failed, ignoring most recent solution." << std::endl;
          break;
        }

      // The norm of the rhs vector is the current residual norm.
      const Real norm_resid = system.rhs->l2_norm();
      libMesh::out << l << ": Nonlinear residual norm = " << norm_resid << std::endl;

      // Update the solution vector.  Add -1 times the current
      // solution to the "last_nonlinear_soln" vector, since we
      // actually solved for -du.
      last_nonlinear_soln->add (-1., *system.solution);
      *system.solution = *last_nonlinear_soln;
      system.solution->close();

      // Make sure current_local_solution is up to date for use in the
      // finite element assembly function.
      system.update();

      if (norm_resid < nonlinear_residual_absolute_tol)
        {
          libMesh::out << "Convergence detected, stopping nonlinear iterations." << std::endl;
          converged = true;
          break;
        }
    }

  // Throw an error if the solver failed to converge.
  if (!converged)
    libmesh_error_msg("Nonlinear iterations failed to converge!");

  // Make a line plot of the solution.
  write_lineplot(system);

  // Compute the error if a reference solution is available.
  if (!jhdata.is_dummy())
    {
      // Compute the error between the computed and reference solution in
      // an actual norm, not just as pointwise values.
      ExactSolution exact_sol(equation_systems);
      exact_sol.attach_exact_value(/*sys_num=*/0, &jhf);
      exact_sol.attach_exact_deriv(/*sys_num=*/0, &jhg);

      // Use higher quadrature order for more accurate error results
      exact_sol.extra_quadrature_order(2);

      // Compute the error.
      exact_sol.compute_error("Jeffery-Hamel", "u");

      // Print out the error values
      libMesh::out << std::scientific << std::setprecision(16)
                   << "L2-Error,               H1-Error\n"
                   << exact_sol.l2_error("Jeffery-Hamel", "u")
                   << ", "
                   << exact_sol.h1_error("Jeffery-Hamel", "u")
                   << std::endl;
    }

  return 0;
}

// --------------------------------------------------------------------------------
// Helper function definitions.

// Function used setting the initial guess.
Number initial_value (const Point & p,
                      const Parameters & /*parameters*/,
                      const std::string &,
                      const std::string &)
{
  return 1. - p(0);
}

// Function used setting the initial gradient guess.
Gradient initial_gradient (const Point & /*p*/,
                           const Parameters & /*parameters*/,
                           const std::string &,
                           const std::string &)
{
  return Gradient(-1, 0, 0);
}

// Helper function for constraining Jacobian and RHS vector for imposing BCs.
void constrain_dof (DenseMatrix<Number> & Ke,
                    DenseVector<Number> & Fe,
                    unsigned local_dof,
                    Real current_value,
                    Real desired_value)
{
  // Set 1 on diagonal.
  Ke(local_dof, local_dof) = 1.0;

  // Set u-g in residual vector.
  Fe(local_dof) = current_value - desired_value;

  // Zero other entries in row.
  for (unsigned int j=0; j<Ke.n(); j++)
    if (j != local_dof)
      Ke(local_dof, j) = 0.0;
}

// Helper function which writes a Matplotlib script
void write_lineplot(System & system)
{
  // Simply write solution values to a file that can be plotted with Matplotlib.
  // We also need x-values, so it's better to loop over the Mesh for this...
  std::ofstream out("jeffery_hamel.py");

  out << "#!/usr/bin/env python" << std::endl;
  out << "import matplotlib.pyplot as plt" << std::endl;
  out << "import numpy as np" << std::endl;
  out << "fig = plt.figure()" << std::endl;
  out << "ax1 = fig.add_subplot(111)" << std::endl;
  out << "data = [" << std::endl;

  // Write values at pre-determined locations.  This vector
  // initialization requires C++11.
  std::vector<Real> output_x = {0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.};
  for (unsigned ctr=0; ctr < output_x.size(); ++ctr)
    {
      Point x(output_x[ctr], 0, 0);
      Number val = system.point_value(/*var=*/0, x);
      Gradient dfdx = system.point_gradient(/*var=*/0, x);

      std::ostringstream line;
      line << std::scientific << std::setprecision(16)
           << output_x[ctr]
           << ", "
           << val
           << ", "
           << dfdx(0)
           << ","
           << std::endl;

      out << line.str();

      // For extra feedback, also write the results to the screen.
      libMesh::out << line.str();
    }

  out << "]" << std::endl;
  out << "xdata = data[0::3]" << std::endl;
  out << "ydata = data[1::3]" << std::endl;
  out << "ax1.plot(xdata, ydata, color='darkblue', marker='o', linestyle='-', linewidth=2)" << std::endl;
  out << "plt.xlim((0,1))" << std::endl;
  out << "ax1.set_xlabel(r'$\\eta$')" << std::endl;
  out << "ax1.set_ylabel(r'$f(\\eta)$')" << std::endl;
  out << "plt.savefig('jeffery_hamel.pdf', format='pdf')" << std::endl;
}

// Define the matrix assembly function for the 1D PDE we are solving
void assemble_1D(EquationSystems & es,
                 const std::string & libmesh_dbg_var(system_name))
{
  // It is a good idea to check we are solving the correct system
  libmesh_assert_equal_to (system_name, "Jeffery-Hamel");

  // Get a reference to the mesh object
  const MeshBase & mesh = es.get_mesh();

  // The dimension we are using, i.e. dim==1
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the system we are solving
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Jeffery-Hamel");

  // Get a reference to the DofMap object for this system.
  const DofMap & dof_map = system.get_dof_map();

  // Get a constant reference to the finite element for variable 0.
  FEType fe_type = dof_map.variable_type(0);

  // Build a finite element object of the specified type.
  UniquePtr<FEBase> fe(FEBase::build(dim, fe_type));

  // This problem has a term "f' * f * test" which is a
  // polynomial of order:
  // (p-1) + p + p = 3*p - 1,
  // or
  // .) 8th order (for cubic Hermites)
  // .) 11th order (for quartic Hermites)
  // However, using EIGHTH vs. ELEVENTH (5 pt vs. 6 pt) quadrature did
  // not seem to make much of a difference in the results, so this is
  // hard-coded to 8th.
  Order required_quadrature_order = EIGHTH;
  QGauss qrule(dim, required_quadrature_order);
  fe->attach_quadrature_rule (&qrule);

  // FE stuff for boundary integration.
  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  QGauss qface(dim-1, required_quadrature_order);
  fe_face->attach_quadrature_rule (&qface);

  // Request cell-specific data to be computed.

  // The element Jacobian * quadrature weight at each integration point.
  const std::vector<Real> & JxW = fe->get_JxW();

  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real> > & phi = fe->get_phi();

  // The element shape function gradients evaluated at the quadrature points.
  const std::vector<std::vector<RealGradient> > & dphi = fe->get_dphi();

  // The element shape function second derivatives (Hessians)
  // evaluated at the quadrature points.
  const std::vector<std::vector<RealTensor> > & d2phi = fe->get_d2phi();

  // The value of the shape function derivatives at the quadrature points.
  const std::vector<std::vector<RealGradient> > & dphi_face = fe_face->get_dphi();

  // Declare a dense matrix and dense vector to hold the element matrix
  // and right-hand-side contribution
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  // This vector will hold the degree of freedom indices for the element.
  // These define where in the global system the element degrees of freedom
  // get mapped.
  std::vector<dof_id_type> dof_indices;

  // Get constants used in the equation.
  const Real alpha = es.parameters.get<Real> ("alpha_degrees") * libMesh::pi / 180;
  const Real Re = es.parameters.get<Real> ("Re");

  // Loop over all the active elements, computing the matrix and
  // right-hand-side contribution from each.
  MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
  const MeshBase::const_element_iterator el_end = mesh.active_local_elements_end();

  for ( ; el != el_end; ++el)
    {
      // It is convenient to store a pointer to the current element
      const Elem * elem = *el;

      // Get the degree of freedom indices for the current element.
      dof_map.dof_indices(elem, dof_indices);

      // Compute the element-specific data for the current element.
      fe->reinit(elem);

      // Store the number of local degrees of freedom contained in this element
      const unsigned int n_dofs = dof_indices.size();

      // Zero out the element stiffness matrix and rhs in preparation for filling.
      Ke.resize(n_dofs, n_dofs);
      Fe.resize(n_dofs);

      // Loop over quadrature points to handle numerical integration.
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          // Reconstruct the value and gradient at the current qp.
          Number u = 0.;
          Gradient grad_u;

          for (unsigned int i=0; i<n_dofs; ++i)
            {
              u += phi[i][qp] * system.current_solution (dof_indices[i]);
              grad_u.add_scaled (dphi[i][qp], system.current_solution (dof_indices[i]));
            }

          // Short cut notation
          const Number dudx = grad_u(0);

          // Build the element matrix and right-hand-side using loops
          // to integrate the test functions (i) against the trial
          // functions (j).
          for (unsigned int i=0; i<phi.size(); ++i)
            {
              // Short cut notation
              const Number phi_i_xx = d2phi[i][qp](0,0);

              Fe(i) += JxW[qp] * dudx * (phi_i_xx +
                                         2 * alpha * Re * u * phi[i][qp] +
                                         4 * alpha * alpha * phi[i][qp]);

              for (unsigned int j=0; j<phi.size(); ++j)
                {
                  Ke(i,j) += JxW[qp] * (dudx * (2 * alpha * Re * phi[j][qp] * phi[i][qp]) +
                                        dphi[j][qp](0) * (phi_i_xx + 2 * alpha * Re * u * phi[i][qp] + 4 * alpha * alpha * phi[i][qp]));
                }
            }
        }

      // Loop over the sides of this element. For a 1D element, the "sides"
      // are defined as the nodes on each edge of the element, i.e. 1D elements
      // have 2 sides.
      for (unsigned int s=0; s<elem->n_sides(); s++)
        if (elem->neighbor_ptr(s) == libmesh_nullptr)
          {
            // Compute the shape function values on the element
            // face.
            fe_face->reinit(elem, s);

            // BCs:
            // On the left (s==0) side, we set:
            // u(0)  = 1
            // u'(0) = 0
            //
            // On the right hand (s==1) side, we set:
            // u(1)  = 0

            // Imposed solution value.
            Real value = (s==0 ? 1. : 0.);

            // On the left, constrain the first and second dofs as
            // they correspond to the value and first derivative on
            // Hermite elements.
            if (s == 0)
              {
                // Value Dirichlet BC
                constrain_dof (Ke, Fe,
                               /*local_dof=*/0,
                               /*current_value=*/system.current_solution(dof_indices[0]),
                               /*desired_value=*/value);

                // Gradient Dirichlet BC
                constrain_dof (Ke, Fe,
                               /*local_dof=*/1,
                               /*current_value=*/system.current_solution(dof_indices[1]),
                               /*desired_value=*/0.0);
              }

            if (s == 1)
              {
                // Value Dirichlet BC
                constrain_dof (Ke, Fe,
                               /*local_dof=*/2,
                               /*current_value=*/system.current_solution(dof_indices[2]),
                               /*desired_value=*/value);

                // Handle the -dudx(1) * dphidx(1) contribution.  We do
                // this as a quadrature loop, although in 1D there is only
                // 1 qp...
                for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
                    // Reconstruct gradient at the current qp.
                    Gradient grad_u;

                    for (unsigned int i=0; i<dphi_face.size(); ++i)
                      grad_u.add_scaled (dphi_face[i][qp], system.current_solution (dof_indices[i]));

                    for (unsigned int i=0; i<dphi_face.size(); i++)
                      {
                        Fe(i) -= grad_u(0) * dphi_face[i][qp](0);

                        for (unsigned int j=0; j<dphi_face.size(); j++)
                          Ke(i,j) -= dphi_face[j][qp](0) * dphi_face[i][qp](0);
                      }
                  }
              }
          }

      // Constrain the element vector and rhs.  This is required if
      // any constraints are active, just as with adaptive mesh
      // refinement or DirichletBoundary objects.
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

      // Add Ke and Fe to the global matrix and right-hand-side.
      system.matrix->add_matrix(Ke, dof_indices);
      system.rhs->add_vector(Fe, dof_indices);
    }
}


// --------------------------------------------------------------------------------
// JefferyHamelData member functions.

// Constructor reads in data file and sets up data arrays.
JefferyHamelData::JefferyHamelData(const std::string & filename) :
  _is_dummy(false)
{
  std::ifstream in(filename.c_str());

  if (!in.good())
    {
      libMesh::out << "\n" << std::endl;
      libMesh::out << "Warning: input file " << filename << " not found." << std::endl;
      libMesh::out << "Building dummy error computation object." << std::endl;
      libMesh::out << "The solver will run but the error will not be computed." << std::endl;
      libMesh::out << "\n" << std::endl;
      _is_dummy = true;
      return;
    }

  // Read CSV file line by line.
  std::string line;
  while (true)
    {
      // Try to read something.  This may set EOF!
      std::getline(in, line);

      if (in)
        {
          // Remove any whitespace characters from line.  This
          // should leave only numbers with commas separating them.
          line.erase(std::remove_if(line.begin(), line.end(), isspace), line.end());

          // Make a stream out of the modified line so we can stream values
          // from it in the usual way.
          std::stringstream ss(line);
          Real x, f, fprime;
          char c;
          ss >> x >> c >> f >> c >> fprime;

          // Store values in our internal vectors
          _x.push_back(x);
          _f.push_back(f);
          _fprime.push_back(fprime);

          // Done processing the line, go to the next while loop iteration
          continue;
        } // if (in)

      // If !in, check to see if EOF was set.  If so, break out
      // of while loop.
      if (in.eof())
        break;

      // If !in and !in.eof(), stream is in a bad state!
      libmesh_error_msg("Stream is bad! Perhaps the file does not exist?");
    } // end while(true)
}


Real
JefferyHamelData::interpolate_value(Real x, const std::vector<Real> & data_vector) const
{
  // This might be a dummy object with no data.
  if (_is_dummy)
    return 0.;

  // Returns an iterator to the first element that is *not* less than value.
  std::vector<Real>::const_iterator lb = std::lower_bound(_x.begin(), _x.end(), x);

  // If x<0, lower_bound returns begin()
  // If x>1, lower_bound returns end().
  if (lb == _x.begin() || lb == _x.end())
    libmesh_error_msg("Requested point was not in the range!");

  // Compute index so we can access other arrays in the same position.
  unsigned idx = std::distance(_x.begin(), lb);

  // Linear slope approximation
  Real slope = (data_vector[idx] - data_vector[idx-1]) / (_x[idx] - _x[idx-1]);

  // dx > 0
  Real dx = x - _x[idx-1];

  Real interpolated_value = data_vector[idx-1] + slope*dx;

  return interpolated_value;
}