time_solver_test_common.h

Go to the documentation of this file.

#include <libmesh/dof_map.h>
#include <libmesh/enum_solver_type.h>
#include <libmesh/enum_preconditioner_type.h>
#include <libmesh/fem_system.h>
#include <libmesh/newton_solver.h>
#include <libmesh/numeric_vector.h>
#include <libmesh/parallel.h>
 
#include "test_comm.h"
#include "libmesh_cppunit.h"
 
#include <memory>
 
using namespace libMesh;
 
template<typename TimeSolverType>
class TimeSolverTestImplementation
{
protected:
 
  // Any specialized initialization that's needed for the test
  virtual void aux_time_solver_init( TimeSolverType & /*time_solver*/ ){}
 
  // Implementation for solving ODE of SystemType
  // Note this test assumes that the time integrator gets the *exact* solution
  // to within floating point tolerance.
  template<typename SystemType>
  void run_test_with_exact_soln(Real deltat, unsigned int n_timesteps)
  {
    Mesh mesh(*TestCommWorld);
    MeshTools::Generation::build_point(mesh);
    EquationSystems es(mesh);
    SystemType & system = es.add_system<SystemType>("ScalarSystem");
 
    system.time_solver = std::make_unique<TimeSolverType>(system);
 
    es.init();
 
    DiffSolver & solver = *(system.time_solver->diff_solver().get());
    solver.relative_step_tolerance = std::numeric_limits<Real>::epsilon()*10;
    solver.relative_residual_tolerance = std::numeric_limits<Real>::epsilon()*10;
    solver.absolute_residual_tolerance = std::numeric_limits<Real>::epsilon()*10;
 
    NewtonSolver & newton = cast_ref<NewtonSolver &>(solver);
 
    // LASPACK GMRES + ILU defaults don't like these problems, so
    // we'll use a sophisticated "just divide the scalars" solver instead.
    newton.get_linear_solver().set_solver_type(JACOBI);
    newton.get_linear_solver().set_preconditioner_type(IDENTITY_PRECOND);
 
    system.deltat = deltat;
 
    TimeSolverType * time_solver = cast_ptr<TimeSolverType *>(system.time_solver.get());
    this->aux_time_solver_init(*time_solver);
 
    // We're going to want to check our solution, and when we run
    // "make check" with LIBMESH_RUN='mpirun -np N" for N>1 then we'll
    // need to keep that check in sync with the processors that are just
    // twiddling their thumbs, not owning our mesh point.
    std::vector<dof_id_type> solution_index;
    solution_index.push_back(0);
    const bool has_solution = system.get_dof_map().all_semilocal_indices(solution_index);
 
    for (unsigned int t_step=0; t_step != n_timesteps; ++t_step)
      {
        system.solve();
        system.time_solver->advance_timestep();
 
        Real rel_error = 0;
 
        if (has_solution)
          {
            Number exact_soln = system.u(system.time);
            rel_error =  std::abs((exact_soln - (*system.solution)(0))/exact_soln);
          }
        system.comm().max(rel_error);
 
        // Using relative error for comparison, so "exact" is 0
        LIBMESH_ASSERT_FP_EQUAL( rel_error,
                                 0.0,
                                 std::numeric_limits<Real>::epsilon()*10 );
      }
  }
 
};
 
 
class FirstOrderScalarSystemBase : public FEMSystem
{
public:
  FirstOrderScalarSystemBase(EquationSystems & es,
                             const std::string & name_in,
                             const unsigned int number_in)
    : FEMSystem(es, name_in, number_in)
  {}
 
 
  virtual Number F( FEMContext & context, unsigned int qp ) =0;
 
  virtual Number M( FEMContext & context, unsigned int qp ) =0;
 
  virtual Number u( Real t ) =0;
 
  virtual void init_data () override
  {
    _u_var = this->add_variable ("u", FIRST, LAGRANGE);
    this->time_evolving(_u_var,1);
    FEMSystem::init_data();
  }
 
  virtual bool element_time_derivative (bool request_jacobian,
                                        DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
    DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Number Fval = this->F(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          Fu(i) += Fval;
      }
 
    return request_jacobian;
  }
 
  virtual bool mass_residual (bool request_jacobian,
                              DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
    DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
    DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(_u_var, _u_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Number udot;
        c.interior_rate(_u_var,qp,udot);
 
        Number Mval = this->M(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fu(i) -= Mval*udot;
 
            if (request_jacobian)
              for (unsigned int j=0; j != n_u_dofs; j++)
                Kuu(i,j) -= Mval*context.get_elem_solution_rate_derivative();
          }
      }
 
    return request_jacobian;
  }
 
protected:
  unsigned int _u_var;
};
 
 
class SecondOrderScalarSystemSecondOrderTimeSolverBase : public FirstOrderScalarSystemBase
{
public:
  SecondOrderScalarSystemSecondOrderTimeSolverBase(EquationSystems & es,
                                                   const std::string & name_in,
                                                   const unsigned int number_in)
    : FirstOrderScalarSystemBase(es, name_in, number_in)
  {}
 
  virtual void init_data () override
  {
    _u_var = this->add_variable ("u", FIRST, LAGRANGE);
    this->time_evolving(_u_var,2);
    FEMSystem::init_data();
  }
 
  virtual Number C( FEMContext & context, unsigned int qp ) =0;
 
  virtual bool damping_residual (bool request_jacobian,
                                 DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
    DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
    DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(_u_var, _u_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Number udot;
        c.interior_rate(_u_var,qp,udot);
 
        Number Cval = this->C(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fu(i) += Cval*udot;
 
            if (request_jacobian)
              for (unsigned int j=0; j != n_u_dofs; j++)
                Kuu(i,j) += Cval*context.get_elem_solution_rate_derivative();
          }
      }
 
    return request_jacobian;
  }
 
  virtual bool mass_residual (bool request_jacobian,
                              DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
    DenseSubVector<Number> & Fu = c.get_elem_residual(_u_var);
    DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(_u_var, _u_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Number uddot;
        c.interior_accel(_u_var,qp,uddot);
 
        Number Mval = this->M(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fu(i) += Mval*uddot;
 
            if (request_jacobian)
              for (unsigned int j=0; j != n_u_dofs; j++)
                Kuu(i,j) += Mval*context.get_elem_solution_accel_derivative();
          }
      }
 
    return request_jacobian;
  }
};
 
 
 
class SecondOrderScalarSystemFirstOrderTimeSolverBase : public SecondOrderScalarSystemSecondOrderTimeSolverBase
{
public:
  SecondOrderScalarSystemFirstOrderTimeSolverBase(EquationSystems & es,
                                                  const std::string & name_in,
                                                  const unsigned int number_in)
    : SecondOrderScalarSystemSecondOrderTimeSolverBase(es, name_in, number_in)
  {}
 
  virtual bool element_time_derivative (bool request_jacobian,
                                        DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
 
    unsigned int v_var = this->get_second_order_dot_var(_u_var);
 
    DenseSubVector<Number> & Fv = c.get_elem_residual(v_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        Number Fval = this->F(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          Fv(i) += Fval;
      }
 
    return request_jacobian;
  }
 
  virtual bool damping_residual (bool request_jacobian,
                                 DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
 
    unsigned int v_var = this->get_second_order_dot_var(_u_var);
 
    DenseSubVector<Number> & Fv = c.get_elem_residual(v_var);
 
    DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Number udot;
        c.interior_rate(v_var,qp,udot);
 
        Number Cval = this->C(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fv(i) += Cval*udot;
 
            if (request_jacobian)
              for (unsigned int j=0; j != n_u_dofs; j++)
                Kvv(i,j) += Cval*context.get_elem_solution_rate_derivative();
          }
      }
 
    return request_jacobian;
  }
 
  virtual bool mass_residual (bool request_jacobian,
                              DiffContext & context) override
  {
    FEMContext & c = cast_ref<FEMContext &>(context);
 
    unsigned int v_var = this->get_second_order_dot_var(_u_var);
 
    DenseSubVector<Number> & Fv = c.get_elem_residual(v_var);
    DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
 
    const unsigned int n_u_dofs =
      cast_int<unsigned int>(c.get_dof_indices(_u_var).size());
    unsigned int n_qpoints = c.get_element_qrule().n_points();
 
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        libMesh::Number uddot;
        c.interior_accel(v_var,qp,uddot);
 
        Number Mval = this->M(c,qp);
 
        for (unsigned int i=0; i != n_u_dofs; i++)
          {
            Fv(i) += Mval*uddot;
 
            if (request_jacobian)
              for (unsigned int j=0; j != n_u_dofs; j++)
                Kvv(i,j) += Mval*context.get_elem_solution_accel_derivative();
          }
      }
 
    return request_jacobian;
  }
};