Merge branch 'libmesh' into 'master'

Libmesh

added libMesh example

See merge request !75

CMakeLists.txt

@@ -23,6 +23,7 @@ set(DWARF_PREFIX 8_unstructured) # The prefix of the name of the binaries produc
 
 # Add the examples
 add_subdirectory(halo_exchange)
+add_subdirectory(libmesh)
 
 # ==================================================================================================
 

libmesh/.clang-format

+Language: Cpp
+DerivePointerAlignment: false
+PointerAlignment: Left
+IndentWidth: 4
+TabWidth: 4
+UseTab: ForIndentation
+ColumnLimit: 80
+AllowShortFunctionsOnASingleLine: true
+AllowShortIfStatementsOnASingleLine: true
+AllowShortLoopsOnASingleLine: true
+AllowShortBlocksOnASingleLine: true
\ No newline at end of file

libmesh/CMakeLists.txt

+# Packages are optional: if they are not present, certain code samples are not compiled
+cmake_minimum_required(VERSION 2.8.10 FATAL_ERROR)
+
+find_package(MPI)      # Built-in in CMake
+find_package(Threads)
+
+include(${CMAKE_CURRENT_SOURCE_DIR}/../../cmake/common.cmake)
+
+set(LIBMESH_INCLUDE_PATH /usr/local/include CACHE PATH "path to libMesh include directory")
+set(LIBMESH_LIBRARY_PATH /usr/local/lib CACHE PATH "path to libMesh library directory")
+set(LIBMESH_LIBRARY mesh_opt CACHE FILEPATH "name of libmesh library" )
+
+# ==================================================================================================
+
+if ("${DWARF_PREFIX}" STREQUAL "")
+  set(DWARF_PREFIX 8_unstructured)
+endif()
+
+set(NAME ${DWARF_PREFIX}_libmesh)
+
+#TODO: check for libmesh too
+if (MPI_FOUND)
+	enable_language(CXX)
+    include_directories(${MPI_INCLUDE_PATH} ${LIBMESH_INCLUDE_PATH})
+    link_directories(${LIBMESH_LIBRARY_PATH}) 
+    add_executable(${NAME} fem_system_ex4.cpp heatsystem.cpp)
+    set(CMAKE_BUILD_TYPE RelWithDebInfo)
+    if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
+    	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
+    elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Intel")
+    	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -xHost -std=c++14")
+    endif()
+    set_target_properties(${NAME} PROPERTIES CXX_STANDARD 14 CXX_STANDARD_REQUIRED YES)
+    target_link_libraries(${NAME} ${LIBMESH_LIBRARY} ${MPI_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
+    install(TARGETS ${NAME} DESTINATION bin)
+    message("** Enabling '${NAME}': with MPI")
+else()
+    message("## Skipping '${NAME}': MPI support missing")
+#    dummy_install(${NAME} "MPI")
+endif()
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${C_FLAGS}")
+
+unset(NAME)
+# ==================================================================================================

libmesh/README.md

+# README - libMesh Example
+
+## Description
+
+The libMesh library provides a framework for the numerical simulation of partial differential equations using arbitrary unstructured discretizations on serial and parallel platforms. (http://libmesh.github.io/index.html)
+
+This example is a slightly modified version of the official libMesh example ["fem_system_ex4"] (https://github.com/libMesh/libmesh/tree/master/examples/fem_system/fem_system_ex4). In this example, a heat transfer equation is solved with an FEM System.
+
+The example demonstrates:
+ * Initialization of libMesh
+ * Definition of a equation system
+ * Usage of adaptive mesh refinement
+ * IO: reading and exporting meshes in different file formats
+
+The example is structured as followed:
+ * `fem_system_ex4.C`: The main program, Initialization, IO, refinement
+ * `heatsystem.C|h`: Heat system with laplace heat equation
+ * `meshes/bridge.e|xda`: FEM-Mesh of a bridge
+ 
+ Screenshot of the result:
+ ![Screenshot of Result](hpc_kernel_samples/unstructured_grids/libmesh/meshes/bridge_screenshot.PNG)
+ 
+ Further reading: [PRACE Summer School 2013 Slides](http://www.training.prace-ri.eu/uploads/tx_pracetmo/libmesh.pdf)
+ 
+## Release Date
+ 
+ 2016-09-08
+ 
+## Version History
+
+ *  2016-09-08 Initial Release on PRACE CodeVault repository
+ 
+## Contributors
+
+ * Thomas Steinreiter - [thomas.steinreiter@risc-software.at](mailto:thomas.steinreiter@risc-software.at)
+
+## Copyright
+
+This code is available under LGPL, Version 2.1 - see also the license file in the CodeVault root directory.
+
+## Languages
+
+This sample is written in C++ 11.
+
+
+## Parallelisation
+
+This sample uses PETSCs internal parallelisation.
+
+
+## Level of the code sample complexity
+
+Advanced
+
+
+## Compiling
+
+Follow the compilation instructions given in the main directory of the kernel samples directory (`/hpc_kernel_samples/README.md`).
+
+Note: libMesh needs to be built against the correct version of MPI. PETSC is needed for parallelism.
+
+## Running
+
+To run the program, use something similar to
+
+	mpiexec -n [nprocs] ./8_unstructured_libmesh
+
+either on the command line or in your batch script, where `nprocs` specifies the number of processes used. Node: if `nprocs` > 1, libMesh must be built with PETSC enabled.
+ 
+### Example
+ 
+If you run
+	
+	mpiexec -n 8 ./8_unstructured_libmesh
+
+the output should look similar to
+
+	 Mesh Information:
+	  elem_dimensions()={2, 3}
+	  spatial_dimension()=3
+	  n_nodes()=25227
+	    n_local_nodes()=3322
+	  n_elem()=127818
+	    n_local_elem()=15972
+	    n_active_elem()=127818
+	  n_subdomains()=2
+	  n_partitions()=8
+	  n_processors()=8
+	  n_threads()=1
+	  processor_id()=0
+	
+	 EquationSystems
+	  n_systems()=1
+	   System #0, "Heat"
+	    Type "Implicit"
+	    Variables="T"
+	    Finite Element Types="LAGRANGE"
+	    Approximation Orders="FIRST"
+	    n_dofs()=25227
+	    n_local_dofs()=3322
+	    n_constrained_dofs()=1887
+	    n_local_constrained_dofs()=122
+	    n_vectors()=1
+	    n_matrices()=1
+	    DofMap Sparsity
+	      Average  On-Processor Bandwidth <= 13.6686
+	      Average Off-Processor Bandwidth <= 0.533635
+	      Maximum  On-Processor Bandwidth <= 27
+	      Maximum Off-Processor Bandwidth <= 15
+	    DofMap Constraints
+	      Number of DoF Constraints = 1887
+	      Average DoF Constraint Length= 0
+	
+	Assembling the System
+	*** Warning, This code is deprecated, and likely to be removed in future library versions! /usr/local/include/libmesh/libmesh_common.h, line 497, compiled Sep  1 2016 at 15:28:09 ***
+	Nonlinear Residual: 8495.66
+	Linear solve starting, tolerance 0.001
+	Linear solve finished, step 55, residual 3.10027
+	Trying full Newton step
+	  Current Residual: 14.2871
+	  Nonlinear step: |du|/|u| = 1, |du| = 510394
+	Assembling the System
+	Nonlinear Residual: 14.2871
+	Linear solve starting, tolerance 0.001
+	Linear solve finished, step 58, residual 0.00263574
+	Trying full Newton step
+	  Current Residual: 0.012002
+	  Nonlinear step: |du|/|u| = 0.000564056, |du| = 287.935
+	Assembling the System
+	Nonlinear Residual: 0.012002
+	Linear solve starting, tolerance 1.2002e-05
+	Linear solve finished, step 94, residual 2.75721e-08
+	Trying full Newton step
+	  Current Residual: 1.29858e-07
+	  Nonlinear solver converged, step 2, residual reduction 1.52852e-11 < 1e-07
+	  Nonlinear solver relative step size 5.66919e-07 > 1e-07
+	L2-Error is: 59155
\ No newline at end of file

libmesh/fem_system_ex4.cpp

+// 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
+
+// <h1>FEMSystem Example 4 - Mixed-dimension heat transfer equation
+// with FEMSystem</h1>
+// \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 <iomanip>
+
+// 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<HeatSystem>("Heat");
+
+	// Solve this as a steady system
+	system.time_solver = libMesh::UniquePtr<libMesh::TimeSolver>(
+	    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<libMesh::ErrorEstimator> 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);
+
+			auto 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<libMesh::Number> 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;
+}

libmesh/fem_system_ex4.in

+# The global FEM error tolerance at each timestep
+# Make this nonzero to solve to a specified tolerance
+# This will probably break with KellyErrorIndicator
+# const Real global_tolerance = 1.e-3;
+global_tolerance = 0.0
+
+# The desired number of active mesh elements
+# Make this nonzero to solve to a specified mesh size
+n_elements = 1200
+
+# The coarse grid size from which to start adaptivity
+coarsegridsize = 20
+
+# The maximum number of adaptive steps per timestep
+max_adaptivesteps = 2
+
+# Turn this off to see the NewtonSolver chatter
+solver_quiet = false
+
+# Solve the 2D or 3D problem
+dimension = 3
+
+# Nonlinear solver tolerance
+relative_step_tolerance = 1e-7
+
+relative_residual_tolerance = 1.e-7

libmesh/heatsystem.cpp

+#include "heatsystem.hpp"
+
+#include "libmesh/dirichlet_boundaries.h"
+#include "libmesh/dof_map.h"
+#include "libmesh/elem.h"
+#include "libmesh/fe_base.h"
+#include "libmesh/fe_interface.h"
+#include "libmesh/fem_context.h"
+#include "libmesh/getpot.h"
+#include "libmesh/mesh.h"
+#include "libmesh/quadrature.h"
+#include "libmesh/string_to_enum.h"
+#include "libmesh/zero_function.h"
+
+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() {
+	T_var = this->add_variable("T", static_cast<Order>(_fe_order),
+	                           Utility::string_to_enum<FEFamily>(_fe_family));
+
+	ZeroFunction<Number> zero;
+
+	this->get_dof_map().add_dirichlet_boundary(
+	    DirichletBoundary({0, 1, 2}, {T_var}, &zero));
+
+	// Do the parent's initialization after variables are defined
+	FEMSystem::init_data();
+
+	this->time_evolving(0);
+}
+
+void HeatSystem::init_context(DiffContext& context) {
+	auto& c = libmesh_cast_ref<FEMContext&>(context);
+
+	const auto& elem_dims = c.elem_dimensions();
+
+	for (const auto& dim : elem_dims) {
+		FEBase* fe = nullptr;
+		c.get_element_fe(T_var, fe, dim);
+
+		fe->get_JxW();  // For integration
+		fe->get_dphi(); // For bilinear form
+		fe->get_xyz();  // For forcing
+		fe->get_phi();  // For forcing
+	}
+
+	FEMSystem::init_context(context);
+}
+
+bool HeatSystem::element_time_derivative(bool request_jacobian,
+                                         DiffContext& context) {
+	auto& c = libmesh_cast_ref<FEMContext&>(context);
+
+	const auto& mesh_dim = c.get_system().get_mesh().mesh_dimension();
+
+	// First we get some references to cell-specific data that
+	// will be used to assemble the linear system.
+	const auto& dim = c.get_elem().dim();
+	FEBase* fe = nullptr;
+	c.get_element_fe(T_var, fe, dim);
+
+	// Element Jacobian * quadrature weights for interior integration
+	const auto& JxW = fe->get_JxW();
+
+	const auto& xyz = fe->get_xyz();
+
+	const auto& phi = fe->get_phi();
+
+	const auto& dphi = fe->get_dphi();
+
+	// The number of local degrees of freedom in each variable
+	const auto& n_T_dofs = c.get_dof_indices(T_var).size();
+
+	// The subvectors and submatrices we need to fill:
+	auto& K = c.get_elem_jacobian(T_var, T_var);
+	auto& F = c.get_elem_residual(T_var);
+
+	// Now we will build the element Jacobian and residual.
+	// Constructing the residual requires the solution and its
+	// gradient from the previous timestep.  This must be
+	// calculated at each quadrature point by summing the
+	// solution degree-of-freedom values by the appropriate
+	// weight functions.
+	auto n_qpoints = c.get_element_qrule().n_points();
+
+	for (std::size_t qp{0}; qp != n_qpoints; qp++) {
+		// Compute the solution gradient at the Newton iterate
+		Gradient grad_T = c.interior_gradient(T_var, qp);
+
+		const auto& k = _k[dim];
+
+		const auto& u_exact = +1.0_r;
+
+		// Only apply forcing to interior elements
+		const auto& forcing =
+		    (dim == mesh_dim) ? -k * u_exact * (dim * libMesh::pi * libMesh::pi)
+		                      : 0;
+
+		const auto& JxWxNK = JxW[qp] * -k;
+
+		for (std::size_t i{0}; i != n_T_dofs; i++)
+			F(i) += JxWxNK * (grad_T * dphi[i][qp] + forcing * phi[i][qp]);
+		if (request_jacobian) {
+			const auto& JxWxNKxD =
+			    JxWxNK * context.get_elem_solution_derivative();
+
+			for (std::size_t i{0}; i != n_T_dofs; i++)
+				for (std::size_t j{0}; j != n_T_dofs; ++j)
+					K(i, j) += JxWxNKxD * (dphi[i][qp] * dphi[j][qp]);
+		}
+	} // end of the quadrature point qp-loop
+
+	return request_jacobian;