The source file adjoint_initial.h with comments:

        #include "libmesh/parameters.h"
        #include "libmesh/point.h"
        #include "libmesh/vector_value.h"
        
        using namespace libMesh;
        
        void adjoint_read_initial_parameters();
        void adjoint_finish_initialization();
        
        Number adjoint_initial_value(const Point& p,
                             const Parameters&,
                             const std::string&,
                             const std::string&);
        
        Gradient adjoint_initial_grad(const Point& p,
                              const Parameters&,
                              const std::string&,
                              const std::string&);



The source file femparameters.h with comments:

        #ifndef __fem_parameters_h__
        #define __fem_parameters_h__
        
        #include "libmesh/libmesh_common.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/enum_norm_type.h"
        #include "libmesh/function_base.h"
        #include "libmesh/getpot.h"
        #include "libmesh/id_types.h"
        #include "libmesh/periodic_boundaries.h"
        
        #include <limits>
        #include <map>
        #include <string>
        #include <vector>
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        class FEMParameters
        {
        public:
            FEMParameters();
        
            ~FEMParameters();
        
            void read(GetPot &input);
        
Parameters applicable to entire EquationSystems:

            unsigned int initial_timestep, n_timesteps;
            bool transient;
            unsigned int deltat_reductions;
            std::string timesolver_core;
            Real end_time, deltat, timesolver_theta,
        	 timesolver_maxgrowth, timesolver_tolerance,
        	 timesolver_upper_tolerance, steadystate_tolerance;
            std::vector<FEMNormType> timesolver_norm;
        
Mesh generation

            unsigned int dimension;
            std::string domaintype, domainfile, elementtype;
            Real elementorder;
            Real domain_xmin, domain_ymin, domain_zmin;
            Real domain_edge_width, domain_edge_length, domain_edge_height;
            unsigned int coarsegridx, coarsegridy, coarsegridz;
            unsigned int coarserefinements, extrarefinements;
        
Mesh refinement

            unsigned int nelem_target;
            Real global_tolerance;
            Real refine_fraction, coarsen_fraction, coarsen_threshold;
            unsigned int max_adaptivesteps;
            unsigned int initial_adaptivesteps;
        
Output

            unsigned int write_interval;
            bool write_gmv_error, write_tecplot_error, output_xda, output_xdr,
        	 output_bz2, output_gz, output_gmv, output_tecplot;
        
Types of Systems to create

            std::vector<std::string> system_types;
        
Parameters applicable to each system:

Boundary and initial conditions

std::vector periodic_boundaries;

            std::map<subdomain_id_type, FunctionBase<Number> *> initial_conditions;
            std::map<boundary_id_type, FunctionBase<Number> *>  dirichlet_conditions,
                                                                neumann_conditions;
            std::map<boundary_id_type, std::vector<unsigned int> > dirichlet_condition_variables,
                                                                   neumann_condition_variables;
            std::map<int, std::map<subdomain_id_type, FunctionBase<Number> *> >
              other_interior_functions;
            std::map<int, std::map<boundary_id_type, FunctionBase<Number> *> >
              other_boundary_functions;
        
Execution type

            bool run_simulation, run_postprocess;
        
Approximation type

            std::vector<std::string> fe_family;
            std::vector<unsigned int> fe_order;
            int extra_quadrature_order;
        
Assembly options

            bool analytic_jacobians;
            Real verify_analytic_jacobians;
            Real numerical_jacobian_h;
        
            bool print_solution_norms, print_solutions,
                 print_residual_norms, print_residuals,
                 print_jacobian_norms, print_jacobians,
                 print_element_jacobians;
        
Solver options

            bool use_petsc_snes;
            bool time_solver_quiet, solver_quiet, solver_verbose, 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;
        
Initialization

            unsigned int initial_sobolev_order;
            unsigned int initial_extra_quadrature;
        
Error indicators

            std::string indicator_type;
            unsigned int sobolev_order;
        
System-specific parameters file

            std::string system_config_file;
        };
        
        #endif // __fem_parameters_h__



The source file heatsystem.h with comments:

        #include "libmesh/enum_fe_family.h"
        #include "libmesh/fem_system.h"
        #include "libmesh/parameter_vector.h"
        
FEMSystem, TimeSolver and NewtonSolver will handle most tasks, but we must specify element residuals
        class HeatSystem : public FEMSystem
        {
        public:
Constructor
          HeatSystem(EquationSystems& es,
                       const std::string& name_in,
                       const unsigned int number_in)
          : FEMSystem(es, name_in, number_in),
            _k(1.0),
            _fe_family("LAGRANGE"), _fe_order(1),
            _analytic_jacobians(true), R_plus_dp(0.0), R_minus_dp(0.0), dp(1.e-6) { qoi.resize(1); }
        
          std::string & fe_family() { return _fe_family;  }
          unsigned int & fe_order() { return _fe_order;  }
          Real & k() { return _k;  }
          bool & analytic_jacobians() { return _analytic_jacobians; }
        
A function to compute and accumulate residuals
          void perturb_accumulate_residuals(const ParameterVector& parameters);
        
Sensitivity Calculation
          Number & compute_final_sensitivity()
            {
              final_sensitivity = -(R_plus_dp - R_minus_dp)/(2*dp);
        
              return final_sensitivity;
            }
        
            void set_tf(Real val)
          {
            tf = val;
          }
        
            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;
            }
        
          Number &get_QoI_value(unsigned int QoI_index)
            {
              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);

RHS for adjoint problem
          virtual void element_qoi_derivative (DiffContext &context,
        				       const QoISet & /* qois */);
        
virtual void element_qoi (DiffContext &context, const QoISet & qois);

Parameters associated with the system
          std::vector<Number> parameters;
        
The ParameterVector object that will contain pointers to the system parameters
          ParameterVector parameter_vector;
        
The parameters to solve for
          Real _k;
        
The final time parameter
          Real tf;
        
Variables to hold the computed QoIs
          Number computed_QoI[1];
        
The FE type to use
          std::string _fe_family;
          unsigned int _fe_order;
        
Index for T variable
          unsigned int T_var;
        
Calculate Jacobians analytically or not?
          bool _analytic_jacobians;
        
Variables to hold the perturbed residuals
          Number R_plus_dp;
          Number R_minus_dp;
        
Perturbation parameter
          Real dp;
        
The final computed sensitivity
          Number final_sensitivity;
        
        };



The source file initial.h with comments:

        #include "libmesh/parameters.h"
        #include "libmesh/point.h"
        #include "libmesh/vector_value.h"
        
        using namespace libMesh;
        
        void read_initial_parameters();
        void finish_initialization();
        
        Number initial_value(const Point& p,
                             const Parameters&,
                             const std::string&,
                             const std::string&);
        
        Gradient initial_grad(const Point& p,
                              const Parameters&,
                              const std::string&,
                              const std::string&);



The source file adjoint_initial.C with comments:

        #include "adjoint_initial.h"
        
        using namespace libMesh;
        
        void adjoint_read_initial_parameters()
        {
        }
        
        void adjoint_finish_initialization()
        {
        }
        
        
        
Initial conditions
        Number adjoint_initial_value(const Point& p,
                             const Parameters&,
                             const std::string&,
                             const std::string&)
        {
          Real x = p(0), y = p(1);
        
          Number val = 0.;
        
          val = sin(M_PI * x) * sin(M_PI * y);
        
          return val;
        }
        
        
        
        Gradient adjoint_initial_grad(const Point& p,
                              const Parameters&,
                              const std::string&,
                              const std::string&)
        {
          Real x = p(0), y = p(1);
        
          return Gradient(M_PI*cos(M_PI * x) * sin(M_PI * y),
          M_PI*sin(M_PI * x) * cos(M_PI * y));
        }



The source file adjoints_ex5.C with comments:

partial(T)/partial(t) - K Laplacian(T) = 0 with initial condition: T(x,y;0) = sin(pi*x) sin(pi*y) and boundary conditions: T(boundary;t) = 0

For these initial and boundary conditions, the exact solution u = exp(-K pi^2 t) * sin(pi*x) * sin(pi*y)

We specify our Quantity of Interest (QoI) as Q(u) = int_{domain} u(x,y;1) sin(pi*x) sin(pi*y) dx dy, and are interested in computing the sensitivity dQ/dK

For this QoI, the continuous adjoint problem reads, -partial(z)/partial(t) - K Laplacian(z) = 0 with initial condition: T(x,y;1) = sin(pi*x) sin(pi*y) and boundary condition: T(boundary;t) = 0

which has the exact solution, z = exp(-K pi^2 (1 - t)) * sin(pi*x) * sin(pi*y) which is the mirror image in time of the forward solution

For an adjoint consistent space-time formulation, the discrete adjoint can be obtained by marching backwards from the adjoint initial condition and solving the transpose of the discrete primal problem at the last nonlinear solve of the corresponding primal timestep. This necessitates the storage of the primal solution at all timesteps, which is accomplished here using a MemorySolutionHistory object. As the name suggests, this object simply stores the primal solution (and other vectors we may choose to save), so that we can retrieve them later, whenever necessary.

The discrete adjoint system for implicit time steppers requires the localization of vectors other than system.solution, which is accomplished using the localize_vectors method. In this particular example, we use the localized adjoint solution to assemble the residual contribution for the current adjoint timestep from the last computed adjoint timestep.

Finally, The adjoint_advance_timestep method, the backwards time analog of advance_timestep prepares the time solver for solving the adjoint system, while the retrieve_timestep method retrieves the saved solutions at the current system.time, so that the adjoint sensitivity contribution for the current time can be computed.

Local includes
        #include "initial.h"
        #include "adjoint_initial.h"
        #include "femparameters.h"
        #include "heatsystem.h"
        
Libmesh includes
        #include "libmesh/equation_systems.h"
        #include "libmesh/dirichlet_boundaries.h"
        #include "libmesh/system_norm.h"
        #include "libmesh/numeric_vector.h"
        
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/mesh_refinement.h"
        
        #include "libmesh/petsc_diff_solver.h"
        #include "libmesh/steady_solver.h"
        #include "libmesh/euler_solver.h"
        #include "libmesh/euler2_solver.h"
        #include "libmesh/newton_solver.h"
        #include "libmesh/twostep_time_solver.h"
        
        #include "libmesh/getpot.h"
        #include "libmesh/tecplot_io.h"
        #include "libmesh/gmv_io.h"
        
SolutionHistory Includes
        #include "libmesh/solution_history.h"
        #include "libmesh/memory_solution_history.h"
        
C++ includes
        #include <iostream>
        #include <sys/time.h>
        #include <iomanip>
        
        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(HeatSystem &system, FEMParameters ¶m)
        {
Use the prescribed FE type
          system.fe_family() = param.fe_family[0];
          system.fe_order() = param.fe_order[0];
        
Use analytical jacobians?
          system.analytic_jacobians() = param.analytic_jacobians;
        
Verify analytic jacobians against numerical ones?
          system.verify_analytic_jacobians = param.verify_analytic_jacobians;
          system.numerical_jacobian_h = param.numerical_jacobian_h;
        
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;
          system.print_element_jacobians = param.print_element_jacobians;
        
Solve this as a time-dependent or steady system
          if (param.transient)
            {
              UnsteadySolver *innersolver;
               if (param.timesolver_core == "euler")
                {
                  EulerSolver *eulersolver =
                    new EulerSolver(system);
        
                  eulersolver->theta = param.timesolver_theta;
                  innersolver = eulersolver;
                }
              else
                {
                  std::cerr << "This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods."
                            <<  std::endl;
                  libmesh_error();
                }
               system.time_solver =
                  AutoPtr<TimeSolver>(innersolver);
            }
          else
            system.time_solver =
              AutoPtr<TimeSolver>(new SteadySolver(system));
        
The Memory Solution History object we will set the system SolutionHistory object to
          MemorySolutionHistory heatsystem_solution_history(system);
          system.time_solver->set_solution_history(heatsystem_solution_history);
        
          system.time_solver->reduce_deltat_on_diffsolver_failure =
                                                param.deltat_reductions;
          system.time_solver->quiet           = param.time_solver_quiet;
        
Create any Dirichlet boundary conditions
          typedef
            std::map<boundary_id_type, FunctionBase<Number> *>::
            const_iterator Iter;
        
          for (Iter i = param.dirichlet_conditions.begin();
               i != param.dirichlet_conditions.end(); ++i)
            {
              boundary_id_type b = i->first;
              FunctionBase<Number> *f = i->second;
              std::set<boundary_id_type> bdys; bdys.insert(b);
              system.get_dof_map().add_dirichlet_boundary
                (DirichletBoundary
                   (bdys, param.dirichlet_condition_variables[b], f));
              std::cout << "Added Dirichlet boundary " << b << " for variables ";
              for (unsigned int vi=0; vi !=
                   param.dirichlet_condition_variables[b].size(); ++vi)
                std::cout << param.dirichlet_condition_variables[b][vi];
              std::cout << std::endl;
            }
        
Set the time stepping options
          system.deltat = param.deltat;
        
And the integration options
          system.extra_quadrature_order = param.extra_quadrature_order;
        
And the nonlinear solver options
          if (param.use_petsc_snes)
            {
        #ifdef LIBMESH_HAVE_PETSC
              PetscDiffSolver *solver = new PetscDiffSolver(system);
              system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
        #else
              libmesh_error();
        #endif
            }
          else
            {
              NewtonSolver *solver = new NewtonSolver(system);
              system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
        
              solver->quiet                       = param.solver_quiet;
              solver->verbose                     = param.solver_verbose;
              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;
            }
        }
        
The main program.
        int main (int argc, char** argv)
        {
Skip adaptive examples on a non-adaptive libMesh build
        #ifndef LIBMESH_ENABLE_AMR
          libmesh_example_assert(false, "--enable-amr");
        #else
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
          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);
        
Create a mesh with the given dimension, distributed across the default MPI communicator.
          Mesh mesh(init.comm(), param.dimension);
        
And an object to refine it
          AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
        
And an EquationSystems to run on it
          EquationSystems equation_systems (mesh);
        
          std::cout << "Building mesh" << std::endl;
        
Build a unit square
          ElemType elemtype;
        
          if (param.elementtype == "tri" ||
              param.elementtype == "unstructured")
            elemtype = TRI3;
          else
            elemtype = QUAD4;
        
          MeshTools::Generation::build_square
            (mesh, param.coarsegridx, param.coarsegridy,
             param.domain_xmin, param.domain_xmin + param.domain_edge_width,
             param.domain_ymin, param.domain_ymin + param.domain_edge_length,
             elemtype);
        
          std::cout << "Building system" << std::endl;
        
          HeatSystem &system = equation_systems.add_system<HeatSystem> ("HeatSystem");
        
          set_system_parameters(system, param);
        
          std::cout << "Initializing systems" << std::endl;
        
Initialize the system
          equation_systems.init ();
        
Refine the grid again if requested
          for (unsigned int i=0; i != param.extrarefinements; ++i)
            {
              mesh_refinement->uniformly_refine(1);
              equation_systems.reinit();
            }
        
          std::cout<<"Setting primal initial conditions"<<std::endl;
        
          read_initial_parameters();
        
          system.project_solution(initial_value, initial_grad,
                                          equation_systems.parameters);
        
Output the H1 norm of the initial conditions
          libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl<<std::endl;
        
Add an adjoint vector, this will be computed after the forward time stepping is complete

Tell the library not to save adjoint solutions during the forward solve

Tell the library not to project this vector, and hence, memory solution history to not save it.

Make this vector ghosted so we can localize it to each element later.
          const std::string & adjoint_solution_name = "adjoint_solution0";
          system.add_vector("adjoint_solution0", false, GHOSTED);
        
Close up any resources initial.C needed
          finish_initialization();
        
Plot the initial conditions
          write_output(equation_systems, 0, "primal");
        
Print information about the mesh and system to the screen.
          mesh.print_info();
          equation_systems.print_info();
        
In optimized mode we catch any solver errors, so that we can write the proper footers before closing. In debug mode we just let the exception throw so that gdb can grab it.
        #ifdef NDEBUG
          try
          {
        #endif
Now we begin the timestep loop to compute the time-accurate solution of the equations.
            for (unsigned int t_step=param.initial_timestep;
                 t_step != param.initial_timestep + param.n_timesteps; ++t_step)
              {
A pretty update message
                std::cout << " Solving time step " << t_step << ", time = "
                          << system.time << std::endl;
        
Solve the forward problem at time t, to obtain the solution at time t + dt
                system.solve();
        
Output the H1 norm of the computed solution
                libMesh::out << "|U(" <<system.time + system.deltat<< ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
        
Advance to the next timestep in a transient problem
                std::cout<<"Advancing timestep"<<std::endl<<std::endl;
                system.time_solver->advance_timestep();
        
Write out this timestep
                  write_output(equation_systems, t_step+1, "primal");
              }
End timestep loop

/////////////// Now for the Adjoint Solution //////////////////////////////////////

Now we will solve the backwards in time adjoint problem
          std::cout << std::endl << "Solving the adjoint problem" << std::endl;
        
We need to tell the library that it needs to project the adjoint, so MemorySolutionHistory knows it has to save it

Tell the library to project the adjoint vector, and hence, memory solution history to save it
          system.set_vector_preservation(adjoint_solution_name, true);
        
          std::cout << "Setting adjoint initial conditions Z("<<system.time<<")"<<std::endl;
        
Need to call adjoint_advance_timestep once for the initial condition setup
          std::cout<<"Retrieving solutions at time t="<<system.time<<std::endl;
          system.time_solver->adjoint_advance_timestep();
        
Output the H1 norm of the retrieved solutions (u^i and u^i+1)
          libMesh::out << "|U(" <<system.time + system.deltat<< ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
        
          libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
        
The first thing we have to do is to apply the adjoint initial condition. The user should supply these. Here they are specified in the functions adjoint_initial_value and adjoint_initial_gradient
          system.project_vector(adjoint_initial_value, adjoint_initial_grad, equation_systems.parameters, system.get_adjoint_solution(0));
        
Since we have specified an adjoint solution for the current time (T), set the adjoint_already_solved boolean to true, so we dont solve unneccesarily in the adjoint sensitivity method
          system.set_adjoint_already_solved(true);
        
          libMesh::out << "|Z(" <<system.time<< ")|= " << system.calculate_norm(system.get_adjoint_solution(), 0, H1) << std::endl<<std::endl;
        
          write_output(equation_systems, param.n_timesteps, "dual");
        
Now that the adjoint initial condition is set, we will start the backwards in time adjoint integration

For loop stepping backwards in time
          for (unsigned int t_step=param.initial_timestep;
               t_step != param.initial_timestep + param.n_timesteps; ++t_step)
            {
A pretty update message
              std::cout << " Solving adjoint time step " << t_step << ", time = "
        		<< system.time << std::endl;
        
The adjoint_advance_timestep function calls the retrieve function of the memory_solution_history class via the memory_solution_history object we declared earlier. The retrieve function sets the system primal vectors to their values at the current timestep
              std::cout<<"Retrieving solutions at time t="<<system.time<<std::endl;
              system.time_solver->adjoint_advance_timestep();
        
Output the H1 norm of the retrieved solution
              libMesh::out << "|U(" <<system.time + system.deltat << ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
        
              libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
        
              system.set_adjoint_already_solved(false);
        
              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);
        
              libMesh::out << "|Z(" <<system.time<< ")|= "<< system.calculate_norm(system.get_adjoint_solution(), 0, H1) << std::endl << std::endl;
        
Get a pointer to the primal solution vector
              NumericVector<Number> &primal_solution = *system.solution;
        
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, param.n_timesteps - (t_step + 1), "dual");
        
Swap back
              primal_solution.swap(dual_solution_0);
            }
End adjoint timestep loop

Now that we have computed both the primal and adjoint solutions, we compute the sensitivties to the parameter p dQ/dp = partialQ/partialp - partialR/partialp partialQ/partialp = (Q(p+dp) - Q(p-dp))/(2*dp), this is not supported by the library yet partialR/partialp = (R(u,z;p+dp) - R(u,z;p-dp))/(2*dp), where R(u,z;p+dp) = int_{0}^{T} f(z;p+dp) - (p+dp) - (p+dp) To do this we need to step forward in time, and compute the perturbed R at each time step and accumulate it Then once all time steps are over, we can compute (R(u,z;p+dp) - R(u,z;p-dp))/(2*dp)

Now we begin the timestep loop to compute the time-accurate adjoint sensitivities
            for (unsigned int t_step=param.initial_timestep;
                 t_step != param.initial_timestep + param.n_timesteps; ++t_step)
              {
A pretty update message
                std::cout << "Retrieving " << t_step << ", time = "
                          << system.time << std::endl;
        
Retrieve the primal and adjoint solutions at the current timestep
                system.time_solver->retrieve_timestep();
        
        	libMesh::out << "|U(" <<system.time + system.deltat << ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
        
              libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
        
              libMesh::out << "|Z(" <<system.time<< ")|= "<< system.calculate_norm(system.get_adjoint_solution(0), 0, H1) << std::endl << std::endl;
        
Call the postprocess function which we have overloaded to compute accumulate the perturbed residuals
                (dynamic_cast<HeatSystem&>(system)).perturb_accumulate_residuals(dynamic_cast<HeatSystem&>(system).get_parameter_vector());
        
Move the system time forward (retrieve_timestep does not do this)
                system.time += system.deltat;
              }
        
A pretty update message
            std::cout << "Retrieving " << " final time = "
        	      << system.time << std::endl;
        
Retrieve the primal and adjoint solutions at the current timestep
            system.time_solver->retrieve_timestep();
        
            libMesh::out << "|U(" <<system.time + system.deltat << ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
        
              libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
        
              libMesh::out << "|Z(" <<system.time<< ")|= "<< system.calculate_norm(system.get_adjoint_solution(0), 0, H1) << std::endl<<std::endl;
        
Call the postprocess function which we have overloaded to compute accumulate the perturbed residuals
            (dynamic_cast<HeatSystem&>(system)).perturb_accumulate_residuals(dynamic_cast<HeatSystem&>(system).get_parameter_vector());
        
Now that we computed the accumulated, perturbed residuals, we can compute the approximate sensitivity
            Number sensitivity_0_0 = (dynamic_cast<HeatSystem&>(system)).compute_final_sensitivity();
        
Print it out
            std::cout<<"Sensitivity of QoI 0 w.r.t parameter 0 is: " << sensitivity_0_0 << std::endl;
        
        #ifdef NDEBUG
        }
          catch (...)
          {
            std::cerr << '[' << libMesh::processor_id()
                      << "] Caught exception; exiting early." << std::endl;
          }
        #endif
        
          std::cerr << '[' << libMesh::processor_id()
                    << "] Completing output." << std::endl;
        
All done.
          return 0;
        
        #endif // LIBMESH_ENABLE_AMR
        }



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 "heatsystem.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 HeatSystem::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
          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 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];
        
A reference to the system context is built with
          const System & sys = c.get_system();
        
Get a pointer to the adjoint solution vector
          NumericVector<Number> &adjoint_solution = const_cast<System &>(sys).get_adjoint_solution(0);
        
Get the previous adjoint solution values at all the qps

          std::vector<Number> old_adjoint (n_qpoints, 0);
        
          c.interior_values<Number>(0, adjoint_solution, old_adjoint);
        
Our QoI depends solely on the final time, so there are no QoI contributions. However, there is a contribution from the adjoint solution timestep, for the time part of the residual of the adjoint problem Loop over the qps
          for (unsigned int qp=0; qp != n_qpoints; qp++)
            {
              for (unsigned int i=0; i != n_T_dofs; i++)
        	{
        	  Q(i) += -JxW[qp] * old_adjoint[qp] * phi[i][qp] ;
              	}
        
            } // end of the quadrature point qp-loop
        }



The source file factoryfunction.C with comments:

        #include "libmesh/dense_vector.h"
        #include "libmesh/factory.h"
        #include "libmesh/function_base.h"
        
        using namespace libMesh;
        
An example of a hand-coded function

        class ExampleOneFunction : public FunctionBase<Number>
        {
          virtual Number operator() (const Point&  /*p*/,
                                     const Real /*time*/)
            {
              return 1;
            }
        
          virtual void operator() (const Point&  /*p*/,
                                   const Real /*time*/,
                                   DenseVector<Number>& output)
            {
              for (unsigned int i=0; i != output.size(); ++i)
                output(i) = 1;
            }
        
          virtual void init() {}
          virtual void clear() {}
          virtual AutoPtr<FunctionBase<Number> > clone() const {
            return AutoPtr<FunctionBase<Number> >
              (new ExampleOneFunction());
          }
        };
        
        #ifdef LIBMESH_USE_SEPARATE_NAMESPACE
        namespace libMesh {
        #endif
        
------------------------------------------------- Full specialization for the Factory > So we can look up hand-coded functions by name string
        template<>
        std::map<std::string, Factory<FunctionBase<Number> >*>&
        Factory<FunctionBase<Number> >::factory_map()
        {
          static std::map<std::string, Factory<FunctionBase<Number> >*> _map;
          return _map;
        }
        
        FactoryImp<ExampleOneFunction, FunctionBase<Number> > example_one_factory ("example_one");
        
        #ifdef LIBMESH_USE_SEPARATE_NAMESPACE
        } // namespace libMesh
        #endif



The source file femparameters.C with comments:

        #include "femparameters.h"
        
        #include "libmesh/factory.h"
        #include "libmesh/parsed_function.h"
        #include "libmesh/zero_function.h"
        #include "libmesh/periodic_boundaries.h"
        #include "libmesh/vector_value.h" // RealVectorValue
        
        #define GETPOT_INPUT(A) { A = input(#A, A);\
          const 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);\
          const std::string stringval = input(#A, std::string());\
          variable_names.push_back(std::string(#A "=") + stringval); }
        
        #define GETPOT_FUNCTION_INPUT(A) {\
          const std::string type  = input(#A "_type", "zero");\
          const std::string value = input(#A "_value", "");\
          A = new_function_base(type, value); }
        #define GETPOT_REGISTER(A) { \
          std::string stringval = input(#A, std::string());\
          variable_names.push_back(std::string(#A "=") + stringval); }
        
        FEMParameters::FEMParameters() :
              initial_timestep(0), n_timesteps(100),
              transient(true),
              deltat_reductions(0),
              timesolver_core("euler"),
              end_time(std::numeric_limits<Real>::max()),
              deltat(0.0001), timesolver_theta(0.5),
              timesolver_maxgrowth(0.), timesolver_tolerance(0.),
              timesolver_upper_tolerance(0.),
              steadystate_tolerance(0.),
              timesolver_norm(0,L2),
        
              dimension(2),
              domaintype("square"), domainfile("mesh.xda"), elementtype("quad"),
              elementorder(2),
              domain_xmin(0.0), domain_ymin(0.0), domain_zmin(0.0),
              domain_edge_width(1.0), domain_edge_length(1.0), domain_edge_height(1.0),
              coarsegridx(1), coarsegridy(1), coarsegridz(1),
              coarserefinements(4), extrarefinements(0),
        
              nelem_target(8000), global_tolerance(0.0),
              refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
              max_adaptivesteps(1),
              initial_adaptivesteps(0),
        
              write_interval(10),
              write_gmv_error(false), write_tecplot_error(false),
              output_xda(false), output_xdr(false),
              output_bz2(true), output_gz(true),
              output_gmv(false), output_tecplot(false),
        
              system_types(0),
        
periodic_boundaries(0),

              run_simulation(true), run_postprocess(false),
        
              fe_family(1, "LAGRANGE"), fe_order(1, 1),
              extra_quadrature_order(0),
        
              analytic_jacobians(true), verify_analytic_jacobians(0.0),
              numerical_jacobian_h(TOLERANCE),
              print_solution_norms(false), print_solutions(false),
              print_residual_norms(false), print_residuals(false),
              print_jacobian_norms(false), print_jacobians(false),
              print_element_jacobians(false),
        
              use_petsc_snes(false),
              time_solver_quiet(true), solver_quiet(true), solver_verbose(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),
        
              initial_sobolev_order(1),
              initial_extra_quadrature(0),
              indicator_type("kelly"), sobolev_order(1)
        {
        }
        
        FEMParameters::~FEMParameters()
        {
          for (std::map<subdomain_id_type, FunctionBase<Number> *>::iterator
               i = initial_conditions.begin(); i != initial_conditions.end();
               ++i)
            delete i->second;
        
          for (std::map<boundary_id_type, FunctionBase<Number> *>::iterator
               i = dirichlet_conditions.begin(); i != dirichlet_conditions.end();
               ++i)
            delete i->second;
        
          for (std::map<boundary_id_type, FunctionBase<Number> *>::iterator
               i = neumann_conditions.begin(); i != neumann_conditions.end();
               ++i)
            delete i->second;
        
          for (std::map<int,
        	 std::map<boundary_id_type, FunctionBase<Number> *>
               >::iterator
               i = other_boundary_functions.begin(); i != other_boundary_functions.end();
               ++i)
            for (std::map<boundary_id_type, FunctionBase<Number> *>::iterator
                 j = i->second.begin(); j != i->second.end();
                 ++j)
              delete j->second;
        
          for (std::map<int,
        	 std::map<subdomain_id_type, FunctionBase<Number> *>
               >::iterator
               i = other_interior_functions.begin(); i != other_interior_functions.end();
               ++i)
            for (std::map<subdomain_id_type, FunctionBase<Number> *>::iterator
                 j = i->second.begin(); j != i->second.end();
                 ++j)
              delete j->second;
        }
        
        
        AutoPtr<FunctionBase<Number> > new_function_base(const std::string& func_type,
                                                const std::string& func_value)
        {
          if (func_type == "parsed")
            return AutoPtr<FunctionBase<Number> >
              (new ParsedFunction<Number>(func_value));
          else if (func_type == "factory")
            return AutoPtr<FunctionBase<Number> >
              (Factory<FunctionBase<Number> >::build(func_value).release());
          else if (func_type == "zero")
            return AutoPtr<FunctionBase<Number> >
              (new ZeroFunction<Number>);
          else
            libmesh_not_implemented();
        
          return AutoPtr<FunctionBase<Number> >();
        }
        
        
        void FEMParameters::read(GetPot &input)
        {
            std::vector<std::string> variable_names;
        
            GETPOT_INT_INPUT(initial_timestep);
            GETPOT_INT_INPUT(n_timesteps);
            GETPOT_INPUT(transient);
            GETPOT_INT_INPUT(deltat_reductions);
            GETPOT_INPUT(timesolver_core);
            GETPOT_INPUT(end_time);
            GETPOT_INPUT(deltat);
            GETPOT_INPUT(timesolver_theta);
            GETPOT_INPUT(timesolver_maxgrowth);
            GETPOT_INPUT(timesolver_tolerance);
            GETPOT_INPUT(timesolver_upper_tolerance);
            GETPOT_INPUT(steadystate_tolerance);
        
            GETPOT_REGISTER(timesolver_norm);
            const unsigned int n_timesolver_norm = input.vector_variable_size("timesolver_norm");
            timesolver_norm.resize(n_timesolver_norm, L2);
            for (unsigned int i=0; i != n_timesolver_norm; ++i)
              {
                int current_norm = 0; // L2
                if (timesolver_norm[i] == H1)
                  current_norm = 1;
                if (timesolver_norm[i] == H2)
                  current_norm = 2;
                current_norm = input("timesolver_norm", current_norm, i);
                if (current_norm == 0)
                  timesolver_norm[i] = L2;
                else if (current_norm == 1)
                  timesolver_norm[i] = H1;
                else if (current_norm == 2)
                  timesolver_norm[i] = H2;
                else
                  timesolver_norm[i] = DISCRETE_L2;
              }
        
        
            GETPOT_INT_INPUT(dimension);
            GETPOT_INPUT(domaintype);
            GETPOT_INPUT(domainfile);
            GETPOT_INPUT(elementtype);
            GETPOT_INPUT(elementorder);
            GETPOT_INPUT(domain_xmin);
            GETPOT_INPUT(domain_ymin);
            GETPOT_INPUT(domain_zmin);
            GETPOT_INPUT(domain_edge_width);
            GETPOT_INPUT(domain_edge_length);
            GETPOT_INPUT(domain_edge_height);
            GETPOT_INT_INPUT(coarsegridx);
            GETPOT_INT_INPUT(coarsegridy);
            GETPOT_INT_INPUT(coarsegridz);
            GETPOT_INT_INPUT(coarserefinements);
            GETPOT_INT_INPUT(extrarefinements);
        
        
            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_INT_INPUT(initial_adaptivesteps);
        
        
            GETPOT_INT_INPUT(write_interval);
            GETPOT_INPUT(output_xda);
            GETPOT_INPUT(output_xdr);
            GETPOT_INPUT(output_gz);
        #ifndef LIBMESH_HAVE_GZSTREAM
            output_gz                   = false;
        #endif
            GETPOT_INPUT(output_bz2);
        #ifndef LIBMESH_HAVE_BZ2
            output_bz2                  = false;
        #endif
            GETPOT_INPUT(output_gmv);
            GETPOT_INPUT(write_gmv_error);
        #ifndef LIBMESH_HAVE_GMV
            output_gmv                  = false;
            write_gmv_error             = false;
        #endif
            GETPOT_INPUT(output_tecplot);
            GETPOT_INPUT(write_tecplot_error);
        #ifndef LIBMESH_HAVE_TECPLOT_API
            output_tecplot              = false;
            write_tecplot_error         = false;
        #endif
        
        
            GETPOT_REGISTER(system_types);
            const unsigned int n_system_types =
              input.vector_variable_size("system_types");
            if (n_system_types)
              {
                system_types.resize(n_system_types, "");
                for (unsigned int i=0; i != n_system_types; ++i)
                  {
                    system_types[i] = input("system_types", system_types[i], i);
                  }
              }
        
        
            GETPOT_REGISTER(periodic_boundaries);
            const unsigned int n_periodic_bcs =
              input.vector_variable_size("periodic_boundaries");
        
            if (n_periodic_bcs)
              {
                if (domaintype != "square" &&
                    domaintype != "cylinder" &&
                    domaintype != "file" &&
                    domaintype != "od2")
                  {
                    libMesh::out << "Periodic boundaries need rectilinear domains" << std::endl;;
                    libmesh_error();
                  }
                for (unsigned int i=0; i != n_periodic_bcs; ++i)
                  {
                    unsigned int myboundary =
                      input("periodic_boundaries", -1, i);
unsigned int pairedboundary = 0;
                    RealVectorValue translation_vector;
                    if (dimension == 2)
                      switch (myboundary)
                      {
                      case 0:
pairedboundary = 2;
                        translation_vector = RealVectorValue(0., domain_edge_length);
                        break;
                      case 1:
pairedboundary = 3;
                        translation_vector = RealVectorValue(-domain_edge_width, 0);
                        break;
                      default:
                        libMesh::out << "Unrecognized periodic boundary id " <<
                                        myboundary << std::endl;;
                        libmesh_error();
                      }
                    else if (dimension == 3)
                      switch (myboundary)
                      {
                      case 0:
pairedboundary = 5;
                        translation_vector = RealVectorValue((Real) 0., (Real) 0., domain_edge_height);
                        break;
                      case 1:
pairedboundary = 3;
                        translation_vector = RealVectorValue((Real) 0., domain_edge_length, (Real) 0.);
                        break;
                      case 2:
pairedboundary = 4;
                        translation_vector = RealVectorValue(-domain_edge_width, (Real) 0., (Real) 0.);
                        break;
                      default:
                        libMesh::out << "Unrecognized periodic boundary id " <<
                                        myboundary << std::endl;;
                        libmesh_error();
                      }
periodic_boundaries.push_back(PeriodicBoundary(translation_vector)); periodic_boundaries[i].myboundary = myboundary; periodic_boundaries[i].pairedboundary = pairedboundary;
                  }
              }
        
Use std::string inputs so GetPot doesn't have to make a bunch of internal C string copies
            std::string zero_string = "zero";
            std::string empty_string = "";
        
            GETPOT_REGISTER(dirichlet_condition_types);
            GETPOT_REGISTER(dirichlet_condition_values);
            GETPOT_REGISTER(dirichlet_condition_boundaries);
            GETPOT_REGISTER(dirichlet_condition_variables);
        
            const unsigned int n_dirichlet_conditions=
              input.vector_variable_size("dirichlet_condition_types");
        
            if (n_dirichlet_conditions !=
                input.vector_variable_size("dirichlet_condition_values"))
              {
                libMesh::out << "Error: " << n_dirichlet_conditions
                             << " Dirichlet condition types does not match "
                             << input.vector_variable_size("dirichlet_condition_values")
                             << " Dirichlet condition values." << std::endl;
        
                libmesh_error();
              }
        
            if (n_dirichlet_conditions !=
                input.vector_variable_size("dirichlet_condition_boundaries"))
              {
                libMesh::out << "Error: " << n_dirichlet_conditions
                             << " Dirichlet condition types does not match "
                             << input.vector_variable_size("dirichlet_condition_boundaries")
                             << " Dirichlet condition boundaries." << std::endl;
        
                libmesh_error();
              }
        
            if (n_dirichlet_conditions !=
                input.vector_variable_size("dirichlet_condition_variables"))
              {
                libMesh::out << "Error: " << n_dirichlet_conditions
                             << " Dirichlet condition types does not match "
                             << input.vector_variable_size("dirichlet_condition_variables")
                             << " Dirichlet condition variables sets." << std::endl;
        
                libmesh_error();
              }
        
            for (unsigned int i=0; i != n_dirichlet_conditions; ++i)
              {
                const std::string func_type =
                  input("dirichlet_condition_types", zero_string, i);
        
                const std::string func_value =
                  input("dirichlet_condition_values", empty_string, i);
        
                const boundary_id_type func_boundary =
                  input("dirichlet_condition_boundaries", boundary_id_type(0), i);
        
                dirichlet_conditions[func_boundary] =
                  (new_function_base(func_type, func_value).release());
        
                const std::string variable_set =
                  input("dirichlet_condition_variables", empty_string, i);
        
                for (unsigned int j=0; j != variable_set.size(); ++j)
                  {
                    if (variable_set[j] == '1')
                      dirichlet_condition_variables[func_boundary].push_back(j);
                    else if (variable_set[j] != '0')
                      {
                        libMesh::out << "Unable to understand Dirichlet variable set"
                                     << variable_set << std::endl;
                        libmesh_error();
                      }
                  }
              }
        
            GETPOT_REGISTER(neumann_condition_types);
            GETPOT_REGISTER(neumann_condition_values);
            GETPOT_REGISTER(neumann_condition_boundaries);
            GETPOT_REGISTER(neumann_condition_variables);
        
            const unsigned int n_neumann_conditions=
              input.vector_variable_size("neumann_condition_types");
        
            if (n_neumann_conditions !=
                input.vector_variable_size("neumann_condition_values"))
              {
                libMesh::out << "Error: " << n_neumann_conditions
                             << " Neumann condition types does not match "
                             << input.vector_variable_size("neumann_condition_values")
                             << " Neumann condition values." << std::endl;
        
                libmesh_error();
              }
        
            if (n_neumann_conditions !=
                input.vector_variable_size("neumann_condition_boundaries"))
              {
                libMesh::out << "Error: " << n_neumann_conditions
                             << " Neumann condition types does not match "
                             << input.vector_variable_size("neumann_condition_boundaries")
                             << " Neumann condition boundaries." << std::endl;
        
                libmesh_error();
              }
        
            if (n_neumann_conditions !=
                input.vector_variable_size("neumann_condition_variables"))
              {
                libMesh::out << "Error: " << n_neumann_conditions
                             << " Neumann condition types does not match "
                             << input.vector_variable_size("neumann_condition_variables")
                             << " Neumann condition variables sets." << std::endl;
        
                libmesh_error();
              }
        
            for (unsigned int i=0; i != n_neumann_conditions; ++i)
              {
                const std::string func_type =
                  input("neumann_condition_types", zero_string, i);
        
                const std::string func_value =
                  input("neumann_condition_values", empty_string, i);
        
                const boundary_id_type func_boundary =
                  input("neumann_condition_boundaries", boundary_id_type(0), i);
        
                neumann_conditions[func_boundary] =
                  (new_function_base(func_type, func_value).release());
        
                const std::string variable_set =
                  input("neumann_condition_variables", empty_string, i);
        
                for (unsigned int j=0; j != variable_set.size(); ++j)
                  {
                    if (variable_set[j] == '1')
                      neumann_condition_variables[func_boundary].push_back(j);
                    else if (variable_set[j] != '0')
                      {
                        libMesh::out << "Unable to understand Neumann variable set"
                                     << variable_set << std::endl;
                        libmesh_error();
                      }
                  }
              }
        
            GETPOT_REGISTER(initial_condition_types);
            GETPOT_REGISTER(initial_condition_values);
            GETPOT_REGISTER(initial_condition_subdomains);
        
            const unsigned int n_initial_conditions=
              input.vector_variable_size("initial_condition_types");
        
            if (n_initial_conditions !=
                input.vector_variable_size("initial_condition_values"))
              {
                libMesh::out << "Error: " << n_initial_conditions
                             << " initial condition types does not match "
                             << input.vector_variable_size("initial_condition_values")
                             << " initial condition values." << std::endl;
        
                libmesh_error();
              }
        
            if (n_initial_conditions !=
                input.vector_variable_size("initial_condition_subdomains"))
              {
                libMesh::out << "Error: " << n_initial_conditions
                             << " initial condition types does not match "
                             << input.vector_variable_size("initial_condition_subdomains")
                             << " initial condition subdomains." << std::endl;
        
                libmesh_error();
              }
        
            for (unsigned int i=0; i != n_initial_conditions; ++i)
              {
                const std::string func_type =
                  input("initial_condition_types", zero_string, i);
        
                const std::string func_value =
                  input("initial_condition_values", empty_string, i);
        
                const subdomain_id_type func_subdomain =
                  input("initial_condition_subdomains", subdomain_id_type(0), i);
        
                initial_conditions[func_subdomain] =
                  (new_function_base(func_type, func_value).release());
              }
        
            GETPOT_REGISTER(other_interior_function_types);
            GETPOT_REGISTER(other_interior_function_values);
            GETPOT_REGISTER(other_interior_function_subdomains);
            GETPOT_REGISTER(other_interior_function_ids);
        
            const unsigned int n_other_interior_functions =
              input.vector_variable_size("other_interior_function_types");
        
            if (n_other_interior_functions !=
                input.vector_variable_size("other_interior_function_values"))
              {
                libMesh::out << "Error: " << n_other_interior_functions
                             << " other interior function types does not match "
                             << input.vector_variable_size("other_interior_function_values")
                             << " other interior function values." << std::endl;
        
                libmesh_error();
              }
        
            if (n_other_interior_functions !=
                input.vector_variable_size("other_interior_function_subdomains"))
              {
                libMesh::out << "Error: " << n_other_interior_functions
                             << " other interior function types does not match "
                             << input.vector_variable_size("other_interior_function_subdomains")
                             << " other interior function subdomains." << std::endl;
        
                libmesh_error();
              }
        
            if (n_other_interior_functions !=
                input.vector_variable_size("other_interior_function_ids"))
              {
                libMesh::out << "Error: " << n_other_interior_functions
                             << " other interior function types does not match "
                             << input.vector_variable_size("other_interior_function_ids")
                             << " other interior function ids." << std::endl;
        
                libmesh_error();
              }
        
            for (unsigned int i=0; i != n_other_interior_functions; ++i)
              {
                const std::string func_type =
                  input("other_interior_function_types", zero_string, i);
        
                const std::string func_value =
                  input("other_interior_function_values", empty_string, i);
        
                const subdomain_id_type func_subdomain =
                  input("other_interior_condition_subdomains", subdomain_id_type(0), i);
        
                const int func_id =
                  input("other_interior_condition_ids", int(0), i);
        
                other_interior_functions[func_id][func_subdomain] =
                  (new_function_base(func_type, func_value).release());
              }
        
            GETPOT_REGISTER(other_boundary_function_types);
            GETPOT_REGISTER(other_boundary_function_values);
            GETPOT_REGISTER(other_boundary_function_boundaries);
            GETPOT_REGISTER(other_boundary_function_ids);
        
            const unsigned int n_other_boundary_functions =
              input.vector_variable_size("other_boundary_function_types");
        
            if (n_other_boundary_functions !=
                input.vector_variable_size("other_boundary_function_values"))
              {
                libMesh::out << "Error: " << n_other_boundary_functions
                             << " other boundary function types does not match "
                             << input.vector_variable_size("other_boundary_function_values")
                             << " other boundary function values." << std::endl;
        
                libmesh_error();
              }
        
            if (n_other_boundary_functions !=
                input.vector_variable_size("other_boundary_function_boundaries"))
              {
                libMesh::out << "Error: " << n_other_boundary_functions
                             << " other boundary function types does not match "
                             << input.vector_variable_size("other_boundary_function_boundaries")
                             << " other boundary function boundaries." << std::endl;
        
                libmesh_error();
              }
        
            if (n_other_boundary_functions !=
                input.vector_variable_size("other_boundary_function_ids"))
              {
                libMesh::out << "Error: " << n_other_boundary_functions
                             << " other boundary function types does not match "
                             << input.vector_variable_size("other_boundary_function_ids")
                             << " other boundary function ids." << std::endl;
        
                libmesh_error();
              }
        
            for (unsigned int i=0; i != n_other_boundary_functions; ++i)
              {
                const std::string func_type =
                  input("other_boundary_function_types", zero_string, i);
        
                const std::string func_value =
                  input("other_boundary_function_values", empty_string, i);
        
                const boundary_id_type func_boundary =
                  input("other_boundary_function_boundaries", boundary_id_type(0), i);
        
                const int func_id =
                  input("other_boundary_function_ids", int(0), i);
        
                other_boundary_functions[func_id][func_boundary] =
                  (new_function_base(func_type, func_value).release());
              }
        
            GETPOT_INPUT(run_simulation);
            GETPOT_INPUT(run_postprocess);
        
        
            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(extra_quadrature_order);
        
        
            GETPOT_INPUT(analytic_jacobians);
            GETPOT_INPUT(verify_analytic_jacobians);
            GETPOT_INPUT(numerical_jacobian_h);
            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);
            GETPOT_INPUT(print_element_jacobians);
        
        
            GETPOT_INPUT(use_petsc_snes);
            GETPOT_INPUT(time_solver_quiet);
            GETPOT_INPUT(solver_quiet);
            GETPOT_INPUT(solver_verbose);
            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(initial_sobolev_order);
            GETPOT_INT_INPUT(initial_extra_quadrature);
            GETPOT_INPUT(indicator_type);
            GETPOT_INT_INPUT(sobolev_order);
        
            GETPOT_INPUT(system_config_file);
        
          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 heatsystem.C with comments:

        #include "libmesh/getpot.h"
        
        #include "heatsystem.h"
        
        #include "libmesh/fe_base.h"
        #include "libmesh/fe_interface.h"
        #include "libmesh/fem_context.h"
        #include "libmesh/mesh.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/system.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/zero_function.h"
        #include "libmesh/boundary_info.h"
        #include "libmesh/dirichlet_boundaries.h"
        #include "libmesh/dof_map.h"
        
        void HeatSystem::init_data ()
        {
          T_var = this->add_variable ("T", static_cast<Order>(_fe_order),
                              Utility::string_to_enum<FEFamily>(_fe_family));
        
Make sure the input file heat.in exists, and parse it.
          {
            std::ifstream i("heat.in");
            if (!i)
              {
                std::cerr << '[' << libMesh::processor_id()
                          << "] Can't find heat.in; exiting early."
                          << std::endl;
                libmesh_error();
              }
          }
          GetPot infile("heat.in");
          _k = infile("k", 1.0);
          _analytic_jacobians = infile("analytic_jacobians", true);
        
          parameters.push_back(_k);
        
Set equation system parameters _k and theta, so they can be read by the exact solution
          this->get_equation_systems().parameters.set<Real>("_k") = _k;
        
          this->time_evolving(T_var);
        
          const boundary_id_type all_ids[4] = {0, 1, 2, 3};
          std::set<boundary_id_type> all_bdys(all_ids, all_ids+(2*2));
        
          std::vector<unsigned int> T_only(1, T_var);
        
          ZeroFunction<Number> zero;
        
          this->get_dof_map().add_dirichlet_boundary
                (DirichletBoundary (all_bdys, T_only, &zero));
        
          FEMSystem::init_data();
        }
        
        
        
        void HeatSystem::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_dphi();
        
We'll have a more automatic solution to preparing adjoint solutions for time integration, eventually...
          if (c.is_adjoint())
            {
A reference to the system context is built with
              const System & sys = c.get_system();
        
Get a pointer to the adjoint solution vector
              NumericVector<Number> &adjoint_solution =
                const_cast<System &>(sys).get_adjoint_solution(0);
        
Add this adjoint solution to the vectors that diff context should localize
              c.add_localized_vector(adjoint_solution, sys);
            }
        
          FEMSystem::init_context(context);
        }
        
        #define optassert(X) {if (!(X)) libmesh_error();}
#define optassert(X) libmesh_assert(X);

        bool HeatSystem::element_time_derivative (bool request_jacobian,
                                                  DiffContext &context)
        {
          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
          const std::vector<Real> &JxW = c.element_fe_var[0]->get_JxW();
        
const std::vector > &phi = c.element_fe_var[0]->get_phi();
          const std::vector<std::vector<RealGradient> > &dphi = c.element_fe_var[0]->get_dphi();
        
Workaround for weird FC6 bug
          optassert(c.dof_indices_var.size() > 0);
        
The number of local degrees of freedom in each variable
          const unsigned int n_u_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.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);
        
              for (unsigned int i=0; i != n_u_dofs; i++)
                F(i) += JxW[qp] * -parameters[0] * (grad_T * dphi[i][qp]);
              if (compute_jacobian)
                for (unsigned int i=0; i != n_u_dofs; i++)
                  for (unsigned int j=0; j != n_u_dofs; ++j)
                    K(i,j) += JxW[qp] * -parameters[0] * (dphi[i][qp] * dphi[j][qp]);
            } // end of the quadrature point qp-loop
        
          return compute_jacobian;
        }
        
Perturb and accumulate dual weighted residuals
        void HeatSystem::perturb_accumulate_residuals(const ParameterVector& parameters_in)
        {
          const unsigned int Np = parameters_in.size();
        
          for (unsigned int j=0; j != Np; ++j)
            {
              Number old_parameter = *parameters_in[j];
        
              *parameters_in[j] = old_parameter - dp;
        
              this->assembly(true, false);
        
              this->rhs->close();
        
              AutoPtr<NumericVector<Number> > R_minus = this->rhs->clone();
        
The contribution at a single time step would be [f(z;p+dp) - (p+dp) - (p+dp)] * dt But since we compute the residual already scaled by dt, there is no need for the * dt
              R_minus_dp += -R_minus->dot(this->get_adjoint_solution(0));
        
              *parameters_in[j] = old_parameter + dp;
        
              this->assembly(true, false);
        
              this->rhs->close();
        
              AutoPtr<NumericVector<Number> > R_plus = this->rhs->clone();
        
              R_plus_dp += -R_plus->dot(this->get_adjoint_solution(0));
        
              *parameters_in[j] = old_parameter;
        
            }
        }



The source file initial.C with comments:

        #include "initial.h"
        
        using namespace libMesh;
        
        void read_initial_parameters()
        {
        }
        
        void finish_initialization()
        {
        }
        
        
        
Initial conditions
        Number initial_value(const Point& p,
                             const Parameters&,
                             const std::string&,
                             const std::string&)
        {
          Real x = p(0), y = p(1);
        
          return sin(M_PI * x) * sin(M_PI * y);
        }
        
        
        
        Gradient initial_grad(const Point& p,
                              const Parameters&,
                              const std::string&,
                              const std::string&)
        {
          Real x = p(0), y = p(1);
        
          return Gradient(M_PI*cos(M_PI * x) * sin(M_PI * y),
          M_PI*sin(M_PI * x) * cos(M_PI * y));
        }



The source file adjoint_initial.h without comments:

 
  #include "libmesh/parameters.h"
  #include "libmesh/point.h"
  #include "libmesh/vector_value.h"
  
  using namespace libMesh;
  
  void adjoint_read_initial_parameters();
  void adjoint_finish_initialization();
  
  Number adjoint_initial_value(const Point& p,
                       const Parameters&,
                       const std::string&,
                       const std::string&);
  
  Gradient adjoint_initial_grad(const Point& p,
                        const Parameters&,
                        const std::string&,
                        const std::string&);



The source file femparameters.h without comments:

 
  #ifndef __fem_parameters_h__
  #define __fem_parameters_h__
  
  #include "libmesh/libmesh_common.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/enum_norm_type.h"
  #include "libmesh/function_base.h"
  #include "libmesh/getpot.h"
  #include "libmesh/id_types.h"
  #include "libmesh/periodic_boundaries.h"
  
  #include <limits>
  #include <map>
  #include <string>
  #include <vector>
  
  using namespace libMesh;
  
  class FEMParameters
  {
  public:
      FEMParameters();
  
      ~FEMParameters();
  
      void read(GetPot &input);
  
  
      unsigned int initial_timestep, n_timesteps;
      bool transient;
      unsigned int deltat_reductions;
      std::string timesolver_core;
      Real end_time, deltat, timesolver_theta,
  	 timesolver_maxgrowth, timesolver_tolerance,
  	 timesolver_upper_tolerance, steadystate_tolerance;
      std::vector<FEMNormType> timesolver_norm;
  
  
      unsigned int dimension;
      std::string domaintype, domainfile, elementtype;
      Real elementorder;
      Real domain_xmin, domain_ymin, domain_zmin;
      Real domain_edge_width, domain_edge_length, domain_edge_height;
      unsigned int coarsegridx, coarsegridy, coarsegridz;
      unsigned int coarserefinements, extrarefinements;
  
  
      unsigned int nelem_target;
      Real global_tolerance;
      Real refine_fraction, coarsen_fraction, coarsen_threshold;
      unsigned int max_adaptivesteps;
      unsigned int initial_adaptivesteps;
  
  
      unsigned int write_interval;
      bool write_gmv_error, write_tecplot_error, output_xda, output_xdr,
  	 output_bz2, output_gz, output_gmv, output_tecplot;
  
  
      std::vector<std::string> system_types;
  
  
  
  
      std::map<subdomain_id_type, FunctionBase<Number> *> initial_conditions;
      std::map<boundary_id_type, FunctionBase<Number> *>  dirichlet_conditions,
                                                          neumann_conditions;
      std::map<boundary_id_type, std::vector<unsigned int> > dirichlet_condition_variables,
                                                             neumann_condition_variables;
      std::map<int, std::map<subdomain_id_type, FunctionBase<Number> *> >
        other_interior_functions;
      std::map<int, std::map<boundary_id_type, FunctionBase<Number> *> >
        other_boundary_functions;
  
  
      bool run_simulation, run_postprocess;
  
  
      std::vector<std::string> fe_family;
      std::vector<unsigned int> fe_order;
      int extra_quadrature_order;
  
  
      bool analytic_jacobians;
      Real verify_analytic_jacobians;
      Real numerical_jacobian_h;
  
      bool print_solution_norms, print_solutions,
           print_residual_norms, print_residuals,
           print_jacobian_norms, print_jacobians,
           print_element_jacobians;
  
  
      bool use_petsc_snes;
      bool time_solver_quiet, solver_quiet, solver_verbose, 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 initial_sobolev_order;
      unsigned int initial_extra_quadrature;
  
  
      std::string indicator_type;
      unsigned int sobolev_order;
  
  
      std::string system_config_file;
  };
  
  #endif // __fem_parameters_h__



The source file heatsystem.h without comments:

 
  #include "libmesh/enum_fe_family.h"
  #include "libmesh/fem_system.h"
  #include "libmesh/parameter_vector.h"
  
  class HeatSystem : public FEMSystem
  {
  public:
    HeatSystem(EquationSystems& es,
                 const std::string& name_in,
                 const unsigned int number_in)
    : FEMSystem(es, name_in, number_in),
      _k(1.0),
      _fe_family("LAGRANGE"), _fe_order(1),
      _analytic_jacobians(true), R_plus_dp(0.0), R_minus_dp(0.0), dp(1.e-6) { qoi.resize(1); }
  
    std::string & fe_family() { return _fe_family;  }
    unsigned int & fe_order() { return _fe_order;  }
    Real & k() { return _k;  }
    bool & analytic_jacobians() { return _analytic_jacobians; }
  
    void perturb_accumulate_residuals(const ParameterVector& parameters);
  
    Number & compute_final_sensitivity()
      {
        final_sensitivity = -(R_plus_dp - R_minus_dp)/(2*dp);
  
        return final_sensitivity;
      }
  
      void set_tf(Real val)
    {
      tf = val;
    }
  
      ParameterVector &get_parameter_vector()
      {
        parameter_vector.resize(parameters.size());
        for(unsigned int i = 0; i != parameters.size(); ++i)
  	{
  	  parameter_vector[i] = &parameters[i];
  	}
  
        return parameter_vector;
      }
  
    Number &get_QoI_value(unsigned int QoI_index)
      {
        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 void element_qoi_derivative (DiffContext &context,
  				       const QoISet & /* qois */);
  
  
    std::vector<Number> parameters;
  
    ParameterVector parameter_vector;
  
    Real _k;
  
    Real tf;
  
    Number computed_QoI[1];
  
    std::string _fe_family;
    unsigned int _fe_order;
  
    unsigned int T_var;
  
    bool _analytic_jacobians;
  
    Number R_plus_dp;
    Number R_minus_dp;
  
    Real dp;
  
    Number final_sensitivity;
  
  };



The source file initial.h without comments:

 
  #include "libmesh/parameters.h"
  #include "libmesh/point.h"
  #include "libmesh/vector_value.h"
  
  using namespace libMesh;
  
  void read_initial_parameters();
  void finish_initialization();
  
  Number initial_value(const Point& p,
                       const Parameters&,
                       const std::string&,
                       const std::string&);
  
  Gradient initial_grad(const Point& p,
                        const Parameters&,
                        const std::string&,
                        const std::string&);



The source file adjoint_initial.C without comments:

 
  #include "adjoint_initial.h"
  
  using namespace libMesh;
  
  void adjoint_read_initial_parameters()
  {
  }
  
  void adjoint_finish_initialization()
  {
  }
  
  
  
  Number adjoint_initial_value(const Point& p,
                       const Parameters&,
                       const std::string&,
                       const std::string&)
  {
    Real x = p(0), y = p(1);
  
    Number val = 0.;
  
    val = sin(M_PI * x) * sin(M_PI * y);
  
    return val;
  }
  
  
  
  Gradient adjoint_initial_grad(const Point& p,
                        const Parameters&,
                        const std::string&,
                        const std::string&)
  {
    Real x = p(0), y = p(1);
  
    return Gradient(M_PI*cos(M_PI * x) * sin(M_PI * y),
    M_PI*sin(M_PI * x) * cos(M_PI * y));
  }



The source file adjoints_ex5.C without comments:

 
  
  
  
  
  
  
  
  
  #include "initial.h"
  #include "adjoint_initial.h"
  #include "femparameters.h"
  #include "heatsystem.h"
  
  #include "libmesh/equation_systems.h"
  #include "libmesh/dirichlet_boundaries.h"
  #include "libmesh/system_norm.h"
  #include "libmesh/numeric_vector.h"
  
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/mesh_refinement.h"
  
  #include "libmesh/petsc_diff_solver.h"
  #include "libmesh/steady_solver.h"
  #include "libmesh/euler_solver.h"
  #include "libmesh/euler2_solver.h"
  #include "libmesh/newton_solver.h"
  #include "libmesh/twostep_time_solver.h"
  
  #include "libmesh/getpot.h"
  #include "libmesh/tecplot_io.h"
  #include "libmesh/gmv_io.h"
  
  #include "libmesh/solution_history.h"
  #include "libmesh/memory_solution_history.h"
  
  #include <iostream>
  #include <sys/time.h>
  #include <iomanip>
  
  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(HeatSystem &system, FEMParameters &param)
  {
    system.fe_family() = param.fe_family[0];
    system.fe_order() = param.fe_order[0];
  
    system.analytic_jacobians() = param.analytic_jacobians;
  
    system.verify_analytic_jacobians = param.verify_analytic_jacobians;
    system.numerical_jacobian_h = param.numerical_jacobian_h;
  
    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.print_element_jacobians = param.print_element_jacobians;
  
    if (param.transient)
      {
        UnsteadySolver *innersolver;
         if (param.timesolver_core == "euler")
          {
            EulerSolver *eulersolver =
              new EulerSolver(system);
  
            eulersolver->theta = param.timesolver_theta;
            innersolver = eulersolver;
          }
        else
          {
            std::cerr << "This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods."
                      <<  std::endl;
            libmesh_error();
          }
         system.time_solver =
            AutoPtr<TimeSolver>(innersolver);
      }
    else
      system.time_solver =
        AutoPtr<TimeSolver>(new SteadySolver(system));
  
    MemorySolutionHistory heatsystem_solution_history(system);
    system.time_solver->set_solution_history(heatsystem_solution_history);
  
    system.time_solver->reduce_deltat_on_diffsolver_failure =
                                          param.deltat_reductions;
    system.time_solver->quiet           = param.time_solver_quiet;
  
    typedef
      std::map<boundary_id_type, FunctionBase<Number> *>::
      const_iterator Iter;
  
    for (Iter i = param.dirichlet_conditions.begin();
         i != param.dirichlet_conditions.end(); ++i)
      {
        boundary_id_type b = i->first;
        FunctionBase<Number> *f = i->second;
        std::set<boundary_id_type> bdys; bdys.insert(b);
        system.get_dof_map().add_dirichlet_boundary
          (DirichletBoundary
             (bdys, param.dirichlet_condition_variables[b], f));
        std::cout << "Added Dirichlet boundary " << b << " for variables ";
        for (unsigned int vi=0; vi !=
             param.dirichlet_condition_variables[b].size(); ++vi)
          std::cout << param.dirichlet_condition_variables[b][vi];
        std::cout << std::endl;
      }
  
    system.deltat = param.deltat;
  
    system.extra_quadrature_order = param.extra_quadrature_order;
  
    if (param.use_petsc_snes)
      {
  #ifdef LIBMESH_HAVE_PETSC
        PetscDiffSolver *solver = new PetscDiffSolver(system);
        system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
  #else
        libmesh_error();
  #endif
      }
    else
      {
        NewtonSolver *solver = new NewtonSolver(system);
        system.time_solver->diff_solver() = AutoPtr<DiffSolver>(solver);
  
        solver->quiet                       = param.solver_quiet;
        solver->verbose                     = param.solver_verbose;
        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;
      }
  }
  
  int main (int argc, char** argv)
  {
  #ifndef LIBMESH_ENABLE_AMR
    libmesh_example_assert(false, "--enable-amr");
  #else
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    LibMeshInit init (argc, argv);
  
    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);
  
    Mesh mesh(init.comm(), param.dimension);
  
    AutoPtr<MeshRefinement> mesh_refinement(new MeshRefinement(mesh));
  
    EquationSystems equation_systems (mesh);
  
    std::cout << "Building mesh" << std::endl;
  
    ElemType elemtype;
  
    if (param.elementtype == "tri" ||
        param.elementtype == "unstructured")
      elemtype = TRI3;
    else
      elemtype = QUAD4;
  
    MeshTools::Generation::build_square
      (mesh, param.coarsegridx, param.coarsegridy,
       param.domain_xmin, param.domain_xmin + param.domain_edge_width,
       param.domain_ymin, param.domain_ymin + param.domain_edge_length,
       elemtype);
  
    std::cout << "Building system" << std::endl;
  
    HeatSystem &system = equation_systems.add_system<HeatSystem> ("HeatSystem");
  
    set_system_parameters(system, param);
  
    std::cout << "Initializing systems" << std::endl;
  
    equation_systems.init ();
  
    for (unsigned int i=0; i != param.extrarefinements; ++i)
      {
        mesh_refinement->uniformly_refine(1);
        equation_systems.reinit();
      }
  
    std::cout<<"Setting primal initial conditions"<<std::endl;
  
    read_initial_parameters();
  
    system.project_solution(initial_value, initial_grad,
                                    equation_systems.parameters);
  
    libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl<<std::endl;
  
    const std::string & adjoint_solution_name = "adjoint_solution0";
    system.add_vector("adjoint_solution0", false, GHOSTED);
  
    finish_initialization();
  
    write_output(equation_systems, 0, "primal");
  
    mesh.print_info();
    equation_systems.print_info();
  
  #ifdef NDEBUG
    try
    {
  #endif
      for (unsigned int t_step=param.initial_timestep;
           t_step != param.initial_timestep + param.n_timesteps; ++t_step)
        {
  	std::cout << " Solving time step " << t_step << ", time = "
                    << system.time << std::endl;
  
  	system.solve();
  
  	libMesh::out << "|U(" <<system.time + system.deltat<< ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
  
  	std::cout<<"Advancing timestep"<<std::endl<<std::endl;
          system.time_solver->advance_timestep();
  
            write_output(equation_systems, t_step+1, "primal");
        }
  
  
    std::cout << std::endl << "Solving the adjoint problem" << std::endl;
  
  
    system.set_vector_preservation(adjoint_solution_name, true);
  
    std::cout << "Setting adjoint initial conditions Z("<<system.time<<")"<<std::endl;
  
    std::cout<<"Retrieving solutions at time t="<<system.time<<std::endl;
    system.time_solver->adjoint_advance_timestep();
  
    libMesh::out << "|U(" <<system.time + system.deltat<< ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
  
    libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
  
    system.project_vector(adjoint_initial_value, adjoint_initial_grad, equation_systems.parameters, system.get_adjoint_solution(0));
  
    system.set_adjoint_already_solved(true);
  
    libMesh::out << "|Z(" <<system.time<< ")|= " << system.calculate_norm(system.get_adjoint_solution(), 0, H1) << std::endl<<std::endl;
  
    write_output(equation_systems, param.n_timesteps, "dual");
  
  
    for (unsigned int t_step=param.initial_timestep;
         t_step != param.initial_timestep + param.n_timesteps; ++t_step)
      {
        std::cout << " Solving adjoint time step " << t_step << ", time = "
  		<< system.time << std::endl;
  
        std::cout<<"Retrieving solutions at time t="<<system.time<<std::endl;
        system.time_solver->adjoint_advance_timestep();
  
        libMesh::out << "|U(" <<system.time + system.deltat << ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
  
        libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
  
        system.set_adjoint_already_solved(false);
  
        system.adjoint_solve();
  
        system.set_adjoint_already_solved(true);
  
        libMesh::out << "|Z(" <<system.time<< ")|= "<< system.calculate_norm(system.get_adjoint_solution(), 0, H1) << std::endl << std::endl;
  
        NumericVector<Number> &primal_solution = *system.solution;
  
        NumericVector<Number> &dual_solution_0 = system.get_adjoint_solution(0);
  
        primal_solution.swap(dual_solution_0);
  
        write_output(equation_systems, param.n_timesteps - (t_step + 1), "dual");
  
        primal_solution.swap(dual_solution_0);
      }
  
  
      for (unsigned int t_step=param.initial_timestep;
           t_step != param.initial_timestep + param.n_timesteps; ++t_step)
        {
  	std::cout << "Retrieving " << t_step << ", time = "
                    << system.time << std::endl;
  
  	system.time_solver->retrieve_timestep();
  
  	libMesh::out << "|U(" <<system.time + system.deltat << ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
  
        libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
  
        libMesh::out << "|Z(" <<system.time<< ")|= "<< system.calculate_norm(system.get_adjoint_solution(0), 0, H1) << std::endl << std::endl;
  
  	(dynamic_cast<HeatSystem&>(system)).perturb_accumulate_residuals(dynamic_cast<HeatSystem&>(system).get_parameter_vector());
  
  	system.time += system.deltat;
        }
  
      std::cout << "Retrieving " << " final time = "
  	      << system.time << std::endl;
  
      system.time_solver->retrieve_timestep();
  
      libMesh::out << "|U(" <<system.time + system.deltat << ")|= " << system.calculate_norm(*system.solution, 0, H1) << std::endl;
  
        libMesh::out << "|U(" <<system.time<< ")|= " << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1) << std::endl;
  
        libMesh::out << "|Z(" <<system.time<< ")|= "<< system.calculate_norm(system.get_adjoint_solution(0), 0, H1) << std::endl<<std::endl;
  
      (dynamic_cast<HeatSystem&>(system)).perturb_accumulate_residuals(dynamic_cast<HeatSystem&>(system).get_parameter_vector());
  
      Number sensitivity_0_0 = (dynamic_cast<HeatSystem&>(system)).compute_final_sensitivity();
  
      std::cout<<"Sensitivity of QoI 0 w.r.t parameter 0 is: " << sensitivity_0_0 << std::endl;
  
  #ifdef NDEBUG
  }
    catch (...)
    {
      std::cerr << '[' << libMesh::processor_id()
                << "] Caught exception; exiting early." << std::endl;
    }
  #endif
  
    std::cerr << '[' << libMesh::processor_id()
              << "] Completing output." << std::endl;
  
    return 0;
  
  #endif // LIBMESH_ENABLE_AMR
  }



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 "heatsystem.h"
  
  using namespace libMesh;
  
  void HeatSystem::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 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];
  
    const System & sys = c.get_system();
  
    NumericVector<Number> &adjoint_solution = const_cast<System &>(sys).get_adjoint_solution(0);
  
  
    std::vector<Number> old_adjoint (n_qpoints, 0);
  
    c.interior_values<Number>(0, adjoint_solution, old_adjoint);
  
    for (unsigned int qp=0; qp != n_qpoints; qp++)
      {
        for (unsigned int i=0; i != n_T_dofs; i++)
  	{
  	  Q(i) += -JxW[qp] * old_adjoint[qp] * phi[i][qp] ;
        	}
  
      } // end of the quadrature point qp-loop
  }



The source file factoryfunction.C without comments:

 
  #include "libmesh/dense_vector.h"
  #include "libmesh/factory.h"
  #include "libmesh/function_base.h"
  
  using namespace libMesh;
  
  
  class ExampleOneFunction : public FunctionBase<Number>
  {
    virtual Number operator() (const Point&  /*p*/,
                               const Real /*time*/)
      {
        return 1;
      }
  
    virtual void operator() (const Point&  /*p*/,
                             const Real /*time*/,
                             DenseVector<Number>& output)
      {
        for (unsigned int i=0; i != output.size(); ++i)
          output(i) = 1;
      }
  
    virtual void init() {}
    virtual void clear() {}
    virtual AutoPtr<FunctionBase<Number> > clone() const {
      return AutoPtr<FunctionBase<Number> >
        (new ExampleOneFunction());
    }
  };
  
  #ifdef LIBMESH_USE_SEPARATE_NAMESPACE
  namespace libMesh {
  #endif
  
  template<>
  std::map<std::string, Factory<FunctionBase<Number> >*>&
  Factory<FunctionBase<Number> >::factory_map()
  {
    static std::map<std::string, Factory<FunctionBase<Number> >*> _map;
    return _map;
  }
  
  FactoryImp<ExampleOneFunction, FunctionBase<Number> > example_one_factory ("example_one");
  
  #ifdef LIBMESH_USE_SEPARATE_NAMESPACE
  } // namespace libMesh
  #endif



The source file femparameters.C without comments:

 
  #include "femparameters.h"
  
  #include "libmesh/factory.h"
  #include "libmesh/parsed_function.h"
  #include "libmesh/zero_function.h"
  #include "libmesh/periodic_boundaries.h"
  #include "libmesh/vector_value.h" // RealVectorValue
  
  #define GETPOT_INPUT(A) { A = input(#A, A);\
    const 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);\
    const std::string stringval = input(#A, std::string());\
    variable_names.push_back(std::string(#A "=") + stringval); }
  
  #define GETPOT_FUNCTION_INPUT(A) {\
    const std::string type  = input(#A "_type", "zero");\
    const std::string value = input(#A "_value", "");\
    A = new_function_base(type, value); }
  #define GETPOT_REGISTER(A) { \
    std::string stringval = input(#A, std::string());\
    variable_names.push_back(std::string(#A "=") + stringval); }
  
  FEMParameters::FEMParameters() :
        initial_timestep(0), n_timesteps(100),
        transient(true),
        deltat_reductions(0),
        timesolver_core("euler"),
        end_time(std::numeric_limits<Real>::max()),
        deltat(0.0001), timesolver_theta(0.5),
        timesolver_maxgrowth(0.), timesolver_tolerance(0.),
        timesolver_upper_tolerance(0.),
        steadystate_tolerance(0.),
        timesolver_norm(0,L2),
  
        dimension(2),
        domaintype("square"), domainfile("mesh.xda"), elementtype("quad"),
        elementorder(2),
        domain_xmin(0.0), domain_ymin(0.0), domain_zmin(0.0),
        domain_edge_width(1.0), domain_edge_length(1.0), domain_edge_height(1.0),
        coarsegridx(1), coarsegridy(1), coarsegridz(1),
        coarserefinements(4), extrarefinements(0),
  
        nelem_target(8000), global_tolerance(0.0),
        refine_fraction(0.3), coarsen_fraction(0.3), coarsen_threshold(10),
        max_adaptivesteps(1),
        initial_adaptivesteps(0),
  
        write_interval(10),
        write_gmv_error(false), write_tecplot_error(false),
        output_xda(false), output_xdr(false),
        output_bz2(true), output_gz(true),
        output_gmv(false), output_tecplot(false),
  
        system_types(0),
  
  
        run_simulation(true), run_postprocess(false),
  
        fe_family(1, "LAGRANGE"), fe_order(1, 1),
        extra_quadrature_order(0),
  
        analytic_jacobians(true), verify_analytic_jacobians(0.0),
        numerical_jacobian_h(TOLERANCE),
        print_solution_norms(false), print_solutions(false),
        print_residual_norms(false), print_residuals(false),
        print_jacobian_norms(false), print_jacobians(false),
        print_element_jacobians(false),
  
        use_petsc_snes(false),
        time_solver_quiet(true), solver_quiet(true), solver_verbose(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),
  
        initial_sobolev_order(1),
        initial_extra_quadrature(0),
        indicator_type("kelly"), sobolev_order(1)
  {
  }
  
  FEMParameters::~FEMParameters()
  {
    for (std::map<subdomain_id_type, FunctionBase<Number> *>::iterator
         i = initial_conditions.begin(); i != initial_conditions.end();
         ++i)
      delete i->second;
  
    for (std::map<boundary_id_type, FunctionBase<Number> *>::iterator
         i = dirichlet_conditions.begin(); i != dirichlet_conditions.end();
         ++i)
      delete i->second;
  
    for (std::map<boundary_id_type, FunctionBase<Number> *>::iterator
         i = neumann_conditions.begin(); i != neumann_conditions.end();
         ++i)
      delete i->second;
  
    for (std::map<int,
  	 std::map<boundary_id_type, FunctionBase<Number> *>
         >::iterator
         i = other_boundary_functions.begin(); i != other_boundary_functions.end();
         ++i)
      for (std::map<boundary_id_type, FunctionBase<Number> *>::iterator
           j = i->second.begin(); j != i->second.end();
           ++j)
        delete j->second;
  
    for (std::map<int,
  	 std::map<subdomain_id_type, FunctionBase<Number> *>
         >::iterator
         i = other_interior_functions.begin(); i != other_interior_functions.end();
         ++i)
      for (std::map<subdomain_id_type, FunctionBase<Number> *>::iterator
           j = i->second.begin(); j != i->second.end();
           ++j)
        delete j->second;
  }
  
  
  AutoPtr<FunctionBase<Number> > new_function_base(const std::string& func_type,
                                          const std::string& func_value)
  {
    if (func_type == "parsed")
      return AutoPtr<FunctionBase<Number> >
        (new ParsedFunction<Number>(func_value));
    else if (func_type == "factory")
      return AutoPtr<FunctionBase<Number> >
        (Factory<FunctionBase<Number> >::build(func_value).release());
    else if (func_type == "zero")
      return AutoPtr<FunctionBase<Number> >
        (new ZeroFunction<Number>);
    else
      libmesh_not_implemented();
  
    return AutoPtr<FunctionBase<Number> >();
  }
  
  
  void FEMParameters::read(GetPot &input)
  {
      std::vector<std::string> variable_names;
  
      GETPOT_INT_INPUT(initial_timestep);
      GETPOT_INT_INPUT(n_timesteps);
      GETPOT_INPUT(transient);
      GETPOT_INT_INPUT(deltat_reductions);
      GETPOT_INPUT(timesolver_core);
      GETPOT_INPUT(end_time);
      GETPOT_INPUT(deltat);
      GETPOT_INPUT(timesolver_theta);
      GETPOT_INPUT(timesolver_maxgrowth);
      GETPOT_INPUT(timesolver_tolerance);
      GETPOT_INPUT(timesolver_upper_tolerance);
      GETPOT_INPUT(steadystate_tolerance);
  
      GETPOT_REGISTER(timesolver_norm);
      const unsigned int n_timesolver_norm = input.vector_variable_size("timesolver_norm");
      timesolver_norm.resize(n_timesolver_norm, L2);
      for (unsigned int i=0; i != n_timesolver_norm; ++i)
        {
          int current_norm = 0; // L2
          if (timesolver_norm[i] == H1)
            current_norm = 1;
          if (timesolver_norm[i] == H2)
            current_norm = 2;
          current_norm = input("timesolver_norm", current_norm, i);
          if (current_norm == 0)
            timesolver_norm[i] = L2;
          else if (current_norm == 1)
            timesolver_norm[i] = H1;
          else if (current_norm == 2)
            timesolver_norm[i] = H2;
          else
            timesolver_norm[i] = DISCRETE_L2;
        }
  
  
      GETPOT_INT_INPUT(dimension);
      GETPOT_INPUT(domaintype);
      GETPOT_INPUT(domainfile);
      GETPOT_INPUT(elementtype);
      GETPOT_INPUT(elementorder);
      GETPOT_INPUT(domain_xmin);
      GETPOT_INPUT(domain_ymin);
      GETPOT_INPUT(domain_zmin);
      GETPOT_INPUT(domain_edge_width);
      GETPOT_INPUT(domain_edge_length);
      GETPOT_INPUT(domain_edge_height);
      GETPOT_INT_INPUT(coarsegridx);
      GETPOT_INT_INPUT(coarsegridy);
      GETPOT_INT_INPUT(coarsegridz);
      GETPOT_INT_INPUT(coarserefinements);
      GETPOT_INT_INPUT(extrarefinements);
  
  
      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_INT_INPUT(initial_adaptivesteps);
  
  
      GETPOT_INT_INPUT(write_interval);
      GETPOT_INPUT(output_xda);
      GETPOT_INPUT(output_xdr);
      GETPOT_INPUT(output_gz);
  #ifndef LIBMESH_HAVE_GZSTREAM
      output_gz                   = false;
  #endif
      GETPOT_INPUT(output_bz2);
  #ifndef LIBMESH_HAVE_BZ2
      output_bz2                  = false;
  #endif
      GETPOT_INPUT(output_gmv);
      GETPOT_INPUT(write_gmv_error);
  #ifndef LIBMESH_HAVE_GMV
      output_gmv                  = false;
      write_gmv_error             = false;
  #endif
      GETPOT_INPUT(output_tecplot);
      GETPOT_INPUT(write_tecplot_error);
  #ifndef LIBMESH_HAVE_TECPLOT_API
      output_tecplot              = false;
      write_tecplot_error         = false;
  #endif
  
  
      GETPOT_REGISTER(system_types);
      const unsigned int n_system_types =
        input.vector_variable_size("system_types");
      if (n_system_types)
        {
          system_types.resize(n_system_types, "");
          for (unsigned int i=0; i != n_system_types; ++i)
            {
              system_types[i] = input("system_types", system_types[i], i);
            }
        }
  
  
      GETPOT_REGISTER(periodic_boundaries);
      const unsigned int n_periodic_bcs =
        input.vector_variable_size("periodic_boundaries");
  
      if (n_periodic_bcs)
        {
          if (domaintype != "square" &&
              domaintype != "cylinder" &&
              domaintype != "file" &&
              domaintype != "od2")
            {
              libMesh::out << "Periodic boundaries need rectilinear domains" << std::endl;;
              libmesh_error();
            }
          for (unsigned int i=0; i != n_periodic_bcs; ++i)
            {
              unsigned int myboundary =
                input("periodic_boundaries", -1, i);
              RealVectorValue translation_vector;
              if (dimension == 2)
                switch (myboundary)
                {
                case 0:
                  translation_vector = RealVectorValue(0., domain_edge_length);
                  break;
                case 1:
                  translation_vector = RealVectorValue(-domain_edge_width, 0);
                  break;
                default:
                  libMesh::out << "Unrecognized periodic boundary id " <<
                                  myboundary << std::endl;;
                  libmesh_error();
                }
              else if (dimension == 3)
                switch (myboundary)
                {
                case 0:
                  translation_vector = RealVectorValue((Real) 0., (Real) 0., domain_edge_height);
                  break;
                case 1:
                  translation_vector = RealVectorValue((Real) 0., domain_edge_length, (Real) 0.);
                  break;
                case 2:
                  translation_vector = RealVectorValue(-domain_edge_width, (Real) 0., (Real) 0.);
                  break;
                default:
                  libMesh::out << "Unrecognized periodic boundary id " <<
                                  myboundary << std::endl;;
                  libmesh_error();
                }
            }
        }
  
      std::string zero_string = "zero";
      std::string empty_string = "";
  
      GETPOT_REGISTER(dirichlet_condition_types);
      GETPOT_REGISTER(dirichlet_condition_values);
      GETPOT_REGISTER(dirichlet_condition_boundaries);
      GETPOT_REGISTER(dirichlet_condition_variables);
  
      const unsigned int n_dirichlet_conditions=
        input.vector_variable_size("dirichlet_condition_types");
  
      if (n_dirichlet_conditions !=
          input.vector_variable_size("dirichlet_condition_values"))
        {
          libMesh::out << "Error: " << n_dirichlet_conditions
                       << " Dirichlet condition types does not match "
                       << input.vector_variable_size("dirichlet_condition_values")
                       << " Dirichlet condition values." << std::endl;
  
          libmesh_error();
        }
  
      if (n_dirichlet_conditions !=
          input.vector_variable_size("dirichlet_condition_boundaries"))
        {
          libMesh::out << "Error: " << n_dirichlet_conditions
                       << " Dirichlet condition types does not match "
                       << input.vector_variable_size("dirichlet_condition_boundaries")
                       << " Dirichlet condition boundaries." << std::endl;
  
          libmesh_error();
        }
  
      if (n_dirichlet_conditions !=
          input.vector_variable_size("dirichlet_condition_variables"))
        {
          libMesh::out << "Error: " << n_dirichlet_conditions
                       << " Dirichlet condition types does not match "
                       << input.vector_variable_size("dirichlet_condition_variables")
                       << " Dirichlet condition variables sets." << std::endl;
  
          libmesh_error();
        }
  
      for (unsigned int i=0; i != n_dirichlet_conditions; ++i)
        {
          const std::string func_type =
            input("dirichlet_condition_types", zero_string, i);
  
          const std::string func_value =
            input("dirichlet_condition_values", empty_string, i);
  
          const boundary_id_type func_boundary =
            input("dirichlet_condition_boundaries", boundary_id_type(0), i);
  
          dirichlet_conditions[func_boundary] =
            (new_function_base(func_type, func_value).release());
  
          const std::string variable_set =
            input("dirichlet_condition_variables", empty_string, i);
  
          for (unsigned int j=0; j != variable_set.size(); ++j)
            {
              if (variable_set[j] == '1')
                dirichlet_condition_variables[func_boundary].push_back(j);
              else if (variable_set[j] != '0')
                {
                  libMesh::out << "Unable to understand Dirichlet variable set"
                               << variable_set << std::endl;
                  libmesh_error();
                }
            }
        }
  
      GETPOT_REGISTER(neumann_condition_types);
      GETPOT_REGISTER(neumann_condition_values);
      GETPOT_REGISTER(neumann_condition_boundaries);
      GETPOT_REGISTER(neumann_condition_variables);
  
      const unsigned int n_neumann_conditions=
        input.vector_variable_size("neumann_condition_types");
  
      if (n_neumann_conditions !=
          input.vector_variable_size("neumann_condition_values"))
        {
          libMesh::out << "Error: " << n_neumann_conditions
                       << " Neumann condition types does not match "
                       << input.vector_variable_size("neumann_condition_values")
                       << " Neumann condition values." << std::endl;
  
          libmesh_error();
        }
  
      if (n_neumann_conditions !=
          input.vector_variable_size("neumann_condition_boundaries"))
        {
          libMesh::out << "Error: " << n_neumann_conditions
                       << " Neumann condition types does not match "
                       << input.vector_variable_size("neumann_condition_boundaries")
                       << " Neumann condition boundaries." << std::endl;
  
          libmesh_error();
        }
  
      if (n_neumann_conditions !=
          input.vector_variable_size("neumann_condition_variables"))
        {
          libMesh::out << "Error: " << n_neumann_conditions
                       << " Neumann condition types does not match "
                       << input.vector_variable_size("neumann_condition_variables")
                       << " Neumann condition variables sets." << std::endl;
  
          libmesh_error();
        }
  
      for (unsigned int i=0; i != n_neumann_conditions; ++i)
        {
          const std::string func_type =
            input("neumann_condition_types", zero_string, i);
  
          const std::string func_value =
            input("neumann_condition_values", empty_string, i);
  
          const boundary_id_type func_boundary =
            input("neumann_condition_boundaries", boundary_id_type(0), i);
  
          neumann_conditions[func_boundary] =
            (new_function_base(func_type, func_value).release());
  
          const std::string variable_set =
            input("neumann_condition_variables", empty_string, i);
  
          for (unsigned int j=0; j != variable_set.size(); ++j)
            {
              if (variable_set[j] == '1')
                neumann_condition_variables[func_boundary].push_back(j);
              else if (variable_set[j] != '0')
                {
                  libMesh::out << "Unable to understand Neumann variable set"
                               << variable_set << std::endl;
                  libmesh_error();
                }
            }
        }
  
      GETPOT_REGISTER(initial_condition_types);
      GETPOT_REGISTER(initial_condition_values);
      GETPOT_REGISTER(initial_condition_subdomains);
  
      const unsigned int n_initial_conditions=
        input.vector_variable_size("initial_condition_types");
  
      if (n_initial_conditions !=
          input.vector_variable_size("initial_condition_values"))
        {
          libMesh::out << "Error: " << n_initial_conditions
                       << " initial condition types does not match "
                       << input.vector_variable_size("initial_condition_values")
                       << " initial condition values." << std::endl;
  
          libmesh_error();
        }
  
      if (n_initial_conditions !=
          input.vector_variable_size("initial_condition_subdomains"))
        {
          libMesh::out << "Error: " << n_initial_conditions
                       << " initial condition types does not match "
                       << input.vector_variable_size("initial_condition_subdomains")
                       << " initial condition subdomains." << std::endl;
  
          libmesh_error();
        }
  
      for (unsigned int i=0; i != n_initial_conditions; ++i)
        {
          const std::string func_type =
            input("initial_condition_types", zero_string, i);
  
          const std::string func_value =
            input("initial_condition_values", empty_string, i);
  
          const subdomain_id_type func_subdomain =
            input("initial_condition_subdomains", subdomain_id_type(0), i);
  
          initial_conditions[func_subdomain] =
            (new_function_base(func_type, func_value).release());
        }
  
      GETPOT_REGISTER(other_interior_function_types);
      GETPOT_REGISTER(other_interior_function_values);
      GETPOT_REGISTER(other_interior_function_subdomains);
      GETPOT_REGISTER(other_interior_function_ids);
  
      const unsigned int n_other_interior_functions =
        input.vector_variable_size("other_interior_function_types");
  
      if (n_other_interior_functions !=
          input.vector_variable_size("other_interior_function_values"))
        {
          libMesh::out << "Error: " << n_other_interior_functions
                       << " other interior function types does not match "
                       << input.vector_variable_size("other_interior_function_values")
                       << " other interior function values." << std::endl;
  
          libmesh_error();
        }
  
      if (n_other_interior_functions !=
          input.vector_variable_size("other_interior_function_subdomains"))
        {
          libMesh::out << "Error: " << n_other_interior_functions
                       << " other interior function types does not match "
                       << input.vector_variable_size("other_interior_function_subdomains")
                       << " other interior function subdomains." << std::endl;
  
          libmesh_error();
        }
  
      if (n_other_interior_functions !=
          input.vector_variable_size("other_interior_function_ids"))
        {
          libMesh::out << "Error: " << n_other_interior_functions
                       << " other interior function types does not match "
                       << input.vector_variable_size("other_interior_function_ids")
                       << " other interior function ids." << std::endl;
  
          libmesh_error();
        }
  
      for (unsigned int i=0; i != n_other_interior_functions; ++i)
        {
          const std::string func_type =
            input("other_interior_function_types", zero_string, i);
  
          const std::string func_value =
            input("other_interior_function_values", empty_string, i);
  
          const subdomain_id_type func_subdomain =
            input("other_interior_condition_subdomains", subdomain_id_type(0), i);
  
          const int func_id =
            input("other_interior_condition_ids", int(0), i);
  
          other_interior_functions[func_id][func_subdomain] =
            (new_function_base(func_type, func_value).release());
        }
  
      GETPOT_REGISTER(other_boundary_function_types);
      GETPOT_REGISTER(other_boundary_function_values);
      GETPOT_REGISTER(other_boundary_function_boundaries);
      GETPOT_REGISTER(other_boundary_function_ids);
  
      const unsigned int n_other_boundary_functions =
        input.vector_variable_size("other_boundary_function_types");
  
      if (n_other_boundary_functions !=
          input.vector_variable_size("other_boundary_function_values"))
        {
          libMesh::out << "Error: " << n_other_boundary_functions
                       << " other boundary function types does not match "
                       << input.vector_variable_size("other_boundary_function_values")
                       << " other boundary function values." << std::endl;
  
          libmesh_error();
        }
  
      if (n_other_boundary_functions !=
          input.vector_variable_size("other_boundary_function_boundaries"))
        {
          libMesh::out << "Error: " << n_other_boundary_functions
                       << " other boundary function types does not match "
                       << input.vector_variable_size("other_boundary_function_boundaries")
                       << " other boundary function boundaries." << std::endl;
  
          libmesh_error();
        }
  
      if (n_other_boundary_functions !=
          input.vector_variable_size("other_boundary_function_ids"))
        {
          libMesh::out << "Error: " << n_other_boundary_functions
                       << " other boundary function types does not match "
                       << input.vector_variable_size("other_boundary_function_ids")
                       << " other boundary function ids." << std::endl;
  
          libmesh_error();
        }
  
      for (unsigned int i=0; i != n_other_boundary_functions; ++i)
        {
          const std::string func_type =
            input("other_boundary_function_types", zero_string, i);
  
          const std::string func_value =
            input("other_boundary_function_values", empty_string, i);
  
          const boundary_id_type func_boundary =
            input("other_boundary_function_boundaries", boundary_id_type(0), i);
  
          const int func_id =
            input("other_boundary_function_ids", int(0), i);
  
          other_boundary_functions[func_id][func_boundary] =
            (new_function_base(func_type, func_value).release());
        }
  
      GETPOT_INPUT(run_simulation);
      GETPOT_INPUT(run_postprocess);
  
  
      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(extra_quadrature_order);
  
  
      GETPOT_INPUT(analytic_jacobians);
      GETPOT_INPUT(verify_analytic_jacobians);
      GETPOT_INPUT(numerical_jacobian_h);
      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);
      GETPOT_INPUT(print_element_jacobians);
  
  
      GETPOT_INPUT(use_petsc_snes);
      GETPOT_INPUT(time_solver_quiet);
      GETPOT_INPUT(solver_quiet);
      GETPOT_INPUT(solver_verbose);
      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(initial_sobolev_order);
      GETPOT_INT_INPUT(initial_extra_quadrature);
      GETPOT_INPUT(indicator_type);
      GETPOT_INT_INPUT(sobolev_order);
  
      GETPOT_INPUT(system_config_file);
  
    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 heatsystem.C without comments:

 
  #include "libmesh/getpot.h"
  
  #include "heatsystem.h"
  
  #include "libmesh/fe_base.h"
  #include "libmesh/fe_interface.h"
  #include "libmesh/fem_context.h"
  #include "libmesh/mesh.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/system.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/zero_function.h"
  #include "libmesh/boundary_info.h"
  #include "libmesh/dirichlet_boundaries.h"
  #include "libmesh/dof_map.h"
  
  void HeatSystem::init_data ()
  {
    T_var = this->add_variable ("T", static_cast<Order>(_fe_order),
                        Utility::string_to_enum<FEFamily>(_fe_family));
  
    {
      std::ifstream i("heat.in");
      if (!i)
        {
          std::cerr << '[' << libMesh::processor_id()
                    << "] Can't find heat.in; exiting early."
                    << std::endl;
          libmesh_error();
        }
    }
    GetPot infile("heat.in");
    _k = infile("k", 1.0);
    _analytic_jacobians = infile("analytic_jacobians", true);
  
    parameters.push_back(_k);
  
    this->get_equation_systems().parameters.set<Real>("_k") = _k;
  
    this->time_evolving(T_var);
  
    const boundary_id_type all_ids[4] = {0, 1, 2, 3};
    std::set<boundary_id_type> all_bdys(all_ids, all_ids+(2*2));
  
    std::vector<unsigned int> T_only(1, T_var);
  
    ZeroFunction<Number> zero;
  
    this->get_dof_map().add_dirichlet_boundary
          (DirichletBoundary (all_bdys, T_only, &zero));
  
    FEMSystem::init_data();
  }
  
  
  
  void HeatSystem::init_context(DiffContext &context)
  {
    FEMContext &c = libmesh_cast_ref<FEMContext&>(context);
  
    c.element_fe_var[0]->get_JxW();
    c.element_fe_var[0]->get_dphi();
  
    if (c.is_adjoint())
      {
        const System & sys = c.get_system();
  
        NumericVector<Number> &adjoint_solution =
          const_cast<System &>(sys).get_adjoint_solution(0);
  
        c.add_localized_vector(adjoint_solution, sys);
      }
  
    FEMSystem::init_context(context);
  }
  
  #define optassert(X) {if (!(X)) libmesh_error();}
  
  bool HeatSystem::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();
  
    optassert(c.dof_indices_var.size() > 0);
  
    const unsigned int n_u_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.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_u_dofs; i++)
          F(i) += JxW[qp] * -parameters[0] * (grad_T * dphi[i][qp]);
        if (compute_jacobian)
          for (unsigned int i=0; i != n_u_dofs; i++)
            for (unsigned int j=0; j != n_u_dofs; ++j)
              K(i,j) += JxW[qp] * -parameters[0] * (dphi[i][qp] * dphi[j][qp]);
      } // end of the quadrature point qp-loop
  
    return compute_jacobian;
  }
  
  void HeatSystem::perturb_accumulate_residuals(const ParameterVector& parameters_in)
  {
    const unsigned int Np = parameters_in.size();
  
    for (unsigned int j=0; j != Np; ++j)
      {
        Number old_parameter = *parameters_in[j];
  
        *parameters_in[j] = old_parameter - dp;
  
        this->assembly(true, false);
  
        this->rhs->close();
  
        AutoPtr<NumericVector<Number> > R_minus = this->rhs->clone();
  
        R_minus_dp += -R_minus->dot(this->get_adjoint_solution(0));
  
        *parameters_in[j] = old_parameter + dp;
  
        this->assembly(true, false);
  
        this->rhs->close();
  
        AutoPtr<NumericVector<Number> > R_plus = this->rhs->clone();
  
        R_plus_dp += -R_plus->dot(this->get_adjoint_solution(0));
  
        *parameters_in[j] = old_parameter;
  
      }
  }



The source file initial.C without comments:

 
  #include "initial.h"
  
  using namespace libMesh;
  
  void read_initial_parameters()
  {
  }
  
  void finish_initialization()
  {
  }
  
  
  
  Number initial_value(const Point& p,
                       const Parameters&,
                       const std::string&,
                       const std::string&)
  {
    Real x = p(0), y = p(1);
  
    return sin(M_PI * x) * sin(M_PI * y);
  }
  
  
  
  Gradient initial_grad(const Point& p,
                        const Parameters&,
                        const std::string&,
                        const std::string&)
  {
    Real x = p(0), y = p(1);
  
    return Gradient(M_PI*cos(M_PI * x) * sin(M_PI * y),
    M_PI*sin(M_PI * x) * cos(M_PI * y));
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex5'
***************************************************************
* Running Example adjoints_ex5:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
Started /net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex5/.libs/lt-example-devel
Building mesh
Building system
*** Warning, This code is untested, experimental, or likely to see future API changes: ../../../include/libmesh/memory_solution_history.h, line 41, compiled Apr 19 2013 at 11:48:31 ***
Initializing systems
Setting primal initial conditions
|U(0)|= 2.27632

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=4225
    n_local_nodes()=1099
  n_elem()=5456
    n_local_elem()=1395
    n_active_elem()=4096
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "HeatSystem"
    Type "Implicit"
    Variables="T" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=4225
    n_local_dofs()=1099
    n_constrained_dofs()=256
    n_local_constrained_dofs()=63
    n_vectors()=3
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 8.71976
      Average Off-Processor Bandwidth <= 0.205917
      Maximum  On-Processor Bandwidth <= 11
      Maximum Off-Processor Bandwidth <= 6
    DofMap Constraints
      Number of DoF Constraints = 256
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0

 Solving time step 0, time = 0
  Nonlinear solver converged, step 0, residual reduction 9.43741e-13 < 1e-09
|U(0.1)|= 2.27183
Advancing timestep

 Solving time step 1, time = 0.1
  Nonlinear solver converged, step 0, residual reduction 9.41233e-13 < 1e-09
|U(0.2)|= 2.26736
Advancing timestep

 Solving time step 2, time = 0.2
  Nonlinear solver converged, step 0, residual reduction 9.36696e-13 < 1e-09
|U(0.3)|= 2.26289
Advancing timestep

 Solving time step 3, time = 0.3
  Nonlinear solver converged, step 0, residual reduction 9.36474e-13 < 1e-09
|U(0.4)|= 2.25843
Advancing timestep

 Solving time step 4, time = 0.4
  Nonlinear solver converged, step 0, residual reduction 9.35902e-13 < 1e-09
|U(0.5)|= 2.25398
Advancing timestep

 Solving time step 5, time = 0.5
  Nonlinear solver converged, step 0, residual reduction 9.33977e-13 < 1e-09
|U(0.6)|= 2.24954
Advancing timestep

 Solving time step 6, time = 0.6
  Nonlinear solver converged, step 0, residual reduction 9.31019e-13 < 1e-09
|U(0.7)|= 2.24511
Advancing timestep

 Solving time step 7, time = 0.7
  Nonlinear solver converged, step 0, residual reduction 9.3018e-13 < 1e-09
|U(0.8)|= 2.24068
Advancing timestep

 Solving time step 8, time = 0.8
  Nonlinear solver converged, step 0, residual reduction 9.2767e-13 < 1e-09
|U(0.9)|= 2.23627
Advancing timestep

 Solving time step 9, time = 0.9
  Nonlinear solver converged, step 0, residual reduction 9.2706e-13 < 1e-09
|U(1)|= 2.23186
Advancing timestep


Solving the adjoint problem
Setting adjoint initial conditions Z(1)
Retrieving solutions at time t=1
|U(1.1)|= 2.23186
|U(1)|= 2.23186
|Z(1)|= 2.27632

 Solving adjoint time step 0, time = 1
Retrieving solutions at time t=1
|U(1)|= 2.23186
|U(0.9)|= 2.23627
|Z(0.9)|= 2.27183

 Solving adjoint time step 1, time = 0.9
Retrieving solutions at time t=0.9
|U(0.9)|= 2.23627
|U(0.8)|= 2.24068
|Z(0.8)|= 2.26736

 Solving adjoint time step 2, time = 0.8
Retrieving solutions at time t=0.8
|U(0.8)|= 2.24068
|U(0.7)|= 2.24511
|Z(0.7)|= 2.26289

 Solving adjoint time step 3, time = 0.7
Retrieving solutions at time t=0.7
|U(0.7)|= 2.24511
|U(0.6)|= 2.24954
|Z(0.6)|= 2.25843

 Solving adjoint time step 4, time = 0.6
Retrieving solutions at time t=0.6
|U(0.6)|= 2.24954
|U(0.5)|= 2.25398
|Z(0.5)|= 2.25398

 Solving adjoint time step 5, time = 0.5
Retrieving solutions at time t=0.5
|U(0.5)|= 2.25398
|U(0.4)|= 2.25843
|Z(0.4)|= 2.24954

 Solving adjoint time step 6, time = 0.4
Retrieving solutions at time t=0.4
|U(0.4)|= 2.25843
|U(0.3)|= 2.26289
|Z(0.3)|= 2.24511

 Solving adjoint time step 7, time = 0.3
Retrieving solutions at time t=0.3
|U(0.3)|= 2.26289
|U(0.2)|= 2.26736
|Z(0.2)|= 2.24068

 Solving adjoint time step 8, time = 0.2
Retrieving solutions at time t=0.2
|U(0.2)|= 2.26736
|U(0.1)|= 2.27183
|Z(0.1)|= 2.23627

 Solving adjoint time step 9, time = 0.1
Retrieving solutions at time t=0.1
|U(0.1)|= 2.27183
|U(2.77556e-17)|= 2.27632
|Z(2.77556e-17)|= 2.23186

Retrieving 0, time = 2.77556e-17
|U(0.1)|= 2.27183
|U(2.77556e-17)|= 2.27632
|Z(2.77556e-17)|= 2.23186

Retrieving 1, time = 0.1
|U(0.2)|= 2.26736
|U(0.1)|= 2.27183
|Z(0.1)|= 2.23627

Retrieving 2, time = 0.2
|U(0.3)|= 2.26289
|U(0.2)|= 2.26736
|Z(0.2)|= 2.24068

Retrieving 3, time = 0.3
|U(0.4)|= 2.25843
|U(0.3)|= 2.26289
|Z(0.3)|= 2.24511

Retrieving 4, time = 0.4
|U(0.5)|= 2.25398
|U(0.4)|= 2.25843
|Z(0.4)|= 2.24954

Retrieving 5, time = 0.5
|U(0.6)|= 2.24954
|U(0.5)|= 2.25398
|Z(0.5)|= 2.25398

Retrieving 6, time = 0.6
|U(0.7)|= 2.24511
|U(0.6)|= 2.24954
|Z(0.6)|= 2.25843

Retrieving 7, time = 0.7
|U(0.8)|= 2.24068
|U(0.7)|= 2.24511
|Z(0.7)|= 2.26289

Retrieving 8, time = 0.8
|U(0.9)|= 2.23627
|U(0.8)|= 2.24068
|Z(0.8)|= 2.26736

Retrieving 9, time = 0.9
|U(1)|= 2.23186
|U(0.9)|= 2.23627
|Z(0.9)|= 2.27183

Retrieving  final time = 1
|U(1.1)|= 2.23186
|U(1)|= 2.23186
|Z(1)|= 2.27632

Sensitivity of QoI 0 w.r.t parameter 0 is: -5.30953

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:49:06 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-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/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=9.21305, Active time=7.81359                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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()     5         0.0050      0.001002    0.0063      0.001253    0.06     0.08     |
|   build_constraint_matrix()        63488     0.0654      0.000001    0.0654      0.000001    0.84     0.84     |
|   build_sparsity()                 5         0.0039      0.000784    0.0141      0.002822    0.05     0.18     |
|   cnstrn_elem_mat_vec()            10240     0.0155      0.000002    0.0155      0.000002    0.20     0.20     |
|   constrain_elem_matrix()          10240     0.0108      0.000001    0.0108      0.000001    0.14     0.14     |
|   constrain_elem_vector()          43008     0.0398      0.000001    0.0398      0.000001    0.51     0.51     |
|   create_dof_constraints()         5         0.0192      0.003849    0.0454      0.009088    0.25     0.58     |
|   distribute_dofs()                5         0.0108      0.002164    0.0412      0.008248    0.14     0.53     |
|   dof_indices()                    285053    0.8995      0.000003    0.8995      0.000003    11.51    11.51    |
|   enforce_constraints_exactly()    50        0.0119      0.000239    0.0119      0.000239    0.15     0.15     |
|   old_dof_indices()                5440      0.0190      0.000003    0.0190      0.000003    0.24     0.24     |
|   prepare_send_list()              5         0.0000      0.000006    0.0000      0.000006    0.00     0.00     |
|   reinit()                         5         0.0164      0.003283    0.0164      0.003283    0.21     0.21     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          22        0.0438      0.001993    0.1992      0.009056    0.56     2.55     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        425264    1.1779      0.000003    1.1779      0.000003    15.07    15.07    |
|   init_shape_functions()           6649      0.0265      0.000004    0.0265      0.000004    0.34     0.34     |
|   inverse_map()                    15890     0.0402      0.000003    0.0402      0.000003    0.51     0.51     |
|                                                                                                                |
| FEMSystem                                                                                                      |
|   assemble_qoi_derivative()        10        0.6396      0.063960    0.8595      0.085948    8.19     11.00    |
|   assembly()                       10        0.2145      0.021448    0.6289      0.062888    2.74     8.05     |
|   assembly(get_jacobian)           10        0.2191      0.021910    0.6592      0.065919    2.80     8.44     |
|   assembly(get_residual)           32        0.5898      0.018432    1.9361      0.060502    7.55     24.78    |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             425264    0.6627      0.000002    0.6627      0.000002    8.48     8.48     |
|   compute_face_map()               6448      0.0264      0.000004    0.0601      0.000009    0.34     0.77     |
|   init_face_shape_functions()      104       0.0003      0.000003    0.0003      0.000003    0.00     0.00     |
|   init_reference_to_physical_map() 6649      0.0206      0.000003    0.0206      0.000003    0.26     0.26     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               22        0.3610      0.016410    0.3627      0.016485    4.62     4.64     |
|                                                                                                                |
| ImplicitSystem                                                                                                 |
|   adjoint_solve()                  10        0.0006      0.000057    1.8579      0.185793    0.01     23.78    |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           16320     0.0202      0.000001    0.0202      0.000001    0.26     0.26     |
|   init()                           8         0.0025      0.000312    0.0025      0.000312    0.03     0.03     |
|                                                                                                                |
| Mesh                                                                                                           |
|   contract()                       4         0.0008      0.000209    0.0016      0.000400    0.01     0.02     |
|   find_neighbors()                 5         0.0281      0.005616    0.0285      0.005697    0.36     0.36     |
|   renumber_nodes_and_elem()        14        0.0027      0.000194    0.0027      0.000194    0.03     0.03     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        6         0.0228      0.003807    0.0228      0.003807    0.29     0.29     |
|   find_global_indices()            6         0.0034      0.000563    0.0284      0.004729    0.04     0.36     |
|   parallel_sort()                  6         0.0009      0.000146    0.0014      0.000225    0.01     0.02     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         22        0.0006      0.000029    0.5641      0.025641    0.01     7.22     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _coarsen_elements()              4         0.0007      0.000182    0.0011      0.000285    0.01     0.01     |
|   _refine_elements()               8         0.0298      0.003724    0.0904      0.011300    0.38     1.16     |
|   add_point()                      16320     0.0305      0.000002    0.0526      0.000003    0.39     0.67     |
|   make_coarsening_compatible()     8         0.0098      0.001230    0.0098      0.001230    0.13     0.13     |
|   make_flags_parallel_consistent() 8         0.0070      0.000880    0.0132      0.001645    0.09     0.17     |
|   make_refinement_compatible()     8         0.0000      0.000002    0.0000      0.000006    0.00     0.00     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0001      0.000091    0.0001      0.000091    0.00     0.00     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      5         0.0669      0.013372    0.0951      0.019025    0.86     1.22     |
|                                                                                                                |
| NewtonSolver                                                                                                   |
|   solve()                          10        0.9284      0.092839    2.2089      0.220890    11.88    28.27    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      29        0.0098      0.000338    0.0103      0.000356    0.13     0.13     |
|   max(bool)                        21        0.0077      0.000367    0.0077      0.000367    0.10     0.10     |
|   max(scalar)                      1736      0.0075      0.000004    0.0075      0.000004    0.10     0.10     |
|   max(vector)                      414       0.0027      0.000006    0.0078      0.000019    0.03     0.10     |
|   min(bool)                        2154      0.0150      0.000007    0.0150      0.000007    0.19     0.19     |
|   min(scalar)                      1722      0.1287      0.000075    0.1287      0.000075    1.65     1.65     |
|   min(vector)                      414       0.0031      0.000008    0.0083      0.000020    0.04     0.11     |
|   probe()                          252       0.0073      0.000029    0.0073      0.000029    0.09     0.09     |
|   receive()                        252       0.0008      0.000003    0.0081      0.000032    0.01     0.10     |
|   send()                           252       0.0005      0.000002    0.0005      0.000002    0.01     0.01     |
|   send_receive()                   264       0.0011      0.000004    0.0101      0.000038    0.01     0.13     |
|   sum()                            222       0.5955      0.002683    0.6735      0.003034    7.62     8.62     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           252       0.0003      0.000001    0.0003      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         5         0.0046      0.000926    0.0151      0.003029    0.06     0.19     |
|   set_parent_processor_ids()       5         0.0022      0.000431    0.0022      0.000431    0.03     0.03     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          20        0.3934      0.019670    0.3934      0.019670    5.03     5.03     |
|                                                                                                                |
| ProjectVector                                                                                                  |
|   operator()                       8         0.0097      0.001212    0.0375      0.004688    0.12     0.48     |
|                                                                                                                |
| System                                                                                                         |
|   calculate_norm()                 77        0.2823      0.003666    1.4223      0.018472    3.61     18.20    |
|   project_vector()                 10        0.0457      0.004569    0.1026      0.010262    0.58     1.31     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            1344805   7.8136                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example adjoints_ex5:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/adjoints/adjoints_ex5'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo