// The libMesh Finite Element Library.
// Copyright (C) 2002-2016 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
FEMSystem Example 4 - Mixed-dimension heat transfer equation
// with FEMSystem
// \author Roy Stogner
// \date 2015
//
// This example shows how elements of multiple dimensions can be
// linked and computed upon simultaneously within the
// DifferentiableSystem class framework
// C++ includes
#include
// Basic include files
#include "libmesh/elem.h"
#include "libmesh/equation_systems.h"
#include "libmesh/error_vector.h"
#include "libmesh/exact_solution.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/parsed_function.h"
#include "libmesh/uniform_refinement_estimator.h"
// The systems and solvers we may use
#include "heatsystem.hpp"
#include "libmesh/diff_solver.h"
#include "libmesh/euler_solver.h"
#include "libmesh/steady_solver.h"
// user defined literal for libmesh's Real
constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }
// The main program.
int main(int argc, char** argv) {
// Initialize libMesh.
libMesh::LibMeshInit init{argc, argv};
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_requires(false, "--enable-amr");
#else
// This doesn't converge with Eigen BICGSTAB for some reason...
libmesh_example_requires(libMesh::default_solver_package() !=
libMesh::EIGEN_SOLVERS,
"--enable-petsc");
// This doesn't converge without at least double precision
libmesh_example_requires(sizeof(libMesh::Real) > 4,
"--disable-singleprecision");
// Parse the input file
GetPot infile{"../fem_system_ex4.in"};
// Read in parameters from the input file
const auto& global_tolerance{infile("global_tolerance", 0._r)};
const auto& nelem_target{infile("n_elements", 400u)};
const auto& deltat{infile("deltat", 0.005_r)};
const auto& coarsegridsize{infile("coarsegridsize", 20u)};
const auto& coarserefinements{infile("coarserefinements", 0u)};
const auto& max_adaptivesteps{infile("max_adaptivesteps", 10u)};
const auto& dim{infile("dimension", 2u)};
// Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
// We have only defined 2 and 3 dimensional problems
libmesh_assert(dim == 2 || dim == 3);
// Create a mesh, with dimension to be overridden later, distributed
// across the default MPI communicator.
libMesh::Mesh mesh{init.comm()};
// And an object to refine it
libMesh::MeshRefinement mesh_refinement{mesh};
mesh_refinement.coarsen_by_parents() = true;
mesh_refinement.absolute_global_tolerance() = global_tolerance;
mesh_refinement.nelem_target() = nelem_target;
mesh_refinement.refine_fraction() = 0.3;
mesh_refinement.coarsen_fraction() = 0.3;
mesh_refinement.coarsen_threshold() = 0.1;
// Use the MeshTools::Generation mesh generator to create a uniform
// grid on the square or cube. We crop the domain at y=2/3 to allow
// for a homogeneous Neumann BC in our benchmark there.
libMesh::boundary_id_type bcid{3}; // +y in 3D
mesh.read("../meshes/bridge.xda");
// Add boundary elements corresponding to the +y boundary of our
// volume mesh
mesh.get_boundary_info().add_elements({bcid}, mesh);
mesh.prepare_for_use();
// To work around ExodusII file format limitations, we need elements
// of different dimensionality to belong to different subdomains.
// Our interior elements defaulted to subdomain id 0, so we'll set
// boundary elements to subdomain 2.
std::for_each(mesh.elements_begin(), mesh.elements_end(), [&](auto elem) {
if (elem->dim() < dim) elem->subdomain_id() = 2;
});
mesh_refinement.uniformly_refine(coarserefinements);
// Print information about the mesh to the screen.
mesh.print_info();
// Create an equation systems object.
libMesh::EquationSystems equation_systems{mesh};
// Declare the system "Heat" and its variables.
auto& system = equation_systems.add_system("Heat");
// Solve this as a steady system
system.time_solver = libMesh::UniquePtr(
new libMesh::SteadySolver(system));
// Initialize the system
equation_systems.init();
// Set the time stepping options
system.deltat = deltat;
// And the nonlinear solver options
auto& solver = *(system.time_solver->diff_solver().get());
solver.quiet = infile("solver_quiet", true);
solver.verbose = !solver.quiet;
solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15u);
solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3_r);
solver.relative_residual_tolerance =
infile("relative_residual_tolerance", 0.0_r);
solver.absolute_residual_tolerance =
infile("absolute_residual_tolerance", 0.0_r);
// And the linear solver options
solver.max_linear_iterations = infile("max_linear_iterations", 50000u);
solver.initial_linear_tolerance =
infile("initial_linear_tolerance", 1.e-3_r);
// Print information about the system to the screen.
equation_systems.print_info();
// Adaptively solve the steady solution
auto a_step = max_adaptivesteps;
for (; a_step != max_adaptivesteps; ++a_step) {
system.solve();
system.postprocess();
libMesh::UniquePtr error_estimator;
// To solve to a tolerance in this problem we
// need a better estimator than Kelly
if (global_tolerance != 0.) {
// We can't adapt to both a tolerance and a mesh
// size at once
libmesh_assert_equal_to(nelem_target, 0);
libMesh::UniformRefinementEstimator* u{
new libMesh::UniformRefinementEstimator};
// The lid-driven cavity problem isn't in H1, so
// lets estimate L2 error
u->error_norm = libMesh::L2;
error_estimator.reset(u);
} else {
// If we aren't adapting to a tolerance we need a
// target mesh size
libmesh_assert_greater(nelem_target, 0);
// Kelly is a lousy estimator to use for a problem
// not in H1 - if we were doing more than a few
// timesteps we'd need to turn off or limit the
// maximum level of our adaptivity eventually
error_estimator.reset(new libMesh::KellyErrorEstimator);
}
libMesh::ErrorVector error;
error_estimator->estimate_error(system, error);
// Print out status at each adaptive step.
libMesh::out << "Adaptive step " << a_step << ": " << '\n';
const auto& global_error = error.l2_norm();
if (global_tolerance != 0.)
libMesh::out << "Global_error = " << global_error << '\n';
if (global_tolerance != 0.)
libMesh::out << "Worst element error = " << error.maximum()
<< ", mean = " << error.mean() << '\n';
if (global_tolerance != 0.) {
// If we've reached our desired tolerance, we
// don't need any more adaptive steps
if (global_error < global_tolerance) break;
mesh_refinement.flag_elements_by_error_tolerance(error);
} else {
// If flag_elements_by_nelem_target returns true, this
// should be our last adaptive step.
if (mesh_refinement.flag_elements_by_nelem_target(error)) {
mesh_refinement.refine_and_coarsen_elements();
equation_systems.reinit();
a_step = max_adaptivesteps;
break;
}
}
// Carry out the adaptive mesh refinement/coarsening
mesh_refinement.refine_and_coarsen_elements();
equation_systems.reinit();
libMesh::out << "Refined mesh to " << mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs() << " active dofs.\n";
}
// Do one last solve if necessary
if (a_step == max_adaptivesteps) {
system.solve();
system.postprocess();
}
#ifdef LIBMESH_HAVE_EXODUS_API
libMesh::ExodusII_IO(mesh).write_equation_systems("out.e",
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
#ifdef LIBMESH_HAVE_GMV
libMesh::GMVIO(mesh).write_equation_systems("out.gmv", equation_systems);
#endif // #ifdef LIBMESH_HAVE_GMV
#ifdef LIBMESH_HAVE_FPARSER
// Check that we got close to the analytic solution
libMesh::ExactSolution exact_sol(equation_systems);
const std::string exact_str =
(dim == 2) ? "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
libMesh::ParsedFunction exact_func(exact_str);
exact_sol.attach_exact_value(0, &exact_func);
exact_sol.compute_error("Heat", "T");
libMesh::Number err{exact_sol.l2_error("Heat", "T")};
// Print out the error value
libMesh::out << "L2-Error is: " << err << '\n';
libmesh_assert_less(libmesh_real(err), 2e-3);
#endif // #ifdef LIBMESH_HAVE_FPARSER
#endif // #ifndef LIBMESH_ENABLE_AMR
// All done.
return 0;
}