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_ex4.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"
Adjoint Related includes
#include "libmesh/qoi_set.h"
#include "libmesh/adjoint_refinement_estimator.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;
}
This is where declare the adjoint refined error estimator. This estimator builds an error bound
for Q(u) - Q(u_h), by solving the adjoint problem on a finer Finite Element space. For more details
see the description of the Adjoint Refinement Error Estimator in adjoint_refinement_error_estimator.C
AutoPtr<AdjointRefinementEstimator> build_adjoint_refinement_error_estimator(QoISet &qois)
{
AutoPtr<AdjointRefinementEstimator> error_estimator;
std::cout<<"Computing the error estimate using the Adjoint Refinement Error Estimator"<<std::endl<<std::endl;
AdjointRefinementEstimator *adjoint_refinement_estimator = new AdjointRefinementEstimator;
error_estimator.reset (adjoint_refinement_estimator);
adjoint_refinement_estimator->qoi_set() = qois;
We enrich the FE space for the dual problem by doing 2 uniform h refinements
adjoint_refinement_estimator->number_h_refinements = 2;
return error_estimator;
}
#endif // LIBMESH_ENABLE_AMR
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;
We will declare an error vector for passing to the adjoint refinement error estimator
ErrorVector QoI_elementwise_error;
Build an adjoint refinement error estimator object
AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
build_adjoint_refinement_error_estimator(qois);
Estimate the error in each element using the Adjoint Refinement estimator
adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
Print out the computed error estimate, note that we access the global error estimates
using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_exact) << std::endl << std::endl;
Also print out effecitivity indices (estimated error/true error)
std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_computed - QoI_1_exact) << std::endl << std::endl;
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 (QoI_elementwise_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 (QoI_elementwise_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;
We will declare an error vector for passing to the adjoint refinement error estimator
Right now, only the first entry of this vector will be filled (with the global QoI error estimate)
Later, each entry of the vector will contain elementwise error that the user can sum to get the total error
ErrorVector QoI_elementwise_error;
Build an adjoint refinement error estimator object
AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
build_adjoint_refinement_error_estimator(qois);
Estimate the error in each element using the Adjoint Refinement estimator
adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
Print out the computed error estimate, note that we access the global error estimates
using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_exact) << std::endl << std::endl;
Also print out effecitivity indices (estimated error/true error)
std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_computed - 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();
CommWorld.sum(computed_QoI[0]);
CommWorld.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] * x * (x - 1.)) ;
}
} // 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();
Basis Functions
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> &Q1 = *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++)
{
Q1(i) += JxW[qp] * (dphi[i][qp] * face_normals[qp]) * x * (x - 1.);
}
}
} // 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_ex4.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/qoi_set.h"
#include "libmesh/adjoint_refinement_estimator.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;
}
AutoPtr<AdjointRefinementEstimator> build_adjoint_refinement_error_estimator(QoISet &qois)
{
AutoPtr<AdjointRefinementEstimator> error_estimator;
std::cout<<"Computing the error estimate using the Adjoint Refinement Error Estimator"<<std::endl<<std::endl;
AdjointRefinementEstimator *adjoint_refinement_estimator = new AdjointRefinementEstimator;
error_estimator.reset (adjoint_refinement_estimator);
adjoint_refinement_estimator->qoi_set() = qois;
adjoint_refinement_estimator->number_h_refinements = 2;
return error_estimator;
}
#endif // LIBMESH_ENABLE_AMR
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 QoI_elementwise_error;
AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
build_adjoint_refinement_error_estimator(qois);
adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_exact) << std::endl << std::endl;
std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_computed - QoI_1_exact) << std::endl << std::endl;
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 (QoI_elementwise_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 (QoI_elementwise_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;
ErrorVector QoI_elementwise_error;
AutoPtr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
build_adjoint_refinement_error_estimator(qois);
adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
std::cout<< "The computed relative error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_exact) << std::endl;
std::cout<< "The computed relative error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_exact) << std::endl << std::endl;
std::cout<< "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
std::cout<< "The effectivity index for the computed error in QoI 1 is " << std::setprecision(17)
<< std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) /
std::abs(QoI_1_computed - 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();
CommWorld.sum(computed_QoI[0]);
CommWorld.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] * x * (x - 1.)) ;
}
} // 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> &Q1 = *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++)
{
Q1(i) += JxW[qp] * (dphi[i][qp] * face_normals[qp]) * x * (x - 1.);
}
}
} // end of the quadrature point qp-loop
}
The console output of the program:
***************************************************************
* Running Example adjoints_ex4:
* 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_ex4/.libs/lt-example-devel
Reading in and building the mesh
Building system
Initializing systems
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=65
n_local_nodes()=37
n_elem()=15
n_local_elem()=8
n_active_elem()=12
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()=65
n_local_dofs()=37
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 11.9077
Average Off-Processor Bandwidth <= 1.29231
Maximum On-Processor Bandwidth <= 26
Maximum Off-Processor Bandwidth <= 10
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Nonlinear solver converged, step 0, residual reduction 9.00697e-11 < 1e-09
Adaptive step 0, we have 12 active elements and 65 active dofs.
Postprocessing:
The relative error in QoI 0 is 1.7448570691365428
The relative error in QoI 1 is 0.0044573233492977817
Computing the error estimate using the Adjoint Refinement Error Estimator
The computed relative error in QoI 0 is 0.00017026143445813402
The computed relative error in QoI 1 is 0.0045778053548259532
The effectivity index for the computed error in QoI 0 is 9.7579015192567794e-05
The effectivity index for the computed error in QoI 1 is 1.0270301246031777
Refined mesh to 48 active elements and 225 active dofs.
Nonlinear solver converged, step 0, residual reduction 5.6491314599024831e-11 < 1.0000000000000001e-09
Adaptive step 1, we have 48 active elements and 225 active dofs.
Postprocessing:
The relative error in QoI 0 is 0.0002799284426766521
The relative error in QoI 1 is 0.00047776972749558597
Computing the error estimate using the Adjoint Refinement Error Estimator
The computed relative error in QoI 0 is 0.00022981246783585551
The computed relative error in QoI 1 is 0.00059422664700472855
The effectivity index for the computed error in QoI 0 is 0.82096862197498799
The effectivity index for the computed error in QoI 1 is 1.2437511479004675
Refined mesh to 192 active elements and 833 active dofs.
Nonlinear solver converged, step 0, residual reduction 4.9476448167753118e-11 < 1.0000000000000001e-09
Adaptive step 2, we have 192 active elements and 833 active dofs.
Postprocessing:
The relative error in QoI 0 is 0.00012462745998891614
The relative error in QoI 1 is 0.0001204819517696618
Computing the error estimate using the Adjoint Refinement Error Estimator
The computed relative error in QoI 0 is 0.00010469999414423847
The computed relative error in QoI 1 is 5.8180075906450974e-05
The effectivity index for the computed error in QoI 0 is 0.8401037311805124
The effectivity index for the computed error in QoI 1 is 0.48289453359520634
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adjoints/adjoints_ex4/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:39:59 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): 7.303e+00 1.00000 7.303e+00
Objects: 8.850e+02 1.00000 8.850e+02
Flops: 2.454e+08 1.02556 2.423e+08 4.847e+08
Flops/sec: 3.360e+07 1.02556 3.318e+07 6.636e+07
MPI Messages: 6.825e+02 1.00000 6.825e+02 1.365e+03
MPI Message Lengths: 3.991e+05 1.00000 5.848e+02 7.982e+05
MPI Reductions: 2.241e+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: 7.3032e+00 100.0% 4.8468e+08 100.0% 1.365e+03 100.0% 5.848e+02 100.0% 2.240e+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 333 1.0 1.5001e-02 1.2 3.66e+07 1.0 0.0e+00 0.0e+00 3.3e+02 0 15 0 0 15 0 15 0 0 15 4879
KSPSetUp 24 1.0 2.1482e-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
KSPSolve 3 1.0 8.5311e-03 1.1 4.67e+06 1.0 1.3e+02 2.3e+02 1.5e+02 0 2 9 4 7 0 2 9 4 7 1093
PCSetUp 18 1.0 9.3838e-02 1.0 4.27e+07 1.1 0.0e+00 0.0e+00 6.7e+01 1 17 0 0 3 1 17 0 0 3 882
PCSetUpOnBlocks 3 1.0 3.4328e-03 1.0 1.36e+06 1.1 0.0e+00 0.0e+00 1.7e+01 0 1 0 0 1 0 1 0 0 1 750
PCApply 365 1.0 1.9943e-01 1.0 1.81e+08 1.0 0.0e+00 0.0e+00 3.2e+01 3 73 0 0 1 3 73 0 0 1 1786
VecDot 6 1.0 3.7193e-05 1.0 3.32e+04 1.0 0.0e+00 0.0e+00 6.0e+00 0 0 0 0 0 0 0 0 0 0 1783
VecMDot 333 1.0 9.6173e-03 1.3 1.83e+07 1.0 0.0e+00 0.0e+00 3.3e+02 0 8 0 0 15 0 8 0 0 15 3804
VecNorm 377 1.0 1.6322e-02 1.2 1.49e+06 1.0 0.0e+00 0.0e+00 3.8e+02 0 1 0 0 17 0 1 0 0 17 182
VecScale 350 1.0 3.5548e-04 1.1 7.24e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 4066
VecCopy 156 1.0 1.5545e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecSet 1239 1.0 1.1516e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 37 1.0 8.2731e-05 1.1 1.25e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 3006
VecMAXPY 350 1.0 5.3318e-03 1.0 1.97e+07 1.0 0.0e+00 0.0e+00 0.0e+00 0 8 0 0 0 0 8 0 0 0 7383
VecAssemblyBegin 312 1.0 5.8199e-02 2.7 0.00e+00 0.0 1.8e+02 1.5e+02 7.2e+02 1 0 13 3 32 1 0 13 3 32 0
VecAssemblyEnd 312 1.0 1.9073e-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 456 1.0 9.0909e-04 1.0 0.00e+00 0.0 9.0e+02 7.6e+02 0.0e+00 0 0 66 86 0 0 0 66 86 0 0
VecScatterEnd 456 1.0 9.8631e-0313.4 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecNormalize 350 1.0 6.0112e-03 2.1 2.17e+06 1.0 0.0e+00 0.0e+00 3.5e+02 0 1 0 0 16 0 1 0 0 16 721
MatMult 64 1.0 2.0459e-03 4.1 4.59e+05 1.1 1.3e+02 2.3e+02 0.0e+00 0 0 9 4 0 0 0 9 4 0 427
MatMultTranspose 286 1.0 1.8388e-02 1.0 2.22e+07 1.0 5.7e+02 9.5e+02 0.0e+00 0 9 42 68 0 0 9 42 68 0 2413
MatSolve 67 1.0 8.0252e-04 1.0 2.04e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 1 0 0 0 0 1 0 0 0 5077
MatSolveTranspos 298 1.0 1.0757e-01 1.0 1.38e+08 1.0 0.0e+00 0.0e+00 0.0e+00 1 56 0 0 0 1 56 0 0 0 2528
MatLUFactorNum 9 1.0 2.8982e-02 1.0 4.27e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 17 0 0 0 0 17 0 0 0 2855
MatILUFactorSym 9 1.0 6.3424e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.7e+01 1 0 0 0 1 1 0 0 0 1 0
MatAssemblyBegin 30 1.0 7.9618e-03 1.5 0.00e+00 0.0 3.6e+01 1.6e+03 6.0e+01 0 0 3 7 3 0 0 3 7 3 0
MatAssemblyEnd 30 1.0 2.0347e-03 1.2 0.00e+00 0.0 2.4e+01 1.6e+02 4.8e+01 0 0 2 0 2 0 0 2 0 2 0
MatGetRowIJ 9 1.0 4.0531e-06 1.3 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
MatGetOrdering 9 1.0 2.3818e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.2e+01 0 0 0 0 1 0 0 0 0 1 0
MatZeroEntries 39 1.0 5.2643e-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
------------------------------------------------------------------------------------------------------------------------
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 27 27 185328 0
Preconditioner 27 27 23472 0
Vector 565 565 5969312 0
Vector Scatter 66 66 68376 0
Index Set 124 124 140836 0
IS L to G Mapping 21 21 11844 0
Matrix 54 54 13187476 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 2.00272e-06
Average time for zero size MPI_Send(): 1.84774e-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:59 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=7.31416, Active time=7.2113 |
----------------------------------------------------------------------------------------------------------------
| 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() 15 0.1608 0.010719 0.1867 0.012444 2.23 2.59 |
| build_sparsity() 15 0.1124 0.007491 0.3033 0.020222 1.56 4.21 |
| create_dof_constraints() 15 0.0172 0.001145 0.0172 0.001145 0.24 0.24 |
| distribute_dofs() 15 0.2383 0.015886 0.5891 0.039273 3.30 8.17 |
| dof_indices() 30667 1.7709 0.000058 1.7709 0.000058 24.56 24.56 |
| old_dof_indices() 22140 1.2649 0.000057 1.2649 0.000057 17.54 17.54 |
| prepare_send_list() 15 0.0007 0.000050 0.0007 0.000050 0.01 0.01 |
| reinit() 15 0.3462 0.023077 0.3462 0.023077 4.80 4.80 |
| |
| EquationSystems |
| build_solution_vector() 9 0.0024 0.000267 0.0272 0.003017 0.03 0.38 |
| |
| FE |
| compute_shape_functions() 15328 0.7951 0.000052 0.7951 0.000052 11.03 11.03 |
| init_shape_functions() 2020 0.0441 0.000022 0.0441 0.000022 0.61 0.61 |
| inverse_map() 39264 0.3020 0.000008 0.3020 0.000008 4.19 4.19 |
| |
| FEMSystem |
| assemble_qoi_derivative() 6 0.1242 0.020705 0.7209 0.120149 1.72 10.00 |
| assembly() 3 0.0105 0.003489 0.0894 0.029814 0.15 1.24 |
| assembly(get_jacobian) 6 0.1014 0.016907 1.1508 0.191806 1.41 15.96 |
| assembly(get_residual) 6 0.0912 0.015194 0.6926 0.115441 1.26 9.61 |
| numerical_elem_jacobian() 2268 0.4529 0.000200 0.4529 0.000200 6.28 6.28 |
| numerical_side_jacobian() 348 0.0222 0.000064 0.0222 0.000064 0.31 0.31 |
| postprocess() 3 0.0058 0.001926 0.0512 0.017055 0.08 0.71 |
| |
| FEMap |
| compute_affine_map() 15328 0.1451 0.000009 0.1451 0.000009 2.01 2.01 |
| compute_face_map() 1972 0.0334 0.000017 0.0816 0.000041 0.46 1.13 |
| init_face_shape_functions() 48 0.0005 0.000010 0.0005 0.000010 0.01 0.01 |
| init_reference_to_physical_map() 2020 0.0430 0.000021 0.0430 0.000021 0.60 0.60 |
| |
| GMVIO |
| write_nodal_data() 9 0.0106 0.001174 0.0114 0.001271 0.15 0.16 |
| |
| ImplicitSystem |
| adjoint_solve() 6 0.0007 0.000122 2.1153 0.352551 0.01 29.33 |
| |
| LocationMap |
| find() 26460 0.0962 0.000004 0.0962 0.000004 1.33 1.33 |
| init() 23 0.0278 0.001208 0.0278 0.001208 0.39 0.39 |
| |
| Mesh |
| all_second_order() 1 0.0003 0.000332 0.0003 0.000332 0.00 0.00 |
| contract() 14 0.0128 0.000915 0.0481 0.003434 0.18 0.67 |
| find_neighbors() 17 0.1704 0.010025 0.1714 0.010084 2.36 2.38 |
| renumber_nodes_and_elem() 14 0.0353 0.002519 0.0353 0.002519 0.49 0.49 |
| |
| MeshCommunication |
| broadcast() 1 0.0009 0.000883 0.0012 0.001248 0.01 0.02 |
| compute_hilbert_indices() 6 0.0013 0.000215 0.0013 0.000215 0.02 0.02 |
| find_global_indices() 6 0.0010 0.000171 0.0043 0.000715 0.01 0.06 |
| parallel_sort() 6 0.0011 0.000187 0.0013 0.000225 0.02 0.02 |
| |
| MeshOutput |
| write_equation_systems() 9 0.0003 0.000030 0.0397 0.004412 0.00 0.55 |
| |
| MeshRefinement |
| _coarsen_elements() 20 0.0116 0.000582 0.0117 0.000587 0.16 0.16 |
| _refine_elements() 23 0.1056 0.004593 0.2984 0.012975 1.46 4.14 |
| add_point() 26460 0.0878 0.000003 0.1874 0.000007 1.22 2.60 |
| make_coarsening_compatible() 28 0.0375 0.001340 0.0375 0.001340 0.52 0.52 |
| make_refinement_compatible() 28 0.0001 0.000002 0.0001 0.000004 0.00 0.00 |
| |
| MetisPartitioner |
| partition() 5 0.0045 0.000892 0.0079 0.001586 0.06 0.11 |
| |
| NewtonSolver |
| solve() 3 0.0156 0.005215 0.1762 0.058745 0.22 2.44 |
| |
| Parallel |
| allgather() 59 0.0005 0.000009 0.0006 0.000010 0.01 0.01 |
| broadcast() 9 0.0003 0.000034 0.0002 0.000027 0.00 0.00 |
| max(bool) 73 0.0019 0.000025 0.0019 0.000025 0.03 0.03 |
| max(scalar) 2118 0.0066 0.000003 0.0066 0.000003 0.09 0.09 |
| max(vector) 510 0.0034 0.000007 0.0080 0.000016 0.05 0.11 |
| min(bool) 2629 0.0075 0.000003 0.0075 0.000003 0.10 0.10 |
| min(scalar) 2069 0.0122 0.000006 0.0122 0.000006 0.17 0.17 |
| min(vector) 510 0.0037 0.000007 0.0084 0.000017 0.05 0.12 |
| probe() 112 0.0016 0.000014 0.0016 0.000014 0.02 0.02 |
| receive() 112 0.0008 0.000007 0.0024 0.000021 0.01 0.03 |
| send() 112 0.0003 0.000003 0.0003 0.000003 0.00 0.00 |
| send_receive() 124 0.0011 0.000009 0.0040 0.000033 0.01 0.06 |
| sum() 124 0.0019 0.000015 0.0030 0.000024 0.03 0.04 |
| |
| Parallel::Request |
| wait() 112 0.0002 0.000002 0.0002 0.000002 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 5 0.0012 0.000245 0.0017 0.000338 0.02 0.02 |
| set_parent_processor_ids() 5 0.0004 0.000083 0.0004 0.000083 0.01 0.01 |
| |
| PetscLinearSolver |
| solve() 15 0.2521 0.016809 0.2521 0.016809 3.50 3.50 |
| |
| ProjectVector |
| operator() 36 0.1580 0.004388 1.4734 0.040927 2.19 20.43 |
| |
| System |
| project_vector() 36 0.0513 0.001424 2.2666 0.062962 0.71 31.43 |
----------------------------------------------------------------------------------------------------------------
| Totals: 193435 7.2113 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adjoints_ex4:
* 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
***************************************************************