finite_volume_ex1.C File Reference

Go to the source code of this file.

Functions
void	assemble_advection_diffusion (EquationSystems &es, const std::string &system_name)
Real	centroid_l2_error (const System &system)
int	main (int argc, char **argv)

Function Documentation

assemble_advection_diffusion()

void assemble_advection_diffusion	(	EquationSystems &	es,
		const std::string &	system_name
	)

Definition at line 231 of file finite_volume_ex1.C.

{
  const MeshBase & mesh = es.get_mesh();
  auto & system = es.get_system<LinearImplicitSystem>(system_name);
  const unsigned int sys_num = system.number();
 
  const Real diffusivity = es.parameters.get<Real>("diffusivity");
  const RealVectorValue velocity = es.parameters.get<RealVectorValue>("velocity");
 
  for (const Elem * elem : mesh.active_local_element_ptr_range())
    {
      const dof_id_type elem_dof = elem->dof_number(sys_num, 0, 0);
      const Point x_elem = elem->true_centroid();
      const GradientStencil grad_elem = gradient_stencil(*elem, sys_num);
 
      // The net outward flux through the cell boundary
      ScalarStencil net_flux;
 
      for (const auto s : elem->side_index_range())
        {
          const FaceGeometry face = face_geometry(*elem, s);
          const Elem * neighbor = elem->neighbor_ptr(s);
 
          // On interior faces, the neighbor centroid and gradient; on
          // boundary faces, the face centroid, where the Dirichlet
          // value is known, and the gradient of this cell.
          const Point x_across = neighbor ? neighbor->true_centroid() : face.centroid;
          const GradientStencil grad_neighbor =
            neighbor ? gradient_stencil(*neighbor, sys_num) : GradientStencil();
          const GradientStencil grad_face = neighbor ? average(grad_elem, grad_neighbor) : grad_elem;
 
          // Diffusion: with d the offset from x_elem to x_across,
          //   n.grad(u) ~= (u_across - u_elem) / (d.n) + (n - d / (d.n)).grad(u)_face,
          // which is exact for linear u. The second term vanishes when d
          // is parallel to n.
          const RealVectorValue d = x_across - x_elem;
          const Real d_dot_n = d * face.normal;
          const Real diffusive_scale = -diffusivity * face.area;
 
          net_flux.coefficients[elem_dof] -= diffusive_scale / d_dot_n;
          if (neighbor)
            net_flux.coefficients[neighbor->dof_number(sys_num, 0, 0)] += diffusive_scale / d_dot_n;
          else
            net_flux.constant += diffusive_scale / d_dot_n * exact_solution(face.centroid);
          add_dot(net_flux, grad_face, face.normal - d / d_dot_n, diffusive_scale);
 
          // Advection: the upwind cell value, extrapolated to the face
          // centroid with that cell's gradient. Inflow boundary faces
          // take the Dirichlet value.
          const Real mass_flux = (velocity * face.normal) * face.area;
          if (mass_flux >= 0)
            {
              net_flux.coefficients[elem_dof] += mass_flux;
              add_dot(net_flux, grad_elem, face.centroid - x_elem, mass_flux);
            }
          else if (neighbor)
            {
              net_flux.coefficients[neighbor->dof_number(sys_num, 0, 0)] += mass_flux;
              add_dot(net_flux, grad_neighbor, face.centroid - x_across, mass_flux);
            }
          else
            net_flux.constant += mass_flux * exact_solution(face.centroid);
        }
 
      std::vector<dof_id_type> row_dofs = {elem_dof};
      std::vector<dof_id_type> column_dofs;
      DenseMatrix<Number> Ke(1, net_flux.coefficients.size());
      for (const auto & [dof, coef] : net_flux.coefficients)
        {
          Ke(0, column_dofs.size()) = coef;
          column_dofs.push_back(dof);
        }
 
      DenseVector<Number> Fe(1);
      Fe(0) = forcing(x_elem, diffusivity, velocity) * elem->volume() - net_flux.constant;
 
      system.matrix->add_matrix(Ke, row_dofs, column_dofs);
      system.rhs->add_vector(Fe, row_dofs);
    }
}

References libMesh::DofObject::dof_number(), exact_solution, forcing(), libMesh::Parameters::get(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::Elem::neighbor_ptr(), libMesh::System::number(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::Elem::side_index_range(), libMesh::Elem::true_centroid(), and libMesh::Elem::volume().

Referenced by main().

centroid_l2_error()

Real centroid_l2_error

(

const System &

system

)

Definition at line 315 of file finite_volume_ex1.C.

{
  const unsigned int sys_num = system.number();
 
  Real error_sq = 0;
  for (const Elem * elem : system.get_mesh().active_local_element_ptr_range())
    {
      const Real u_h = libmesh_real(system.current_solution(elem->dof_number(sys_num, 0, 0)));
      const Real diff = u_h - exact_solution(elem->true_centroid());
      error_sq += elem->volume() * diff * diff;
    }
  system.comm().sum(error_sq);
 
  return std::sqrt(error_sq);
}

References libMesh::ParallelObject::comm(), libMesh::System::current_solution(), libMesh::DofObject::dof_number(), exact_solution, libMesh::System::get_mesh(), libMesh::libmesh_real(), libMesh::System::number(), libMesh::Real, libMesh::Parallel::Communicator::sum(), libMesh::Elem::true_centroid(), and libMesh::Elem::volume().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 332 of file finite_volume_ex1.C.

{
  // Initialize libMesh.
  LibMeshInit init(argc, argv);
 
  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");
 
  // The mesh is two-dimensional.
  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
 
  // Parse the input file, allowing the command line to override it.
  GetPot infile("finite_volume_ex1.in");
  infile.parse_command_line(argc, argv);
 
  const unsigned int coarse_grid_size = infile("coarse_grid_size", 4);
  const unsigned int max_refinements = infile("max_refinements", 6);
  const Real diffusivity = infile("diffusivity", 1.);
  const RealVectorValue velocity(infile("velocity_x", 1.), infile("velocity_y", 0.5));
 
  // The scheme is second order, and a measured rate within this
  // distance of 2 is taken as confirmation of that order. The bound is
  // loose enough to absorb the pre-asymptotic behavior of the error on
  // the coarsest grids and the round-off of any supported Real
  // precision, and tight enough to catch a loss of second-order
  // accuracy. Dropping the advection reconstruction reduces the scheme
  // to first order, and dropping the non-orthogonal diffusion
  // correction stops the error from converging at all.
  const Real expected_rate = 2;
  const Real rate_tolerance = 0.1;
 
  std::vector<Real> errors;
  bool converged = false;
  for (const auto level : make_range(max_refinements + 1))
    {
      const unsigned int grid_size = coarse_grid_size << level;
 
      Mesh mesh(init.comm());
      EquationSystems equation_systems(mesh);
      auto & system = equation_systems.add_system<LinearImplicitSystem>("AdvectionDiffusion");
      system.add_variable("u", CONSTANT, MONOMIAL);
      system.attach_assemble_function(assemble_advection_diffusion);
 
      // Couple each cell to two layers of side neighbors. The DofMap
      // registers its coupling functor with the mesh, so building the
      // mesh after this call also retains the corresponding ghost
      // elements on a distributed mesh.
      system.get_dof_map().default_coupling().set_n_levels(2);
 
      MeshTools::Generation::build_square(mesh, grid_size, grid_size, 0., 1., 0., 1., C0POLYGON);
 
      equation_systems.parameters.set<Real>("diffusivity") = diffusivity;
      equation_systems.parameters.set<RealVectorValue>("velocity") = velocity;
      equation_systems.init();
 
      system.solve();
 
      const Real error = centroid_l2_error(system);
      errors.push_back(error);
      libMesh::out << "Grid size " << grid_size << ": " << mesh.n_active_elem()
                   << " cells, centroid L2 error " << error;
      if (level)
        libMesh::out << ", convergence rate " << std::log2(errors[level - 1] / error);
 
      // Fit a line to log2(error) against log2(h) over the last three
      // grids. The grid spacing h halves at each level, so the sample
      // points are equally spaced and the least-squares slope reduces
      // to the difference between the first and last points divided by
      // their separation of two levels.
      if (level >= 2)
        {
          const Real fit_rate = std::log2(errors[level - 2] / error) / 2;
          libMesh::out << ", three-grid fit rate " << fit_rate;
          converged = std::abs(fit_rate - expected_rate) < rate_tolerance;
        }
      libMesh::out << std::endl;
 
      if (converged)
        {
#ifdef LIBMESH_HAVE_EXODUS_API
          ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
#endif
          break;
        }
    }
 
  libmesh_error_msg_if(!converged,
                       "The three-grid fit of the L2 error convergence rate did not reach "
                           << expected_rate << " +/- " << rate_tolerance << " within "
                           << max_refinements << " refinements");
 
  return 0;
}

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_advection_diffusion(), libMesh::MeshTools::Generation::build_square(), libMesh::C0POLYGON, centroid_l2_error(), libMesh::CONSTANT, libMesh::default_solver_package(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, main(), libMesh::make_range(), mesh, libMesh::MONOMIAL, libMesh::MeshBase::n_active_elem(), libMesh::out, libMesh::EquationSystems::parameters, libMesh::Real, libMesh::Parameters::set(), and libMesh::ExodusII_IO::write_equation_systems().