Commit 0adbd844 authored by Thomas Steinreiter's avatar Thomas Steinreiter
Browse files

refactoring

parent 0ebca407
...@@ -28,9 +28,9 @@ if (MPI_FOUND) ...@@ -28,9 +28,9 @@ if (MPI_FOUND)
if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel") elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xHost -std=c++11") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xHost -std=c++14")
endif() endif()
set_target_properties(${NAME} PROPERTIES CXX_STANDARD 11 CXX_STANDARD_REQUIRED YES) set_target_properties(${NAME} PROPERTIES CXX_STANDARD 14 CXX_STANDARD_REQUIRED YES)
target_link_libraries(${NAME} ${LIBMESH_LIBRARY} ${MPI_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT}) target_link_libraries(${NAME} ${LIBMESH_LIBRARY} ${MPI_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
install(TARGETS ${NAME} DESTINATION bin) install(TARGETS ${NAME} DESTINATION bin)
message("** Enabling '${NAME}': with MPI") message("** Enabling '${NAME}': with MPI")
......
...@@ -48,36 +48,38 @@ ...@@ -48,36 +48,38 @@
#include "libmesh/euler_solver.h" #include "libmesh/euler_solver.h"
#include "libmesh/steady_solver.h" #include "libmesh/steady_solver.h"
// Bring in everything from the libMesh namespace // user defined literal for libmesh's Real
using namespace libMesh; constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }
// The main program. // The main program.
int main(int argc, char** argv) { int main(int argc, char** argv) {
// Initialize libMesh. // Initialize libMesh.
LibMeshInit init(argc, argv); libMesh::LibMeshInit init{argc, argv};
#ifndef LIBMESH_ENABLE_AMR #ifndef LIBMESH_ENABLE_AMR
libmesh_example_requires(false, "--enable-amr"); libmesh_example_requires(false, "--enable-amr");
#else #else
// This doesn't converge with Eigen BICGSTAB for some reason... // This doesn't converge with Eigen BICGSTAB for some reason...
libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, libmesh_example_requires(libMesh::default_solver_package() !=
libMesh::EIGEN_SOLVERS,
"--enable-petsc"); "--enable-petsc");
// This doesn't converge without at least double precision // This doesn't converge without at least double precision
libmesh_example_requires(sizeof(Real) > 4, "--disable-singleprecision"); libmesh_example_requires(sizeof(libMesh::Real) > 4,
"--disable-singleprecision");
// Parse the input file // Parse the input file
GetPot infile("../fem_system_ex4.in"); GetPot infile{"../fem_system_ex4.in"};
// Read in parameters from the input file // Read in parameters from the input file
const Real global_tolerance = infile("global_tolerance", 0.); const auto& global_tolerance{infile("global_tolerance", 0._r)};
const unsigned int nelem_target = infile("n_elements", 400); const auto& nelem_target{infile("n_elements", 400u)};
const Real deltat = infile("deltat", 0.005); const auto& deltat{infile("deltat", 0.005_r)};
const unsigned int coarsegridsize = infile("coarsegridsize", 20); const auto& coarsegridsize{infile("coarsegridsize", 20u)};
const unsigned int coarserefinements = infile("coarserefinements", 0); const auto& coarserefinements{infile("coarserefinements", 0u)};
const unsigned int max_adaptivesteps = infile("max_adaptivesteps", 10); const auto& max_adaptivesteps{infile("max_adaptivesteps", 10u)};
const unsigned int dim = infile("dimension", 2); const auto& dim{infile("dimension", 2u)};
// Skip higher-dimensional examples on a lower-dimensional libMesh build // Skip higher-dimensional examples on a lower-dimensional libMesh build
libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support"); libmesh_example_requires(dim <= LIBMESH_DIM, "2D/3D support");
...@@ -87,10 +89,10 @@ int main(int argc, char** argv) { ...@@ -87,10 +89,10 @@ int main(int argc, char** argv) {
// Create a mesh, with dimension to be overridden later, distributed // Create a mesh, with dimension to be overridden later, distributed
// across the default MPI communicator. // across the default MPI communicator.
Mesh mesh(init.comm()); libMesh::Mesh mesh{init.comm()};
// And an object to refine it // And an object to refine it
MeshRefinement mesh_refinement(mesh); libMesh::MeshRefinement mesh_refinement{mesh};
mesh_refinement.coarsen_by_parents() = true; mesh_refinement.coarsen_by_parents() = true;
mesh_refinement.absolute_global_tolerance() = global_tolerance; mesh_refinement.absolute_global_tolerance() = global_tolerance;
mesh_refinement.nelem_target() = nelem_target; mesh_refinement.nelem_target() = nelem_target;
...@@ -101,30 +103,22 @@ int main(int argc, char** argv) { ...@@ -101,30 +103,22 @@ int main(int argc, char** argv) {
// Use the MeshTools::Generation mesh generator to create a uniform // 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 // 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. // for a homogeneous Neumann BC in our benchmark there.
boundary_id_type bcid = 3; // +y in 3D libMesh::boundary_id_type bcid{3}; // +y in 3D
mesh.read("../meshes/bridge.xda"); mesh.read("../meshes/bridge.xda");
{
// Add boundary elements corresponding to the +y boundary of our // Add boundary elements corresponding to the +y boundary of our
// volume mesh // volume mesh
std::set<boundary_id_type> bcids; mesh.get_boundary_info().add_elements({bcid}, mesh);
bcids.insert(bcid);
mesh.get_boundary_info().add_elements(bcids, mesh);
mesh.prepare_for_use(); mesh.prepare_for_use();
}
// To work around ExodusII file format limitations, we need elements // To work around ExodusII file format limitations, we need elements
// of different dimensionality to belong to different subdomains. // of different dimensionality to belong to different subdomains.
// Our interior elements defaulted to subdomain id 0, so we'll set // Our interior elements defaulted to subdomain id 0, so we'll set
// boundary elements to subdomain 2. // boundary elements to subdomain 2.
{ std::for_each(mesh.elements_begin(), mesh.elements_end(), [&](auto elem) {
const MeshBase::element_iterator end_el = mesh.elements_end();
for (MeshBase::element_iterator el = mesh.elements_begin();
el != end_el; ++el) {
Elem* elem = *el;
if (elem->dim() < dim) elem->subdomain_id() = 2; if (elem->dim() < dim) elem->subdomain_id() = 2;
} });
}
mesh_refinement.uniformly_refine(coarserefinements); mesh_refinement.uniformly_refine(coarserefinements);
...@@ -132,13 +126,14 @@ int main(int argc, char** argv) { ...@@ -132,13 +126,14 @@ int main(int argc, char** argv) {
mesh.print_info(); mesh.print_info();
// Create an equation systems object. // Create an equation systems object.
EquationSystems equation_systems(mesh); libMesh::EquationSystems equation_systems{mesh};
// Declare the system "Heat" and its variables. // Declare the system "Heat" and its variables.
HeatSystem& system = equation_systems.add_system<HeatSystem>("Heat"); auto& system = equation_systems.add_system<HeatSystem>("Heat");
// Solve this as a steady system // Solve this as a steady system
system.time_solver = UniquePtr<TimeSolver>(new SteadySolver(system)); system.time_solver = libMesh::UniquePtr<libMesh::TimeSolver>(
new libMesh::SteadySolver(system));
// Initialize the system // Initialize the system
equation_systems.init(); equation_systems.init();
...@@ -147,34 +142,32 @@ int main(int argc, char** argv) { ...@@ -147,34 +142,32 @@ int main(int argc, char** argv) {
system.deltat = deltat; system.deltat = deltat;
// And the nonlinear solver options // And the nonlinear solver options
DiffSolver& solver = *(system.time_solver->diff_solver().get()); auto& solver = *(system.time_solver->diff_solver().get());
solver.quiet = infile("solver_quiet", true); solver.quiet = infile("solver_quiet", true);
solver.verbose = !solver.quiet; solver.verbose = !solver.quiet;
solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15); solver.max_nonlinear_iterations = infile("max_nonlinear_iterations", 15u);
solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3); solver.relative_step_tolerance = infile("relative_step_tolerance", 1.e-3_r);
solver.relative_residual_tolerance = solver.relative_residual_tolerance =
infile("relative_residual_tolerance", 0.0); infile("relative_residual_tolerance", 0.0_r);
solver.absolute_residual_tolerance = solver.absolute_residual_tolerance =
infile("absolute_residual_tolerance", 0.0); infile("absolute_residual_tolerance", 0.0_r);
// And the linear solver options // And the linear solver options
solver.max_linear_iterations = infile("max_linear_iterations", 50000); solver.max_linear_iterations = infile("max_linear_iterations", 50000u);
solver.initial_linear_tolerance = infile("initial_linear_tolerance", 1.e-3); solver.initial_linear_tolerance =
infile("initial_linear_tolerance", 1.e-3_r);
// Print information about the system to the screen. // Print information about the system to the screen.
equation_systems.print_info(); equation_systems.print_info();
// Adaptively solve the steady solution // Adaptively solve the steady solution
unsigned int a_step = max_adaptivesteps; auto a_step = max_adaptivesteps;
for (; a_step != max_adaptivesteps; ++a_step) { for (; a_step != max_adaptivesteps; ++a_step) {
system.solve(); system.solve();
system.postprocess(); system.postprocess();
ErrorVector error; libMesh::UniquePtr<libMesh::ErrorEstimator> error_estimator;
UniquePtr<ErrorEstimator> error_estimator;
// To solve to a tolerance in this problem we // To solve to a tolerance in this problem we
// need a better estimator than Kelly // need a better estimator than Kelly
if (global_tolerance != 0.) { if (global_tolerance != 0.) {
...@@ -182,11 +175,12 @@ int main(int argc, char** argv) { ...@@ -182,11 +175,12 @@ int main(int argc, char** argv) {
// size at once // size at once
libmesh_assert_equal_to(nelem_target, 0); libmesh_assert_equal_to(nelem_target, 0);
UniformRefinementEstimator* u = new UniformRefinementEstimator; libMesh::UniformRefinementEstimator* u{
new libMesh::UniformRefinementEstimator};
// The lid-driven cavity problem isn't in H1, so // The lid-driven cavity problem isn't in H1, so
// lets estimate L2 error // lets estimate L2 error
u->error_norm = L2; u->error_norm = libMesh::L2;
error_estimator.reset(u); error_estimator.reset(u);
} else { } else {
...@@ -198,21 +192,22 @@ int main(int argc, char** argv) { ...@@ -198,21 +192,22 @@ int main(int argc, char** argv) {
// not in H1 - if we were doing more than a few // not in H1 - if we were doing more than a few
// timesteps we'd need to turn off or limit the // timesteps we'd need to turn off or limit the
// maximum level of our adaptivity eventually // maximum level of our adaptivity eventually
error_estimator.reset(new KellyErrorEstimator); error_estimator.reset(new libMesh::KellyErrorEstimator);
} }
libMesh::ErrorVector error;
error_estimator->estimate_error(system, error); error_estimator->estimate_error(system, error);
// Print out status at each adaptive step. // Print out status at each adaptive step.
Real global_error = error.l2_norm(); libMesh::out << "Adaptive step " << a_step << ": " << '\n';
libMesh::out << "Adaptive step " << a_step << ": " << std::endl;
const auto& global_error = error.l2_norm();
if (global_tolerance != 0.) if (global_tolerance != 0.)
libMesh::out << "Global_error = " << global_error << std::endl; libMesh::out << "Global_error = " << global_error << '\n';
if (global_tolerance != 0.) if (global_tolerance != 0.)
libMesh::out << "Worst element error = " << error.maximum() libMesh::out << "Worst element error = " << error.maximum()
<< ", mean = " << error.mean() << std::endl; << ", mean = " << error.mean() << '\n';
if (global_tolerance != 0.) { if (global_tolerance != 0.) {
// If we've reached our desired tolerance, we // If we've reached our desired tolerance, we
...@@ -236,8 +231,7 @@ int main(int argc, char** argv) { ...@@ -236,8 +231,7 @@ int main(int argc, char** argv) {
libMesh::out << "Refined mesh to " << mesh.n_active_elem() libMesh::out << "Refined mesh to " << mesh.n_active_elem()
<< " active elements and " << " active elements and "
<< equation_systems.n_active_dofs() << " active dofs." << equation_systems.n_active_dofs() << " active dofs.\n";
<< std::endl;
} }
// Do one last solve if necessary // Do one last solve if necessary
if (a_step == max_adaptivesteps) { if (a_step == max_adaptivesteps) {
...@@ -247,26 +241,27 @@ int main(int argc, char** argv) { ...@@ -247,26 +241,27 @@ int main(int argc, char** argv) {
} }
#ifdef LIBMESH_HAVE_EXODUS_API #ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems); libMesh::ExodusII_IO(mesh).write_equation_systems("out.e",
equation_systems);
#endif // #ifdef LIBMESH_HAVE_EXODUS_API #endif // #ifdef LIBMESH_HAVE_EXODUS_API
#ifdef LIBMESH_HAVE_GMV #ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems("out.gmv", equation_systems); libMesh::GMVIO(mesh).write_equation_systems("out.gmv", equation_systems);
#endif // #ifdef LIBMESH_HAVE_GMV #endif // #ifdef LIBMESH_HAVE_GMV
#ifdef LIBMESH_HAVE_FPARSER #ifdef LIBMESH_HAVE_FPARSER
// Check that we got close to the analytic solution // Check that we got close to the analytic solution
ExactSolution exact_sol(equation_systems); libMesh::ExactSolution exact_sol(equation_systems);
const std::string exact_str = const std::string exact_str =
(dim == 2) ? "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)"; (dim == 2) ? "sin(pi*x)*sin(pi*y)" : "sin(pi*x)*sin(pi*y)*sin(pi*z)";
ParsedFunction<Number> exact_func(exact_str); libMesh::ParsedFunction<libMesh::Number> exact_func(exact_str);
exact_sol.attach_exact_value(0, &exact_func); exact_sol.attach_exact_value(0, &exact_func);
exact_sol.compute_error("Heat", "T"); exact_sol.compute_error("Heat", "T");
Number err = exact_sol.l2_error("Heat", "T"); libMesh::Number err{exact_sol.l2_error("Heat", "T")};
// Print out the error value // Print out the error value
libMesh::out << "L2-Error is: " << err << std::endl; libMesh::out << "L2-Error is: " << err << '\n';
libmesh_assert_less(libmesh_real(err), 2e-3); libmesh_assert_less(libmesh_real(err), 2e-3);
......
...@@ -13,17 +13,17 @@ ...@@ -13,17 +13,17 @@
#include "libmesh/zero_function.h" #include "libmesh/zero_function.h"
using namespace libMesh; using namespace libMesh;
// user defined literal for libmesh's Real
constexpr auto operator"" _r(long double n) { return libMesh::Real(n); }
void HeatSystem::init_data() { void HeatSystem::init_data() {
T_var = this->add_variable("T", static_cast<Order>(_fe_order), T_var = this->add_variable("T", static_cast<Order>(_fe_order),
Utility::string_to_enum<FEFamily>(_fe_family)); Utility::string_to_enum<FEFamily>(_fe_family));
std::set<boundary_id_type> bdys{0, 1, 2};
std::vector<unsigned int> T_only(1, T_var);
ZeroFunction<Number> zero; ZeroFunction<Number> zero;
this->get_dof_map().add_dirichlet_boundary( this->get_dof_map().add_dirichlet_boundary(
DirichletBoundary(bdys, T_only, &zero)); DirichletBoundary({0, 1, 2}, {T_var}, &zero));
// Do the parent's initialization after variables are defined // Do the parent's initialization after variables are defined
FEMSystem::init_data(); FEMSystem::init_data();
...@@ -32,16 +32,12 @@ void HeatSystem::init_data() { ...@@ -32,16 +32,12 @@ void HeatSystem::init_data() {
} }
void HeatSystem::init_context(DiffContext& context) { void HeatSystem::init_context(DiffContext& context) {
FEMContext& c = libmesh_cast_ref<FEMContext&>(context); auto& c = libmesh_cast_ref<FEMContext&>(context);
const std::set<unsigned char>& elem_dims = c.elem_dimensions(); const auto& elem_dims = c.elem_dimensions();
for (std::set<unsigned char>::const_iterator dim_it = elem_dims.begin();
dim_it != elem_dims.end(); ++dim_it) {
const unsigned char dim = *dim_it;
FEBase* fe = libmesh_nullptr;
for (const auto& dim : elem_dims) {
FEBase* fe = nullptr;
c.get_element_fe(T_var, fe, dim); c.get_element_fe(T_var, fe, dim);
fe->get_JxW(); // For integration fe->get_JxW(); // For integration
...@@ -55,31 +51,31 @@ void HeatSystem::init_context(DiffContext& context) { ...@@ -55,31 +51,31 @@ void HeatSystem::init_context(DiffContext& context) {
bool HeatSystem::element_time_derivative(bool request_jacobian, bool HeatSystem::element_time_derivative(bool request_jacobian,
DiffContext& context) { DiffContext& context) {
FEMContext& c = libmesh_cast_ref<FEMContext&>(context); auto& c = libmesh_cast_ref<FEMContext&>(context);
const unsigned int mesh_dim = c.get_system().get_mesh().mesh_dimension(); const auto& mesh_dim = c.get_system().get_mesh().mesh_dimension();
// First we get some references to cell-specific data that // First we get some references to cell-specific data that
// will be used to assemble the linear system. // will be used to assemble the linear system.
const unsigned int dim = c.get_elem().dim(); const auto& dim = c.get_elem().dim();
FEBase* fe = libmesh_nullptr; FEBase* fe = nullptr;
c.get_element_fe(T_var, fe, dim); c.get_element_fe(T_var, fe, dim);
// Element Jacobian * quadrature weights for interior integration // Element Jacobian * quadrature weights for interior integration
const std::vector<Real>& JxW = fe->get_JxW(); const auto& JxW = fe->get_JxW();
const std::vector<Point>& xyz = fe->get_xyz(); const auto& xyz = fe->get_xyz();
const std::vector<std::vector<Real>>& phi = fe->get_phi(); const auto& phi = fe->get_phi();
const std::vector<std::vector<RealGradient>>& dphi = fe->get_dphi(); const auto& dphi = fe->get_dphi();
// The number of local degrees of freedom in each variable // The number of local degrees of freedom in each variable
const unsigned int n_T_dofs = c.get_dof_indices(T_var).size(); const auto& n_T_dofs = c.get_dof_indices(T_var).size();
// The subvectors and submatrices we need to fill: // The subvectors and submatrices we need to fill:
DenseSubMatrix<Number>& K = c.get_elem_jacobian(T_var, T_var); auto& K = c.get_elem_jacobian(T_var, T_var);
DenseSubVector<Number>& F = c.get_elem_residual(T_var); auto& F = c.get_elem_residual(T_var);
// Now we will build the element Jacobian and residual. // Now we will build the element Jacobian and residual.
// Constructing the residual requires the solution and its // Constructing the residual requires the solution and its
...@@ -87,33 +83,31 @@ bool HeatSystem::element_time_derivative(bool request_jacobian, ...@@ -87,33 +83,31 @@ bool HeatSystem::element_time_derivative(bool request_jacobian,
// calculated at each quadrature point by summing the // calculated at each quadrature point by summing the
// solution degree-of-freedom values by the appropriate // solution degree-of-freedom values by the appropriate
// weight functions. // weight functions.
unsigned int n_qpoints = c.get_element_qrule().n_points(); auto n_qpoints = c.get_element_qrule().n_points();
for (unsigned int qp = 0; qp != n_qpoints; qp++) { for (std::size_t qp{0}; qp != n_qpoints; qp++) {
// Compute the solution gradient at the Newton iterate // Compute the solution gradient at the Newton iterate
Gradient grad_T = c.interior_gradient(T_var, qp); Gradient grad_T = c.interior_gradient(T_var, qp);
const Number k = _k[dim]; const auto& k = _k[dim];
const Point& p = xyz[qp];
const Number u_exact = +1.0; const auto& u_exact = +1.0_r;
// Only apply forcing to interior elements // Only apply forcing to interior elements
const Number forcing = const auto& forcing =
(dim == mesh_dim) ? -k * u_exact * (dim * libMesh::pi * libMesh::pi) (dim == mesh_dim) ? -k * u_exact * (dim * libMesh::pi * libMesh::pi)
: 0; : 0;
const Number JxWxNK = JxW[qp] * -k; const auto& JxWxNK = JxW[qp] * -k;
for (unsigned int i = 0; i != n_T_dofs; i++) for (std::size_t i{0}; i != n_T_dofs; i++)
F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]); F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
if (request_jacobian) { if (request_jacobian) {
const Number JxWxNKxD = const auto& JxWxNKxD =
JxWxNK * context.get_elem_solution_derivative(); JxWxNK * context.get_elem_solution_derivative();
for (unsigned int i = 0; i != n_T_dofs; i++) for (std::size_t i{0}; i != n_T_dofs; i++)
for (unsigned int j = 0; j != n_T_dofs; ++j) for (std::size_t j{0}; j != n_T_dofs; ++j)
K(i, j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]); K(i, j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
} }
} // end of the quadrature point qp-loop } // end of the quadrature point qp-loop
......
...@@ -37,9 +37,9 @@ class HeatSystem : public libMesh::FEMSystem { ...@@ -37,9 +37,9 @@ class HeatSystem : public libMesh::FEMSystem {
libMesh::Real _k[4]; libMesh::Real _k[4];
// The FE type to use // The FE type to use
std::string _fe_family; std::string _fe_family{};
unsigned int _fe_order; unsigned int _fe_order{};
// The variable index (yes, this will be 0...) // The variable index (yes, this will be 0...)
unsigned int T_var; unsigned int T_var{};
}; };
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment