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;
/**
*
*/
Real exact_solution (const Real x,
const Real y,
const Real t)
{
const Real xo = 0.2;
const Real yo = 0.2;
const Real u = 0.8;
const Real v = 0.8;
const Real num =
pow(x - u*t - xo, 2.) +
pow(y - v*t - yo, 2.);
const Real den =
0.01*(4.*t + 1.);
return exp(-num/den)/(4.*t + 1.);
}
The source file transient_ex1.C with comments:
Transient Example 1 - Solving a Transient Linear System in Parallel
This example shows how a simple, linear transient system can be solved in parallel. The system is simple scalar convection-diffusion with a specified external velocity. The initial condition is given, and the solution is advanced in time with a standard Crank-Nicolson time-stepping strategy.
C++ include files that we need
#include <iostream>
#include <algorithm>
#include <sstream>
#include <math.h>
Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/gmv_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
This example will solve a linear transient system,
so we need to include the \p TransientLinearImplicitSystem definition.
#include "libmesh/linear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/vector_value.h"
The definition of a geometric element
#include "libmesh/elem.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Function prototype. This function will assemble the system
matrix and right-hand-side at each time step. Note that
since the system is linear we technically do not need to
assmeble the matrix at each time step, but we will anyway.
In subsequent examples we will employ adaptive mesh refinement,
and with a changing mesh it will be necessary to rebuild the
system matrix.
void assemble_cd (EquationSystems& es,
const std::string& system_name);
Function prototype. This function will initialize the system.
Initialization functions are optional for systems. They allow
you to specify the initial values of the solution. If an
initialization function is not provided then the default (0)
solution is provided.
void init_cd (EquationSystems& es,
const std::string& system_name);
Exact solution function prototype. This gives the exact
solution as a function of space and time. In this case the
initial condition will be taken as the exact solution at time 0,
as will the Dirichlet boundary conditions at time t.
Real exact_solution (const Real x,
const Real y,
const Real t);
Number exact_value (const Point& p,
const Parameters& parameters,
const std::string&,
const std::string&)
{
return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
}
We can now begin the main program. Note that this
example will fail if you are using complex numbers
since it was designed to be run only with real numbers.
int main (int argc, char** argv)
{
Initialize libMesh.
LibMeshInit init (argc, argv);
This example requires Adaptive Mesh Refinement support - although
it only refines uniformly, the refinement code used is the same
underneath
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
Skip this 2D example if libMesh was compiled as 1D-only.
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Read the mesh from file. This is the coarse mesh that will be used
in example 10 to demonstrate adaptive mesh refinement. Here we will
simply read it in and uniformly refine it 5 times before we compute
with it.
Mesh mesh;
mesh.read ("mesh.xda");
Create a MeshRefinement object to handle refinement of our mesh.
This class handles all the details of mesh refinement and coarsening.
MeshRefinement mesh_refinement (mesh);
Uniformly refine the mesh 5 times. This is the
first time we use the mesh refinement capabilities
of the library.
mesh_refinement.uniformly_refine (5);
Print information about the mesh to the screen.
mesh.print_info();
Create an equation systems object.
EquationSystems equation_systems (mesh);
Add a transient system to the EquationSystems
object named "Convection-Diffusion".
TransientLinearImplicitSystem & system =
equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
Adds the variable "u" to "Convection-Diffusion". "u"
will be approximated using first-order approximation.
system.add_variable ("u", FIRST);
Give the system a pointer to the matrix assembly
and initialization functions.
system.attach_assemble_function (assemble_cd);
system.attach_init_function (init_cd);
Initialize the data structures for the equation system.
equation_systems.init ();
Prints information about the system to the screen.
equation_systems.print_info();
Write out the initial conditions.
GMVIO(mesh).write_equation_systems ("out_000.gmv",
equation_systems);
The Convection-Diffusion system requires that we specify
the flow velocity. We will specify it as a RealVectorValue
data type and then use the Parameters object to pass it to
the assemble function.
equation_systems.parameters.set<RealVectorValue>("velocity") =
RealVectorValue (0.8, 0.8);
Solve the system "Convection-Diffusion". This will be done by
looping over the specified time interval and calling the
solve() member at each time step. This will assemble the
system and call the linear solver.
const Real dt = 0.025;
system.time = 0.;
for (unsigned int t_step = 0; t_step < 50; t_step++)
{
Incremenet the time counter, set the time and the
time step size as parameters in the EquationSystem.
system.time += dt;
equation_systems.parameters.set<Real> ("time") = system.time;
equation_systems.parameters.set<Real> ("dt") = dt;
A pretty update message
std::cout << " Solving time step ";
Since some compilers fail to offer full stream
functionality, libMesh offers a string stream
to work around this. Note that for other compilers,
this is just a set of preprocessor macros and therefore
should cost nothing (compared to a hand-coded string stream).
We use additional curly braces here simply to enforce data
locality.
{
std::ostringstream out;
out << std::setw(2)
<< std::right
<< t_step
<< ", time="
<< std::fixed
<< std::setw(6)
<< std::setprecision(3)
<< std::setfill('0')
<< std::left
<< system.time
<< "...";
std::cout << out.str() << std::endl;
}
At this point we need to update the old
solution vector. The old solution vector
will be the current solution vector from the
previous time step. We will do this by extracting the
system from the \p EquationSystems object and using
vector assignment. Since only \p TransientSystems
(and systems derived from them) contain old solutions
we need to specify the system type when we ask for it.
*system.old_local_solution = *system.current_local_solution;
Assemble & solve the linear system
equation_systems.get_system("Convection-Diffusion").solve();
Output evey 10 timesteps to file.
if ( (t_step+1)%10 == 0)
{
std::ostringstream file_name;
file_name << "out_"
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step+1
<< ".gmv";
GMVIO(mesh).write_equation_systems (file_name.str(),
equation_systems);
}
}
#endif // #ifdef LIBMESH_ENABLE_AMR
All done.
return 0;
}
We now define the function which provides the
initialization routines for the "Convection-Diffusion"
system. This handles things like setting initial
conditions and boundary conditions.
void init_cd (EquationSystems& es,
const std::string& system_name)
{
It is a good idea to make sure we are initializing
the proper system.
libmesh_assert_equal_to (system_name, "Convection-Diffusion");
Get a reference to the Convection-Diffusion system object.
TransientLinearImplicitSystem & system =
es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
Project initial conditions at time 0
es.parameters.set<Real> ("time") = system.time = 0;
system.project_solution(exact_value, NULL, es.parameters);
}
Now we define the assemble function which will be used
by the EquationSystems object at each timestep to assemble
the linear system for solution.
void assemble_cd (EquationSystems& es,
const std::string& system_name)
{
#ifdef LIBMESH_ENABLE_AMR
It is a good idea to make sure we are assembling
the proper system.
libmesh_assert_equal_to (system_name, "Convection-Diffusion");
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 Convection-Diffusion system object.
TransientLinearImplicitSystem & system =
es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
Get a constant reference to the Finite Element type
for the first (and only) variable in the system.
FEType fe_type = system.variable_type(0);
Build a Finite Element object of the specified type. Since the
\p FEBase::build() member dynamically creates memory we will
store the object as an \p AutoPtr. This can be thought
of as a pointer that will clean up after itself.
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
A Gauss quadrature rule for numerical integration.
Let the \p FEType object decide what order rule is appropriate.
QGauss qrule (dim, fe_type.default_quadrature_order());
QGauss qface (dim-1, fe_type.default_quadrature_order());
Tell the finite element object to use our quadrature rule.
fe->attach_quadrature_rule (&qrule);
fe_face->attach_quadrature_rule (&qface);
Here we define some references to cell-specific data that
will be used to assemble the linear system. We will start
with the element Jacobian * quadrature weight at each integration point.
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
The element shape function gradients evaluated at the quadrature
points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
The XY locations of the quadrature points used for face integration
const std::vector<Point>& qface_points = fe_face->get_xyz();
A reference to the \p DofMap object for this system. The \p DofMap
object handles the index translation from node and element numbers
to degree of freedom numbers. We will talk more about the \p DofMap
in future examples.
const DofMap& dof_map = system.get_dof_map();
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".
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;
Here we extract the velocity & parameters that we put in the
EquationSystems object.
const RealVectorValue velocity =
es.parameters.get<RealVectorValue> ("velocity");
const Real dt = es.parameters.get<Real> ("dt");
Now we will loop over all the elements in the mesh that
live on the local processor. We will compute the element
matrix and right-hand-side contribution. Since the mesh
will be refined we want to 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();
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).
Ke.resize (dof_indices.size(),
dof_indices.size());
Fe.resize (dof_indices.size());
Now we will build the element matrix and right-hand-side.
Constructing the RHS requires the solution and its
gradient from the previous timestep. This myst be
calculated at each quadrature point by summing the
solution degree-of-freedom values by the appropriate
weight functions.
for (unsigned int qp=0; qp<qrule.n_points(); qp++)
{
Values to hold the old solution & its gradient.
Number u_old = 0.;
Gradient grad_u_old;
Compute the old solution & its gradient.
for (unsigned int l=0; l<phi.size(); l++)
{
u_old += phi[l][qp]*system.old_solution (dof_indices[l]);
This will work,
grad_u_old += dphi[l][qp]*system.old_solution (dof_indices[l]);
but we can do it without creating a temporary like this:
grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
}
Now compute the element matrix and RHS contributions.
for (unsigned int i=0; i<phi.size(); i++)
{
The RHS contribution
Fe(i) += JxW[qp]*(
Mass matrix term
u_old*phi[i][qp] +
-.5*dt*(
Convection term
(grad_u_old may be complex, so the
order here is important!)
(grad_u_old*velocity)*phi[i][qp] +
Diffusion term
0.01*(grad_u_old*dphi[i][qp]))
);
for (unsigned int j=0; j<phi.size(); j++)
{
The matrix contribution
Ke(i,j) += JxW[qp]*(
Mass-matrix
phi[i][qp]*phi[j][qp] +
.5*dt*(
Convection term
(velocity*dphi[j][qp])*phi[i][qp] +
Diffusion term
0.01*(dphi[i][qp]*dphi[j][qp]))
);
}
}
}
At this point 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 imposed
via the penalty method.
The following loops 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.
The following loops 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.
{
The penalty value.
const Real penalty = 1.e10;
The following loops 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 s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
fe_face->reinit(elem,s);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
const Number value = exact_solution (qface_points[qp](0),
qface_points[qp](1),
system.time);
RHS contribution
for (unsigned int i=0; i<psi.size(); i++)
Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
Matrix contribution
for (unsigned int i=0; i<psi.size(); i++)
for (unsigned int j=0; j<psi.size(); j++)
Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
}
}
}
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 \p SparseMatrix::add_matrix()
and \p NumericVector::add_vector() members do this for us.
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
That concludes the system matrix assembly routine.
#endif // #ifdef LIBMESH_ENABLE_AMR
}
The source file exact_solution.C without comments:
#include <math.h>
#include "libmesh/libmesh_common.h"
using namespace libMesh;
/**
*
*/
Real exact_solution (const Real x,
const Real y,
const Real t)
{
const Real xo = 0.2;
const Real yo = 0.2;
const Real u = 0.8;
const Real v = 0.8;
const Real num =
pow(x - u*t - xo, 2.) +
pow(y - v*t - yo, 2.);
const Real den =
0.01*(4.*t + 1.);
return exp(-num/den)/(4.*t + 1.);
}
The source file transient_ex1.C without comments:
#include <iostream>
#include <algorithm>
#include <sstream>
#include <math.h>
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/gmv_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/fe.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/dof_map.h"
#include "libmesh/sparse_matrix.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_vector.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/vector_value.h"
#include "libmesh/elem.h"
using namespace libMesh;
void assemble_cd (EquationSystems& es,
const std::string& system_name);
void init_cd (EquationSystems& es,
const std::string& system_name);
Real exact_solution (const Real x,
const Real y,
const Real t);
Number exact_value (const Point& p,
const Parameters& parameters,
const std::string&,
const std::string&)
{
return exact_solution(p(0), p(1), parameters.get<Real> ("time"));
}
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Mesh mesh;
mesh.read ("mesh.xda");
MeshRefinement mesh_refinement (mesh);
mesh_refinement.uniformly_refine (5);
mesh.print_info();
EquationSystems equation_systems (mesh);
TransientLinearImplicitSystem & system =
equation_systems.add_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
system.add_variable ("u", FIRST);
system.attach_assemble_function (assemble_cd);
system.attach_init_function (init_cd);
equation_systems.init ();
equation_systems.print_info();
GMVIO(mesh).write_equation_systems ("out_000.gmv",
equation_systems);
equation_systems.parameters.set<RealVectorValue>("velocity") =
RealVectorValue (0.8, 0.8);
const Real dt = 0.025;
system.time = 0.;
for (unsigned int t_step = 0; t_step < 50; t_step++)
{
system.time += dt;
equation_systems.parameters.set<Real> ("time") = system.time;
equation_systems.parameters.set<Real> ("dt") = dt;
std::cout << " Solving time step ";
{
std::ostringstream out;
out << std::setw(2)
<< std::right
<< t_step
<< ", time="
<< std::fixed
<< std::setw(6)
<< std::setprecision(3)
<< std::setfill('0')
<< std::left
<< system.time
<< "...";
std::cout << out.str() << std::endl;
}
*system.old_local_solution = *system.current_local_solution;
equation_systems.get_system("Convection-Diffusion").solve();
if ( (t_step+1)%10 == 0)
{
std::ostringstream file_name;
file_name << "out_"
<< std::setw(3)
<< std::setfill('0')
<< std::right
<< t_step+1
<< ".gmv";
GMVIO(mesh).write_equation_systems (file_name.str(),
equation_systems);
}
}
#endif // #ifdef LIBMESH_ENABLE_AMR
return 0;
}
void init_cd (EquationSystems& es,
const std::string& system_name)
{
libmesh_assert_equal_to (system_name, "Convection-Diffusion");
TransientLinearImplicitSystem & system =
es.get_system<TransientLinearImplicitSystem>("Convection-Diffusion");
es.parameters.set<Real> ("time") = system.time = 0;
system.project_solution(exact_value, NULL, es.parameters);
}
void assemble_cd (EquationSystems& es,
const std::string& system_name)
{
#ifdef LIBMESH_ENABLE_AMR
libmesh_assert_equal_to (system_name, "Convection-Diffusion");
const MeshBase& mesh = es.get_mesh();
const unsigned int dim = mesh.mesh_dimension();
TransientLinearImplicitSystem & system =
es.get_system<TransientLinearImplicitSystem> ("Convection-Diffusion");
FEType fe_type = system.variable_type(0);
AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
QGauss qrule (dim, fe_type.default_quadrature_order());
QGauss qface (dim-1, fe_type.default_quadrature_order());
fe->attach_quadrature_rule (&qrule);
fe_face->attach_quadrature_rule (&qface);
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<Real>& JxW_face = fe_face->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<Real> >& psi = fe_face->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
const std::vector<Point>& qface_points = fe_face->get_xyz();
const DofMap& dof_map = system.get_dof_map();
DenseMatrix<Number> Ke;
DenseVector<Number> Fe;
std::vector<dof_id_type> dof_indices;
const RealVectorValue velocity =
es.parameters.get<RealVectorValue> ("velocity");
const Real dt = es.parameters.get<Real> ("dt");
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++)
{
Number u_old = 0.;
Gradient grad_u_old;
for (unsigned int l=0; l<phi.size(); l++)
{
u_old += phi[l][qp]*system.old_solution (dof_indices[l]);
grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l]));
}
for (unsigned int i=0; i<phi.size(); i++)
{
Fe(i) += JxW[qp]*(
u_old*phi[i][qp] +
-.5*dt*(
(grad_u_old*velocity)*phi[i][qp] +
0.01*(grad_u_old*dphi[i][qp]))
);
for (unsigned int j=0; j<phi.size(); j++)
{
Ke(i,j) += JxW[qp]*(
phi[i][qp]*phi[j][qp] +
.5*dt*(
(velocity*dphi[j][qp])*phi[i][qp] +
0.01*(dphi[i][qp]*dphi[j][qp]))
);
}
}
}
{
const Real penalty = 1.e10;
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
fe_face->reinit(elem,s);
for (unsigned int qp=0; qp<qface.n_points(); qp++)
{
const Number value = exact_solution (qface_points[qp](0),
qface_points[qp](1),
system.time);
for (unsigned int i=0; i<psi.size(); i++)
Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp];
for (unsigned int i=0; i<psi.size(); i++)
for (unsigned int j=0; j<psi.size(); j++)
Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp];
}
}
}
dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
system.matrix->add_matrix (Ke, dof_indices);
system.rhs->add_vector (Fe, dof_indices);
}
#endif // #ifdef LIBMESH_ENABLE_AMR
}
The console output of the program:
***************************************************************
* Running Example transient_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
***************************************************************
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=6273
n_local_nodes()=3183
n_elem()=13650
n_local_elem()=6872
n_active_elem()=10240
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "Convection-Diffusion"
Type "TransientLinearImplicit"
Variables="u"
Finite Element Types="LAGRANGE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="FIRST", "THIRD"
n_dofs()=6273
n_local_dofs()=3183
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=3
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 7.54918
Average Off-Processor Bandwidth <= 0.0596206
Maximum On-Processor Bandwidth <= 11
Maximum Off-Processor Bandwidth <= 5
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Solving time step 0, time=0.0250...
Solving time step 1, time=0.0500...
Solving time step 2, time=0.0750...
Solving time step 3, time=0.1000...
Solving time step 4, time=0.1250...
Solving time step 5, time=0.1500...
Solving time step 6, time=0.1750...
Solving time step 7, time=0.2000...
Solving time step 8, time=0.2250...
Solving time step 9, time=0.2500...
Solving time step 10, time=0.2750...
Solving time step 11, time=0.3000...
Solving time step 12, time=0.3250...
Solving time step 13, time=0.3500...
Solving time step 14, time=0.3750...
Solving time step 15, time=0.4000...
Solving time step 16, time=0.4250...
Solving time step 17, time=0.4500...
Solving time step 18, time=0.4750...
Solving time step 19, time=0.5000...
Solving time step 20, time=0.5250...
Solving time step 21, time=0.5500...
Solving time step 22, time=0.5750...
Solving time step 23, time=0.6000...
Solving time step 24, time=0.6250...
Solving time step 25, time=0.6500...
Solving time step 26, time=0.6750...
Solving time step 27, time=0.7000...
Solving time step 28, time=0.7250...
Solving time step 29, time=0.7500...
Solving time step 30, time=0.7750...
Solving time step 31, time=0.8000...
Solving time step 32, time=0.8250...
Solving time step 33, time=0.8500...
Solving time step 34, time=0.8750...
Solving time step 35, time=0.9000...
Solving time step 36, time=0.9250...
Solving time step 37, time=0.9500...
Solving time step 38, time=0.9750...
Solving time step 39, time=1.0000...
Solving time step 40, time=1.0250...
Solving time step 41, time=1.0500...
Solving time step 42, time=1.0750...
Solving time step 43, time=1.1000...
Solving time step 44, time=1.1250...
Solving time step 45, time=1.1500...
Solving time step 46, time=1.1750...
Solving time step 47, time=1.2000...
Solving time step 48, time=1.2250...
Solving time step 49, time=1.2500...
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/transient/transient_ex1/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:52:10 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): 1.507e+01 1.00000 1.507e+01
Objects: 3.120e+02 1.00000 3.120e+02
Flops: 2.106e+08 1.01451 2.091e+08 4.182e+08
Flops/sec: 1.397e+07 1.01451 1.387e+07 2.775e+07
MPI Messages: 7.830e+02 1.00000 7.830e+02 1.566e+03
MPI Message Lengths: 7.890e+05 1.00000 1.008e+03 1.578e+06
MPI Reductions: 2.205e+03 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: 1.5073e+01 100.0% 4.1825e+08 100.0% 1.566e+03 100.0% 1.008e+03 100.0% 2.204e+03 100.0%
------------------------------------------------------------------------------------------------------------------------
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 495 1.0 5.9774e-03 1.2 1.90e+07 1.0 0.0e+00 0.0e+00 5.0e+02 0 9 0 0 22 0 9 0 0 22 6266
VecNorm 595 1.0 3.2074e-03 1.5 3.79e+06 1.0 0.0e+00 0.0e+00 6.0e+02 0 2 0 0 27 0 2 0 0 27 2327
VecScale 545 1.0 5.6815e-04 1.0 1.73e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 1 0 0 0 0 1 0 0 0 6017
VecCopy 151 1.0 3.4285e-04 1.0 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 703 1.0 1.1609e-03 1.0 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 100 1.0 6.9690e-03 1.0 6.37e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 180
VecMAXPY 545 1.0 6.0866e-03 1.0 2.22e+07 1.0 0.0e+00 0.0e+00 0.0e+00 0 10 0 0 0 0 10 0 0 0 7175
VecAssemblyBegin 202 1.0 1.1765e-0159.8 0.00e+00 0.0 1.0e+02 1.2e+03 6.1e+02 0 0 6 8 27 0 0 6 8 27 0
VecAssemblyEnd 202 1.0 1.7405e-04 1.3 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 646 1.0 1.0023e-03 1.0 0.00e+00 0.0 1.3e+03 7.0e+02 0.0e+00 0 0 83 57 0 0 0 83 57 0 0
VecScatterEnd 646 1.0 1.1274e-0211.7 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 545 1.0 3.7513e-03 1.4 5.20e+06 1.0 0.0e+00 0.0e+00 5.4e+02 0 2 0 0 25 0 2 0 0 25 2734
MatMult 545 1.0 2.9635e-02 1.6 2.46e+07 1.0 1.1e+03 6.5e+02 0.0e+00 0 12 70 45 0 0 12 70 45 0 1631
MatSolve 595 1.0 4.3775e-02 1.0 8.98e+07 1.0 0.0e+00 0.0e+00 0.0e+00 0 43 0 0 0 0 43 0 0 0 4069
MatLUFactorNum 50 1.0 6.6998e-02 1.1 4.95e+07 1.0 0.0e+00 0.0e+00 0.0e+00 0 24 0 0 0 0 24 0 0 0 1470
MatILUFactorSym 50 1.0 1.5487e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 1.5e+02 1 0 0 0 7 1 0 0 0 7 0
MatAssemblyBegin 100 1.0 1.0072e-0163.6 0.00e+00 0.0 1.5e+02 3.6e+03 2.0e+02 0 0 10 35 9 0 0 10 35 9 0
MatAssemblyEnd 100 1.0 6.3522e-03 1.3 0.00e+00 0.0 4.0e+00 1.6e+02 8.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetRowIJ 50 1.0 1.1921e-05 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
MatGetOrdering 50 1.0 1.2894e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 1.0e+02 0 0 0 0 5 0 0 0 0 5 0
MatZeroEntries 52 1.0 4.6635e-04 1.0 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 495 1.0 1.1497e-02 1.1 3.80e+07 1.0 0.0e+00 0.0e+00 5.0e+02 0 18 0 0 22 0 18 0 0 22 6516
KSPSetUp 100 1.0 3.4571e-04 1.0 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 50 1.0 3.2147e-01 1.0 2.11e+08 1.0 1.1e+03 6.5e+02 1.3e+03 2100 70 45 61 2100 70 45 61 1301
PCSetUp 100 1.0 2.3050e-01 1.1 4.95e+07 1.0 0.0e+00 0.0e+00 2.5e+02 1 24 0 0 11 1 24 0 0 11 427
PCSetUpOnBlocks 50 1.0 2.3004e-01 1.1 4.95e+07 1.0 0.0e+00 0.0e+00 2.5e+02 1 24 0 0 11 1 24 0 0 11 428
PCApply 595 1.0 4.8267e-02 1.0 8.98e+07 1.0 0.0e+00 0.0e+00 0.0e+00 0 43 0 0 0 0 43 0 0 0 3690
------------------------------------------------------------------------------------------------------------------------
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 81 81 1957960 0
Vector Scatter 6 6 6216 0
Index Set 162 162 759024 0
IS L to G Mapping 5 5 2820 0
Matrix 53 53 46757060 0
Krylov Solver 2 2 19360 0
Preconditioner 2 2 1784 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 1.19209e-07
Average time for MPI_Barrier(): 1.43051e-06
Average time for zero size MPI_Send(): 1.29938e-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:52:10 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=15.0881, Active time=14.8845 |
----------------------------------------------------------------------------------------------------------------
| 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.0927 0.092685 0.0974 0.097403 0.62 0.65 |
| build_sparsity() 1 0.0537 0.053685 0.1702 0.170213 0.36 1.14 |
| create_dof_constraints() 1 0.0230 0.022992 0.0230 0.022992 0.15 0.15 |
| distribute_dofs() 1 0.0812 0.081248 0.2229 0.222853 0.55 1.50 |
| dof_indices() 296932 6.7468 0.000023 6.7468 0.000023 45.33 45.33 |
| prepare_send_list() 1 0.0001 0.000070 0.0001 0.000070 0.00 0.00 |
| reinit() 1 0.1386 0.138646 0.1386 0.138646 0.93 0.93 |
| |
| EquationSystems |
| build_solution_vector() 6 0.1092 0.018196 0.7874 0.131226 0.73 5.29 |
| |
| FE |
| compute_shape_functions() 262400 1.4125 0.000005 1.4125 0.000005 9.49 9.49 |
| init_shape_functions() 6700 0.0161 0.000002 0.0161 0.000002 0.11 0.11 |
| inverse_map() 13200 0.0533 0.000004 0.0533 0.000004 0.36 0.36 |
| |
| FEMap |
| compute_affine_map() 262400 1.2738 0.000005 1.2738 0.000005 8.56 8.56 |
| compute_face_map() 6600 0.0563 0.000009 0.1114 0.000017 0.38 0.75 |
| init_face_shape_functions() 100 0.0005 0.000005 0.0005 0.000005 0.00 0.00 |
| init_reference_to_physical_map() 6700 0.0410 0.000006 0.0410 0.000006 0.28 0.28 |
| |
| GMVIO |
| write_nodal_data() 6 0.1583 0.026382 0.1586 0.026430 1.06 1.07 |
| |
| LocationMap |
| find() 32736 0.0690 0.000002 0.0690 0.000002 0.46 0.46 |
| init() 5 0.0016 0.000318 0.0016 0.000318 0.01 0.01 |
| |
| Mesh |
| find_neighbors() 2 0.1417 0.070859 0.1418 0.070917 0.95 0.95 |
| |
| MeshCommunication |
| broadcast() 1 0.0012 0.001205 0.0020 0.001969 0.01 0.01 |
| compute_hilbert_indices() 3 0.0404 0.013465 0.0404 0.013465 0.27 0.27 |
| find_global_indices() 3 0.0167 0.005561 0.0628 0.020938 0.11 0.42 |
| parallel_sort() 3 0.0031 0.001050 0.0046 0.001541 0.02 0.03 |
| |
| MeshOutput |
| write_equation_systems() 6 0.0003 0.000057 0.9466 0.157763 0.00 6.36 |
| |
| MeshRefinement |
| _refine_elements() 5 0.1188 0.023767 0.2695 0.053909 0.80 1.81 |
| add_point() 32736 0.0713 0.000002 0.1445 0.000004 0.48 0.97 |
| |
| MetisPartitioner |
| partition() 2 0.1662 0.083116 0.2282 0.114099 1.12 1.53 |
| |
| Parallel |
| allgather() 11 0.0023 0.000209 0.0023 0.000211 0.02 0.02 |
| broadcast() 9 0.0006 0.000065 0.0005 0.000056 0.00 0.00 |
| max(bool) 6 0.0019 0.000323 0.0019 0.000323 0.01 0.01 |
| max(scalar) 272 0.0008 0.000003 0.0008 0.000003 0.01 0.01 |
| max(vector) 61 0.0004 0.000007 0.0009 0.000016 0.00 0.01 |
| min(bool) 322 0.0008 0.000002 0.0008 0.000002 0.01 0.01 |
| min(scalar) 264 0.0068 0.000026 0.0068 0.000026 0.05 0.05 |
| min(vector) 61 0.0005 0.000008 0.0010 0.000017 0.00 0.01 |
| probe() 16 0.0005 0.000031 0.0005 0.000031 0.00 0.00 |
| receive() 16 0.0002 0.000014 0.0007 0.000045 0.00 0.00 |
| send() 16 0.0001 0.000005 0.0001 0.000005 0.00 0.00 |
| send_receive() 22 0.0004 0.000016 0.0012 0.000055 0.00 0.01 |
| sum() 39 0.0020 0.000052 0.0024 0.000061 0.01 0.02 |
| |
| Parallel::Request |
| wait() 16 0.0000 0.000003 0.0000 0.000003 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 2 0.0102 0.005116 0.0112 0.005577 0.07 0.07 |
| set_parent_processor_ids() 2 0.0136 0.006811 0.0136 0.006811 0.09 0.09 |
| |
| PetscLinearSolver |
| solve() 50 0.3325 0.006649 0.3325 0.006649 2.23 2.23 |
| |
| System |
| assemble() 50 3.5807 0.071615 12.3876 0.247752 24.06 83.22 |
| project_vector() 1 0.0425 0.042453 0.1573 0.157297 0.29 1.06 |
----------------------------------------------------------------------------------------------------------------
| Totals: 921788 14.8845 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example transient_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