The source file femparameters.h with comments:
#ifndef __fem_parameters_h__
#define __fem_parameters_h__
#include <limits>
#include <string>
#include "libmesh/libmesh_common.h"
#include "libmesh/getpot.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
class FEMParameters
{
public:
FEMParameters() :
domainfile("lshaped.xda"),
coarserefinements(0),
solver_quiet(true),
reuse_preconditioner(false),
require_residual_reduction(true),
min_step_length(1e-5),
max_linear_iterations(200000), max_nonlinear_iterations(20),
relative_step_tolerance(1.e-7), relative_residual_tolerance(1.e-10),
initial_linear_tolerance(1.e-3), minimum_linear_tolerance(TOLERANCE*TOLERANCE),
linear_tolerance_multiplier(1.e-3),
nelem_target(30000), global_tolerance(0.0),
refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
refine_uniformly(false),
max_adaptivesteps(1),
indicator_type("kelly"), patch_reuse(true),
fe_family(1, "LAGRANGE"), fe_order(1, 1),
analytic_jacobians(true), verify_analytic_jacobians(0.0),
print_solution_norms(false), print_solutions(false),
print_residual_norms(false), print_residuals(false),
print_jacobian_norms(false), print_jacobians(false) {}
void read(GetPot &input);
std::string domainfile;
unsigned int coarserefinements;
bool solver_quiet, reuse_preconditioner, require_residual_reduction;
Real min_step_length;
unsigned int max_linear_iterations, max_nonlinear_iterations;
Real relative_step_tolerance, relative_residual_tolerance,
initial_linear_tolerance, minimum_linear_tolerance,
linear_tolerance_multiplier;
unsigned int nelem_target;
Real global_tolerance;
Real refine_fraction, coarsen_fraction, coarsen_threshold;
bool refine_uniformly;
unsigned int max_adaptivesteps;
std::string indicator_type;
bool patch_reuse;
std::vector<std::string> fe_family;
std::vector<unsigned int> fe_order;
bool analytic_jacobians;
Real verify_analytic_jacobians;
bool print_solution_norms, print_solutions,
print_residual_norms, print_residuals,
print_jacobian_norms, print_jacobians;
};
#endif // __fem_parameters_h__
The source file L-shaped.h with comments:
#include "libmesh/enum_fe_family.h"
#include "libmesh/fem_system.h"
#include "libmesh/qoi_set.h"
#include "libmesh/system.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
FEMSystem, TimeSolver and NewtonSolver will handle most tasks,
but we must specify element residuals
class LaplaceSystem : public FEMSystem
{
public:
Constructor
LaplaceSystem(EquationSystems& es,
const std::string& name_in,
const unsigned int number_in)
: FEMSystem(es, name_in, number_in),
_fe_family("LAGRANGE"), _fe_order(1),
_analytic_jacobians(true) { qoi.resize(2); }
std::string & fe_family() { return _fe_family; }
unsigned int & fe_order() { return _fe_order; }
bool & analytic_jacobians() { return _analytic_jacobians; }
Postprocessing function which we are going to override for this application
virtual void postprocess(void);
Number &get_QoI_value(std::string type, unsigned int QoI_index)
{
if(type == "exact")
{
return exact_QoI[QoI_index];
}
else
{
return computed_QoI[QoI_index];
}
}
protected:
System initialization
virtual void init_data ();
Context initialization
virtual void init_context (DiffContext &context);
Element residual and jacobian calculations
Time dependent parts
virtual bool element_time_derivative (bool request_jacobian,
DiffContext &context);
Constraint parts
virtual bool side_constraint (bool request_jacobian,
DiffContext &context);
Overloading the postprocess function
virtual void element_postprocess(DiffContext &context);
virtual void side_postprocess(DiffContext &context);
Overloading the qoi function on elements
virtual void element_qoi_derivative
(DiffContext &context,
const QoISet & qois);
Overloading the qoi function on sides
virtual void side_qoi_derivative
(DiffContext &context,
const QoISet & qois);
Number exact_solution (const Point&);
Variables to hold the computed QoIs
Number computed_QoI[2];
Variables to read in the exact QoIs from l-shaped.in
Number exact_QoI[2];
The FE type to use
std::string _fe_family;
unsigned int _fe_order;
Calculate Jacobians analytically or not?
bool _analytic_jacobians;
};
The source file adjoints_ex1.C with comments:
Adjoints Example 1 - Laplace Equation in the L-Shaped Domain with Adjoint based mesh refinement
This example solves the Laplace equation on the classic "L-shaped" domain with adaptive mesh refinement. The exact solution is u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta). The kelly and adjoint residual error estimators are used to develop error indicators and guide mesh adaptation. Since we use the adjoint capabilities of libMesh in this example, we use the DiffSystem framework. This file (adjoints_ex1.C) contains the declaration of mesh and equation system objects, L-shaped.C contains the assembly of the system, element_qoi_derivative.C and side_qoi_derivative.C contain the RHS for the adjoint systems. Postprocessing to compute the QoIs is done in element_postprocess.C and side_postprocess.C.
The initial mesh contains three QUAD9 elements which represent the standard quadrants I, II, and III of the domain [-1,1]x[-1,1], i.e. Element 0: [-1,0]x[ 0,1] Element 1: [ 0,1]x[ 0,1] Element 2: [-1,0]x[-1,0] The mesh is provided in the standard libMesh ASCII format file named "lshaped.xda". In addition, an input file named "general.in" is provided which allows the user to set several parameters for the solution so that the problem can be re-run without a re-compile. The solution technique employed is to have a refinement loop with a linear (forward and adjoint) solve inside followed by a refinement of the grid and projection of the solution to the new grid In the final loop iteration, there is no additional refinement after the solve. In the input file "general.in", the variable "max_adaptivesteps" controls the number of refinement steps, and "refine_fraction" / "coarsen_fraction" determine the number of elements which will be refined / coarsened at each step.
C++ includes
#include <iostream>
#include <iomanip>
General libMesh includes
#include "libmesh/equation_systems.h"
#include "libmesh/linear_solver.h"
#include "libmesh/error_vector.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/newton_solver.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/steady_solver.h"
#include "libmesh/system_norm.h"
Error Estimator includes
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/patch_recovery_error_estimator.h"
Adjoint Related includes
#include "libmesh/adjoint_residual_error_estimator.h"
#include "libmesh/qoi_set.h"
libMesh I/O includes
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
Local includes
#include "femparameters.h"
#include "L-shaped.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Local function declarations
Number output files, the files are give a prefix of primal or adjoint_i depending on whether the output is the primal solution or the dual solution for the ith QoI
Write gmv output
Number output files, the files are give a prefix of primal or adjoint_i depending on whether the output is the primal solution or the dual solution for the ith QoI
Write gmv output
void write_output(EquationSystems &es,
unsigned int a_step, // The adaptive step count
std::string solution_type) // primal or adjoint solve
{
MeshBase &mesh = es.get_mesh();
#ifdef LIBMESH_HAVE_GMV
std::ostringstream file_name_gmv;
file_name_gmv << solution_type
<< ".out.gmv."
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step;
GMVIO(mesh).write_equation_systems
(file_name_gmv.str(), es);
#endif
}
Set the parameters for the nonlinear and linear solvers to be used during the simulation
void set_system_parameters(LaplaceSystem &system, FEMParameters ¶m)
{
Use analytical jacobians?
system.analytic_jacobians() = param.analytic_jacobians;
Verify analytic jacobians against numerical ones?
system.verify_analytic_jacobians = param.verify_analytic_jacobians;
Use the prescribed FE type
system.fe_family() = param.fe_family[0];
system.fe_order() = param.fe_order[0];
More desperate debugging options
system.print_solution_norms = param.print_solution_norms;
system.print_solutions = param.print_solutions;
system.print_residual_norms = param.print_residual_norms;
system.print_residuals = param.print_residuals;
system.print_jacobian_norms = param.print_jacobian_norms;
system.print_jacobians = param.print_jacobians;
No transient time solver
system.time_solver =
AutoPtr<TimeSolver>(new SteadySolver(system));
Nonlinear solver options
{
NewtonSolver *solver = new NewtonSolver(system);
system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
solver->quiet = param.solver_quiet;
solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
solver->minsteplength = param.min_step_length;
solver->relative_step_tolerance = param.relative_step_tolerance;
solver->relative_residual_tolerance = param.relative_residual_tolerance;
solver->require_residual_reduction = param.require_residual_reduction;
solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
if (system.time_solver->reduce_deltat_on_diffsolver_failure)
{
solver->continue_after_max_iterations = true;
solver->continue_after_backtrack_failure = true;
}
And the linear solver options
solver->max_linear_iterations = param.max_linear_iterations;
solver->initial_linear_tolerance = param.initial_linear_tolerance;
solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
}
}
Build the mesh refinement object and set parameters for refining/coarsening etc
#ifdef LIBMESH_ENABLE_AMR
AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
FEMParameters ¶m)
{
AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
mesh_refinement->coarsen_by_parents() = true;
mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
mesh_refinement->nelem_target() = param.nelem_target;
mesh_refinement->refine_fraction() = param.refine_fraction;
mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
return mesh_refinement;
}
#endif // LIBMESH_ENABLE_AMR
This is where we declare the error estimators to be built and used for
mesh refinement. The adjoint residual estimator needs two estimators.
One for the forward component of the estimate and one for the adjoint
weighting factor. Here we use the Patch Recovery indicator to estimate both the
forward and adjoint weights. The H1 seminorm component of the error is used
as dictated by the weak form the Laplace equation.
AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters ¶m, QoISet &qois)
{
AutoPtr<ErrorEstimator> error_estimator;
if (param.indicator_type == "kelly")
{
std::cout<<"Using Kelly Error Estimator"<<std::endl;
error_estimator.reset(new KellyErrorEstimator);
}
else if (param.indicator_type == "adjoint_residual")
{
std::cout<<"Using Adjoint Residual Error Estimator with Patch Recovery Weights"<<std::endl;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
adjoint_residual_estimator->qoi_set() = qois;
adjoint_residual_estimator->error_plot_suffix = "error.gmv";
PatchRecoveryErrorEstimator *p1 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->primal_error_estimator().reset(p1);
PatchRecoveryErrorEstimator *p2 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->dual_error_estimator().reset(p2);
adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
p1->set_patch_reuse(param.patch_reuse);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
p2->set_patch_reuse(param.patch_reuse);
}
else
{
std::cerr << "Unknown indicator_type" << std::endl;
libmesh_error();
}
return error_estimator;
}
The main program.
int main (int argc, char** argv)
{
Initialize libMesh.
LibMeshInit init (argc, argv);
Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
std::cout << "Started " << argv[0] << std::endl;
Make sure the general input file exists, and parse it
{
std::ifstream i("general.in");
if (!i)
{
std::cerr << '[' << libMesh::processor_id()
<< "] Can't find general.in; exiting early."
<< std::endl;
libmesh_error();
}
}
GetPot infile("general.in");
Read in parameters from the input file
FEMParameters param;
param.read(infile);
Skip this default-2D example if libMesh was compiled as 1D-only.
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Create a mesh.
Mesh mesh;
And an object to refine it
AutoPtr<MeshRefinement> mesh_refinement =
build_mesh_refinement(mesh, param);
And an EquationSystems to run on it
EquationSystems equation_systems (mesh);
std::cout << "Reading in and building the mesh" << std::endl;
Read in the mesh
mesh.read(param.domainfile.c_str());
Make all the elements of the mesh second order so we can compute
with a higher order basis
mesh.all_second_order();
Create a mesh refinement object to do the initial uniform refinements
on the coarse grid read in from lshaped.xda
MeshRefinement initial_uniform_refinements(mesh);
initial_uniform_refinements.uniformly_refine(param.coarserefinements);
std::cout << "Building system" << std::endl;
Build the FEMSystem
LaplaceSystem &system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
Set its parameters
set_system_parameters(system, param);
std::cout << "Initializing systems" << std::endl;
equation_systems.init ();
Print information about the mesh and system to the screen.
mesh.print_info();
equation_systems.print_info();
LinearSolver<Number> *linear_solver = system.get_linear_solver();
{
Adaptively solve the timestep
unsigned int a_step = 0;
for (; a_step != param.max_adaptivesteps; ++a_step)
{
We can't adapt to both a tolerance and a
target mesh size
if (param.global_tolerance != 0.)
libmesh_assert_equal_to (param.nelem_target, 0);
If we aren't adapting to a tolerance we need a
target mesh size
else
libmesh_assert_greater (param.nelem_target, 0);
linear_solver->reuse_preconditioner(false);
Solve the forward problem
system.solve();
Write out the computed primal solution
write_output(equation_systems, a_step, "primal");
Get a pointer to the primal solution vector
NumericVector<Number> &primal_solution = *system.solution;
Declare a QoISet object, we need this object to set weights for our QoI error contributions
QoISet qois;
Declare a qoi_indices vector, each index will correspond to a QoI
std::vector<unsigned int> qoi_indices;
qoi_indices.push_back(0);
qoi_indices.push_back(1);
qois.add_indices(qoi_indices);
Set weights for each index, these will weight the contribution of each QoI in the final error
estimate to be used for flagging elements for refinement
qois.set_weight(0, 0.5);
qois.set_weight(1, 0.5);
Make sure we get the contributions to the adjoint RHS from the sides
system.assemble_qoi_sides = true;
We are about to solve the adjoint system, but before we do this we see the same preconditioner
flag to reuse the preconditioner from the forward solver
linear_solver->reuse_preconditioner(param.reuse_preconditioner);
Solve the adjoint system. This takes the transpose of the stiffness matrix and then
solves the resulting system
system.adjoint_solve();
Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unneccesarily in the error estimator
system.set_adjoint_already_solved(true);
Get a pointer to the solution vector of the adjoint problem for QoI 0
NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
Swap the primal and dual solutions so we can write out the adjoint solution
primal_solution.swap(dual_solution_0);
write_output(equation_systems, a_step, "adjoint_0");
Swap back
primal_solution.swap(dual_solution_0);
Get a pointer to the solution vector of the adjoint problem for QoI 0
NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
Swap again
primal_solution.swap(dual_solution_1);
write_output(equation_systems, a_step, "adjoint_1");
Swap back again
primal_solution.swap(dual_solution_1);
std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs()
<< " active dofs." << std::endl ;
Postprocess, compute the approximate QoIs and write them out to the console
std::cout << "Postprocessing: " << std::endl;
system.postprocess_sides = true;
system.postprocess();
Number QoI_0_computed = system.get_QoI_value("computed", 0);
Number QoI_0_exact = system.get_QoI_value("exact", 0);
Number QoI_1_computed = system.get_QoI_value("computed", 1);
Number QoI_1_exact = system.get_QoI_value("exact", 1);
std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(QoI_0_computed - QoI_0_exact) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(QoI_1_computed - QoI_1_exact) /
std::abs(QoI_1_exact) << std::endl << std::endl;
Now we construct the data structures for the mesh refinement process
ErrorVector error;
Build an error estimator object
AutoPtr<ErrorEstimator> error_estimator =
build_error_estimator(param, qois);
Estimate the error in each element using the Adjoint Residual or Kelly error estimator
error_estimator->estimate_error(system, error);
We have to refine either based on reaching an error tolerance or
a number of elements target, which should be verified above
Otherwise we flag elements by error tolerance or nelem target
Uniform refinement
Uniform refinement
if(param.refine_uniformly)
{
mesh_refinement->uniformly_refine(1);
}
Adaptively refine based on reaching an error tolerance
else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
{
mesh_refinement->flag_elements_by_error_tolerance (error);
mesh_refinement->refine_and_coarsen_elements();
}
Adaptively refine based on reaching a target number of elements
else
{
if (mesh.n_active_elem() >= param.nelem_target)
{
std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
break;
}
mesh_refinement->flag_elements_by_nelem_target (error);
mesh_refinement->refine_and_coarsen_elements();
}
Dont forget to reinit the system after each adaptive refinement !
equation_systems.reinit();
std::cout << "Refined mesh to "
<< mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs()
<< " active dofs." << std::endl;
}
Do one last solve if necessary
if (a_step == param.max_adaptivesteps)
{
linear_solver->reuse_preconditioner(false);
system.solve();
write_output(equation_systems, a_step, "primal");
NumericVector<Number> &primal_solution = *system.solution;
QoISet qois;
std::vector<unsigned int> qoi_indices;
qoi_indices.push_back(0);
qoi_indices.push_back(1);
qois.add_indices(qoi_indices);
qois.set_weight(0, 0.5);
qois.set_weight(1, 0.5);
system.assemble_qoi_sides = true;
linear_solver->reuse_preconditioner(param.reuse_preconditioner);
system.adjoint_solve();
Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unneccesarily in the error estimator
system.set_adjoint_already_solved(true);
NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
primal_solution.swap(dual_solution_0);
write_output(equation_systems, a_step, "adjoint_0");
primal_solution.swap(dual_solution_0);
NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
primal_solution.swap(dual_solution_1);
write_output(equation_systems, a_step, "adjoint_1");
primal_solution.swap(dual_solution_1);
std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs()
<< " active dofs." << std::endl ;
std::cout << "Postprocessing: " << std::endl;
system.postprocess_sides = true;
system.postprocess();
Number QoI_0_computed = system.get_QoI_value("computed", 0);
Number QoI_0_exact = system.get_QoI_value("exact", 0);
Number QoI_1_computed = system.get_QoI_value("computed", 1);
Number QoI_1_exact = system.get_QoI_value("exact", 1);
std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(QoI_0_computed - QoI_0_exact) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(QoI_1_computed - QoI_1_exact) /
std::abs(QoI_1_exact) << std::endl << std::endl;
}
}
std::cerr << '[' << libMesh::processor_id()
<< "] Completing output." << std::endl;
#endif // #ifndef LIBMESH_ENABLE_AMR
All done.
return 0;
}
The source file element_postprocess.C with comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
Local includes
#include "L-shaped.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Last modified: Feb 13, 2009
Define the postprocess function to compute QoI 0, the integral of the the solution over a subdomain
Define the postprocess function to compute QoI 0, the integral of the the solution over a subdomain
void LaplaceSystem::element_postprocess (DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
const std::vector<Point> &xyz = c.element_fe_var[0]->get_xyz();
The number of local degrees of freedom in each variable
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
The function R = int_{omega} T dR
omega is a subset of Omega (the whole domain), omega = [0.75, 1.0] x [0.0, 0.25]
Number dQoI_0 = 0.;
Loop over quadrature points
for (unsigned int qp = 0; qp != n_qpoints; qp++)
{
Get co-ordinate locations of the current quadrature point
const Real x = xyz[qp](0);
const Real y = xyz[qp](1);
If in the sub-domain omega, add the contribution to the integral R
if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125)
{
Get the solution value at the quadrature point
Number T = c.interior_value(0, qp);
Update the elemental increment dR for each qp
dQoI_0 += JxW[qp] * T;
}
}
Update the computed value of the global functional R, by adding the contribution from this element
computed_QoI[0] = computed_QoI[0] + dQoI_0;
}
The source file element_qoi_derivative.C with comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
Local includes
#include "L-shaped.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
We only have one QoI, so we don't bother checking the qois argument
to see if it was requested from us
void LaplaceSystem::element_qoi_derivative (DiffContext &context,
const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
The basis functions for the element
const std::vector<std::vector<Real> > &phi = c.element_fe_var[0]->get_phi();
The element quadrature points
const std::vector<Point > &q_point = c.element_fe_var[0]->get_xyz();
The number of local degrees of freedom in each variable
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
Fill the QoI RHS corresponding to this QoI. Since this is the 0th QoI
we fill in the [0][i] subderivatives, i corresponding to the variable index.
Our system has only one variable, so we only have to fill the [0][0] subderivative
DenseSubVector<Number> &Q = *c.elem_qoi_subderivatives[0][0];
Loop over the qps
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
If in the sub-domain over which QoI 0 is supported, add contributions
to the adjoint rhs
if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125)
{
for (unsigned int i=0; i != n_T_dofs; i++)
Q(i) += JxW[qp] *phi[i][qp] ;
}
} // end of the quadrature point qp-loop
}
The source file femparameters.C with comments:
#include "femparameters.h"
#define GETPOT_INPUT(A) { A = input(#A, A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_INT_INPUT(A) { A = input(#A, (int)A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_REGISTER(A) { \
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
void FEMParameters::read(GetPot &input)
{
std::vector<std::string> variable_names;
GETPOT_INT_INPUT(coarserefinements);
GETPOT_INPUT(domainfile);
GETPOT_INPUT(solver_quiet);
GETPOT_INPUT(reuse_preconditioner);
GETPOT_INPUT(require_residual_reduction);
GETPOT_INPUT(min_step_length);
GETPOT_INT_INPUT(max_linear_iterations);
GETPOT_INT_INPUT(max_nonlinear_iterations);
GETPOT_INPUT(relative_step_tolerance);
GETPOT_INPUT(relative_residual_tolerance);
GETPOT_INPUT(initial_linear_tolerance);
GETPOT_INPUT(minimum_linear_tolerance);
GETPOT_INPUT(linear_tolerance_multiplier);
GETPOT_INT_INPUT(nelem_target);
GETPOT_INPUT(global_tolerance);
GETPOT_INPUT(refine_fraction);
GETPOT_INPUT(coarsen_fraction);
GETPOT_INPUT(coarsen_threshold);
GETPOT_INT_INPUT(max_adaptivesteps);
GETPOT_INPUT(refine_uniformly);
GETPOT_INPUT(indicator_type);
GETPOT_INPUT(patch_reuse);
GETPOT_REGISTER(fe_family);
const unsigned int n_fe_family =
std::max(1u, input.vector_variable_size("fe_family"));
fe_family.resize(n_fe_family, "LAGRANGE");
for (unsigned int i=0; i != n_fe_family; ++i)
fe_family[i] = input("fe_family", fe_family[i].c_str(), i);
GETPOT_REGISTER(fe_order);
const unsigned int n_fe_order =
input.vector_variable_size("fe_order");
fe_order.resize(n_fe_order, 1);
for (unsigned int i=0; i != n_fe_order; ++i)
fe_order[i] = input("fe_order", (int)fe_order[i], i);
GETPOT_INPUT(analytic_jacobians);
GETPOT_INPUT(verify_analytic_jacobians);
GETPOT_INPUT(print_solution_norms);
GETPOT_INPUT(print_solutions);
GETPOT_INPUT(print_residual_norms);
GETPOT_INPUT(print_residuals);
GETPOT_INPUT(print_jacobian_norms);
GETPOT_INPUT(print_jacobians);
std::vector<std::string> bad_variables =
input.unidentified_arguments(variable_names);
if (libMesh::processor_id() == 0 && !bad_variables.empty())
{
std::cerr << "ERROR: Unrecognized variables:" << std::endl;
for (unsigned int i = 0; i != bad_variables.size(); ++i)
std::cerr << bad_variables[i] << std::endl;
std::cerr << "not found among recognized variables." << std::endl;
for (unsigned int i = 0; i != variable_names.size(); ++i)
std::cerr << variable_names[i] << std::endl;
libmesh_error();
}
}
The source file L-shaped.C with comments:
#include "libmesh/getpot.h"
#include "libmesh/fe_base.h"
#include "libmesh/quadrature.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/parallel.h"
#include "libmesh/fem_context.h"
Local includes
#include "L-shaped.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
void LaplaceSystem::init_data ()
{
this->add_variable ("T", static_cast<Order>(_fe_order),
Utility::string_to_enum<FEFamily>(_fe_family));
GetPot infile("l-shaped.in");
exact_QoI[0] = infile("QoI_0", 0.0);
exact_QoI[1] = infile("QoI_1", 0.0);
Do the parent's initialization after variables are defined
FEMSystem::init_data();
this->time_evolving(0);
}
void LaplaceSystem::init_context(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
Now make sure we have requested all the data
we need to build the linear system.
c.element_fe_var[0]->get_JxW();
c.element_fe_var[0]->get_phi();
c.element_fe_var[0]->get_dphi();
c.side_fe_var[0]->get_JxW();
c.side_fe_var[0]->get_phi();
c.side_fe_var[0]->get_dphi();
}
#define optassert(X) {if (!(X)) libmesh_error();}
Assemble the element contributions to the stiffness matrix
bool LaplaceSystem::element_time_derivative (bool request_jacobian,
DiffContext &context)
{
Are the jacobians specified analytically ?
bool compute_jacobian = request_jacobian && _analytic_jacobians;
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
Element basis functions
const std::vector<std::vector<RealGradient> > &dphi = c.element_fe_var[0]->get_dphi();
The number of local degrees of freedom in each variable
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
The subvectors and submatrices we need to fill:
DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
DenseSubVector<Number> &F = *c.elem_subresiduals[0];
Now we will build the element Jacobian and residual.
Constructing the residual requires the solution and its
gradient from the previous timestep. This must be
calculated at each quadrature point by summing the
solution degree-of-freedom values by the appropriate
weight functions.
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Compute the solution gradient at the Newton iterate
Gradient grad_T = c.interior_gradient(0, qp);
The residual contribution from this element
for (unsigned int i=0; i != n_T_dofs; i++)
F(i) += JxW[qp] * ( grad_T * dphi[i][qp] ) ;
if (compute_jacobian)
for (unsigned int i=0; i != n_T_dofs; i++)
for (unsigned int j=0; j != n_T_dofs; ++j)
The analytic jacobian
K(i,j) += JxW[qp] * ( dphi[i][qp] * dphi[j][qp] );
} // end of the quadrature point qp-loop
return compute_jacobian;
}
Set Dirichlet bcs, side contributions to global stiffness matrix
bool LaplaceSystem::side_constraint (bool request_jacobian,
DiffContext &context)
{
Are the jacobians specified analytically ?
bool compute_jacobian = request_jacobian && _analytic_jacobians;
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
Side basis functions
const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
Side Quadrature points
const std::vector<Point > &qside_point = c.side_fe_var[0]->get_xyz();
The number of local degrees of freedom in each variable
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
The subvectors and submatrices we need to fill:
DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
DenseSubVector<Number> &F = *c.elem_subresiduals[0];
unsigned int n_qpoints = (c.get_side_qrule())->n_points();
const Real penalty = 1./(TOLERANCE*TOLERANCE);
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Compute the solution at the old Newton iterate
Number T = c.side_value(0, qp);
We get the Dirichlet bcs from the exact solution
Number u_dirichlet = exact_solution (qside_point[qp]);
The residual from the boundary terms, penalize non-zero temperature
for (unsigned int i=0; i != n_T_dofs; i++)
F(i) += JxW[qp] * penalty * ( T - u_dirichlet) * phi[i][qp];
if (compute_jacobian)
for (unsigned int i=0; i != n_T_dofs; i++)
for (unsigned int j=0; j != n_T_dofs; ++j)
The analytic jacobian
K(i,j) += JxW[qp] * penalty * phi[i][qp] * phi[j][qp];
} // end of the quadrature point qp-loop
return compute_jacobian;
}
Override the default DiffSystem postprocess function to compute the
approximations to the QoIs
void LaplaceSystem::postprocess()
{
Reset the array holding the computed QoIs
computed_QoI[0] = 0.0;
computed_QoI[1] = 0.0;
FEMSystem::postprocess();
Parallel::sum(computed_QoI[0]);
Parallel::sum(computed_QoI[1]);
}
The exact solution to the singular problem,
u_exact = r^(2/3)*sin(2*theta/3). We use this to set the Dirichlet boundary conditions
Number LaplaceSystem::exact_solution(const Point& p)// xyz location
{
const Real x1 = p(0);
const Real x2 = p(1);
Real theta = atan2(x2,x1);
Make sure 0 <= theta <= 2*pi
if (theta < 0)
theta += 2. * libMesh::pi;
return pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
}
The source file side_postprocess.C with comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
Local includes
#include "L-shaped.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
Define the postprocess function to compute QoI 1, the integral of the the normal
derivative of the solution over part of the boundary
void LaplaceSystem::side_postprocess(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
unsigned int n_qpoints = (c.get_side_qrule())->n_points();
Number dQoI_1 = 0. ;
Loop over qp's, compute the function at each qp and add
to get the QoI
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
const Real TOL = 1.e-5;
If on the bottom horizontal bdry (y = -1)
if(fabs(y - 1.0) <= TOL && x > 0.0)
{
Get the value of the gradient at this point
const Gradient grad_T = c.side_gradient(0,qp);
Add the contribution of this qp to the integral QoI
dQoI_1 += JxW[qp] * (grad_T * face_normals[qp]) ;
}
} // end of the quadrature point qp-loop
computed_QoI[1] = computed_QoI[1] + dQoI_1;
}
The source file side_qoi_derivative.C with comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
Local includes
#include "L-shaped.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
We only have one QoI, so we don't bother checking the qois argument
to see if it was requested from us
void LaplaceSystem::side_qoi_derivative (DiffContext &context,
const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
First we get some references to cell-specific data that
will be used to assemble the linear system.
Element Jacobian * quadrature weights for interior integration
Element Jacobian * quadrature weights for interior integration
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
The basis functions for the side
const std::vector<std::vector<RealGradient> > &dphi = c.side_fe_var[0]->get_dphi();
The side quadrature points
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
Get the normal to the side at each qp
const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
The number of local degrees of freedom in each variable
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
unsigned int n_qpoints = (c.get_side_qrule())->n_points();
Fill the QoI RHS corresponding to this QoI. Since this is QoI 1
we fill in the [1][i] subderivatives, i corresponding to the variable index.
Our system has only one variable, so we only have to fill the [1][0] subderivative
DenseSubVector<Number> &Q = *c.elem_qoi_subderivatives[1][0];
const Real TOL = 1.e-5;
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
If on the sides where the boundary QoI is supported, add contributions
to the adjoint rhs
if(fabs(y - 1.0) <= TOL && x > 0.0)
{
for (unsigned int i=0; i != n_T_dofs; i++)
Q(i) += JxW[qp] * (dphi[i][qp] * face_normals[qp]);
}
} // end of the quadrature point qp-loop
}
The source file femparameters.h without comments:
#ifndef __fem_parameters_h__
#define __fem_parameters_h__
#include <limits>
#include <string>
#include "libmesh/libmesh_common.h"
#include "libmesh/getpot.h"
using namespace libMesh;
class FEMParameters
{
public:
FEMParameters() :
domainfile("lshaped.xda"),
coarserefinements(0),
solver_quiet(true),
reuse_preconditioner(false),
require_residual_reduction(true),
min_step_length(1e-5),
max_linear_iterations(200000), max_nonlinear_iterations(20),
relative_step_tolerance(1.e-7), relative_residual_tolerance(1.e-10),
initial_linear_tolerance(1.e-3), minimum_linear_tolerance(TOLERANCE*TOLERANCE),
linear_tolerance_multiplier(1.e-3),
nelem_target(30000), global_tolerance(0.0),
refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
refine_uniformly(false),
max_adaptivesteps(1),
indicator_type("kelly"), patch_reuse(true),
fe_family(1, "LAGRANGE"), fe_order(1, 1),
analytic_jacobians(true), verify_analytic_jacobians(0.0),
print_solution_norms(false), print_solutions(false),
print_residual_norms(false), print_residuals(false),
print_jacobian_norms(false), print_jacobians(false) {}
void read(GetPot &input);
std::string domainfile;
unsigned int coarserefinements;
bool solver_quiet, reuse_preconditioner, require_residual_reduction;
Real min_step_length;
unsigned int max_linear_iterations, max_nonlinear_iterations;
Real relative_step_tolerance, relative_residual_tolerance,
initial_linear_tolerance, minimum_linear_tolerance,
linear_tolerance_multiplier;
unsigned int nelem_target;
Real global_tolerance;
Real refine_fraction, coarsen_fraction, coarsen_threshold;
bool refine_uniformly;
unsigned int max_adaptivesteps;
std::string indicator_type;
bool patch_reuse;
std::vector<std::string> fe_family;
std::vector<unsigned int> fe_order;
bool analytic_jacobians;
Real verify_analytic_jacobians;
bool print_solution_norms, print_solutions,
print_residual_norms, print_residuals,
print_jacobian_norms, print_jacobians;
};
#endif // __fem_parameters_h__
The source file L-shaped.h without comments:
#include "libmesh/enum_fe_family.h"
#include "libmesh/fem_system.h"
#include "libmesh/qoi_set.h"
#include "libmesh/system.h"
using namespace libMesh;
class LaplaceSystem : public FEMSystem
{
public:
LaplaceSystem(EquationSystems& es,
const std::string& name_in,
const unsigned int number_in)
: FEMSystem(es, name_in, number_in),
_fe_family("LAGRANGE"), _fe_order(1),
_analytic_jacobians(true) { qoi.resize(2); }
std::string & fe_family() { return _fe_family; }
unsigned int & fe_order() { return _fe_order; }
bool & analytic_jacobians() { return _analytic_jacobians; }
virtual void postprocess(void);
Number &get_QoI_value(std::string type, unsigned int QoI_index)
{
if(type == "exact")
{
return exact_QoI[QoI_index];
}
else
{
return computed_QoI[QoI_index];
}
}
protected:
virtual void init_data ();
virtual void init_context (DiffContext &context);
virtual bool element_time_derivative (bool request_jacobian,
DiffContext &context);
virtual bool side_constraint (bool request_jacobian,
DiffContext &context);
virtual void element_postprocess(DiffContext &context);
virtual void side_postprocess(DiffContext &context);
virtual void element_qoi_derivative
(DiffContext &context,
const QoISet & qois);
virtual void side_qoi_derivative
(DiffContext &context,
const QoISet & qois);
Number exact_solution (const Point&);
Number computed_QoI[2];
Number exact_QoI[2];
std::string _fe_family;
unsigned int _fe_order;
bool _analytic_jacobians;
};
The source file adjoints_ex1.C without comments:
#include <iostream>
#include <iomanip>
#include "libmesh/equation_systems.h"
#include "libmesh/linear_solver.h"
#include "libmesh/error_vector.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/newton_solver.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/steady_solver.h"
#include "libmesh/system_norm.h"
#include "libmesh/kelly_error_estimator.h"
#include "libmesh/patch_recovery_error_estimator.h"
#include "libmesh/adjoint_residual_error_estimator.h"
#include "libmesh/qoi_set.h"
#include "libmesh/getpot.h"
#include "libmesh/gmv_io.h"
#include "femparameters.h"
#include "L-shaped.h"
using namespace libMesh;
void write_output(EquationSystems &es,
unsigned int a_step, // The adaptive step count
std::string solution_type) // primal or adjoint solve
{
MeshBase &mesh = es.get_mesh();
#ifdef LIBMESH_HAVE_GMV
std::ostringstream file_name_gmv;
file_name_gmv << solution_type
<< ".out.gmv."
<< std::setw(2)
<< std::setfill('0')
<< std::right
<< a_step;
GMVIO(mesh).write_equation_systems
(file_name_gmv.str(), es);
#endif
}
void set_system_parameters(LaplaceSystem &system, FEMParameters ¶m)
{
system.analytic_jacobians() = param.analytic_jacobians;
system.verify_analytic_jacobians = param.verify_analytic_jacobians;
system.fe_family() = param.fe_family[0];
system.fe_order() = param.fe_order[0];
system.print_solution_norms = param.print_solution_norms;
system.print_solutions = param.print_solutions;
system.print_residual_norms = param.print_residual_norms;
system.print_residuals = param.print_residuals;
system.print_jacobian_norms = param.print_jacobian_norms;
system.print_jacobians = param.print_jacobians;
system.time_solver =
AutoPtr<TimeSolver>(new SteadySolver(system));
{
NewtonSolver *solver = new NewtonSolver(system);
system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
solver->quiet = param.solver_quiet;
solver->max_nonlinear_iterations = param.max_nonlinear_iterations;
solver->minsteplength = param.min_step_length;
solver->relative_step_tolerance = param.relative_step_tolerance;
solver->relative_residual_tolerance = param.relative_residual_tolerance;
solver->require_residual_reduction = param.require_residual_reduction;
solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
if (system.time_solver->reduce_deltat_on_diffsolver_failure)
{
solver->continue_after_max_iterations = true;
solver->continue_after_backtrack_failure = true;
}
solver->max_linear_iterations = param.max_linear_iterations;
solver->initial_linear_tolerance = param.initial_linear_tolerance;
solver->minimum_linear_tolerance = param.minimum_linear_tolerance;
}
}
#ifdef LIBMESH_ENABLE_AMR
AutoPtr<MeshRefinement> build_mesh_refinement(MeshBase &mesh,
FEMParameters ¶m)
{
AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
mesh_refinement->coarsen_by_parents() = true;
mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
mesh_refinement->nelem_target() = param.nelem_target;
mesh_refinement->refine_fraction() = param.refine_fraction;
mesh_refinement->coarsen_fraction() = param.coarsen_fraction;
mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
return mesh_refinement;
}
#endif // LIBMESH_ENABLE_AMR
AutoPtr<ErrorEstimator> build_error_estimator(FEMParameters ¶m, QoISet &qois)
{
AutoPtr<ErrorEstimator> error_estimator;
if (param.indicator_type == "kelly")
{
std::cout<<"Using Kelly Error Estimator"<<std::endl;
error_estimator.reset(new KellyErrorEstimator);
}
else if (param.indicator_type == "adjoint_residual")
{
std::cout<<"Using Adjoint Residual Error Estimator with Patch Recovery Weights"<<std::endl;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
adjoint_residual_estimator->qoi_set() = qois;
adjoint_residual_estimator->error_plot_suffix = "error.gmv";
PatchRecoveryErrorEstimator *p1 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->primal_error_estimator().reset(p1);
PatchRecoveryErrorEstimator *p2 =
new PatchRecoveryErrorEstimator;
adjoint_residual_estimator->dual_error_estimator().reset(p2);
adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
p1->set_patch_reuse(param.patch_reuse);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
p2->set_patch_reuse(param.patch_reuse);
}
else
{
std::cerr << "Unknown indicator_type" << std::endl;
libmesh_error();
}
return error_estimator;
}
int main (int argc, char** argv)
{
LibMeshInit init (argc, argv);
#ifndef LIBMESH_ENABLE_AMR
libmesh_example_assert(false, "--enable-amr");
#else
std::cout << "Started " << argv[0] << std::endl;
{
std::ifstream i("general.in");
if (!i)
{
std::cerr << '[' << libMesh::processor_id()
<< "] Can't find general.in; exiting early."
<< std::endl;
libmesh_error();
}
}
GetPot infile("general.in");
FEMParameters param;
param.read(infile);
libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
Mesh mesh;
AutoPtr<MeshRefinement> mesh_refinement =
build_mesh_refinement(mesh, param);
EquationSystems equation_systems (mesh);
std::cout << "Reading in and building the mesh" << std::endl;
mesh.read(param.domainfile.c_str());
mesh.all_second_order();
MeshRefinement initial_uniform_refinements(mesh);
initial_uniform_refinements.uniformly_refine(param.coarserefinements);
std::cout << "Building system" << std::endl;
LaplaceSystem &system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
set_system_parameters(system, param);
std::cout << "Initializing systems" << std::endl;
equation_systems.init ();
mesh.print_info();
equation_systems.print_info();
LinearSolver<Number> *linear_solver = system.get_linear_solver();
{
unsigned int a_step = 0;
for (; a_step != param.max_adaptivesteps; ++a_step)
{
if (param.global_tolerance != 0.)
libmesh_assert_equal_to (param.nelem_target, 0);
else
libmesh_assert_greater (param.nelem_target, 0);
linear_solver->reuse_preconditioner(false);
system.solve();
write_output(equation_systems, a_step, "primal");
NumericVector<Number> &primal_solution = *system.solution;
QoISet qois;
std::vector<unsigned int> qoi_indices;
qoi_indices.push_back(0);
qoi_indices.push_back(1);
qois.add_indices(qoi_indices);
qois.set_weight(0, 0.5);
qois.set_weight(1, 0.5);
system.assemble_qoi_sides = true;
linear_solver->reuse_preconditioner(param.reuse_preconditioner);
system.adjoint_solve();
system.set_adjoint_already_solved(true);
NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
primal_solution.swap(dual_solution_0);
write_output(equation_systems, a_step, "adjoint_0");
primal_solution.swap(dual_solution_0);
NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
primal_solution.swap(dual_solution_1);
write_output(equation_systems, a_step, "adjoint_1");
primal_solution.swap(dual_solution_1);
std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs()
<< " active dofs." << std::endl ;
std::cout << "Postprocessing: " << std::endl;
system.postprocess_sides = true;
system.postprocess();
Number QoI_0_computed = system.get_QoI_value("computed", 0);
Number QoI_0_exact = system.get_QoI_value("exact", 0);
Number QoI_1_computed = system.get_QoI_value("computed", 1);
Number QoI_1_exact = system.get_QoI_value("exact", 1);
std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(QoI_0_computed - QoI_0_exact) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(QoI_1_computed - QoI_1_exact) /
std::abs(QoI_1_exact) << std::endl << std::endl;
ErrorVector error;
AutoPtr<ErrorEstimator> error_estimator =
build_error_estimator(param, qois);
error_estimator->estimate_error(system, error);
if(param.refine_uniformly)
{
mesh_refinement->uniformly_refine(1);
}
else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
{
mesh_refinement->flag_elements_by_error_tolerance (error);
mesh_refinement->refine_and_coarsen_elements();
}
else
{
if (mesh.n_active_elem() >= param.nelem_target)
{
std::cout<<"We reached the target number of elements."<<std::endl <<std::endl;
break;
}
mesh_refinement->flag_elements_by_nelem_target (error);
mesh_refinement->refine_and_coarsen_elements();
}
equation_systems.reinit();
std::cout << "Refined mesh to "
<< mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs()
<< " active dofs." << std::endl;
}
if (a_step == param.max_adaptivesteps)
{
linear_solver->reuse_preconditioner(false);
system.solve();
write_output(equation_systems, a_step, "primal");
NumericVector<Number> &primal_solution = *system.solution;
QoISet qois;
std::vector<unsigned int> qoi_indices;
qoi_indices.push_back(0);
qoi_indices.push_back(1);
qois.add_indices(qoi_indices);
qois.set_weight(0, 0.5);
qois.set_weight(1, 0.5);
system.assemble_qoi_sides = true;
linear_solver->reuse_preconditioner(param.reuse_preconditioner);
system.adjoint_solve();
system.set_adjoint_already_solved(true);
NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
primal_solution.swap(dual_solution_0);
write_output(equation_systems, a_step, "adjoint_0");
primal_solution.swap(dual_solution_0);
NumericVector<Number> &dual_solution_1 = system.get_adjoint_solution(1);
primal_solution.swap(dual_solution_1);
write_output(equation_systems, a_step, "adjoint_1");
primal_solution.swap(dual_solution_1);
std::cout << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
<< " active elements and "
<< equation_systems.n_active_dofs()
<< " active dofs." << std::endl ;
std::cout << "Postprocessing: " << std::endl;
system.postprocess_sides = true;
system.postprocess();
Number QoI_0_computed = system.get_QoI_value("computed", 0);
Number QoI_0_exact = system.get_QoI_value("exact", 0);
Number QoI_1_computed = system.get_QoI_value("computed", 1);
Number QoI_1_exact = system.get_QoI_value("exact", 1);
std::cout<< "The relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(QoI_0_computed - QoI_0_exact) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(QoI_1_computed - QoI_1_exact) /
std::abs(QoI_1_exact) << std::endl << std::endl;
}
}
std::cerr << '[' << libMesh::processor_id()
<< "] Completing output." << std::endl;
#endif // #ifndef LIBMESH_ENABLE_AMR
return 0;
}
The source file element_postprocess.C without comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
#include "L-shaped.h"
using namespace libMesh;
void LaplaceSystem::element_postprocess (DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
const std::vector<Point> &xyz = c.element_fe_var[0]->get_xyz();
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
Number dQoI_0 = 0.;
for (unsigned int qp = 0; qp != n_qpoints; qp++)
{
const Real x = xyz[qp](0);
const Real y = xyz[qp](1);
if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125)
{
Number T = c.interior_value(0, qp);
dQoI_0 += JxW[qp] * T;
}
}
computed_QoI[0] = computed_QoI[0] + dQoI_0;
}
The source file element_qoi_derivative.C without comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
#include "L-shaped.h"
using namespace libMesh;
void LaplaceSystem::element_qoi_derivative (DiffContext &context,
const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
const std::vector<std::vector<Real> > &phi = c.element_fe_var[0]->get_phi();
const std::vector<Point > &q_point = c.element_fe_var[0]->get_xyz();
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
DenseSubVector<Number> &Q = *c.elem_qoi_subderivatives[0][0];
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
if(fabs(x - 0.875) <= 0.125 && fabs(y - 0.125) <= 0.125)
{
for (unsigned int i=0; i != n_T_dofs; i++)
Q(i) += JxW[qp] *phi[i][qp] ;
}
} // end of the quadrature point qp-loop
}
The source file femparameters.C without comments:
#include "femparameters.h"
#define GETPOT_INPUT(A) { A = input(#A, A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_INT_INPUT(A) { A = input(#A, (int)A);\
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
#define GETPOT_REGISTER(A) { \
std::string stringval = input(#A, std::string());\
variable_names.push_back(std::string(#A "=") + stringval); };
void FEMParameters::read(GetPot &input)
{
std::vector<std::string> variable_names;
GETPOT_INT_INPUT(coarserefinements);
GETPOT_INPUT(domainfile);
GETPOT_INPUT(solver_quiet);
GETPOT_INPUT(reuse_preconditioner);
GETPOT_INPUT(require_residual_reduction);
GETPOT_INPUT(min_step_length);
GETPOT_INT_INPUT(max_linear_iterations);
GETPOT_INT_INPUT(max_nonlinear_iterations);
GETPOT_INPUT(relative_step_tolerance);
GETPOT_INPUT(relative_residual_tolerance);
GETPOT_INPUT(initial_linear_tolerance);
GETPOT_INPUT(minimum_linear_tolerance);
GETPOT_INPUT(linear_tolerance_multiplier);
GETPOT_INT_INPUT(nelem_target);
GETPOT_INPUT(global_tolerance);
GETPOT_INPUT(refine_fraction);
GETPOT_INPUT(coarsen_fraction);
GETPOT_INPUT(coarsen_threshold);
GETPOT_INT_INPUT(max_adaptivesteps);
GETPOT_INPUT(refine_uniformly);
GETPOT_INPUT(indicator_type);
GETPOT_INPUT(patch_reuse);
GETPOT_REGISTER(fe_family);
const unsigned int n_fe_family =
std::max(1u, input.vector_variable_size("fe_family"));
fe_family.resize(n_fe_family, "LAGRANGE");
for (unsigned int i=0; i != n_fe_family; ++i)
fe_family[i] = input("fe_family", fe_family[i].c_str(), i);
GETPOT_REGISTER(fe_order);
const unsigned int n_fe_order =
input.vector_variable_size("fe_order");
fe_order.resize(n_fe_order, 1);
for (unsigned int i=0; i != n_fe_order; ++i)
fe_order[i] = input("fe_order", (int)fe_order[i], i);
GETPOT_INPUT(analytic_jacobians);
GETPOT_INPUT(verify_analytic_jacobians);
GETPOT_INPUT(print_solution_norms);
GETPOT_INPUT(print_solutions);
GETPOT_INPUT(print_residual_norms);
GETPOT_INPUT(print_residuals);
GETPOT_INPUT(print_jacobian_norms);
GETPOT_INPUT(print_jacobians);
std::vector<std::string> bad_variables =
input.unidentified_arguments(variable_names);
if (libMesh::processor_id() == 0 && !bad_variables.empty())
{
std::cerr << "ERROR: Unrecognized variables:" << std::endl;
for (unsigned int i = 0; i != bad_variables.size(); ++i)
std::cerr << bad_variables[i] << std::endl;
std::cerr << "not found among recognized variables." << std::endl;
for (unsigned int i = 0; i != variable_names.size(); ++i)
std::cerr << variable_names[i] << std::endl;
libmesh_error();
}
}
The source file L-shaped.C without comments:
#include "libmesh/getpot.h"
#include "libmesh/fe_base.h"
#include "libmesh/quadrature.h"
#include "libmesh/string_to_enum.h"
#include "libmesh/parallel.h"
#include "libmesh/fem_context.h"
#include "L-shaped.h"
using namespace libMesh;
void LaplaceSystem::init_data ()
{
this->add_variable ("T", static_cast<Order>(_fe_order),
Utility::string_to_enum<FEFamily>(_fe_family));
GetPot infile("l-shaped.in");
exact_QoI[0] = infile("QoI_0", 0.0);
exact_QoI[1] = infile("QoI_1", 0.0);
FEMSystem::init_data();
this->time_evolving(0);
}
void LaplaceSystem::init_context(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
c.element_fe_var[0]->get_JxW();
c.element_fe_var[0]->get_phi();
c.element_fe_var[0]->get_dphi();
c.side_fe_var[0]->get_JxW();
c.side_fe_var[0]->get_phi();
c.side_fe_var[0]->get_dphi();
}
#define optassert(X) {if (!(X)) libmesh_error();}
bool LaplaceSystem::element_time_derivative (bool request_jacobian,
DiffContext &context)
{
bool compute_jacobian = request_jacobian && _analytic_jacobians;
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
const std::vector<std::vector<RealGradient> > &dphi = c.element_fe_var[0]->get_dphi();
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
DenseSubVector<Number> &F = *c.elem_subresiduals[0];
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Gradient grad_T = c.interior_gradient(0, qp);
for (unsigned int i=0; i != n_T_dofs; i++)
F(i) += JxW[qp] * ( grad_T * dphi[i][qp] ) ;
if (compute_jacobian)
for (unsigned int i=0; i != n_T_dofs; i++)
for (unsigned int j=0; j != n_T_dofs; ++j)
K(i,j) += JxW[qp] * ( dphi[i][qp] * dphi[j][qp] );
} // end of the quadrature point qp-loop
return compute_jacobian;
}
bool LaplaceSystem::side_constraint (bool request_jacobian,
DiffContext &context)
{
bool compute_jacobian = request_jacobian && _analytic_jacobians;
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<std::vector<Real> > &phi = c.side_fe_var[0]->get_phi();
const std::vector<Point > &qside_point = c.side_fe_var[0]->get_xyz();
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
DenseSubMatrix<Number> &K = *c.elem_subjacobians[0][0];
DenseSubVector<Number> &F = *c.elem_subresiduals[0];
unsigned int n_qpoints = (c.get_side_qrule())->n_points();
const Real penalty = 1./(TOLERANCE*TOLERANCE);
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
Number T = c.side_value(0, qp);
Number u_dirichlet = exact_solution (qside_point[qp]);
for (unsigned int i=0; i != n_T_dofs; i++)
F(i) += JxW[qp] * penalty * ( T - u_dirichlet) * phi[i][qp];
if (compute_jacobian)
for (unsigned int i=0; i != n_T_dofs; i++)
for (unsigned int j=0; j != n_T_dofs; ++j)
K(i,j) += JxW[qp] * penalty * phi[i][qp] * phi[j][qp];
} // end of the quadrature point qp-loop
return compute_jacobian;
}
void LaplaceSystem::postprocess()
{
computed_QoI[0] = 0.0;
computed_QoI[1] = 0.0;
FEMSystem::postprocess();
Parallel::sum(computed_QoI[0]);
Parallel::sum(computed_QoI[1]);
}
Number LaplaceSystem::exact_solution(const Point& p)// xyz location
{
const Real x1 = p(0);
const Real x2 = p(1);
Real theta = atan2(x2,x1);
if (theta < 0)
theta += 2. * libMesh::pi;
return pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
}
The source file side_postprocess.C without comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
#include "L-shaped.h"
using namespace libMesh;
void LaplaceSystem::side_postprocess(DiffContext &context)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
unsigned int n_qpoints = (c.get_side_qrule())->n_points();
Number dQoI_1 = 0. ;
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
const Real TOL = 1.e-5;
if(fabs(y - 1.0) <= TOL && x > 0.0)
{
const Gradient grad_T = c.side_gradient(0,qp);
dQoI_1 += JxW[qp] * (grad_T * face_normals[qp]) ;
}
} // end of the quadrature point qp-loop
computed_QoI[1] = computed_QoI[1] + dQoI_1;
}
The source file side_qoi_derivative.C without comments:
#include "libmesh/libmesh_common.h"
#include "libmesh/elem.h"
#include "libmesh/fe_base.h"
#include "libmesh/fem_context.h"
#include "libmesh/point.h"
#include "libmesh/quadrature.h"
#include "L-shaped.h"
using namespace libMesh;
void LaplaceSystem::side_qoi_derivative (DiffContext &context,
const QoISet & /* qois */)
{
FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
const std::vector<Real> &JxW = c.side_fe_var[0]->get_JxW();
const std::vector<std::vector<RealGradient> > &dphi = c.side_fe_var[0]->get_dphi();
const std::vector<Point > &q_point = c.side_fe_var[0]->get_xyz();
const std::vector<Point> &face_normals = c.side_fe_var[0]->get_normals();
const unsigned int n_T_dofs = c.dof_indices_var[0].size();
unsigned int n_qpoints = (c.get_side_qrule())->n_points();
DenseSubVector<Number> &Q = *c.elem_qoi_subderivatives[1][0];
const Real TOL = 1.e-5;
for (unsigned int qp=0; qp != n_qpoints; qp++)
{
const Real x = q_point[qp](0);
const Real y = q_point[qp](1);
if(fabs(y - 1.0) <= TOL && x > 0.0)
{
for (unsigned int i=0; i != n_T_dofs; i++)
Q(i) += JxW[qp] * (dphi[i][qp] * face_normals[qp]);
}
} // end of the quadrature point qp-loop
}
The console output of the program:
***************************************************************
* Running Example adjoints_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
***************************************************************
Started /workspace/libmesh/examples/adjoints/adjoints_ex1/.libs/lt-example-devel
Reading in and building the mesh
Building system
Initializing systems
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=833
n_local_nodes()=433
n_elem()=255
n_local_elem()=131
n_active_elem()=192
n_subdomains()=1
n_partitions()=2
n_processors()=2
n_threads()=1
processor_id()=0
EquationSystems
n_systems()=1
System #0, "LaplaceSystem"
Type "Implicit"
Variables="T"
Finite Element Types="LAGRANGE", "JACOBI_20_00"
Infinite Element Mapping="CARTESIAN"
Approximation Orders="SECOND", "THIRD"
n_dofs()=833
n_local_dofs()=433
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 14.7203
Average Off-Processor Bandwidth <= 0.465786
Maximum On-Processor Bandwidth <= 26
Maximum Off-Processor Bandwidth <= 16
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Nonlinear solver converged, step 0, residual reduction 9.34193e-11 < 1e-09
Adaptive step 0, we have 192 active elements and 833 active dofs.
Postprocessing:
The relative error in QoI 0 is 0.00012462745923404317
The relative error in QoI 1 is 0.0001425593081568735
Using Adjoint Residual Error Estimator with Patch Recovery Weights
Refined mesh to 222 active elements and 947 active dofs.
Nonlinear solver converged, step 0, residual reduction 1.3689655455026611e-10 < 1.0000000000000001e-09
Adaptive step 1, we have 222 active elements and 947 active dofs.
Postprocessing:
The relative error in QoI 0 is 4.993438647470717e-05
The relative error in QoI 1 is 9.9674077986582966e-05
Using Adjoint Residual Error Estimator with Patch Recovery Weights
Refined mesh to 261 active elements and 1095 active dofs.
Nonlinear solver converged, step 0, residual reduction 1.5311732376561689e-10 < 1.0000000000000001e-09
Adaptive step 2, we have 261 active elements and 1095 active dofs.
Postprocessing:
The relative error in QoI 0 is 3.2764387216297939e-05
The relative error in QoI 1 is 9.1338016542211393e-05
Using Adjoint Residual Error Estimator with Patch Recovery Weights
Refined mesh to 306 active elements and 1261 active dofs.
Nonlinear solver converged, step 0, residual reduction 2.3826982360793464e-10 < 1.0000000000000001e-09
Adaptive step 3, we have 306 active elements and 1261 active dofs.
Postprocessing:
The relative error in QoI 0 is 1.3209665544372335e-05
The relative error in QoI 1 is 2.3940418269928809e-05
Using Adjoint Residual Error Estimator with Patch Recovery Weights
Refined mesh to 366 active elements and 1489 active dofs.
Nonlinear solver reached maximum step 1, latest evaluated residual 2.3112150692446607e-05
Continuing...
Adaptive step 4, we have 366 active elements and 1489 active dofs.
Postprocessing:
The relative error in QoI 0 is 4.959183651586029e-06
The relative error in QoI 1 is 2.3713339116644846e-06
Using Adjoint Residual Error Estimator with Patch Recovery Weights
Refined mesh to 429 active elements and 1727 active dofs.
Nonlinear solver reached maximum step 1, latest evaluated residual 2.3539816263250479e-05
Continuing...
Adaptive step 5, we have 429 active elements and 1727 active dofs.
Postprocessing:
The relative error in QoI 0 is 1.5870073697169233e-06
The relative error in QoI 1 is 7.483567240798096e-06
Using Adjoint Residual Error Estimator with Patch Recovery Weights
Refined mesh to 501 active elements and 1987 active dofs.
Nonlinear solver reached maximum step 1, latest evaluated residual 1.9232872080560362e-05
Continuing...
Adaptive step 6, we have 501 active elements and 1987 active dofs.
Postprocessing:
The relative error in QoI 0 is 3.2482930211046812e-07
The relative error in QoI 1 is 9.3753169864299525e-06
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adjoints/adjoints_ex1/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:39:25 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): 6.438e+00 1.00000 6.438e+00
Objects: 1.058e+03 1.00000 1.058e+03
Flops: 2.488e+08 1.15433 2.322e+08 4.643e+08
Flops/sec: 3.865e+07 1.15433 3.607e+07 7.213e+07
MPI Messages: 8.675e+02 1.00000 8.675e+02 1.735e+03
MPI Message Lengths: 7.802e+05 1.00000 8.993e+02 1.560e+06
MPI Reductions: 2.650e+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: 6.4375e+00 100.0% 4.6434e+08 100.0% 1.735e+03 100.0% 8.993e+02 100.0% 2.649e+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
KSPGMRESOrthog 442 1.0 1.0920e-02 1.5 1.52e+07 1.0 0.0e+00 0.0e+00 4.4e+02 0 6 0 0 17 0 6 0 0 17 2751
KSPSetUp 35 1.0 2.3222e-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 7 1.0 1.4677e-01 1.0 1.18e+08 1.2 4.0e+02 4.0e+02 4.5e+02 2 47 23 10 17 2 47 23 10 17 1496
PCSetUp 28 1.0 2.3085e-01 1.1 1.18e+08 1.2 0.0e+00 0.0e+00 9.8e+01 3 46 0 0 4 3 46 0 0 4 920
PCSetUpOnBlocks 7 1.0 1.1402e-01 1.1 5.88e+07 1.2 0.0e+00 0.0e+00 3.5e+01 2 23 0 0 1 2 23 0 0 1 931
PCApply 484 1.0 1.6287e-01 1.1 1.63e+08 1.2 0.0e+00 0.0e+00 3.5e+01 2 65 0 0 1 2 65 0 0 1 1865
VecMDot 442 1.0 8.5526e-03 1.7 7.60e+06 1.0 0.0e+00 0.0e+00 4.4e+02 0 3 0 0 17 0 3 0 0 17 1756
VecNorm 512 1.0 2.8467e-02 2.2 7.53e+05 1.0 0.0e+00 0.0e+00 5.1e+02 0 0 0 0 19 0 0 0 0 19 52
VecScale 463 1.0 2.3913e-04 1.0 3.40e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 2814
VecCopy 124 1.0 1.0920e-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 1281 1.0 4.5085e-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 49 1.0 6.6519e-05 1.1 7.25e+04 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 2156
VecMAXPY 463 1.0 2.2376e-03 1.0 8.25e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 4 0 0 0 0 4 0 0 0 7290
VecAssemblyBegin 343 1.0 6.0505e-02 1.0 0.00e+00 0.0 1.0e+02 2.1e+02 7.8e+02 1 0 6 1 29 1 0 6 1 29 0
VecAssemblyEnd 343 1.0 1.6141e-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 614 1.0 8.5473e-04 1.0 0.00e+00 0.0 1.2e+03 9.1e+02 0.0e+00 0 0 71 72 0 0 0 71 72 0 0
VecScatterEnd 614 1.0 1.8744e-02 2.8 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 463 1.0 2.2444e-02 3.1 1.02e+06 1.0 0.0e+00 0.0e+00 4.6e+02 0 0 0 0 17 0 0 0 0 17 90
MatMult 202 1.0 2.1081e-02 2.3 4.58e+06 1.0 4.0e+02 4.0e+02 0.0e+00 0 2 23 10 0 0 2 23 10 0 426
MatMultTranspose 261 1.0 4.9677e-03 1.0 6.02e+06 1.0 5.2e+02 4.0e+02 0.0e+00 0 3 30 14 0 0 3 30 14 0 2377
MatSolve 209 1.0 1.6316e-02 1.1 4.57e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 19 0 0 0 0 19 0 0 0 5299
MatSolveTranspos 275 1.0 2.9164e-02 1.1 5.88e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 24 0 0 0 0 24 0 0 0 3812
MatLUFactorNum 14 1.0 7.0183e-02 1.2 1.18e+08 1.2 0.0e+00 0.0e+00 0.0e+00 1 46 0 0 0 1 46 0 0 0 3025
MatILUFactorSym 14 1.0 1.5859e-01 1.1 0.00e+00 0.0 0.0e+00 0.0e+00 4.2e+01 2 0 0 0 2 2 0 0 0 2 0
MatAssemblyBegin 42 1.0 1.3707e-02 2.5 0.00e+00 0.0 4.8e+01 1.9e+03 8.4e+01 0 0 3 6 3 0 0 3 6 3 0
MatAssemblyEnd 42 1.0 1.6310e-03 1.1 0.00e+00 0.0 2.8e+01 1.0e+02 5.6e+01 0 0 2 0 2 0 0 2 0 2 0
MatGetRowIJ 14 1.0 4.0531e-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 14 1.0 2.4891e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.8e+01 0 0 0 0 1 0 0 0 0 1 0
MatZeroEntries 28 1.0 1.4687e-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
------------------------------------------------------------------------------------------------------------------------
Memory usage is given in bytes:
Object Type Creations Destructions Memory Descendants' Mem.
Reports information only for process 0.
--- Event Stage 0: Main Stage
Krylov Solver 28 28 271040 0
Preconditioner 28 28 24976 0
Vector 695 695 5137904 0
Vector Scatter 92 92 95312 0
Index Set 160 160 164700 0
IS L to G Mapping 19 19 10716 0
Matrix 35 35 17664712 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 0
Average time for MPI_Barrier(): 1.38283e-06
Average time for zero size MPI_Send(): 1.2517e-05
#PETSc Option Table entries:
-ksp_right_pc
-log_summary
-pc_type bjacobi
-sub_pc_factor_levels 4
-sub_pc_factor_zeropivot 0
-sub_pc_type ilu
#End of PETSc Option Table entries
Compiled without FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 8 sizeof(PetscInt) 4
Configure run at: Thu Nov 8 11:21:02 2012
Configure options: --with-debugging=false --COPTFLAGS=-O3 --CXXOPTFLAGS=-O3 --FOPTFLAGS=-O3 --with-clanguage=C++ --with-shared-libraries=1 --with-mpi-dir=/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1 --with-mumps=true --download-mumps=1 --with-metis=true --download-metis=1 --with-parmetis=true --download-parmetis=1 --with-superlu=true --download-superlu=1 --with-superludir=true --download-superlu_dist=1 --with-blacs=true --download-blacs=1 --with-scalapack=true --download-scalapack=1 --with-hypre=true --download-hypre=1 --with-blas-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_intel_lp64.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_sequential.so,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_core.so]" --with-lapack-lib="[/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64/libmkl_lapack95_lp64.a]"
-----------------------------------------
Libraries compiled on Thu Nov 8 11:21:02 2012 on daedalus.ices.utexas.edu
Machine characteristics: Linux-2.6.32-279.1.1.el6.x86_64-x86_64-with-redhat-6.3-Carbon
Using PETSc directory: /opt/apps/ossw/libraries/petsc/petsc-3.3-p2
Using PETSc arch: intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt
-----------------------------------------
Using C compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx -wd1572 -O3 -fPIC ${COPTFLAGS} ${CFLAGS}
Using Fortran compiler: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90 -fPIC -O3 ${FOPTFLAGS} ${FFLAGS}
-----------------------------------------
Using include paths: -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/include -I/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/include -I/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/include
-----------------------------------------
Using C linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpicxx
Using Fortran linker: /opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/bin/mpif90
Using libraries: -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lpetsc -lX11 -Wl,-rpath,/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -L/opt/apps/ossw/libraries/petsc/petsc-3.3-p2/intel-12.1-mkl-intel-10.3.12.361-mpich2-1.4.1p1-cxx-opt/lib -lcmumps -ldmumps -lsmumps -lzmumps -lmumps_common -lpord -lHYPRE -lpthread -lsuperlu_dist_3.0 -lparmetis -lmetis -lscalapack -lblacs -lsuperlu_4.3 -Wl,-rpath,/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -L/opt/apps/sysnet/intel/12.1/mkl/10.3.12.361/lib/intel64 -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,-rpath,/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -L/opt/apps/ossw/libraries/mpich2/mpich2-1.4.1p1/sl6/intel-12.1/lib -Wl,-rpath,/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -L/opt/apps/sysnet/intel/12.1/composer_xe_2011_sp1.7.256/compiler/lib/intel64 -Wl,-rpath,/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -L/usr/lib/gcc/x86_64-redhat-linux/4.4.6 -lmpichf90 -lifport -lifcore -lm -lm -lmpichcxx -ldl -lmpich -lopa -lmpl -lrt -lpthread -limf -lsvml -lipgo -ldecimal -lcilkrts -lstdc++ -lgcc_s -lirc -lirc_s -ldl
-----------------------------------------
----------------------------------------------------------------------------------------------------------------------
| Processor id: 0 |
| Num Processors: 2 |
| Time: Tue Feb 5 13:39:25 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=6.4496, Active time=6.37659 |
---------------------------------------------------------------------------------------------------------------------
| 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 |
|---------------------------------------------------------------------------------------------------------------------|
| |
| |
| AdjointResidualErrorEstimator |
| estimate_error() 6 0.0951 0.015852 1.7692 0.294872 1.49 27.75 |
| |
| DofMap |
| add_neighbors_to_send_list() 19 0.0756 0.003980 0.0907 0.004771 1.19 1.42 |
| build_constraint_matrix() 5205 0.0533 0.000010 0.0533 0.000010 0.84 0.84 |
| build_sparsity() 7 0.0471 0.006725 0.1161 0.016585 0.74 1.82 |
| cnstrn_elem_mat_vec() 1041 0.0065 0.000006 0.0065 0.000006 0.10 0.10 |
| constrain_elem_matrix() 1041 0.0057 0.000006 0.0057 0.000006 0.09 0.09 |
| constrain_elem_vector() 3123 0.0095 0.000003 0.0095 0.000003 0.15 0.15 |
| create_dof_constraints() 19 0.0584 0.003076 0.1141 0.006008 0.92 1.79 |
| distribute_dofs() 19 0.1107 0.005826 0.2770 0.014577 1.74 4.34 |
| dof_indices() 35049 1.7153 0.000049 1.7153 0.000049 26.90 26.90 |
| enforce_constraints_exactly() 48 0.0100 0.000208 0.0100 0.000208 0.16 0.16 |
| old_dof_indices() 6246 0.3602 0.000058 0.3602 0.000058 5.65 5.65 |
| prepare_send_list() 19 0.0004 0.000021 0.0004 0.000021 0.01 0.01 |
| reinit() 19 0.1609 0.008466 0.1609 0.008466 2.52 2.52 |
| |
| EquationSystems |
| build_discontinuous_solution_vector() 12 0.0082 0.000681 0.0439 0.003660 0.13 0.69 |
| build_solution_vector() 21 0.0168 0.000801 0.2175 0.010355 0.26 3.41 |
| |
| FE |
| compute_shape_functions() 22904 0.9155 0.000040 0.9155 0.000040 14.36 14.36 |
| init_shape_functions() 3498 0.0884 0.000025 0.0884 0.000025 1.39 1.39 |
| inverse_map() 18324 0.1406 0.000008 0.1406 0.000008 2.21 2.21 |
| |
| FEMSystem |
| assemble_qoi_derivative() 7 0.0819 0.011700 0.4728 0.067536 1.28 7.41 |
| assembly() 7 0.0692 0.009883 0.6840 0.097716 1.08 10.73 |
| assembly(get_jacobian) 7 0.0661 0.009450 0.6765 0.096646 1.04 10.61 |
| assembly(get_residual) 7 0.0597 0.008521 0.4393 0.062762 0.94 6.89 |
| numerical_elem_jacobian() 2274 0.4141 0.000182 0.4141 0.000182 6.49 6.49 |
| numerical_side_jacobian() 640 0.0418 0.000065 0.0418 0.000065 0.66 0.66 |
| postprocess() 7 0.0378 0.005393 0.4008 0.057261 0.59 6.29 |
| |
| FEMap |
| compute_affine_map() 22904 0.2143 0.000009 0.2143 0.000009 3.36 3.36 |
| compute_face_map() 3200 0.0635 0.000020 0.1628 0.000051 1.00 2.55 |
| init_face_shape_functions() 70 0.0006 0.000009 0.0006 0.000009 0.01 0.01 |
| init_reference_to_physical_map() 3498 0.0852 0.000024 0.0852 0.000024 1.34 1.34 |
| |
| GMVIO |
| write_nodal_data() 21 0.0819 0.003900 0.0829 0.003946 1.28 1.30 |
| |
| ImplicitSystem |
| adjoint_solve() 7 0.0009 0.000134 1.3308 0.190118 0.01 20.87 |
| |
| LocationMap |
| find() 3320 0.0130 0.000004 0.0130 0.000004 0.20 0.20 |
| init() 15 0.0159 0.001063 0.0159 0.001063 0.25 0.25 |
| |
| Mesh |
| all_first_order() 12 0.0342 0.002854 0.0342 0.002854 0.54 0.54 |
| all_second_order() 1 0.0002 0.000233 0.0002 0.000233 0.00 0.00 |
| contract() 6 0.0010 0.000171 0.0071 0.001178 0.02 0.11 |
| find_neighbors() 33 0.1465 0.004438 0.1484 0.004497 2.30 2.33 |
| renumber_nodes_and_elem() 6 0.0060 0.001007 0.0060 0.001007 0.09 0.09 |
| |
| MeshCommunication |
| broadcast() 1 0.0007 0.000689 0.0010 0.000978 0.01 0.02 |
| compute_hilbert_indices() 22 0.0241 0.001096 0.0241 0.001096 0.38 0.38 |
| find_global_indices() 22 0.0104 0.000474 0.0408 0.001854 0.16 0.64 |
| parallel_sort() 22 0.0037 0.000166 0.0047 0.000213 0.06 0.07 |
| |
| MeshOutput |
| write_equation_systems() 21 0.0006 0.000031 0.3020 0.014382 0.01 4.74 |
| |
| MeshRefinement |
| _coarsen_elements() 12 0.0007 0.000057 0.0007 0.000061 0.01 0.01 |
| _refine_elements() 15 0.0140 0.000934 0.0408 0.002723 0.22 0.64 |
| add_point() 3320 0.0111 0.000003 0.0245 0.000007 0.17 0.38 |
| make_coarsening_compatible() 29 0.0257 0.000885 0.0257 0.000885 0.40 0.40 |
| make_refinement_compatible() 29 0.0020 0.000069 0.0021 0.000072 0.03 0.03 |
| |
| MetisPartitioner |
| partition() 21 0.0665 0.003166 0.1081 0.005145 1.04 1.69 |
| |
| NewtonSolver |
| solve() 7 0.0397 0.005669 1.3179 0.188276 0.62 20.67 |
| |
| Parallel |
| allgather() 115 0.0024 0.000021 0.0026 0.000023 0.04 0.04 |
| broadcast() 9 0.0002 0.000026 0.0002 0.000020 0.00 0.00 |
| max(bool) 70 0.0022 0.000031 0.0022 0.000031 0.03 0.03 |
| max(scalar) 2964 0.0075 0.000003 0.0075 0.000003 0.12 0.12 |
| max(vector) 697 0.0043 0.000006 0.0094 0.000013 0.07 0.15 |
| max(vector) 12 0.0006 0.000051 0.0008 0.000063 0.01 0.01 |
| min(bool) 3628 0.0087 0.000002 0.0087 0.000002 0.14 0.14 |
| min(scalar) 2895 0.0419 0.000014 0.0419 0.000014 0.66 0.66 |
| min(vector) 703 0.0048 0.000007 0.0105 0.000015 0.08 0.16 |
| probe() 200 0.0013 0.000007 0.0013 0.000007 0.02 0.02 |
| receive() 200 0.0010 0.000005 0.0023 0.000012 0.02 0.04 |
| send() 200 0.0005 0.000002 0.0005 0.000002 0.01 0.01 |
| send_receive() 244 0.0014 0.000006 0.0046 0.000019 0.02 0.07 |
| sum() 244 0.0107 0.000044 0.0120 0.000049 0.17 0.19 |
| |
| Parallel::Request |
| wait() 200 0.0003 0.000002 0.0003 0.000002 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 33 0.0194 0.000589 0.0232 0.000704 0.30 0.36 |
| set_parent_processor_ids() 21 0.0076 0.000360 0.0076 0.000360 0.12 0.12 |
| |
| PatchRecoveryErrorEstimator |
| estimate_error() 18 0.3735 0.020748 1.2206 0.067813 5.86 19.14 |
| |
| PetscLinearSolver |
| solve() 21 0.3256 0.015505 0.3256 0.015505 5.11 5.11 |
| |
| ProjectVector |
| operator() 18 0.0292 0.001625 0.4191 0.023285 0.46 6.57 |
| |
| System |
| project_vector() 18 0.0220 0.001224 0.6256 0.034755 0.35 9.81 |
---------------------------------------------------------------------------------------------------------------------
| Totals: 148693 6.3766 100.00 |
---------------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adjoints_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
***************************************************************