fem_system_ex6.C File Reference

Go to the source code of this file.

Classes
class	NonlinearDiffusionSystem
	An FEMSystem for -div((1 + u^2) grad u) = f whose element Jacobian is computed by forward-mode automatic differentiation of the element residual. More...

Typedefs
typedef MetaPhysicL::DualNumber< Number, MetaPhysicL::DynamicSparseNumberArray< Number, unsigned int > >	ADNumber

Functions
Number	exact_value (const Point &p, const Parameters &, const std::string &, const std::string &)
Gradient	exact_gradient (const Point &p, const Parameters &, const std::string &, const std::string &)
Number	forcing (const Point &p)
int	main (int argc, char **argv)

Typedef Documentation

ADNumber

typedef MetaPhysicL::DualNumber<Number, MetaPhysicL::DynamicSparseNumberArray<Number, unsigned int> > ADNumber

Definition at line 81 of file fem_system_ex6.C.

Function Documentation

exact_gradient()

Gradient exact_gradient	(	const Point &	p,
		const Parameters &	,
		const std::string &	,
		const std::string &
	)

Definition at line 91 of file fem_system_ex6.C.

{
  const Real pi = libMesh::pi;
  return RealGradient(pi * std::cos(pi * p(0)) * std::sin(pi * p(1)),
                      pi * std::sin(pi * p(0)) * std::cos(pi * p(1)));
}

References libMesh::pi, and libMesh::Real.

Referenced by forcing(), and main().

exact_value()

Number exact_value	(	const Point &	p,
		const Parameters &	,
		const std::string &	,
		const std::string &
	)

Definition at line 85 of file fem_system_ex6.C.

{
  return std::sin(libMesh::pi * p(0)) * std::sin(libMesh::pi * p(1));
}

References libMesh::pi.

Referenced by forcing(), and main().

forcing()

Number forcing

(

const Point &

p

)

Definition at line 101 of file fem_system_ex6.C.

{
  const Number u = exact_value(p, Parameters(), "", "");
  const Gradient grad_u = exact_gradient(p, Parameters(), "", "");
  return 2 * libMesh::pi * libMesh::pi * (1. + u * u) * u - 2. * u * grad_u.norm_sq();
}

References exact_gradient(), exact_value(), libMesh::TypeVector< T >::norm_sq(), and libMesh::pi.

Referenced by assemble_advection_diffusion(), HeatSystem::element_time_derivative(), and NonlinearDiffusionSystem::element_time_derivative().

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 236 of file fem_system_ex6.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");
 
  // Dual numbers are provided by MetaPhysicL.
#ifndef LIBMESH_HAVE_METAPHYSICL
  libmesh_example_requires(false, "--enable-metaphysicl");
#else
 
  // The mesh is two-dimensional.
  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
 
  // We use Dirichlet boundary conditions here.
#ifndef LIBMESH_ENABLE_DIRICHLET
  libmesh_example_requires(false, "--enable-dirichlet");
#endif
 
  // Parse the input file, allowing the command line to override it.
  GetPot infile("fem_system_ex6.in");
  infile.parse_command_line(argc, argv);
 
  const unsigned int grid_size = infile("grid_size", 16);
 
  Mesh mesh(init.comm());
  MeshTools::Generation::build_square(mesh, grid_size, grid_size, 0., 1., 0., 1., QUAD9);
  mesh.print_info();
 
  EquationSystems equation_systems(mesh);
  auto & system = equation_systems.add_system<NonlinearDiffusionSystem>("NonlinearDiffusion");
  system.time_solver = std::make_unique<SteadySolver>(system);
 
  // A positive tolerance checks each element Jacobian against a
  // finite-difference approximation.
  system.verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
 
  equation_systems.init();
 
  DiffSolver & solver = *(system.time_solver->diff_solver());
  solver.quiet = false;
  solver.verbose = true;
  solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 10);
  solver.relative_residual_tolerance = infile("relative_residual_tolerance", 1.e-12);
  solver.relative_step_tolerance = infile("relative_step_tolerance", 0.);
  solver.absolute_residual_tolerance = infile("absolute_residual_tolerance", 0.);
 
  equation_systems.print_info();
 
  system.solve();
 
  ExactSolution exact_sol(equation_systems);
  exact_sol.attach_exact_value(exact_value);
  exact_sol.attach_exact_deriv(exact_gradient);
  exact_sol.compute_error("NonlinearDiffusion", "u");
  libMesh::out << "L2 error: " << exact_sol.l2_error("NonlinearDiffusion", "u") << std::endl
               << "H1 error: " << exact_sol.h1_error("NonlinearDiffusion", "u") << std::endl;
 
#endif // LIBMESH_HAVE_METAPHYSICL
 
  return 0;
}

References libMesh::DiffSolver::absolute_residual_tolerance, libMesh::EquationSystems::add_system(), libMesh::ExactSolution::attach_exact_deriv(), libMesh::ExactSolution::attach_exact_value(), libMesh::MeshTools::Generation::build_square(), libMesh::ExactSolution::compute_error(), libMesh::default_solver_package(), exact_gradient(), exact_value(), libMesh::ExactSolution::h1_error(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), main(), libMesh::DiffSolver::max_nonlinear_iterations, mesh, libMesh::out, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::DiffSolver::quiet, libMesh::DiffSolver::relative_residual_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::DifferentiableSystem::time_solver, and libMesh::DiffSolver::verbose.