The source file exact_solution.C with comments:
#include <math.h>
Mesh library includes
#include "libmesh/libmesh_common.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
/**
* This is the exact solution that
* we are trying to obtain. We will solve
*
* - (u_xx + u_yy) = f
*
* and take a finite difference approximation using this
* function to get f. This is the well-known "method of
* manufactured solutions".
*/
Real exact_solution (const int component,
const Real x,
const Real y,
const Real z = 0.)
{
static const Real pi = acos(-1.);
switch( component )
{
case 0:
return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z);
case 1:
return sin(.5*pi*x)*cos(.5*pi*y)*cos(.5*pi*z);
case 2:
return sin(.5*pi*x)*cos(.5*pi*y)*cos(.5*pi*z)*cos(.5*pi*x*y*z);
default:
libmesh_error();
}
dummy
return 0.0;
}
The source file vector_fe_ex1.C with comments:
Vector Finite Element Example 1 - Solving an uncoupled Poisson Problem
This is the first vector FE example program. It builds on the introduction_ex3 example program by showing how to solve a simple uncoupled Poisson system using vector Lagrange elements.
C++ include files that we need
#include <iostream>
#include <algorithm>
#include <math.h>
Basic include files needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/gmv_io.h"
Define the Finite Element object.
#include "libmesh/fe.h"
Define Gauss quadrature rules.
#include "libmesh/quadrature_gauss.h"
Define useful datatypes for finite element
matrix and vector components.
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/elem.h"
Define the DofMap, which handles degree of freedom
indexing.
#include "libmesh/dof_map.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Function prototype. This is the function that will assemble
the linear system for our Poisson problem. Note that the
function will take the EquationSystems object and the
name of the system we are assembling as input. From the
EquationSystems object we have access to the Mesh and
other objects we might need.
void assemble_poisson(EquationSystems& es,
const std::string& system_name);
Function prototype for the exact solution.
Real exact_solution (const int component,
const Real x,
const Real y,
const Real z = 0.);
int main (int argc, char** argv)
{
Initialize libraries.
LibMeshInit init (argc, argv);
Brief message to the user regarding the program name
and command line arguments.
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
Skip this 2D example if libMesh was compiled as 1D-only.
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Mesh mesh;
Use the MeshTools::Generation mesh generator to create a uniform
2D grid on the square [-1,1]^2. We instruct the mesh generator
to build a mesh of 15x15 QUAD9 elements.
MeshTools::Generation::build_square (mesh,
15, 15,
-1., 1.,
-1., 1.,
QUAD9);
Print information about the mesh to the screen.
mesh.print_info();
Create an equation systems object.
EquationSystems equation_systems (mesh);
Declare the Poisson system and its variables.
The Poisson system is another example of a steady system.
equation_systems.add_system<LinearImplicitSystem> ("Poisson");
Adds the variable "u" to "Poisson". "u"
will be approximated using second-order approximation
using vector Lagrange elements. Since the mesh is 2-D, "u" will
have two components.
equation_systems.get_system("Poisson").add_variable("u", SECOND, LAGRANGE_VEC);
Give the system a pointer to the matrix assembly
function. This will be called when needed by the
library.
equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
Initialize the data structures for the equation system.
equation_systems.init();
Prints information about the system to the screen.
equation_systems.print_info();
Solve the system "Poisson". Note that calling this
member will assemble the linear system and invoke
the default numerical solver. With PETSc the solver can be
controlled from the command line. For example,
you can invoke conjugate gradient with:
./vector_fe_ex1 -ksp_type cg
You can also get a nice X-window that monitors the solver convergence with:
./vector_fe_ex1 -ksp_xmonitor
if you linked against the appropriate X libraries when you built PETSc.
./vector_fe_ex1 -ksp_type cg
You can also get a nice X-window that monitors the solver convergence with:
./vector_fe_ex1 -ksp_xmonitor
if you linked against the appropriate X libraries when you built PETSc.
equation_systems.get_system("Poisson").solve();
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems( "out.e", equation_systems);
#endif
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems( "out.gmv", equation_systems);
#endif
All done.
return 0;
}
We now define the matrix assembly function for the
Poisson system. We need to first compute element
matrices and right-hand sides, and then take into
account the boundary conditions, which will be handled
via a penalty method.
void assemble_poisson(EquationSystems& es,
const std::string& system_name)
{
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "Poisson");
Get a constant reference to the mesh object.
const MeshBase& mesh = es.get_mesh();
The dimension that we are running
const unsigned int dim = mesh.mesh_dimension();
Get a reference to the LinearImplicitSystem we are solving
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson");
A reference to the DofMap object for this system. The DofMap
object handles the index translation from node and element numbers
to degree of freedom numbers. We will talk more about the DofMap
in future examples.
const DofMap& dof_map = system.get_dof_map();
Get a constant reference to the Finite Element type
for the first (and only) variable in the system.
FEType fe_type = dof_map.variable_type(0);
Build a Finite Element object of the specified type.
Note that FEVectorBase is a typedef for the templated FE
class.
AutoPtr<FEVectorBase> fe (FEVectorBase::build(dim, fe_type));
A 5th order Gauss quadrature rule for numerical integration.
QGauss qrule (dim, FIFTH);
Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
Declare a special finite element object for
boundary integration.
AutoPtr<FEVectorBase> fe_face (FEVectorBase::build(dim, fe_type));
Boundary integration requires one quadraure rule,
with dimensionality one less than the dimensionality
of the element.
QGauss qface(dim-1, FIFTH);
Tell the finite element object to use our
quadrature rule.
fe_face->attach_quadrature_rule (&qface);
Here we define some references to cell-specific data that
will be used to assemble the linear system.
The element Jacobian * quadrature weight at each integration point.
The element Jacobian * quadrature weight at each integration point.
const std::vector<Real>& JxW = fe->get_JxW();
The physical XY locations of the quadrature points on the element.
These might be useful for evaluating spatially varying material
properties at the quadrature points.
const std::vector<Point>& q_point = fe->get_xyz();
The element shape functions evaluated at the quadrature points.
Notice the shape functions are a vector rather than a scalar.
const std::vector<std::vector<RealGradient> >& phi = fe->get_phi();
The element shape function gradients evaluated at the quadrature
points. Notice that the shape function gradients are a tensor.
const std::vector<std::vector<RealTensor> >& dphi = fe->get_dphi();
Define data structures to contain the element matrix
and right-hand-side vector contribution. Following
basic finite element terminology we will denote these
"Ke" and "Fe". These datatypes are templated on
Number, which allows the same code to work for real
or complex numbers.
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
This vector will hold the degree of freedom indices for
the element. These define where in the global system
the element degrees of freedom get mapped.
std::vector<dof_id_type> dof_indices;
Now we will loop over all the elements in the mesh.
We will compute the element matrix and right-hand-side
contribution.
Element iterators are a nice way to iterate through all the elements, or all the elements that have some property. The iterator el will iterate from the first to the last element on the local processor. The iterator end_el tells us when to stop. It is smart to make this one const so that we don't accidentally mess it up! In case users later modify this program to include refinement, we will be safe and will only consider the active elements; hence we use a variant of the \p active_elem_iterator.
Element iterators are a nice way to iterate through all the elements, or all the elements that have some property. The iterator el will iterate from the first to the last element on the local processor. The iterator end_el tells us when to stop. It is smart to make this one const so that we don't accidentally mess it up! In case users later modify this program to include refinement, we will be safe and will only consider the active elements; hence we use a variant of the \p active_elem_iterator.
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
Loop over the elements. Note that ++el is preferred to
el++ since the latter requires an unnecessary temporary
object.
for ( ; el != end_el ; ++el)
{
Store a pointer to the element we are currently
working on. This allows for nicer syntax later.
const Elem* elem = *el;
Get the degree of freedom indices for the
current element. These define where in the global
matrix and right-hand-side this element will
contribute to.
dof_map.dof_indices (elem, dof_indices);
Compute the element-specific data for the current
element. This involves computing the location of the
quadrature points (q_point) and the shape functions
(phi, dphi) for the current element.
fe->reinit (elem);
Zero the element matrix and right-hand side before
summing them. We use the resize member here because
the number of degrees of freedom might have changed from
the last element. Note that this will be the case if the
element type is different (i.e. the last element was a
triangle, now we are on a quadrilateral).
The DenseMatrix::resize() and the DenseVector::resize() members will automatically zero out the matrix and vector.
The DenseMatrix::resize() and the DenseVector::resize() members will automatically zero out the matrix and vector.
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
Now loop over the quadrature points. This handles
the numeric integration.
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Now we will build the element matrix. This involves
a double loop to integrate the test funcions (i) against
the trial functions (j).
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*( dphi[i][qp].contract(dphi[j][qp]) );
}
This is the end of the matrix summation loop
Now we build the element right-hand-side contribution.
This involves a single loop in which we integrate the
"forcing function" in the PDE against the test functions.
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
const Real eps = 1.e-3;
"f" is the forcing function for the Poisson equation.
In this case we set f to be a finite difference
Laplacian approximation to the (known) exact solution.
We will use the second-order accurate FD Laplacian approximation, which in 2D is
u_xx + u_yy = (u(i,j-1) + u(i,j+1) + u(i-1,j) + u(i+1,j) + -4*u(i,j))/h^2
Since the value of the forcing function depends only on the location of the quadrature point (q_point[qp]) we will compute it here, outside of the i-loop
We will use the second-order accurate FD Laplacian approximation, which in 2D is
u_xx + u_yy = (u(i,j-1) + u(i,j+1) + u(i-1,j) + u(i+1,j) + -4*u(i,j))/h^2
Since the value of the forcing function depends only on the location of the quadrature point (q_point[qp]) we will compute it here, outside of the i-loop
const Real fx = -(exact_solution(0,x,y-eps) +
exact_solution(0,x,y+eps) +
exact_solution(0,x-eps,y) +
exact_solution(0,x+eps,y) -
4.*exact_solution(0,x,y))/eps/eps;
const Real fy = -(exact_solution(1,x,y-eps) +
exact_solution(1,x,y+eps) +
exact_solution(1,x-eps,y) +
exact_solution(1,x+eps,y) -
4.*exact_solution(1,x,y))/eps/eps;
const RealGradient f( fx, fy );
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*f*phi[i][qp];
}
}
We have now reached the end of the RHS summation,
and the end of quadrature point loop, so
the interior element integration has
been completed. However, we have not yet addressed
boundary conditions. For this example we will only
consider simple Dirichlet boundary conditions.
There are several ways Dirichlet boundary conditions can be imposed. A simple approach, which works for interpolary bases like the standard Lagrange polynomials, is to assign function values to the degrees of freedom living on the domain boundary. This works well for interpolary bases, but is more difficult when non-interpolary (e.g Legendre or Hierarchic) bases are used.
Dirichlet boundary conditions can also be imposed with a "penalty" method. In this case essentially the L2 projection of the boundary values are added to the matrix. The projection is multiplied by some large factor so that, in floating point arithmetic, the existing (smaller) entries in the matrix and right-hand-side are effectively ignored.
This amounts to adding a term of the form (in latex notation)
\frac{1}{\epsilon} \int_{\delta \Omega} \phi_i \phi_j = \frac{1}{\epsilon} \int_{\delta \Omega} u \phi_i
where
\frac{1}{\epsilon} is the penalty parameter, defined such that \epsilon << 1
There are several ways Dirichlet boundary conditions can be imposed. A simple approach, which works for interpolary bases like the standard Lagrange polynomials, is to assign function values to the degrees of freedom living on the domain boundary. This works well for interpolary bases, but is more difficult when non-interpolary (e.g Legendre or Hierarchic) bases are used.
Dirichlet boundary conditions can also be imposed with a "penalty" method. In this case essentially the L2 projection of the boundary values are added to the matrix. The projection is multiplied by some large factor so that, in floating point arithmetic, the existing (smaller) entries in the matrix and right-hand-side are effectively ignored.
This amounts to adding a term of the form (in latex notation)
\frac{1}{\epsilon} \int_{\delta \Omega} \phi_i \phi_j = \frac{1}{\epsilon} \int_{\delta \Omega} u \phi_i
where
\frac{1}{\epsilon} is the penalty parameter, defined such that \epsilon << 1
{
The following loop is over the sides of the element.
If the element has no neighbor on a side then that
side MUST live on a boundary of the domain.
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor(side) == NULL)
{
The value of the shape functions at the quadrature
points.
const std::vector<std::vector<RealGradient> >& phi_face = fe_face->get_phi();
The Jacobian * Quadrature Weight at the quadrature
points on the face.
const std::vector<Real>& JxW_face = fe_face->get_JxW();
The XYZ locations (in physical space) of the
quadrature points on the face. This is where
we will interpolate the boundary value function.
const std::vector<Point>& qface_point = fe_face->get_xyz();
Compute the shape function values on the element
face.
fe_face->reinit(elem, side);
Loop over the face quadrature points for integration.
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
The location on the boundary of the current
face quadrature point.
const Real xf = qface_point[qp](0);
const Real yf = qface_point[qp](1);
The penalty value. \frac{1}{\epsilon}
in the discussion above.
const Real penalty = 1.e10;
The boundary values.
const RealGradient f( exact_solution(0, xf, yf),
exact_solution(1, xf, yf) );
Matrix contribution of the L2 projection.
for (unsigned int i=0; i<phi_face.size(); i++)
for (unsigned int j=0; j<phi_face.size(); j++)
Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
Right-hand-side contribution of the L2
projection.
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) += JxW_face[qp]*penalty*f*phi_face[i][qp];
}
}
}
We have now finished the quadrature point loop,
and have therefore applied all the boundary conditions.
If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
The element matrix and right-hand-side are now built for this element. Add them to the global matrix and right-hand-side vector. The SparseMatrix::add_matrix() and NumericVector::add_vector() members do this for us.
If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
The element matrix and right-hand-side are now built for this element. Add them to the global matrix and right-hand-side vector. The SparseMatrix::add_matrix() and NumericVector::add_vector() members do this for us.
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
All done!
}
The source file exact_solution.C without comments:
#include <math.h>
#include "libmesh/libmesh_common.h"
using namespace libMesh;
/**
* This is the exact solution that
* we are trying to obtain. We will solve
*
* - (u_xx + u_yy) = f
*
* and take a finite difference approximation using this
* function to get f. This is the well-known "method of
* manufactured solutions".
*/
Real exact_solution (const int component,
const Real x,
const Real y,
const Real z = 0.)
{
static const Real pi = acos(-1.);
switch( component )
{
case 0:
return cos(.5*pi*x)*sin(.5*pi*y)*cos(.5*pi*z);
case 1:
return sin(.5*pi*x)*cos(.5*pi*y)*cos(.5*pi*z);
case 2:
return sin(.5*pi*x)*cos(.5*pi*y)*cos(.5*pi*z)*cos(.5*pi*x*y*z);
default:
libmesh_error();
}
return 0.0;
}
The source file vector_fe_ex1.C without comments:
#include <iostream>
#include <algorithm>
#include <math.h>
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/equation_systems.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/gmv_io.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/elem.h"
#include "libmesh/dof_map.h"
using namespace libMesh;
void assemble_poisson(EquationSystems& es,
const std::string& system_name);
Real exact_solution (const int component,
const Real x,
const Real y,
const Real z = 0.);
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
std::cout << "Running " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Mesh mesh;
MeshTools::Generation::build_square (mesh,
15, 15,
-1., 1.,
-1., 1.,
QUAD9);
mesh.print_info();
EquationSystems equation_systems (mesh);
equation_systems.add_system<LinearImplicitSystem> ("Poisson");
equation_systems.get_system("Poisson").add_variable("u", SECOND, LAGRANGE_VEC);
equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);
equation_systems.init();
equation_systems.print_info();
equation_systems.get_system("Poisson").solve();
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems( "out.e", equation_systems);
#endif
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems( "out.gmv", equation_systems);
#endif
return 0;
}
void assemble_poisson(EquationSystems& es,
const std::string& system_name)
{
libmesh_assert_equal_to (system_name, "Poisson");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson");
const DofMap& dof_map = system.get_dof_map();
FEType fe_type = dof_map.variable_type(0);
AutoPtr<FEVectorBase> fe (FEVectorBase::build(dim, fe_type));
QGauss qrule (dim, FIFTH);
fe->attach_quadrature_rule (&qrule);
AutoPtr<FEVectorBase> fe_face (FEVectorBase::build(dim, fe_type));
QGauss qface(dim-1, FIFTH);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Point>& q_point = fe->get_xyz();
const std::vector<std::vector<RealGradient> >& phi = fe->get_phi();
const std::vector<std::vector<RealTensor> >& dphi = fe->get_dphi();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
std::vector<dof_id_type> dof_indices;
MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
for ( ; el != end_el ; ++el)
{
const Elem* elem = *el;
dof_map.dof_indices (elem, dof_indices);
fe->reinit (elem);
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
for (unsigned int i=0; i<phi.size(); i++)
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*( dphi[i][qp].contract(dphi[j][qp]) );
}
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
const Real eps = 1.e-3;
const Real fx = -(exact_solution(0,x,y-eps) +
exact_solution(0,x,y+eps) +
exact_solution(0,x-eps,y) +
exact_solution(0,x+eps,y) -
4.*exact_solution(0,x,y))/eps/eps;
const Real fy = -(exact_solution(1,x,y-eps) +
exact_solution(1,x,y+eps) +
exact_solution(1,x-eps,y) +
exact_solution(1,x+eps,y) -
4.*exact_solution(1,x,y))/eps/eps;
const RealGradient f( fx, fy );
for (unsigned int i=0; i<phi.size(); i++)
Fe(i) += JxW[qp]*f*phi[i][qp];
}
}
{
for (unsigned int side=0; side<elem->n_sides(); side++)
if (elem->neighbor(side) == NULL)
{
const std::vector<std::vector<RealGradient> >& phi_face = fe_face->get_phi();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<Point>& qface_point = fe_face->get_xyz();
fe_face->reinit(elem, side);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
const Real xf = qface_point[qp](0);
const Real yf = qface_point[qp](1);
const Real penalty = 1.e10;
const RealGradient f( exact_solution(0, xf, yf),
exact_solution(1, xf, yf) );
for (unsigned int i=0; i<phi_face.size(); i++)
for (unsigned int j=0; j<phi_face.size(); j++)
Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
for (unsigned int i=0; i<phi_face.size(); i++)
Fe(i) += JxW_face[qp]*penalty*f*phi_face[i][qp];
}
}
}
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
}
The console output of the program:
***************************************************************
* Running Example vector_fe_ex1:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Running /workspace/libmesh/examples/vector_fe/vector_fe_ex1/.libs/lt-example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=961
n_local_nodes()=495
n_elem()=225
n_local_elem()=112
n_active_elem()=225
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Poisson"
Type "LinearImplicit"
Variables="u"
Finite Element Types="LAGRANGE_VEC", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="SECOND", "THIRD"
n_dofs()=1922
n_local_dofs()=990
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 29.6982
Average Off-Processor Bandwidth <= 1.04058
Maximum On-Processor Bandwidth <= 52
Maximum Off-Processor Bandwidth <= 28
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/vector_fe/vector_fe_ex1/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:53:36 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012
Max Max/Min Avg Total
Time (sec): 2.012e-01 1.00005 2.012e-01
Objects: 6.400e+01 1.00000 6.400e+01
Flops: 2.073e+07 1.09946 1.979e+07 3.958e+07
Flops/sec: 1.030e+08 1.09952 9.837e+07 1.967e+08
MPI Messages: 4.150e+01 1.00000 4.150e+01 8.300e+01
MPI Message Lengths: 4.371e+04 1.00000 1.053e+03 8.743e+04
MPI Reductions: 1.140e+02 1.00000
Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N flops
and VecAXPY() for complex vectors of length N --> 8N flops
Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions --
Avg %Total Avg %Total counts %Total Avg %Total counts %Total
0: Main Stage: 2.0111e-01 100.0% 3.9579e+07 100.0% 8.300e+01 100.0% 1.053e+03 100.0% 1.130e+02 99.1%
------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
Count: number of times phase was executed
Time and Flops: Max - maximum over all processors
Ratio - ratio of maximum to minimum over all processors
Mess: number of messages sent
Avg. len: average message length
Reduct: number of global reductions
Global: entire computation
Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
%T - percent time in this phase %f - percent flops in this phase
%M - percent messages in this phase %L - percent message lengths in this phase
%R - percent reductions in this phase
Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event Count Time (sec) Flops --- Global --- --- Stage --- Total
Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T %f %M %L %R %T %f %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------
--- Event Stage 0: Main Stage
VecMDot 32 1.0 5.3120e-04 1.8 9.26e+05 1.1 0.0e+00 0.0e+00 3.2e+01 0 5 0 0 28 0 5 0 0 28 3385
VecNorm 35 1.0 1.2898e-04 1.2 6.93e+04 1.1 0.0e+00 0.0e+00 3.5e+01 0 0 0 0 31 0 0 0 0 31 1043
VecScale 34 1.0 3.3855e-05 1.2 3.37e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1930
VecCopy 3 1.0 4.0531e-06 1.4 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecSet 42 1.0 2.2411e-05 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 4 1.0 6.5200e-03 1.0 7.92e+03 1.1 0.0e+00 0.0e+00 0.0e+00 3 0 0 0 0 3 0 0 0 0 2
VecMAXPY 34 1.0 2.8110e-04 1.1 9.90e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 5 0 0 0 0 5 0 0 0 6838
VecAssemblyBegin 3 1.0 5.9843e-05 1.1 0.00e+00 0.0 2.0e+00 5.8e+02 9.0e+00 0 0 2 1 8 0 0 2 1 8 0
VecAssemblyEnd 3 1.0 1.9550e-05 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecScatterBegin 35 1.0 8.0585e-05 1.1 0.00e+00 0.0 7.0e+01 8.1e+02 0.0e+00 0 0 84 65 0 0 0 84 65 0 0
VecScatterEnd 35 1.0 1.8604e-0342.6 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecNormalize 34 1.0 2.0266e-04 1.1 1.01e+05 1.1 0.0e+00 0.0e+00 3.4e+01 0 0 0 0 30 0 0 0 0 30 967
MatMult 34 1.0 2.7747e-03 2.8 2.04e+06 1.1 6.8e+01 7.9e+02 0.0e+00 1 10 82 62 0 1 10 82 62 0 1412
MatSolve 35 1.0 2.7785e-03 1.1 8.63e+06 1.1 0.0e+00 0.0e+00 0.0e+00 1 42 0 0 0 1 42 0 0 0 5954
MatLUFactorNum 1 1.0 3.8412e-03 1.1 8.03e+06 1.1 0.0e+00 0.0e+00 0.0e+00 2 38 0 0 0 2 38 0 0 0 3953
MatILUFactorSym 1 1.0 1.3245e-02 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 3.0e+00 6 0 0 0 3 6 0 0 0 3 0
MatAssemblyBegin 2 1.0 9.0718e-04 6.5 0.00e+00 0.0 3.0e+00 9.2e+03 4.0e+00 0 0 4 32 4 0 0 4 32 4 0
MatAssemblyEnd 2 1.0 3.9005e-04 1.1 0.00e+00 0.0 4.0e+00 2.0e+02 8.0e+00 0 0 5 1 7 0 0 5 1 7 0
MatGetRowIJ 1 1.0 5.0068e-06 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetOrdering 1 1.0 9.7990e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 4.0e+00 0 0 0 0 4 0 0 0 0 4 0
MatZeroEntries 3 1.0 9.3937e-05 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPGMRESOrthog 32 1.0 8.2159e-04 1.4 1.85e+06 1.1 0.0e+00 0.0e+00 3.2e+01 0 9 0 0 28 0 9 0 0 28 4378
KSPSetUp 2 1.0 7.7963e-05 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
KSPSolve 1 1.0 2.9279e-02 1.0 2.07e+07 1.1 6.8e+01 7.9e+02 7.6e+01 15100 82 62 67 15100 82 62 67 1352
PCSetUp 2 1.0 1.7541e-02 1.1 8.03e+06 1.1 0.0e+00 0.0e+00 9.0e+00 8 38 0 0 8 8 38 0 0 8 866
PCSetUpOnBlocks 1 1.0 1.7272e-02 1.1 8.03e+06 1.1 0.0e+00 0.0e+00 7.0e+00 8 38 0 0 6 8 38 0 0 6 879
PCApply 35 1.0 3.0367e-03 1.1 8.63e+06 1.1 0.0e+00 0.0e+00 0.0e+00 1 42 0 0 0 1 42 0 0 0 5448
------------------------------------------------------------------------------------------------------------------------
Memory usage is given in bytes:
Object Type Creations Destructions Memory Descendants' Mem.
Reports information only for process 0.
--- Event Stage 0: Main Stage
Vector 43 43 367144 0
Vector Scatter 2 2 2072 0
Index Set 9 9 17528 0
IS L to G Mapping 1 1 564 0
Matrix 4 4 1899988 0
Krylov Solver 2 2 19360 0
Preconditioner 2 2 1784 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 0
Average time for MPI_Barrier(): 1.00136e-06
Average time for zero size MPI_Send(): 1.2517e-05
#PETSc Option Table entries:
-ksp_right_pc
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4
Configure run at: Thu Nov 8 11:21:02 2012
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared-libraries=1 --with-mpi-dir=/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1 --with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_intel_lp64.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_sequential.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_core.so]" --with-lapack-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_lapack95_lp64.a]"
-----------------------------------------
Libraries compiled on Thu Nov 8 11:21:02 2012 on daedalus.ices.utexas.edu
Machine characteristics: Linux-2.6.32-279.1.1.el6.x86_64-x86_64-with-redhat-6.3-Carbon
Using PETSc directory: /opt/apps/ossw/libraries/petsc/petsc-3.3-p2
Using PETSc arch: intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt
-----------------------------------------
Using C compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx -wd1572 -O3 -fPIC ${COPTFLAGS} ${CFLAGS}
Using Fortran compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90 -fPIC -O3 ${FOPTFLAGS} ${FFLAGS}
-----------------------------------------
Using include paths: -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/include
-----------------------------------------
Using C linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx
Using Fortran linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90
Using libraries: -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lHYPRE -lpthread -lsuperlu_dist_3.0 -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.3 -Wl,-rpath,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -L/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,-rpath,/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -L/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -Wl,-rpath,/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -L/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -L/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lirc_s -ldl
-----------------------------------------
----------------------------------------------------------------------------------------------------------------------
| Processor id: 0 |
| Num Processors: 2 |
| Time: Tue Feb 5 13:53:36 2013 |
| OS: Linux |
| HostName: hbar.ices.utexas.edu |
| OS Release: 2.6.32-279.1.1.el6.x86_64 |
| OS Version: #1 SMP Tue Jul 10 11:24:23 CDT 2012 |
| Machine: x86_64 |
| Username: benkirk |
| Configuration: ./configure '--enable-everything' |
| '--prefix=/workspace/libmesh/install' |
| 'CXX=mpicxx' |
| 'CC=mpicc' |
| 'F77=mpif77' |
| 'FC=mpif90' |
| 'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2' |
| 'PETSC_ARCH=intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt' |
| 'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt' |
| 'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/intel-12.1/mpich2-1.4.1p1/mkl-intel-10.3.12.361'|
| 'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/intel-12.1' |
----------------------------------------------------------------------------------------------------------------------
----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.213393, Active time=0.190446 |
----------------------------------------------------------------------------------------------------------------
| Event nCalls Total Time Avg Time Total Time Avg Time % of Active Time |
| w/o Sub w/o Sub With Sub With Sub w/o S With S |
|----------------------------------------------------------------------------------------------------------------|
| |
| |
| DofMap |
| add_neighbors_to_send_list() 1 0.0088 0.008818 0.0124 0.012417 4.63 6.52 |
| build_sparsity() 1 0.0111 0.011117 0.0250 0.025026 5.84 13.14 |
| create_dof_constraints() 1 0.0045 0.004546 0.0045 0.004546 2.39 2.39 |
| distribute_dofs() 1 0.0068 0.006791 0.0164 0.016417 3.57 8.62 |
| dof_indices() 480 0.0544 0.000113 0.0544 0.000113 28.57 28.57 |
| prepare_send_list() 1 0.0001 0.000102 0.0001 0.000102 0.05 0.05 |
| reinit() 1 0.0094 0.009355 0.0094 0.009355 4.91 4.91 |
| |
| EquationSystems |
| build_solution_vector() 2 0.0014 0.000700 0.0275 0.013763 0.73 14.45 |
| |
| ExodusII_IO |
| write_nodal_data() 1 0.0028 0.002818 0.0028 0.002818 1.48 1.48 |
| |
| FE |
| compute_shape_functions() 142 0.0119 0.000084 0.0119 0.000084 6.25 6.25 |
| init_shape_functions() 31 0.0005 0.000015 0.0005 0.000015 0.24 0.24 |
| inverse_map() 90 0.0008 0.000008 0.0008 0.000008 0.40 0.40 |
| |
| FEMap |
| compute_affine_map() 142 0.0013 0.000009 0.0013 0.000009 0.69 0.69 |
| compute_face_map() 30 0.0005 0.000018 0.0013 0.000044 0.28 0.69 |
| init_face_shape_functions() 1 0.0000 0.000012 0.0000 0.000012 0.01 0.01 |
| init_reference_to_physical_map() 31 0.0008 0.000025 0.0008 0.000025 0.41 0.41 |
| |
| GMVIO |
| write_nodal_data() 1 0.0034 0.003416 0.0035 0.003465 1.79 1.82 |
| |
| Mesh |
| find_neighbors() 1 0.0029 0.002941 0.0031 0.003094 1.54 1.62 |
| renumber_nodes_and_elem() 2 0.0003 0.000161 0.0003 0.000161 0.17 0.17 |
| |
| MeshCommunication |
| compute_hilbert_indices() 2 0.0021 0.001053 0.0021 0.001053 1.11 1.11 |
| find_global_indices() 2 0.0009 0.000475 0.0041 0.002064 0.50 2.17 |
| parallel_sort() 2 0.0006 0.000323 0.0008 0.000376 0.34 0.39 |
| |
| MeshOutput |
| write_equation_systems() 2 0.0001 0.000063 0.0340 0.017019 0.07 17.87 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0014 0.001427 0.0014 0.001427 0.75 0.75 |
| |
| MetisPartitioner |
| partition() 1 0.0031 0.003109 0.0049 0.004851 1.63 2.55 |
| |
| Parallel |
| allgather() 9 0.0002 0.000025 0.0002 0.000027 0.12 0.13 |
| max(bool) 1 0.0000 0.000004 0.0000 0.000004 0.00 0.00 |
| max(scalar) 124 0.0004 0.000003 0.0004 0.000003 0.19 0.19 |
| max(vector) 28 0.0002 0.000006 0.0004 0.000013 0.09 0.19 |
| min(bool) 144 0.0003 0.000002 0.0003 0.000002 0.17 0.17 |
| min(scalar) 118 0.0019 0.000016 0.0019 0.000016 0.98 0.98 |
| min(vector) 28 0.0002 0.000008 0.0005 0.000016 0.11 0.24 |
| probe() 12 0.0002 0.000020 0.0002 0.000020 0.13 0.13 |
| receive() 12 0.0001 0.000008 0.0003 0.000028 0.05 0.18 |
| send() 12 0.0000 0.000004 0.0000 0.000004 0.03 0.03 |
| send_receive() 16 0.0002 0.000015 0.0007 0.000042 0.13 0.35 |
| sum() 23 0.0002 0.000009 0.0008 0.000033 0.11 0.39 |
| |
| Parallel::Request |
| wait() 12 0.0000 0.000003 0.0000 0.000003 0.02 0.02 |
| |
| Partitioner |
| set_node_processor_ids() 1 0.0006 0.000650 0.0008 0.000754 0.34 0.40 |
| set_parent_processor_ids() 1 0.0002 0.000208 0.0002 0.000208 0.11 0.11 |
| |
| PetscLinearSolver |
| solve() 1 0.0313 0.031284 0.0313 0.031284 16.43 16.43 |
| |
| System |
| assemble() 1 0.0241 0.024093 0.0528 0.052755 12.65 27.70 |
----------------------------------------------------------------------------------------------------------------
| Totals: 1513 0.1904 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example vector_fe_ex1:
* mpirun -np 2 example-devel -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
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
C++ Includes