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),
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"),
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);
unsigned int dimension;
Real elementorder;
std::string domainfile;
unsigned int coarserefinements;
bool solver_quiet, 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;
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-qoi.h with comments:
#ifndef L_QOI_H
#define L_QOI_H
#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 "libmesh/diff_qoi.h"
Bring in everything from the libMesh namespace
using namespace libMesh;
class LaplaceQoI : public DifferentiableQoI
{
public:
LaplaceQoI(){}
virtual ~LaplaceQoI(){}
virtual void init_qoi( std::vector<Number>& sys_qoi );
virtual void postprocess( ){}
virtual void element_qoi_derivative(DiffContext &context, const QoISet & qois);
virtual void element_qoi (DiffContext &context, const QoISet & qois);
virtual AutoPtr<DifferentiableQoI> clone( ) {
return AutoPtr<DifferentiableQoI> ( new LaplaceQoI(*this) );
}
};
#endif // L_QOI_H
The source file L-shaped.h with comments:
#include "libmesh/enum_fe_family.h"
#include "libmesh/fem_system.h"
#include "libmesh/parameter_vector.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) { }
std::string & fe_family() { return _fe_family; }
unsigned int & fe_order() { return _fe_order; }
bool & analytic_jacobians() { return _analytic_jacobians; }
Number &get_parameter_value(unsigned int parameter_index)
{
return parameters[parameter_index];
}
ParameterVector &get_parameter_vector()
{
parameter_vector.resize(parameters.size());
for(unsigned int i = 0; i != parameters.size(); ++i)
{
parameter_vector[i] = ¶meters[i];
}
return parameter_vector;
}
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);
Number exact_solution (const Point&);
Parameters associated with the system
std::vector<Number> parameters;
The ParameterVector object that will contain pointers to
the system parameters
ParameterVector parameter_vector;
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_ex2.C with comments:
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). We scale this exact solution by a combination of parameters, (\alpha_{1} + 2 *\alpha_{2}) to get u = (\alpha_{1} + 2 *\alpha_{2}) * r^{2/3} * \sin ( (2/3) * \theta), which again satisfies the Laplace Equation. We define the Quantity of Interest in element_qoi.C, and compute the sensitivity of the QoI to \alpha_{1} and \alpha_{2} using the adjoint sensitivity method. Since we use the adjoint capabilities of libMesh in this example, we use the DiffSystem framework. This file (main.C) contains the declaration of mesh and equation system objects, L-shaped.C contains the assembly of the system, element_qoi_derivative.C contains the RHS for the adjoint system. Postprocessing to compute the the QoIs is done in element_qoi.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/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"
Sensitivity Calculation related includes
#include "libmesh/parameter_vector.h"
#include "libmesh/sensitivity_data.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"
#include "L-qoi.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)
{
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<<std::endl;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
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);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
}
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");
QoISet qois;
std::vector<unsigned int> qoi_indices;
qoi_indices.push_back(0);
qois.add_indices(qoi_indices);
qois.set_weight(0, 0.5);
Put some scope here to test that the cloning is working right
{
LaplaceQoI qoi;
system.attach_qoi( &qoi );
}
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();
{
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);
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;
A SensitivityData object to hold the qois and parameters
SensitivityData sensitivities(qois, system, system.get_parameter_vector());
Make sure we get the contributions to the adjoint RHS from the sides
system.assemble_qoi_sides = true;
Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
function, so we have to set the adjoint_already_solved boolean to false
system.set_adjoint_already_solved(false);
Compute the sensitivities
system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
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);
GetPot infile_l_shaped("l-shaped.in");
Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
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<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
std::cout<< "The relative error in sensitivity QoI_0_0 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) /
std::abs(sensitivity_QoI_0_0_exact) << std::endl;
std::cout<< "The relative error in sensitivity QoI_0_1 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) /
std::abs(sensitivity_QoI_0_1_exact) << std::endl << std::endl;
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);
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)
{
std::cout<<"Refining Uniformly"<<std::endl<<std::endl;
mesh_refinement->uniformly_refine(1);
}
Adaptively refine based on reaching an error tolerance
else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
{
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);
Estimate the error in each element using the Adjoint Residual or Kelly error estimator
error_estimator->estimate_error(system, error);
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
{
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);
Estimate the error in each element using the Adjoint Residual or Kelly error estimator
error_estimator->estimate_error(system, error);
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)
{
system.solve();
write_output(equation_systems, a_step, "primal");
NumericVector<Number> &primal_solution = *system.solution;
SensitivityData sensitivities(qois, system, system.get_parameter_vector());
system.assemble_qoi_sides = true;
Here we solve the adjoint problem inside the adjoint_qoi_parameter_sensitivity
function, so we have to set the adjoint_already_solved boolean to false
system.set_adjoint_already_solved(false);
system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
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);
GetPot infile_l_shaped("l-shaped.in");
Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
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<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
std::cout<< "The error in sensitivity QoI_0_0 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact << std::endl;
std::cout<< "The error in sensitivity QoI_0_1 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact << std::endl << std::endl;
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);
}
}
std::cerr << '[' << libMesh::processor_id()
<< "] Completing output." << std::endl;
#endif // #ifndef LIBMESH_ENABLE_AMR
All done.
return 0;
}
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(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_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-qoi.C with comments:
#include "L-qoi.h"
using namespace libMesh;
void LaplaceQoI::init_qoi( std::vector<Number>& sys_qoi )
{
Only 1 qoi to worry about
sys_qoi.resize(1);
return;
}
We only have one QoI, so we don't bother checking the qois argument
to see if it was requested from us
void LaplaceQoI::element_qoi (DiffContext &context,
const QoISet & /* qois */ )
{
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();
unsigned int n_qpoints = (c.get_element_qrule())->n_points();
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 xf = xyz[qp](0);
const Real yf = xyz[qp](1);
If in the sub-domain omega, add the contribution to the integral R
if(fabs(xf - 0.875) <= 0.125 && fabs(yf - 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
c.elem_qoi[0] = c.elem_qoi[0] + dQoI_0;
}
We only have one QoI, so we don't bother checking the qois argument
to see if it was requested from us
void LaplaceQoI::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 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");
parameters.push_back(infile("alpha_1", 1.0));
parameters.push_back(infile("alpha_2", 1.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 function gradients
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] * (parameters[0] + (2.*parameters[1])) * ( 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] * (parameters[0] + (2.*parameters[1])) * ( 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;
}
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 (parameters[0] + (2.*parameters[1])) * pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
}
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),
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"),
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);
unsigned int dimension;
Real elementorder;
std::string domainfile;
unsigned int coarserefinements;
bool solver_quiet, 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;
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-qoi.h without comments:
#ifndef L_QOI_H
#define L_QOI_H
#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 "libmesh/diff_qoi.h"
using namespace libMesh;
class LaplaceQoI : public DifferentiableQoI
{
public:
LaplaceQoI(){}
virtual ~LaplaceQoI(){}
virtual void init_qoi( std::vector<Number>& sys_qoi );
virtual void postprocess( ){}
virtual void element_qoi_derivative(DiffContext &context, const QoISet & qois);
virtual void element_qoi (DiffContext &context, const QoISet & qois);
virtual AutoPtr<DifferentiableQoI> clone( ) {
return AutoPtr<DifferentiableQoI> ( new LaplaceQoI(*this) );
}
};
#endif // L_QOI_H
The source file L-shaped.h without comments:
#include "libmesh/enum_fe_family.h"
#include "libmesh/fem_system.h"
#include "libmesh/parameter_vector.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) { }
std::string & fe_family() { return _fe_family; }
unsigned int & fe_order() { return _fe_order; }
bool & analytic_jacobians() { return _analytic_jacobians; }
Number &get_parameter_value(unsigned int parameter_index)
{
return parameters[parameter_index];
}
ParameterVector &get_parameter_vector()
{
parameter_vector.resize(parameters.size());
for(unsigned int i = 0; i != parameters.size(); ++i)
{
parameter_vector[i] = ¶meters[i];
}
return parameter_vector;
}
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);
Number exact_solution (const Point&);
std::vector<Number> parameters;
ParameterVector parameter_vector;
std::string _fe_family;
unsigned int _fe_order;
bool _analytic_jacobians;
};
The source file adjoints_ex2.C without comments:
#include <iostream>
#include <iomanip>
#include "libmesh/equation_systems.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/parameter_vector.h"
#include "libmesh/sensitivity_data.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"
#include "L-qoi.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)
{
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<<std::endl;
AdjointResidualErrorEstimator *adjoint_residual_estimator = new AdjointResidualErrorEstimator;
error_estimator.reset (adjoint_residual_estimator);
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);
adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0, H1_SEMINORM);
}
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");
QoISet qois;
std::vector<unsigned int> qoi_indices;
qoi_indices.push_back(0);
qois.add_indices(qoi_indices);
qois.set_weight(0, 0.5);
{
LaplaceQoI qoi;
system.attach_qoi( &qoi );
}
set_system_parameters(system, param);
std::cout << "Initializing systems" << std::endl;
equation_systems.init ();
mesh.print_info();
equation_systems.print_info();
{
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);
system.solve();
write_output(equation_systems, a_step, "primal");
NumericVector<Number> &primal_solution = *system.solution;
SensitivityData sensitivities(qois, system, system.get_parameter_vector());
system.assemble_qoi_sides = true;
system.set_adjoint_already_solved(false);
system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
system.set_adjoint_already_solved(true);
GetPot infile_l_shaped("l-shaped.in");
Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
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<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
std::cout<< "The relative error in sensitivity QoI_0_0 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) /
std::abs(sensitivity_QoI_0_0_exact) << std::endl;
std::cout<< "The relative error in sensitivity QoI_0_1 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) /
std::abs(sensitivity_QoI_0_1_exact) << std::endl << std::endl;
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);
if(param.refine_uniformly)
{
std::cout<<"Refining Uniformly"<<std::endl<<std::endl;
mesh_refinement->uniformly_refine(1);
}
else if(param.global_tolerance >= 0. && param.nelem_target == 0.)
{
ErrorVector error;
AutoPtr<ErrorEstimator> error_estimator =
build_error_estimator(param);
error_estimator->estimate_error(system, error);
mesh_refinement->flag_elements_by_error_tolerance (error);
mesh_refinement->refine_and_coarsen_elements();
}
else
{
ErrorVector error;
AutoPtr<ErrorEstimator> error_estimator =
build_error_estimator(param);
error_estimator->estimate_error(system, error);
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)
{
system.solve();
write_output(equation_systems, a_step, "primal");
NumericVector<Number> &primal_solution = *system.solution;
SensitivityData sensitivities(qois, system, system.get_parameter_vector());
system.assemble_qoi_sides = true;
system.set_adjoint_already_solved(false);
system.adjoint_qoi_parameter_sensitivity(qois, system.get_parameter_vector(), sensitivities);
system.set_adjoint_already_solved(true);
GetPot infile_l_shaped("l-shaped.in");
Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
Number sensitivity_QoI_0_0_exact = infile_l_shaped("sensitivity_0_0", 0.0);
Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
Number sensitivity_QoI_0_1_exact = infile_l_shaped("sensitivity_0_1", 0.0);
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<<"Sensitivity of QoI one to Parameter one is "<<sensitivity_QoI_0_0_computed<<std::endl;
std::cout<<"Sensitivity of QoI one to Parameter two is "<<sensitivity_QoI_0_1_computed<<std::endl;
std::cout<< "The error in sensitivity QoI_0_0 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact << std::endl;
std::cout<< "The error in sensitivity QoI_0_1 is "
<< std::setprecision(17)
<< std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact << std::endl << std::endl;
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);
}
}
std::cerr << '[' << libMesh::processor_id()
<< "] Completing output." << std::endl;
#endif // #ifndef LIBMESH_ENABLE_AMR
return 0;
}
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(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_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-qoi.C without comments:
#include "L-qoi.h"
using namespace libMesh;
void LaplaceQoI::init_qoi( std::vector<Number>& sys_qoi )
{
sys_qoi.resize(1);
return;
}
void LaplaceQoI::element_qoi (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<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 xf = xyz[qp](0);
const Real yf = xyz[qp](1);
if(fabs(xf - 0.875) <= 0.125 && fabs(yf - 0.125) <= 0.125)
{
Number T = c.interior_value(0, qp);
dQoI_0 += JxW[qp] * T;
}
}
c.elem_qoi[0] = c.elem_qoi[0] + dQoI_0;
}
void LaplaceQoI::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 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");
parameters.push_back(infile("alpha_1", 1.0));
parameters.push_back(infile("alpha_2", 1.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] * (parameters[0] + (2.*parameters[1])) * ( 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] * (parameters[0] + (2.*parameters[1])) * ( 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;
}
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 (parameters[0] + (2.*parameters[1])) * pow(x1*x1 + x2*x2, 1./3.)*sin(2./3.*theta);
}
The console output of the program:
***************************************************************
* Running Example adjoints_ex2:
* 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_ex2/.libs/lt-example-devel
Reading in and building the mesh
Building system
Initializing systems
Mesh Information:
mesh_dimension()=2
spatial_dimension()=3
n_nodes()=225
n_local_nodes()=117
n_elem()=63
n_local_elem()=32
n_active_elem()=48
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()=225
n_local_dofs()=117
n_constrained_dofs()=0
n_local_constrained_dofs()=0
n_vectors()=1
n_matrices()=1
DofMap Sparsity
Average On-Processor Bandwidth <= 13.52
Average Off-Processor Bandwidth <= 0.924444
Maximum On-Processor Bandwidth <= 25
Maximum Off-Processor Bandwidth <= 16
DofMap Constraints
Number of DoF Constraints = 0
Number of Node Constraints = 0
Nonlinear solver converged, step 0, residual reduction 3.96886e-13 < 1e-09
Adaptive step 0, we have 48 active elements and 225 active dofs.
Sensitivity of QoI one to Parameter one is 0.00543563
Sensitivity of QoI one to Parameter two is 0.0108713
The relative error in sensitivity QoI_0_0 is 0.00027992896221159407
The relative error in sensitivity QoI_0_1 is 0.0002799289573966468
Refining Uniformly
Refined mesh to 192 active elements and 833 active dofs.
Nonlinear solver converged, step 0, residual reduction 4.1557950841412934e-12 < 1.0000000000000001e-09
Adaptive step 1, we have 192 active elements and 833 active dofs.
Sensitivity of QoI one to Parameter one is 0.0054364699173682631
Sensitivity of QoI one to Parameter two is 0.010872939833425197
The relative error in sensitivity QoI_0_0 is 0.00012463358065163853
The relative error in sensitivity QoI_0_1 is 0.00012463370124116011
Refining Uniformly
Refined mesh to 768 active elements and 3201 active dofs.
Nonlinear solver converged, step 0, residual reduction 2.5163960226169602e-11 < 1.0000000000000001e-09
Adaptive step 2, we have 768 active elements and 3201 active dofs.
Sensitivity of QoI one to Parameter one is 0.0054368750032082382
Sensitivity of QoI one to Parameter two is 0.010873750005157855
The relative error in sensitivity QoI_0_0 is 5.013020643404552e-05
The relative error in sensitivity QoI_0_1 is 5.0130322176555172e-05
Refining Uniformly
Refined mesh to 3072 active elements and 12545 active dofs.
Nonlinear solver converged, step 0, residual reduction 2.1399651493003946e-10 < 1.0000000000000001e-09
Adaptive step 3, we have 3072 active elements and 12545 active dofs.
Sensitivity of QoI one to Parameter one is 0.0054370390969481755
Sensitivity of QoI one to Parameter two is 0.010874078194381213
The error in sensitivity QoI_0_0 is 1.9950091241522059e-05
The error in sensitivity QoI_0_1 is 1.9950046653287726e-05
************************************************************************************************************************
*** WIDEN YOUR WINDOW TO 120 CHARACTERS. Use 'enscript -r -fCourier9' to print this document ***
************************************************************************************************************************
---------------------------------------------- PETSc Performance Summary: ----------------------------------------------
/workspace/libmesh/examples/adjoints/adjoints_ex2/.libs/lt-example-devel on a intel-12. named hbar.ices.utexas.edu with 2 processors, by benkirk Tue Feb 5 13:39:38 2013
Using Petsc Release Version 3.3.0, Patch 2, Fri Jul 13 15:42:00 CDT 2012
Max Max/Min Avg Total
Time (sec): 1.255e+01 1.00000 1.255e+01
Objects: 4.510e+02 1.00000 4.510e+02
Flops: 3.758e+08 1.05624 3.658e+08 7.315e+08
Flops/sec: 2.993e+07 1.05624 2.913e+07 5.827e+07
MPI Messages: 4.620e+02 1.00000 4.620e+02 9.240e+02
MPI Message Lengths: 4.267e+05 1.00000 9.237e+02 8.535e+05
MPI Reductions: 1.328e+03 1.00000
Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N flops
and VecAXPY() for complex vectors of length N --> 8N flops
Summary of Stages: ----- Time ------ ----- Flops ----- --- Messages --- -- Message Lengths -- -- Reductions --
Avg %Total Avg %Total counts %Total Avg %Total counts %Total
0: Main Stage: 1.2554e+01 100.0% 7.3152e+08 100.0% 9.240e+02 100.0% 9.237e+02 100.0% 1.327e+03 99.9%
------------------------------------------------------------------------------------------------------------------------
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 289 1.0 2.4266e-02 1.4 5.60e+07 1.0 0.0e+00 0.0e+00 2.9e+02 0 15 0 0 22 0 15 0 0 22 4576
KSPSetUp 16 1.0 1.7500e-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 4 1.0 1.9420e-01 1.0 2.53e+08 1.1 4.1e+02 8.8e+02 4.3e+02 2 67 44 42 33 2 67 44 42 33 2537
PCSetUp 16 1.0 1.8070e-01 1.1 8.44e+07 1.1 0.0e+00 0.0e+00 5.6e+01 1 22 0 0 4 1 22 0 0 4 891
PCSetUpOnBlocks 4 1.0 8.9184e-02 1.1 4.22e+07 1.1 0.0e+00 0.0e+00 2.0e+01 1 11 0 0 2 1 11 0 0 2 902
PCApply 310 1.0 2.0910e-01 1.1 2.41e+08 1.1 0.0e+00 0.0e+00 2.0e+01 2 64 0 0 2 2 64 0 0 2 2241
VecDot 8 1.0 6.6280e-05 1.3 3.40e+04 1.0 0.0e+00 0.0e+00 8.0e+00 0 0 0 0 1 0 0 0 0 1 1014
VecMDot 289 1.0 1.6466e-02 1.8 2.80e+07 1.0 0.0e+00 0.0e+00 2.9e+02 0 8 0 0 22 0 8 0 0 22 3371
VecNorm 326 1.0 1.5117e-02 1.7 2.11e+06 1.0 0.0e+00 0.0e+00 3.3e+02 0 1 0 0 25 0 1 0 0 25 276
VecScale 318 1.0 4.6706e-04 1.0 1.04e+06 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 4397
VecCopy 97 1.0 1.8740e-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 605 1.0 9.3842e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecAXPY 38 1.0 1.2898e-04 1.1 2.27e+05 1.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 3479
VecMAXPY 302 1.0 7.9386e-03 1.0 2.99e+07 1.0 0.0e+00 0.0e+00 0.0e+00 0 8 0 0 0 0 8 0 0 0 7472
VecAssemblyBegin 165 1.0 1.4161e-01 2.5 0.00e+00 0.0 7.6e+01 4.4e+02 4.6e+02 1 0 8 4 35 1 0 8 4 35 0
VecAssemblyEnd 165 1.0 1.4377e-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
VecScatterBegin 378 1.0 9.7036e-04 1.0 0.00e+00 0.0 7.6e+02 8.9e+02 0.0e+00 0 0 82 79 0 0 0 82 79 0 0
VecScatterEnd 378 1.0 7.0043e-03 7.6 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0
VecNormalize 302 1.0 8.5912e-03 4.1 3.01e+06 1.0 0.0e+00 0.0e+00 3.0e+02 0 1 0 0 23 0 1 0 0 23 694
MatMult 205 1.0 2.0582e-02 1.3 2.21e+07 1.0 4.1e+02 8.8e+02 0.0e+00 0 6 44 42 0 0 6 44 42 0 2116
MatMultTranspose 97 1.0 7.3831e-03 1.0 8.81e+06 1.0 1.9e+02 7.6e+02 0.0e+00 0 2 21 17 0 0 2 21 17 0 2356
MatSolve 209 1.0 7.1565e-02 1.1 1.43e+08 1.1 0.0e+00 0.0e+00 0.0e+00 1 38 0 0 0 1 38 0 0 0 3898
MatSolveTranspos 101 1.0 4.4324e-02 1.0 5.60e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 15 0 0 0 0 15 0 0 0 2462
MatLUFactorNum 8 1.0 5.8062e-02 1.2 8.44e+07 1.1 0.0e+00 0.0e+00 0.0e+00 0 22 0 0 0 0 22 0 0 0 2772
MatILUFactorSym 8 1.0 1.2128e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.4e+01 1 0 0 0 2 1 0 0 0 2 0
MatAssemblyBegin 20 1.0 1.8247e-0231.0 0.00e+00 0.0 3.6e+01 3.6e+03 4.0e+01 0 0 4 15 3 0 0 4 15 3 0
MatAssemblyEnd 20 1.0 1.7569e-03 1.2 0.00e+00 0.0 1.6e+01 1.6e+02 3.2e+01 0 0 2 0 2 0 0 2 0 2 0
MatGetRowIJ 8 1.0 3.8147e-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 8 1.0 1.9097e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 1.6e+01 0 0 0 0 1 0 0 0 0 1 0
MatZeroEntries 20 1.0 5.0545e-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 16 16 154880 0
Preconditioner 16 16 14272 0
Vector 334 334 6195456 0
Vector Scatter 14 14 14504 0
Index Set 46 46 105520 0
IS L to G Mapping 4 4 2256 0
Matrix 20 20 21520096 0
Viewer 1 0 0 0
========================================================================================================================
Average time to get PetscTime(): 9.53674e-08
Average time for MPI_Barrier(): 1.38283e-06
Average time for zero size MPI_Send(): 1.34706e-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:38 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=12.5666, Active time=12.3904 |
----------------------------------------------------------------------------------------------------------------
| 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() 4 0.1012 0.025303 0.1128 0.028209 0.82 0.91 |
| build_sparsity() 4 0.0695 0.017377 0.1896 0.047401 0.56 1.53 |
| create_dof_constraints() 4 0.0105 0.002627 0.0105 0.002627 0.08 0.08 |
| distribute_dofs() 4 0.1208 0.030191 0.2835 0.070863 0.97 2.29 |
| dof_indices() 67474 3.9460 0.000058 3.9460 0.000058 31.85 31.85 |
| old_dof_indices() 8064 0.4664 0.000058 0.4664 0.000058 3.76 3.76 |
| prepare_send_list() 4 0.0003 0.000068 0.0003 0.000068 0.00 0.00 |
| reinit() 4 0.1617 0.040413 0.1617 0.040413 1.30 1.30 |
| |
| EquationSystems |
| build_solution_vector() 8 0.0184 0.002299 0.2529 0.031616 0.15 2.04 |
| |
| FE |
| compute_shape_functions() 60896 4.1161 0.000068 4.1161 0.000068 33.22 33.22 |
| init_shape_functions() 3888 0.0854 0.000022 0.0854 0.000022 0.69 0.69 |
| inverse_map() 28094 0.2199 0.000008 0.2199 0.000008 1.77 1.77 |
| |
| FEMSystem |
| assemble_qoi() 20 0.1847 0.009233 3.4961 0.174807 1.49 28.22 |
| assemble_qoi_derivative() 4 0.0692 0.017290 0.7161 0.179019 0.56 5.78 |
| assembly() 8 0.3756 0.046954 1.5139 0.189239 3.03 12.22 |
| assembly(get_jacobian) 4 0.1835 0.045876 0.7500 0.187493 1.48 6.05 |
| assembly(get_residual) 20 0.3889 0.019446 3.2105 0.160525 3.14 25.91 |
| |
| FEMap |
| compute_affine_map() 60896 0.5925 0.000010 0.5925 0.000010 4.78 4.78 |
| compute_face_map() 3776 0.0627 0.000017 0.1559 0.000041 0.51 1.26 |
| init_face_shape_functions() 64 0.0005 0.000009 0.0005 0.000009 0.00 0.00 |
| init_reference_to_physical_map() 3888 0.0855 0.000022 0.0855 0.000022 0.69 0.69 |
| |
| GMVIO |
| write_nodal_data() 8 0.0986 0.012324 0.0990 0.012373 0.80 0.80 |
| |
| ImplicitSystem |
| adjoint_solve() 4 0.0004 0.000112 1.6351 0.408763 0.00 13.20 |
| |
| LocationMap |
| find() 20460 0.0767 0.000004 0.0767 0.000004 0.62 0.62 |
| init() 8 0.0176 0.002194 0.0176 0.002194 0.14 0.14 |
| |
| Mesh |
| all_second_order() 1 0.0002 0.000234 0.0002 0.000234 0.00 0.00 |
| contract() 3 0.0021 0.000688 0.0136 0.004540 0.02 0.11 |
| find_neighbors() 6 0.0678 0.011295 0.0685 0.011412 0.55 0.55 |
| renumber_nodes_and_elem() 3 0.0116 0.003852 0.0116 0.003852 0.09 0.09 |
| |
| MeshCommunication |
| broadcast() 1 0.0007 0.000708 0.0010 0.000992 0.01 0.01 |
| compute_hilbert_indices() 7 0.0170 0.002427 0.0170 0.002427 0.14 0.14 |
| find_global_indices() 7 0.0070 0.000996 0.0267 0.003820 0.06 0.22 |
| parallel_sort() 7 0.0018 0.000260 0.0021 0.000296 0.01 0.02 |
| |
| MeshOutput |
| write_equation_systems() 8 0.0003 0.000036 0.3526 0.044077 0.00 2.85 |
| |
| MeshRefinement |
| _coarsen_elements() 3 0.0006 0.000216 0.0007 0.000222 0.01 0.01 |
| _refine_elements() 8 0.0812 0.010144 0.2317 0.028962 0.65 1.87 |
| add_point() 20460 0.0674 0.000003 0.1466 0.000007 0.54 1.18 |
| make_coarsening_compatible() 6 0.0170 0.002838 0.0170 0.002838 0.14 0.14 |
| make_refinement_compatible() 6 0.0000 0.000003 0.0000 0.000005 0.00 0.00 |
| |
| MetisPartitioner |
| partition() 6 0.0512 0.008540 0.0775 0.012913 0.41 0.63 |
| |
| NewtonSolver |
| solve() 4 0.0461 0.011516 1.6530 0.413239 0.37 13.34 |
| |
| Parallel |
| allgather() 28 0.0004 0.000013 0.0004 0.000015 0.00 0.00 |
| broadcast() 9 0.0002 0.000025 0.0002 0.000020 0.00 0.00 |
| max(bool) 19 0.0013 0.000069 0.0013 0.000069 0.01 0.01 |
| max(scalar) 865 0.0024 0.000003 0.0024 0.000003 0.02 0.02 |
| max(vector) 201 0.0012 0.000006 0.0028 0.000014 0.01 0.02 |
| min(bool) 1056 0.0025 0.000002 0.0025 0.000002 0.02 0.02 |
| min(scalar) 850 0.0773 0.000091 0.0773 0.000091 0.62 0.62 |
| min(vector) 201 0.0014 0.000007 0.0031 0.000015 0.01 0.02 |
| probe() 50 0.0006 0.000013 0.0006 0.000013 0.01 0.01 |
| receive() 50 0.0004 0.000007 0.0010 0.000020 0.00 0.01 |
| send() 50 0.0002 0.000003 0.0002 0.000003 0.00 0.00 |
| send_receive() 64 0.0006 0.000009 0.0019 0.000029 0.00 0.02 |
| sum() 105 0.0013 0.000013 0.0748 0.000712 0.01 0.60 |
| |
| Parallel::Request |
| wait() 50 0.0001 0.000002 0.0001 0.000002 0.00 0.00 |
| |
| Partitioner |
| set_node_processor_ids() 6 0.0101 0.001683 0.0109 0.001823 0.08 0.09 |
| set_parent_processor_ids() 6 0.0056 0.000934 0.0056 0.000934 0.05 0.05 |
| |
| PetscLinearSolver |
| solve() 8 0.3641 0.045517 0.3641 0.045517 2.94 2.94 |
| |
| ProjectVector |
| operator() 6 0.0753 0.012554 0.6792 0.113193 0.61 5.48 |
| |
| System |
| project_vector() 6 0.0209 0.003490 0.9318 0.155306 0.17 7.52 |
----------------------------------------------------------------------------------------------------------------
| Totals: 281774 12.3904 100.00 |
----------------------------------------------------------------------------------------------------------------
***************************************************************
* Done Running Example adjoints_ex2:
* 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
***************************************************************