The source file adaptivity_ex5.C with comments:
#include <iostream>
#include <algorithm>
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath>
Basic include file needed for the mesh functionality.
#include "libmesh/libmesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/gmv_io.h"
#include "libmesh/exodusII_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/periodic_boundaries.h"
#include "libmesh/periodic_boundary.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/parsed_function.h"
#include "libmesh/getpot.h"
This example will solve a linear transient system,
so we need to include the \p TransientLinearImplicitSystem definition.
#include "libmesh/transient_system.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/vector_value.h"
To refine the mesh we need an \p ErrorEstimator
object to figure out which elements to refine.
#include "libmesh/error_vector.h"
#include "libmesh/kelly_error_estimator.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"));
}
With --enable-fparser, the user can also optionally set their own
exact solution equations.
FunctionBase<Number>* parsed_solution = NULL;
Returns a string with 'number' formatted and placed directly
into the string in some way
std::string exodus_filename(unsigned number);
Begin the main program. Note that the first
statement in the program throws an error if
you are in complex number mode, since this
example is only intended to work with real
numbers.
int main (int argc, char** argv)
{
Initialize libMesh.
LibMeshInit init (argc, argv);
#if !defined(LIBMESH_ENABLE_AMR)
libmesh_example_assert(false, "--enable-amr");
#elif !defined(LIBMESH_HAVE_XDR)
We use XDR support in our output here
libmesh_example_assert(false, "--enable-xdr");
#elif !defined(LIBMESH_ENABLE_PERIODIC)
libmesh_example_assert(false, "--enable-periodic");
#else
Our Trilinos interface does not yet support adaptive transient
problems
libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
Brief message to the user regarding the program name
and command line arguments.
Use commandline parameter to specify if we are to read in an initial solution or generate it ourself
Use commandline parameter to specify if we are to read in an initial solution or generate it ourself
std::cout << "Usage:\n"
<<"\t " << argv[0] << " -init_timestep 0\n"
<< "OR\n"
<<"\t " << argv[0] << " -read_solution -init_timestep 26\n"
<< std::endl;
std::cout << "Running: " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
Create a GetPot object to parse the command line
GetPot command_line (argc, argv);
This boolean value is obtained from the command line, it is true
if the flag "-read_solution" is present, false otherwise.
It indicates whether we are going to read in
the mesh and solution files "saved_mesh.xda" and "saved_solution.xda"
or whether we are going to start from scratch by just reading
"mesh.xda"
const bool read_solution = command_line.search("-read_solution");
This value is also obtained from the commandline and it specifies the
initial value for the t_step looping variable. We must
distinguish between the two cases here, whether we read in the
solution or we started from scratch, so that we do not overwrite the
gmv output files.
unsigned int init_timestep = 0;
Search the command line for the "init_timestep" flag and if it is
present, set init_timestep accordingly.
if(command_line.search("-init_timestep"))
init_timestep = command_line.next(0);
else
{
if (libMesh::processor_id() == 0)
std::cerr << "ERROR: Initial timestep not specified\n" << std::endl;
This handy function will print the file name, line number,
and then abort. Currrently the library does not use C++
exception handling.
libmesh_error();
}
This value is also obtained from the command line, and specifies
the number of time steps to take.
unsigned int n_timesteps = 0;
Again do a search on the command line for the argument
if(command_line.search("-n_timesteps"))
n_timesteps = command_line.next(0);
else
{
std::cout << "ERROR: Number of timesteps not specified\n" << std::endl;
libmesh_error();
}
The user can specify a different exact solution on the command
line, if we have an expression parser compiled in
#ifdef LIBMESH_HAVE_FPARSER
const bool have_expression = command_line.search("-exact_solution");
#else
const bool have_expression = false;
#endif
if (have_expression)
parsed_solution = new ParsedFunction<Number>(command_line.next(std::string()));
Skip this 2D example if libMesh was compiled as 1D-only.
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Create a new mesh.
ParallelMesh doesn't yet understand periodic BCs, plus
we still need some work on automatic parallel restarts
SerialMesh mesh;
Create an equation systems object.
EquationSystems equation_systems (mesh);
MeshRefinement mesh_refinement (mesh);
First we process the case where we do not read in the solution
if(!read_solution)
{
MeshTools::Generation::build_square(mesh, 2, 2, 0., 2., 0., 2., QUAD4);
Again do a search on the command line for an argument
unsigned int n_refinements = 5;
if(command_line.search("-n_refinements"))
n_refinements = command_line.next(0);
Uniformly refine the mesh 5 times
if(!read_solution)
mesh_refinement.uniformly_refine (n_refinements);
Print information about the mesh to the screen.
mesh.print_info();
Declare the system and its variables.
Begin by creating a transient system
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 initialization function.
system.attach_init_function (init_cd);
}
Otherwise we read in the solution and mesh
else
{
Read in the mesh stored in "saved_mesh.xda"
mesh.read("saved_mesh.xdr");
Print information about the mesh to the screen.
mesh.print_info();
Read in the solution stored in "saved_solution.xda"
equation_systems.read("saved_solution.xdr", libMeshEnums::DECODE);
}
Get a reference to the system so that we can attach things to it
TransientLinearImplicitSystem & system =
equation_systems.get_system<TransientLinearImplicitSystem>
("Convection-Diffusion");
Give the system a pointer to the assembly function.
system.attach_assemble_function (assemble_cd);
Creating and attaching Periodic Boundaries
DofMap & dof_map = system.get_dof_map();
Create a boundary periodic with one displaced 2.0 in the x
direction
PeriodicBoundary horz(RealVectorValue(2.0, 0., 0.));
Connect boundary ids 3 and 1 with it
horz.myboundary = 3;
horz.pairedboundary = 1;
Add it to the PeriodicBoundaries
dof_map.add_periodic_boundary(horz);
Create a boundary periodic with one displaced 2.0 in the y
direction
PeriodicBoundary vert(RealVectorValue(0., 2.0, 0.));
Connect boundary ids 0 and 2 with it
vert.myboundary = 0;
vert.pairedboundary = 2;
Add it to the PeriodicBoundaries
dof_map.add_periodic_boundary(vert);
Initialize the data structures for the equation system.
if(!read_solution)
equation_systems.init ();
else
equation_systems.reinit ();
Print out the H1 norm of the initialized or saved solution, for
verification purposes:
Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Initial H1 norm = " << H1norm << std::endl << std::endl;
Prints information about the system to the screen.
equation_systems.print_info();
equation_systems.parameters.set<unsigned int>
("linear solver maximum iterations") = 250;
equation_systems.parameters.set<Real>
("linear solver tolerance") = TOLERANCE;
if(!read_solution)
{
Write out the initial condition
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems ("out.gmv.000",
equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems (exodus_filename(0),
equation_systems);
#endif
}
else
{
Write out the solution that was read in
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems ("solution_read_in.gmv",
equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems ("solution_read_in.e",
equation_systems);
#endif
}
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);
The Convection-Diffusion system also requires a specified
diffusivity. We use an isotropic (hence Real) value.
equation_systems.parameters.set<Real>("diffusivity") = 0.01;
Solve the system "Convection-Diffusion". This will be done by
looping over the specified time interval and calling the
\p solve() member at each time step. This will assemble the
system and call the linear solver.
const Real dt = 0.025;
system.time = init_timestep*dt;
Tell the MeshRefinement object about the periodic boundaries so
that it can get heuristics like level-one conformity and
unrefined island elimination right.
mesh_refinement.set_periodic_boundaries_ptr(dof_map.get_periodic_boundaries());
We do 25 timesteps both before and after writing out the
intermediate solution
for(unsigned int t_step=init_timestep;
t_step<(init_timestep+n_timesteps);
t_step++)
{
Increment 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 ";
{
Save flags to avoid polluting cout with custom precision values, etc.
std::ios_base::fmtflags os_flags = std::cout.flags();
std::cout << t_step
<< ", time="
<< std::setw(6)
<< std::setprecision(3)
<< std::setfill('0')
<< std::left
<< system.time
<< "..."
<< std::endl;
Restore flags
std::cout.flags(os_flags);
}
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 TransientLinearImplicitSystems
(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;
The number of refinement steps per time step.
unsigned int max_r_steps = 1;
if(command_line.search("-max_r_steps"))
max_r_steps = command_line.next(0);
A refinement loop.
for (unsigned int r_step=0; r_step<max_r_steps+1; r_step++)
{
Assemble & solve the linear system
system.solve();
Print out the H1 norm, for verification purposes:
H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "H1 norm = " << H1norm << std::endl;
Possibly refine the mesh
if (r_step+1 <= max_r_steps)
{
std::cout << " Refining the mesh..." << std::endl;
The \p ErrorVector is a particular \p StatisticsVector
for computing error information on a finite element mesh.
ErrorVector error;
The \p ErrorEstimator class interrogates a finite element
solution and assigns to each element a positive error value.
This value is used for deciding which elements to refine
and which to coarsen.
ErrorEstimator* error_estimator = new KellyErrorEstimator;
KellyErrorEstimator error_estimator;
Compute the error for each active element using the provided
\p flux_jump indicator. Note in general you will need to
provide an error estimator specifically designed for your
application.
error_estimator.estimate_error (system,
error);
This takes the error in \p error and decides which elements
will be coarsened or refined. Any element within 20% of the
maximum error on any element will be refined, and any
element within 7% of the minimum error on any element might
be coarsened. Note that the elements flagged for refinement
will be refined, but those flagged for coarsening _might_ be
coarsened.
mesh_refinement.refine_fraction() = 0.80;
mesh_refinement.coarsen_fraction() = 0.07;
mesh_refinement.max_h_level() = 5;
mesh_refinement.flag_elements_by_error_fraction (error);
This call actually refines and coarsens the flagged
elements.
mesh_refinement.refine_and_coarsen_elements();
This call reinitializes the \p EquationSystems object for
the newly refined mesh. One of the steps in the
reinitialization is projecting the \p solution,
\p old_solution, etc... vectors from the old mesh to
the current one.
equation_systems.reinit ();
}
}
Again do a search on the command line for an argument
unsigned int output_freq = 10;
if(command_line.search("-output_freq"))
output_freq = command_line.next(0);
Output every 10 timesteps to file.
if ( (t_step+1)%output_freq == 0)
{
OStringStream file_name;
#ifdef LIBMESH_HAVE_GMV
file_name << "out.gmv.";
OSSRealzeroright(file_name,3,0,t_step+1);
GMVIO(mesh).write_equation_systems (file_name.str(), equation_systems);
GMVIO(mesh).write_equation_systems (file_name.str(), equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
So... if paraview is told to open a file called out.e.{N}, it automatically tries to
open out.e.{N-1}, out.e.{N-2}, etc. If we name the file something else, we can work
around that issue, but the right thing to do (for adaptive meshes) is to write a filename
with the adaptation step into a separate file.
ExodusII_IO(mesh).write_equation_systems (exodus_filename(t_step+1),
equation_systems);
#endif
}
}
if(!read_solution)
{
Print out the H1 norm of the saved solution, for verification purposes:
H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Final H1 norm = " << H1norm << std::endl << std::endl;
mesh.write("saved_mesh.xdr");
equation_systems.write("saved_solution.xdr", libMeshEnums::ENCODE);
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems ("saved_solution.gmv",
equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems ("saved_solution.e",
equation_systems);
#endif
}
#endif // #ifndef LIBMESH_ENABLE_AMR
We might have a parser to clean up
delete parsed_solution;
return 0;
}
Here we define the initialization routine for the
Convection-Diffusion system. This routine is
responsible for applying the initial conditions to
the system.
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;
if (parsed_solution)
system.project_solution(parsed_solution, NULL);
else
system.project_solution(exact_value, NULL, es.parameters);
}
This function defines the assembly routine which
will be called at each time step. It is responsible
for computing the proper matrix entries for the
element stiffness matrices and right-hand sides.
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 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();
The element shape functions evaluated at the quadrature points.
const std::vector<std::vector<Real> >& phi = fe->get_phi();
The element shape function gradients evaluated at the quadrature
points.
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
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 diffusivity =
es.parameters.get<Real> ("diffusivity");
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
diffusivity*(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
diffusivity*(dphi[i][qp]*dphi[j][qp]))
);
}
}
}
We have now built the element matrix and RHS vector in terms
of the element degrees of freedom. However, it is possible
that some of the element DOFs are constrained to enforce
solution continuity, i.e. they are not really "free". We need
to constrain those DOFs in terms of non-constrained DOFs to
ensure a continuous solution. The
\p DofMap::constrain_element_matrix_and_vector() method does
just that.
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);
}
Finished computing the sytem matrix and right-hand side.
#endif // #ifdef LIBMESH_ENABLE_AMR
}
std::string exodus_filename(unsigned number)
{
std::ostringstream oss;
oss << "out_";
oss << std::setw(3) << std::setfill('0') << number;
oss << ".e";
return oss.str();
}
The source file exact_solution.C with comments:
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
#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 adaptivity_ex5.C without comments:
#include <iostream>
#include <algorithm>
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath>
#include "libmesh/libmesh.h"
#include "libmesh/serial_mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/gmv_io.h"
#include "libmesh/exodusII_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/periodic_boundaries.h"
#include "libmesh/periodic_boundary.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/parsed_function.h"
#include "libmesh/getpot.h"
#include "libmesh/transient_system.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/vector_value.h"
#include "libmesh/error_vector.h"
#include "libmesh/kelly_error_estimator.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"));
}
FunctionBase<Number>* parsed_solution = NULL;
std::string exodus_filename(unsigned number);
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
#if !defined(LIBMESH_ENABLE_AMR)
libmesh_example_assert(false, "--enable-amr");
#elif !defined(LIBMESH_HAVE_XDR)
libmesh_example_assert(false, "--enable-xdr");
#elif !defined(LIBMESH_ENABLE_PERIODIC)
libmesh_example_assert(false, "--enable-periodic");
#else
libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
std::cout << "Usage:\n"
<<"\t " << argv[0] << " -init_timestep 0\n"
<< "OR\n"
<<"\t " << argv[0] << " -read_solution -init_timestep 26\n"
<< std::endl;
std::cout << "Running: " << argv[0];
for (int i=1; i<argc; i++)
std::cout << " " << argv[i];
std::cout << std::endl << std::endl;
GetPot command_line (argc, argv);
const bool read_solution = command_line.search("-read_solution");
unsigned int init_timestep = 0;
if(command_line.search("-init_timestep"))
init_timestep = command_line.next(0);
else
{
if (libMesh::processor_id() == 0)
std::cerr << "ERROR: Initial timestep not specified\n" << std::endl;
libmesh_error();
}
unsigned int n_timesteps = 0;
if(command_line.search("-n_timesteps"))
n_timesteps = command_line.next(0);
else
{
std::cout << "ERROR: Number of timesteps not specified\n" << std::endl;
libmesh_error();
}
#ifdef LIBMESH_HAVE_FPARSER
const bool have_expression = command_line.search("-exact_solution");
#else
const bool have_expression = false;
#endif
if (have_expression)
parsed_solution = new ParsedFunction<Number>(command_line.next(std::string()));
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
SerialMesh mesh;
EquationSystems equation_systems (mesh);
MeshRefinement mesh_refinement (mesh);
if(!read_solution)
{
MeshTools::Generation::build_square(mesh, 2, 2, 0., 2., 0., 2., QUAD4);
unsigned int n_refinements = 5;
if(command_line.search("-n_refinements"))
n_refinements = command_line.next(0);
if(!read_solution)
mesh_refinement.uniformly_refine (n_refinements);
mesh.print_info();
TransientLinearImplicitSystem & system =
equation_systems.add_system<TransientLinearImplicitSystem>
("Convection-Diffusion");
system.add_variable ("u", FIRST);
system.attach_init_function (init_cd);
}
else
{
mesh.read("saved_mesh.xdr");
mesh.print_info();
equation_systems.read("saved_solution.xdr", libMeshEnums::DECODE);
}
TransientLinearImplicitSystem & system =
equation_systems.get_system<TransientLinearImplicitSystem>
("Convection-Diffusion");
system.attach_assemble_function (assemble_cd);
DofMap & dof_map = system.get_dof_map();
PeriodicBoundary horz(RealVectorValue(2.0, 0., 0.));
horz.myboundary = 3;
horz.pairedboundary = 1;
dof_map.add_periodic_boundary(horz);
PeriodicBoundary vert(RealVectorValue(0., 2.0, 0.));
vert.myboundary = 0;
vert.pairedboundary = 2;
dof_map.add_periodic_boundary(vert);
if(!read_solution)
equation_systems.init ();
else
equation_systems.reinit ();
Real H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Initial H1 norm = " << H1norm << std::endl << std::endl;
equation_systems.print_info();
equation_systems.parameters.set<unsigned int>
("linear solver maximum iterations") = 250;
equation_systems.parameters.set<Real>
("linear solver tolerance") = TOLERANCE;
if(!read_solution)
{
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems ("out.gmv.000",
equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems (exodus_filename(0),
equation_systems);
#endif
}
else
{
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems ("solution_read_in.gmv",
equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems ("solution_read_in.e",
equation_systems);
#endif
}
equation_systems.parameters.set<RealVectorValue>("velocity") =
RealVectorValue (0.8, 0.8);
equation_systems.parameters.set<Real>("diffusivity") = 0.01;
const Real dt = 0.025;
system.time = init_timestep*dt;
mesh_refinement.set_periodic_boundaries_ptr(dof_map.get_periodic_boundaries());
for(unsigned int t_step=init_timestep;
t_step<(init_timestep+n_timesteps);
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::ios_base::fmtflags os_flags = std::cout.flags();
std::cout << t_step
<< ", time="
<< std::setw(6)
<< std::setprecision(3)
<< std::setfill('0')
<< std::left
<< system.time
<< "..."
<< std::endl;
std::cout.flags(os_flags);
}
*system.old_local_solution = *system.current_local_solution;
unsigned int max_r_steps = 1;
if(command_line.search("-max_r_steps"))
max_r_steps = command_line.next(0);
for (unsigned int r_step=0; r_step<max_r_steps+1; r_step++)
{
system.solve();
H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "H1 norm = " << H1norm << std::endl;
if (r_step+1 <= max_r_steps)
{
std::cout << " Refining the mesh..." << std::endl;
ErrorVector error;
KellyErrorEstimator error_estimator;
error_estimator.estimate_error (system,
error);
mesh_refinement.refine_fraction() = 0.80;
mesh_refinement.coarsen_fraction() = 0.07;
mesh_refinement.max_h_level() = 5;
mesh_refinement.flag_elements_by_error_fraction (error);
mesh_refinement.refine_and_coarsen_elements();
equation_systems.reinit ();
}
}
unsigned int output_freq = 10;
if(command_line.search("-output_freq"))
output_freq = command_line.next(0);
if ( (t_step+1)%output_freq == 0)
{
#ifdef LIBMESH_HAVE_GMV
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems (exodus_filename(t_step+1),
equation_systems);
#endif
}
}
if(!read_solution)
{
H1norm = system.calculate_norm(*system.solution, SystemNorm(H1));
std::cout << "Final H1 norm = " << H1norm << std::endl << std::endl;
mesh.write("saved_mesh.xdr");
equation_systems.write("saved_solution.xdr", libMeshEnums::ENCODE);
#ifdef LIBMESH_HAVE_GMV
GMVIO(mesh).write_equation_systems ("saved_solution.gmv",
equation_systems);
#endif
#ifdef LIBMESH_HAVE_EXODUS_API
ExodusII_IO(mesh).write_equation_systems ("saved_solution.e",
equation_systems);
#endif
}
#endif // #ifndef LIBMESH_ENABLE_AMR
delete parsed_solution;
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;
if (parsed_solution)
system.project_solution(parsed_solution, NULL);
else
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<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
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 diffusivity =
es.parameters.get<Real> ("diffusivity");
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] +
diffusivity*(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] +
diffusivity*(dphi[i][qp]*dphi[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
}
std::string exodus_filename(unsigned number)
{
std::ostringstream oss;
oss << "out_";
oss << std::setw(3) << std::setfill('0') << number;
oss << ".e";
return oss.str();
}
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 console output of the program:
***************************************************************
* Running Example adaptivity_ex5:
* mpirun -np 2 example-devel -n_timesteps 25 -n_refinements 5 -output_freq 10 -init_timestep 0 -exact_solution '10*exp(-(pow(x-0.8*t-0.2,2)+pow(y-0.8*t-0.2,2))/(0.01*(4*t+1)))/(4*t+1)' -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Usage:
/workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel -init_timestep 0
OR
/workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel -read_solution -init_timestep 26
Running: /workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel -n_timesteps 25 -n_refinements 5 -output_freq 10 -init_timestep 0 -exact_solution 10*exp(-(pow(x-0.8*t-0.2,2) + pow(y-0.8*t-0.2,2))/(0.01*(4*t+1)))/(4*t+1) -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()=4225
n_local_nodes()=2152
n_elem()=5460
n_local_elem()=2760
n_active_elem()=4096
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
Initial H1 norm = 17.4175
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()=4225
n_local_dofs()=2152
n_constrained_dofs()=129
n_local_constrained_dofs()=32
n_vectors()=3
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 8.98462
Average Off-Processor Bandwidth <= 0.299172
Maximum On-Processor Bandwidth <= 14
Maximum Off-Processor Bandwidth <= 8
DofMap Constraints
Number of DoF Constraints = 129
Average DoF Constraint Length= 1
Number of Node Constraints = 129
Maximum Node Constraint Length= 2
Average Node Constraint Length= 2
Solving time step 0, time=0.0250...
H1 norm = 15.9
Refining the mesh...
H1 norm = 15.9
Solving time step 1, time=0.0500...
H1 norm = 14.6
Refining the mesh...
H1 norm = 14.6
Solving time step 2, time=0.0750...
H1 norm = 13.5
Refining the mesh...
H1 norm = 13.5
Solving time step 3, time=0.1000...
H1 norm = 12.6
Refining the mesh...
H1 norm = 12.6
Solving time step 4, time=0.1250...
H1 norm = 11.7
Refining the mesh...
H1 norm = 11.7
Solving time step 5, time=0.1500...
H1 norm = 11
Refining the mesh...
H1 norm = 11
Solving time step 6, time=0.1750...
H1 norm = 10.4
Refining the mesh...
H1 norm = 10.4
Solving time step 7, time=0.2000...
H1 norm = 9.82
Refining the mesh...
H1 norm = 9.82
Solving time step 8, time=0.2250...
H1 norm = 9.31
Refining the mesh...
H1 norm = 9.31
Solving time step 9, time=0.2500...
H1 norm = 8.85
Refining the mesh...
H1 norm = 8.85
Solving time step 10, time=0.2750...
H1 norm = 8.43
Refining the mesh...
H1 norm = 8.43
Solving time step 11, time=0.3000...
H1 norm = 8.05
Refining the mesh...
H1 norm = 8.05
Solving time step 12, time=0.3250...
H1 norm = 7.71
Refining the mesh...
H1 norm = 7.71
Solving time step 13, time=0.3500...
H1 norm = 7.39
Refining the mesh...
H1 norm = 7.39
Solving time step 14, time=0.3750...
H1 norm = 7.1
Refining the mesh...
H1 norm = 7.1
Solving time step 15, time=0.4000...
H1 norm = 6.83
Refining the mesh...
H1 norm = 6.83
Solving time step 16, time=0.4250...
H1 norm = 6.58
Refining the mesh...
H1 norm = 6.58
Solving time step 17, time=0.4500...
H1 norm = 6.35
Refining the mesh...
H1 norm = 6.35
Solving time step 18, time=0.4750...
H1 norm = 6.13
Refining the mesh...
H1 norm = 6.13
Solving time step 19, time=0.5000...
H1 norm = 5.93
Refining the mesh...
H1 norm = 5.93
Solving time step 20, time=0.5250...
H1 norm = 5.74
Refining the mesh...
H1 norm = 5.74
Solving time step 21, time=0.5500...
H1 norm = 5.56
Refining the mesh...
H1 norm = 5.56
Solving time step 22, time=0.5750...
H1 norm = 5.4
Refining the mesh...
H1 norm = 5.4
Solving time step 23, time=0.6000...
H1 norm = 5.24
Refining the mesh...
H1 norm = 5.24
Solving time step 24, time=0.6250...
H1 norm = 5.09
Refining the mesh...
H1 norm = 5.09
Final H1 norm = 5.09
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:39:03 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.161e+01 1.00000 1.161e+01
Objects: 2.571e+03 1.00000 2.571e+03
Flops: 3.234e+07 1.18994 2.976e+07 5.951e+07
Flops/sec: 2.786e+06 1.18994 2.564e+06 5.127e+06
MPI Messages: 1.946e+03 1.00000 1.946e+03 3.891e+03
MPI Message Lengths: 1.562e+06 1.01944 7.955e+02 3.095e+06
MPI Reductions: 5.803e+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.1607e+01 100.0% 5.9512e+07 100.0% 3.891e+03 100.0% 7.955e+02 100.0% 5.802e+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 246 1.0 2.0816e-03 1.7 5.05e+05 1.1 0.0e+00 0.0e+00 2.5e+02 0 2 0 0 4 0 2 0 0 4 459
VecNorm 346 1.0 1.3165e-03 1.3 2.36e+05 1.1 0.0e+00 0.0e+00 3.5e+02 0 1 0 0 6 0 1 0 0 6 340
VecScale 296 1.0 1.3804e-04 1.1 1.01e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1381
VecCopy 301 1.0 1.5998e-04 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
VecSet 832 1.0 3.3903e-04 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
VecAXPY 100 1.0 9.3937e-05 1.0 6.97e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1407
VecMAXPY 296 1.0 2.2793e-04 1.0 6.72e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 2 0 0 0 0 2 0 0 0 5578
VecAssemblyBegin 804 1.0 2.8564e-02 1.4 0.00e+00 0.0 2.4e+02 3.1e+02 2.2e+03 0 0 6 2 38 0 0 6 2 38 0
VecAssemblyEnd 804 1.0 4.0174e-04 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 975 1.0 1.3030e-03 1.0 0.00e+00 0.0 1.7e+03 8.3e+02 0.0e+00 0 0 45 47 0 0 0 45 47 0 0
VecScatterEnd 975 1.0 2.2590e-02 2.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
VecNormalize 296 1.0 1.2729e-03 1.1 3.02e+05 1.1 0.0e+00 0.0e+00 3.0e+02 0 1 0 0 5 0 1 0 0 5 449
MatMult 296 1.0 2.3772e-02 1.9 1.88e+06 1.1 5.9e+02 5.4e+02 0.0e+00 0 6 15 10 0 0 6 15 10 0 149
MatSolve 346 1.0 4.8864e-03 1.2 1.07e+07 1.2 0.0e+00 0.0e+00 0.0e+00 0 33 0 0 0 0 33 0 0 0 4066
MatLUFactorNum 50 1.0 1.5325e-02 1.2 1.82e+07 1.2 0.0e+00 0.0e+00 0.0e+00 0 56 0 0 0 0 56 0 0 0 2160
MatILUFactorSym 50 1.0 4.7266e-02 1.2 0.00e+00 0.0 0.0e+00 0.0e+00 1.5e+02 0 0 0 0 3 0 0 0 0 3 0
MatAssemblyBegin 100 1.0 2.3938e-02 1.5 0.00e+00 0.0 2.9e+02 1.5e+03 2.0e+02 0 0 7 14 3 0 0 7 14 3 0
MatAssemblyEnd 100 1.0 3.6945e-03 1.1 0.00e+00 0.0 1.0e+02 1.4e+02 2.1e+02 0 0 3 0 4 0 0 3 0 4 0
MatGetRowIJ 50 1.0 8.8215e-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 50 1.0 7.3385e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 1.0e+02 0 0 0 0 2 0 0 0 0 2 0
MatZeroEntries 102 1.0 2.0170e-04 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 246 1.0 2.4226e-03 1.5 1.01e+06 1.1 0.0e+00 0.0e+00 2.5e+02 0 3 0 0 4 0 3 0 0 4 789
KSPSetUp 100 1.0 4.5943e-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 9.2873e-02 1.0 3.23e+07 1.2 5.9e+02 5.4e+02 8.9e+02 1100 15 10 15 1100 15 10 15 641
PCSetUp 100 1.0 6.8237e-02 1.2 1.82e+07 1.2 0.0e+00 0.0e+00 3.0e+02 1 56 0 0 5 1 56 0 0 5 485
PCSetUpOnBlocks 50 1.0 6.6404e-02 1.2 1.82e+07 1.2 0.0e+00 0.0e+00 2.5e+02 1 56 0 0 4 1 56 0 0 4 499
PCApply 346 1.0 7.0581e-03 1.1 1.07e+07 1.2 0.0e+00 0.0e+00 0.0e+00 0 33 0 0 0 0 33 0 0 0 2815
------------------------------------------------------------------------------------------------------------------------
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 1184 1184 5010288 0
Vector Scatter 359 359 371924 0
Index Set 665 665 640652 0
IS L to G Mapping 130 130 73320 0
Matrix 128 128 11362692 0
Krylov Solver 52 52 503360 0
Preconditioner 52 52 46384 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 0
Average time for MPI_Barrier(): 8.10623e-07
Average time for zero size MPI_Send(): 1.44243e-05
#PETSc Option Table entries:
-exact_solution 10*exp(-(pow(x-0.8*t-0.2,2) + pow(y-0.8*t-0.2,2))/(0.01*(4*t+1)))/(4*t+1)
-init_timestep 0
-ksp_right_pc
-log_summary
-n_refinements 5
-n_timesteps 25
-output_freq 10
-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:39:03 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=11.618, Active time=11.3686 |
----------------------------------------------------------------------------------------------------------------
| 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() 26 0.2018 0.007761 0.2518 0.009683 1.77 2.21 |
| build_constraint_matrix() 14760 0.0906 0.000006 0.0906 0.000006 0.80 0.80 |
| build_sparsity() 26 0.1543 0.005935 0.4009 0.015420 1.36 3.53 |
| cnstrn_elem_mat_vec() 14760 0.0569 0.000004 0.0569 0.000004 0.50 0.50 |
| create_dof_constraints() 26 0.7767 0.029874 1.3287 0.051103 6.83 11.69 |
| distribute_dofs() 26 0.2149 0.008265 0.5792 0.022279 1.89 5.10 |
| dof_indices() 113363 2.9987 0.000026 2.9987 0.000026 26.38 26.38 |
| enforce_constraints_exactly() 76 0.0090 0.000119 0.0090 0.000119 0.08 0.08 |
| old_dof_indices() 46863 1.2541 0.000027 1.2541 0.000027 11.03 11.03 |
| prepare_send_list() 26 0.0009 0.000033 0.0009 0.000033 0.01 0.01 |
| reinit() 26 0.3567 0.013721 0.3567 0.013721 3.14 3.14 |
| |
| EquationSystems |
| build_solution_vector() 6 0.0181 0.003013 0.1606 0.026761 0.16 1.41 |
| write() 1 0.0076 0.007648 0.0079 0.007878 0.07 0.07 |
| |
| ExodusII_IO |
| write_nodal_data() 4 0.0240 0.005990 0.0240 0.005990 0.21 0.21 |
| |
| FE |
| compute_shape_functions() 67572 0.3307 0.000005 0.3307 0.000005 2.91 2.91 |
| init_shape_functions() 35774 0.2138 0.000006 0.2138 0.000006 1.88 1.88 |
| inverse_map() 108423 0.4890 0.000005 0.4890 0.000005 4.30 4.30 |
| |
| FEMap |
| compute_affine_map() 67572 0.3412 0.000005 0.3412 0.000005 3.00 3.00 |
| compute_face_map() 17836 0.1747 0.000010 0.3545 0.000020 1.54 3.12 |
| init_face_shape_functions() 1207 0.0089 0.000007 0.0089 0.000007 0.08 0.08 |
| init_reference_to_physical_map() 35774 0.2701 0.000008 0.2701 0.000008 2.38 2.38 |
| |
| GMVIO |
| write_nodal_data() 2 0.0172 0.008609 0.0173 0.008656 0.15 0.15 |
| |
| JumpErrorEstimator |
| estimate_error() 25 0.7401 0.029602 2.9187 0.116747 6.51 25.67 |
| |
| LocationMap |
| find() 19140 0.0452 0.000002 0.0452 0.000002 0.40 0.40 |
| init() 55 0.0316 0.000575 0.0316 0.000575 0.28 0.28 |
| |
| Mesh |
| contract() 25 0.0117 0.000469 0.0197 0.000790 0.10 0.17 |
| find_neighbors() 27 0.3378 0.012510 0.3452 0.012785 2.97 3.04 |
| renumber_nodes_and_elem() 79 0.0199 0.000251 0.0199 0.000251 0.17 0.17 |
| |
| MeshCommunication |
| assign_global_indices() 1 0.0192 0.019180 0.0195 0.019500 0.17 0.17 |
| compute_hilbert_indices() 28 0.0714 0.002549 0.0714 0.002549 0.63 0.63 |
| find_global_indices() 28 0.0307 0.001097 0.1123 0.004012 0.27 0.99 |
| parallel_sort() 28 0.0071 0.000255 0.0082 0.000293 0.06 0.07 |
| |
| MeshOutput |
| write_equation_systems() 6 0.0002 0.000039 0.2024 0.033733 0.00 1.78 |
| |
| MeshRefinement |
| _coarsen_elements() 50 0.0091 0.000183 0.0094 0.000188 0.08 0.08 |
| _refine_elements() 55 0.0773 0.001405 0.1726 0.003139 0.68 1.52 |
| add_point() 19140 0.0451 0.000002 0.0927 0.000005 0.40 0.82 |
| make_coarsening_compatible() 111 0.2256 0.002033 0.2996 0.002700 1.98 2.64 |
| make_refinement_compatible() 111 0.0103 0.000093 0.0112 0.000101 0.09 0.10 |
| |
| MeshTools::Generation |
| build_cube() 1 0.0004 0.000371 0.0004 0.000371 0.00 0.00 |
| |
| MetisPartitioner |
| partition() 27 0.2107 0.007802 0.3244 0.012015 1.85 2.85 |
| |
| Parallel |
| allgather() 143 0.0009 0.000007 0.0012 0.000008 0.01 0.01 |
| barrier() 1 0.0000 0.000015 0.0000 0.000015 0.00 0.00 |
| broadcast() 2 0.0000 0.000009 0.0000 0.000009 0.00 0.00 |
| gather() 14 0.0001 0.000008 0.0001 0.000008 0.00 0.00 |
| max(bool) 267 0.0056 0.000021 0.0056 0.000021 0.05 0.05 |
| max(scalar) 5720 0.0140 0.000002 0.0140 0.000002 0.12 0.12 |
| max(vector) 1448 0.0086 0.000006 0.0192 0.000013 0.08 0.17 |
| min(bool) 7177 0.0165 0.000002 0.0165 0.000002 0.14 0.14 |
| min(scalar) 5661 0.1078 0.000019 0.1078 0.000019 0.95 0.95 |
| min(vector) 1448 0.0095 0.000007 0.0207 0.000014 0.08 0.18 |
| probe() 303 0.0029 0.000009 0.0029 0.000009 0.03 0.03 |
| receive() 297 0.0017 0.000006 0.0046 0.000015 0.01 0.04 |
| send() 284 0.0007 0.000003 0.0007 0.000003 0.01 0.01 |
| send_receive() 340 0.0027 0.000008 0.0083 0.000024 0.02 0.07 |
| sum() 278 0.0100 0.000036 0.0129 0.000046 0.09 0.11 |
| |
| Parallel::Request |
| wait() 286 0.0004 0.000001 0.0004 0.000001 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 27 0.0282 0.001046 0.0306 0.001133 0.25 0.27 |
| set_parent_processor_ids() 27 0.0259 0.000958 0.0259 0.000958 0.23 0.23 |
| |
| PetscLinearSolver |
| solve() 50 0.1246 0.002492 0.1246 0.002492 1.10 1.10 |
| |
| PointLocatorTree |
| init(no master) 50 0.2220 0.004440 0.2294 0.004588 1.95 2.02 |
| operator() 6946 0.1599 0.000023 0.1983 0.000029 1.41 1.74 |
| |
| ProjectVector |
| operator() 75 0.1146 0.001528 1.1704 0.015606 1.01 10.30 |
| |
| System |
| assemble() 50 0.3141 0.006282 1.0812 0.021623 2.76 9.51 |
| calculate_norm() 52 0.1368 0.002632 0.8250 0.015865 1.20 7.26 |
| project_vector() 76 0.1516 0.001995 2.1272 0.027989 1.33 18.71 |
| |
| XdrIO |
| write() 1 0.0059 0.005853 0.0064 0.006438 0.05 0.06 |
----------------------------------------------------------------------------------------------------------------
| Totals: 594109 11.3686 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adaptivity_ex5:
* mpirun -np 2 example-devel -n_timesteps 25 -n_refinements 5 -output_freq 10 -init_timestep 0 -exact_solution '10*exp(-(pow(x-0.8*t-0.2,2)+pow(y-0.8*t-0.2,2))/(0.01*(4*t+1)))/(4*t+1)' -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
***** Finished first 25 steps, now read in saved solution and continue *****
***************************************************************
* Running Example adaptivity_ex5:
* mpirun -np 2 example-devel -read_solution -n_timesteps 25 -output_freq 10 -init_timestep 25 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Usage:
/workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel -init_timestep 0
OR
/workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel -read_solution -init_timestep 26
Running: /workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel -read_solution -n_timesteps 25 -output_freq 10 -init_timestep 25 -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()=737
n_local_nodes()=385
n_elem()=884
n_local_elem()=456
n_active_elem()=664
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
Initial H1 norm = 5.09105
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()=737
n_local_dofs()=385
n_constrained_dofs()=130
n_local_constrained_dofs()=61
n_vectors()=3
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 9.35007
Average Off-Processor Bandwidth <= 0.626866
Maximum On-Processor Bandwidth <= 21
Maximum Off-Processor Bandwidth <= 12
DofMap Constraints
Number of DoF Constraints = 130
Average DoF Constraint Length= 1.92308
Number of Node Constraints = 236
Maximum Node Constraint Length= 5
Average Node Constraint Length= 2.51271
Solving time step 25, time=0.6500...
H1 norm = 4.95
Refining the mesh...
H1 norm = 4.95
Solving time step 26, time=0.6750...
H1 norm = 4.82
Refining the mesh...
H1 norm = 4.82
Solving time step 27, time=0.7000...
H1 norm = 4.69
Refining the mesh...
H1 norm = 4.69
Solving time step 28, time=0.7250...
H1 norm = 4.58
Refining the mesh...
H1 norm = 4.58
Solving time step 29, time=0.7500...
H1 norm = 4.46
Refining the mesh...
H1 norm = 4.46
Solving time step 30, time=0.7750...
H1 norm = 4.36
Refining the mesh...
H1 norm = 4.36
Solving time step 31, time=0.8000...
H1 norm = 4.25
Refining the mesh...
H1 norm = 4.25
Solving time step 32, time=0.8250...
H1 norm = 4.16
Refining the mesh...
H1 norm = 4.16
Solving time step 33, time=0.8500...
H1 norm = 4.06
Refining the mesh...
H1 norm = 4.06
Solving time step 34, time=0.8750...
H1 norm = 3.97
Refining the mesh...
H1 norm = 3.97
Solving time step 35, time=0.9000...
H1 norm = 3.89
Refining the mesh...
H1 norm = 3.89
Solving time step 36, time=0.9250...
H1 norm = 3.81
Refining the mesh...
H1 norm = 3.81
Solving time step 37, time=0.9500...
H1 norm = 3.73
Refining the mesh...
H1 norm = 3.73
Solving time step 38, time=0.9750...
H1 norm = 3.66
Refining the mesh...
H1 norm = 3.66
Solving time step 39, time=100000...
H1 norm = 3.58
Refining the mesh...
H1 norm = 3.58
Solving time step 40, time=1.0300...
H1 norm = 3.51
Refining the mesh...
H1 norm = 3.51
Solving time step 41, time=1.0500...
H1 norm = 3.45
Refining the mesh...
H1 norm = 3.45
Solving time step 42, time=1.0700...
H1 norm = 3.38
Refining the mesh...
H1 norm = 3.38
Solving time step 43, time=1.1000...
H1 norm = 3.32
Refining the mesh...
H1 norm = 3.32
Solving time step 44, time=1.1200...
H1 norm = 3.26
Refining the mesh...
H1 norm = 3.26
Solving time step 45, time=1.1500...
H1 norm = 3.21
Refining the mesh...
H1 norm = 3.21
Solving time step 46, time=1.1700...
H1 norm = 3.15
Refining the mesh...
H1 norm = 3.15
Solving time step 47, time=1.2000...
H1 norm = 3.1
Refining the mesh...
H1 norm = 3.1
Solving time step 48, time=1.2200...
H1 norm = 3.05
Refining the mesh...
H1 norm = 3.05
Solving time step 49, time=1.2500...
H1 norm = 3
Refining the mesh...
H1 norm = 3
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adaptivity/adaptivity_ex5/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:39:18 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.457e+01 1.00000 1.457e+01
Objects: 2.613e+03 1.00000 2.613e+03
Flops: 4.711e+07 1.09112 4.515e+07 9.029e+07
Flops/sec: 3.233e+06 1.09112 3.098e+06 6.196e+06
MPI Messages: 1.970e+03 1.00000 1.970e+03 3.941e+03
MPI Message Lengths: 1.718e+06 1.02315 8.622e+02 3.398e+06
MPI Reductions: 5.878e+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.4573e+01 100.0% 9.0294e+07 100.0% 3.941e+03 100.0% 8.622e+02 100.0% 5.877e+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 226 1.0 1.6725e-03 1.8 6.64e+05 1.1 0.0e+00 0.0e+00 2.3e+02 0 1 0 0 4 0 1 0 0 4 755
VecNorm 326 1.0 1.2486e-03 1.3 3.23e+05 1.1 0.0e+00 0.0e+00 3.3e+02 0 1 0 0 6 0 1 0 0 6 492
VecScale 276 1.0 1.4091e-04 1.1 1.37e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1846
VecCopy 308 1.0 1.6403e-04 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
VecSet 822 1.0 3.0255e-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
VecAXPY 100 1.0 1.0395e-04 1.1 9.91e+04 1.1 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 1815
VecMAXPY 276 1.0 3.0851e-04 1.1 8.88e+05 1.1 0.0e+00 0.0e+00 0.0e+00 0 2 0 0 0 0 2 0 0 0 5480
VecAssemblyBegin 827 1.0 2.6139e-02 2.1 0.00e+00 0.0 2.5e+02 3.6e+02 2.2e+03 0 0 6 3 38 0 0 6 3 38 0
VecAssemblyEnd 827 1.0 3.9387e-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
VecScatterBegin 975 1.0 1.2980e-03 1.0 0.00e+00 0.0 1.7e+03 8.9e+02 0.0e+00 0 0 44 46 0 0 0 44 46 0 0
VecScatterEnd 975 1.0 1.1665e-02 2.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 276 1.0 1.2870e-03 1.2 4.10e+05 1.1 0.0e+00 0.0e+00 2.8e+02 0 1 0 0 5 0 1 0 0 5 606
MatMult 276 1.0 1.3233e-02 2.2 2.49e+06 1.1 5.5e+02 5.9e+02 0.0e+00 0 5 14 10 0 0 5 14 10 0 358
MatSolve 326 1.0 6.8872e-03 1.1 1.50e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 32 0 0 0 0 32 0 0 0 4158
MatLUFactorNum 50 1.0 2.4397e-02 1.1 2.75e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 59 0 0 0 0 59 0 0 0 2168
MatILUFactorSym 50 1.0 6.7325e-02 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 1.5e+02 0 0 0 0 3 0 0 0 0 3 0
MatAssemblyBegin 100 1.0 3.9896e-0221.3 0.00e+00 0.0 2.9e+02 1.6e+03 2.0e+02 0 0 7 14 3 0 0 7 14 3 0
MatAssemblyEnd 100 1.0 4.1165e-03 1.2 0.00e+00 0.0 1.0e+02 1.5e+02 2.1e+02 0 0 3 0 4 0 0 3 0 4 0
MatGetRowIJ 50 1.0 6.1989e-06 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
MatGetOrdering 50 1.0 7.3934e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 1.0e+02 0 0 0 0 2 0 0 0 0 2 0
MatZeroEntries 104 1.0 1.8358e-04 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 226 1.0 2.0599e-03 1.6 1.33e+06 1.1 0.0e+00 0.0e+00 2.3e+02 0 3 0 0 4 0 3 0 0 4 1227
KSPSetUp 100 1.0 4.2391e-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 1.1723e-01 1.0 4.71e+07 1.1 5.5e+02 5.9e+02 8.5e+02 1100 14 10 15 1100 14 10 15 770
PCSetUp 100 1.0 9.7395e-02 1.1 2.75e+07 1.1 0.0e+00 0.0e+00 3.0e+02 1 59 0 0 5 1 59 0 0 5 543
PCSetUpOnBlocks 50 1.0 9.5548e-02 1.1 2.75e+07 1.1 0.0e+00 0.0e+00 2.5e+02 1 59 0 0 4 1 59 0 0 4 554
PCApply 326 1.0 9.0082e-03 1.1 1.50e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 32 0 0 0 0 32 0 0 0 3179
------------------------------------------------------------------------------------------------------------------------
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 1203 1203 5998712 0
Vector Scatter 366 366 379176 0
Index Set 675 675 669908 0
IS L to G Mapping 133 133 75012 0
Matrix 131 131 16270956 0
Krylov Solver 52 52 503360 0
Preconditioner 52 52 46384 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 5.72205e-07
Average time for zero size MPI_Send(): 1.2517e-05
#PETSc Option Table entries:
-init_timestep 25
-ksp_right_pc
-log_summary
-n_timesteps 25
-output_freq 10
-pc_type bjacobi
-read_solution
-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:39:18 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=14.5834, Active time=14.3004 |
----------------------------------------------------------------------------------------------------------------
| 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() 27 0.2656 0.009838 0.3343 0.012381 1.86 2.34 |
| build_constraint_matrix() 21596 0.1252 0.000006 0.1252 0.000006 0.88 0.88 |
| build_sparsity() 27 0.2130 0.007889 0.5412 0.020044 1.49 3.78 |
| cnstrn_elem_mat_vec() 21596 0.0772 0.000004 0.0772 0.000004 0.54 0.54 |
| create_dof_constraints() 27 0.8225 0.030462 1.2771 0.047301 5.75 8.93 |
| distribute_dofs() 27 0.2600 0.009629 0.7062 0.026155 1.82 4.94 |
| dof_indices() 152261 4.0219 0.000026 4.0219 0.000026 28.12 28.12 |
| enforce_constraints_exactly() 78 0.0099 0.000127 0.0099 0.000127 0.07 0.07 |
| old_dof_indices() 68301 1.8147 0.000027 1.8147 0.000027 12.69 12.69 |
| prepare_send_list() 27 0.0011 0.000042 0.0011 0.000042 0.01 0.01 |
| reinit() 27 0.4403 0.016308 0.4403 0.016308 3.08 3.08 |
| |
| EquationSystems |
| build_solution_vector() 5 0.0077 0.001533 0.0630 0.012599 0.05 0.44 |
| read() 1 0.0105 0.010470 0.0907 0.090694 0.07 0.63 |
| update() 1 0.0001 0.000121 0.0001 0.000121 0.00 0.00 |
| |
| ExodusII_IO |
| write_nodal_data() 4 0.0154 0.003849 0.0154 0.003849 0.11 0.11 |
| |
| FE |
| compute_shape_functions() 89062 0.4554 0.000005 0.4554 0.000005 3.18 3.18 |
| init_shape_functions() 45639 0.2788 0.000006 0.2788 0.000006 1.95 1.95 |
| inverse_map() 123282 0.5787 0.000005 0.5787 0.000005 4.05 4.05 |
| |
| FEMap |
| compute_affine_map() 89062 0.4366 0.000005 0.4366 0.000005 3.05 3.05 |
| compute_face_map() 22769 0.2191 0.000010 0.4605 0.000020 1.53 3.22 |
| init_face_shape_functions() 662 0.0049 0.000007 0.0049 0.000007 0.03 0.03 |
| init_reference_to_physical_map() 45639 0.3261 0.000007 0.3261 0.000007 2.28 2.28 |
| |
| GMVIO |
| write_nodal_data() 1 0.0028 0.002794 0.0028 0.002842 0.02 0.02 |
| |
| JumpErrorEstimator |
| estimate_error() 25 0.9934 0.039737 3.9038 0.156150 6.95 27.30 |
| |
| LocationMap |
| find() 3672 0.0091 0.000002 0.0091 0.000002 0.06 0.06 |
| init() 51 0.0378 0.000742 0.0378 0.000742 0.26 0.26 |
| |
| Mesh |
| contract() 26 0.0070 0.000268 0.0152 0.000586 0.05 0.11 |
| find_neighbors() 26 0.3768 0.014491 0.3786 0.014563 2.63 2.65 |
| renumber_nodes_and_elem() 78 0.0210 0.000269 0.0210 0.000269 0.15 0.15 |
| |
| MeshCommunication |
| assign_global_indices() 1 0.0190 0.019044 0.0194 0.019400 0.13 0.14 |
| compute_hilbert_indices() 26 0.0946 0.003639 0.0946 0.003639 0.66 0.66 |
| find_global_indices() 26 0.0386 0.001484 0.1453 0.005590 0.27 1.02 |
| parallel_sort() 26 0.0092 0.000352 0.0101 0.000387 0.06 0.07 |
| |
| MeshOutput |
| write_equation_systems() 5 0.0002 0.000039 0.0817 0.016340 0.00 0.57 |
| |
| MeshRefinement |
| _coarsen_elements() 51 0.0077 0.000151 0.0080 0.000156 0.05 0.06 |
| _refine_elements() 51 0.0324 0.000634 0.0526 0.001031 0.23 0.37 |
| add_point() 3672 0.0099 0.000003 0.0194 0.000005 0.07 0.14 |
| make_coarsening_compatible() 111 0.2780 0.002504 0.2976 0.002681 1.94 2.08 |
| make_refinement_compatible() 111 0.0137 0.000123 0.0146 0.000132 0.10 0.10 |
| |
| MetisPartitioner |
| partition() 26 0.2716 0.010447 0.4191 0.016119 1.90 2.93 |
| |
| Parallel |
| allgather() 140 0.0011 0.000008 0.0014 0.000010 0.01 0.01 |
| broadcast() 46 0.0008 0.000016 0.0006 0.000014 0.01 0.00 |
| max(bool) 264 0.0096 0.000036 0.0096 0.000036 0.07 0.07 |
| max(scalar) 5705 0.0138 0.000002 0.0138 0.000002 0.10 0.10 |
| max(vector) 1446 0.0088 0.000006 0.0191 0.000013 0.06 0.13 |
| min(bool) 7162 0.0170 0.000002 0.0170 0.000002 0.12 0.12 |
| min(scalar) 5648 0.1036 0.000018 0.1036 0.000018 0.72 0.72 |
| min(vector) 1446 0.0097 0.000007 0.0209 0.000014 0.07 0.15 |
| probe() 280 0.0030 0.000011 0.0030 0.000011 0.02 0.02 |
| receive() 278 0.0017 0.000006 0.0047 0.000017 0.01 0.03 |
| send() 278 0.0007 0.000003 0.0007 0.000003 0.00 0.00 |
| send_receive() 328 0.0029 0.000009 0.0088 0.000027 0.02 0.06 |
| sum() 267 0.0145 0.000054 0.0153 0.000057 0.10 0.11 |
| |
| Parallel::Request |
| wait() 278 0.0004 0.000001 0.0004 0.000001 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 27 0.0304 0.001128 0.0342 0.001266 0.21 0.24 |
| set_parent_processor_ids() 26 0.0315 0.001212 0.0315 0.001212 0.22 0.22 |
| |
| PetscLinearSolver |
| solve() 50 0.1654 0.003308 0.1654 0.003308 1.16 1.16 |
| |
| PointLocatorTree |
| init(no master) 51 0.2796 0.005483 0.2872 0.005632 1.96 2.01 |
| operator() 3698 0.0499 0.000013 0.0712 0.000019 0.35 0.50 |
| |
| ProjectVector |
| operator() 78 0.1881 0.002411 2.0003 0.025645 1.32 13.99 |
| |
| System |
| assemble() 50 0.4551 0.009101 1.5632 0.031264 3.18 10.93 |
| calculate_norm() 51 0.1776 0.003483 1.0681 0.020943 1.24 7.47 |
| project_vector() 78 0.1227 0.001573 3.0614 0.039249 0.86 21.41 |
| |
| XdrIO |
| read() 1 0.0050 0.005049 0.0052 0.005192 0.04 0.04 |
----------------------------------------------------------------------------------------------------------------
| Totals: 715781 14.3004 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adaptivity_ex5:
* mpirun -np 2 example-devel -read_solution -n_timesteps 25 -output_freq 10 -init_timestep 25 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc -log_summary
***************************************************************
Adaptivity Example 5 - Periodic Boundary Conditions with Adaptive Mesh Refinement
This example uses the same simple, linear transient system as in example 10; however in this case periodic boundary conditions are applied at the sides of the domain.
This code also contains an example use of ParsedFunction, to allow users to specify an exact solution on the command line.
C++ include files that we need