The source file assembly.h with comments:

        #ifndef __assembly_h__
        #define __assembly_h__
        
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/fe.h"
        #include "libmesh/fe_interface.h"
        #include "libmesh/fe_base.h"
        #include "libmesh/elem_assembly.h"
        #include "libmesh/quadrature_gauss.h"
        
rbOOmit includes
        #include "libmesh/rb_theta.h"
        #include "libmesh/rb_assembly_expansion.h"
        
Bring in bits from the libMesh namespace. Just the bits we're using, since this is a header.
        using libMesh::ElemAssembly;
        using libMesh::FEInterface;
        using libMesh::FEMContext;
        using libMesh::Number;
        using libMesh::Point;
        using libMesh::RBAssemblyExpansion;
        using libMesh::RBParameters;
        using libMesh::RBTheta;
        using libMesh::RBThetaExpansion;
        using libMesh::Real;
        using libMesh::RealGradient;
        
Functors for the parameter-dependent part of the affine decomposition of the PDE The RHS and outputs just require a constant value of 1, so use a default RBTheta object there
        struct ThetaA0 : RBTheta { virtual Number evaluate(const RBParameters& )   { return 0.05;  } };
        struct ThetaA1 : RBTheta { virtual Number evaluate(const RBParameters& mu) { return mu.get_value("x_vel"); } };
        struct ThetaA2 : RBTheta { virtual Number evaluate(const RBParameters& mu) { return mu.get_value("y_vel"); } };
        
        struct A0 : ElemAssembly
        {
Assemble the Laplacian operator
          virtual void interior_assembly(FEMContext &c)
          {
            const unsigned int u_var = 0;
        
            const std::vector<Real> &JxW =
              c.element_fe_var[u_var]->get_JxW();
        
The velocity shape function gradients at interior quadrature points.
            const std::vector<std::vector<RealGradient> >& dphi =
              c.element_fe_var[u_var]->get_dphi();
        
The number of local degrees of freedom in each variable
            const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
Now we will build the affine operator
            unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
            for (unsigned int qp=0; qp != n_qpoints; qp++)
              for (unsigned int i=0; i != n_u_dofs; i++)
                for (unsigned int j=0; j != n_u_dofs; j++)
                  c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
          }
        };
        
        
        struct A1 : ElemAssembly
        {
Convection in the x-direction
          virtual void interior_assembly(FEMContext &c)
          {
            const unsigned int u_var = 0;
        
            const std::vector<Real> &JxW =
              c.element_fe_var[u_var]->get_JxW();
        
            const std::vector<std::vector<Real> >& phi =
              c.element_fe_var[u_var]->get_phi();
        
            const std::vector<std::vector<RealGradient> >& dphi =
              c.element_fe_var[u_var]->get_dphi();
        
The number of local degrees of freedom in each variable
            const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
Now we will build the affine operator
            unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
            for (unsigned int qp=0; qp != n_qpoints; qp++)
              for (unsigned int i=0; i != n_u_dofs; i++)
                for (unsigned int j=0; j != n_u_dofs; j++)
                  c.get_elem_jacobian()(i,j) += JxW[qp] *dphi[j][qp](0)*phi[i][qp];
          }
        };
        
        struct A2 : ElemAssembly
        {
Convection in the y-direction
          virtual void interior_assembly(FEMContext &c)
          {
            const unsigned int u_var = 0;
        
            const std::vector<Real> &JxW =
              c.element_fe_var[u_var]->get_JxW();
        
            const std::vector<std::vector<Real> >& phi =
              c.element_fe_var[u_var]->get_phi();
        
            const std::vector<std::vector<RealGradient> >& dphi =
              c.element_fe_var[u_var]->get_dphi();
        
The number of local degrees of freedom in each variable
            const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
Now we will build the affine operator
            unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
            for (unsigned int qp=0; qp != n_qpoints; qp++)
              for (unsigned int i=0; i != n_u_dofs; i++)
                for (unsigned int j=0; j != n_u_dofs; j++)
                  c.get_elem_jacobian()(i,j) += JxW[qp] *dphi[j][qp](1)*phi[i][qp];
          }
        };
        
        struct F0 : ElemAssembly
        {
Source term, 1 throughout the domain
          virtual void interior_assembly(FEMContext &c)
          {
            const unsigned int u_var = 0;
        
            const std::vector<Real> &JxW =
              c.element_fe_var[u_var]->get_JxW();
        
            const std::vector<std::vector<Real> >& phi =
              c.element_fe_var[u_var]->get_phi();
        
The number of local degrees of freedom in each variable
            const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
Now we will build the affine operator
            unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
            for (unsigned int qp=0; qp != n_qpoints; qp++)
              for (unsigned int i=0; i != n_u_dofs; i++)
                c.get_elem_residual()(i) += JxW[qp] * ( 1.*phi[i][qp] );
          }
        };
        
        struct OutputAssembly : ElemAssembly
        {
          OutputAssembly(Real min_x_in, Real max_x_in,
                         Real min_y_in, Real max_y_in)
                        :
                        min_x(min_x_in),
                        max_x(max_x_in),
                        min_y(min_y_in),
                        max_y(max_y_in)
          {}
        
Output: Average value over the region [min_x,max_x]x[min_y,max_y]
          virtual void interior_assembly(FEMContext &c)
          {
            const unsigned int u_var = 0;
        
            const std::vector<Real> &JxW =
              c.element_fe_var[u_var]->get_JxW();
        
            const std::vector<std::vector<Real> >& phi =
              c.element_fe_var[u_var]->get_phi();
        
The number of local degrees of freedom in each variable
            const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
        
Now we will build the affine operator
            unsigned int n_qpoints = (c.get_element_qrule())->n_points();
        
            Real output_area = (max_x-min_x) * (max_y-min_y);
        
            Point centroid = c.elem->centroid();
            if( (min_x <= centroid(0)) && (centroid(0) <= max_x) &&
                (min_y <= centroid(1)) && (centroid(1) <= max_y) )
              for (unsigned int qp=0; qp != n_qpoints; qp++)
                for (unsigned int i=0; i != n_u_dofs; i++)
                  c.get_elem_residual()(i) += JxW[qp] * ( 1.*phi[i][qp] ) / output_area;
          }
        
Member variables that define the output region in 2D
          Real min_x, max_x, min_y, max_y;
        };
        
Define an RBThetaExpansion class for this PDE
        struct CDRBThetaExpansion : RBThetaExpansion
        {
        
          /**
           * Constructor.
           */
          CDRBThetaExpansion()
          {
set up the RBThetaExpansion object
            attach_A_theta(&theta_a_0);   // Attach the lhs theta
            attach_A_theta(&theta_a_1);
            attach_A_theta(&theta_a_2);
        
            attach_F_theta(&rb_theta);    // Attach the rhs theta
        
            attach_output_theta(&rb_theta); // Attach output 0 theta
            attach_output_theta(&rb_theta); // Attach output 1 theta
            attach_output_theta(&rb_theta); // Attach output 2 theta
            attach_output_theta(&rb_theta); // Attach output 3 theta
          }
        
The RBTheta member variables
          ThetaA0 theta_a_0;
          ThetaA1 theta_a_1;
          ThetaA2 theta_a_2;
          RBTheta rb_theta; // Default RBTheta object, just returns 1.
        };
        
Define an RBAssemblyExpansion class for this PDE
        struct CDRBAssemblyExpansion : RBAssemblyExpansion
        {
        
          /**
           * Constructor.
           */
          CDRBAssemblyExpansion()
            :
            L0(0.7,0.8,0.7,0.8),
            L1(0.2,0.3,0.7,0.8),
            L2(0.2,0.3,0.2,0.3),
            L3(0.7,0.8,0.2,0.3)
          {
And set up the RBAssemblyExpansion object
            attach_A_assembly(&A0_assembly); // Attach the lhs assembly
            attach_A_assembly(&A1_assembly);
            attach_A_assembly(&A2_assembly);
        
            attach_F_assembly(&F0_assembly); // Attach the rhs assembly
        
            attach_output_assembly(&L0);       // Attach output 0 assembly
            attach_output_assembly(&L1);       // Attach output 1 assembly
            attach_output_assembly(&L2);       // Attach output 2 assembly
            attach_output_assembly(&L3);       // Attach output 3 assembly
          }
        
The ElemAssembly objects
          A0 A0_assembly;
          A1 A1_assembly;
          A2 A2_assembly;
          F0 F0_assembly;
          OutputAssembly L0;
          OutputAssembly L1;
          OutputAssembly L2;
          OutputAssembly L3;
        };
        
        #endif



The source file rb_classes.h with comments:

        #ifndef __rb_classes_h__
        #define __rb_classes_h__
        
        #include "libmesh/rb_construction.h"
        #include "libmesh/fe_base.h"
        
local include
        #include "assembly.h"
        
Bring in bits from the libMesh namespace. Just the bits we're using, since this is a header.
        using libMesh::AutoPtr;
        using libMesh::DirichletBoundary;
        using libMesh::EquationSystems;
        using libMesh::FEMContext;
        using libMesh::RBConstruction;
        using libMesh::RBEvaluation;
        using libMesh::Real;
        
        
A simple subclass of RBEvaluation, which just needs to specify (a lower bound for) the coercivity constant for this problem. For this simple convection-diffusion problem, we can set the coercivity constant lower bound to 0.05.
        class SimpleRBEvaluation : public RBEvaluation
        {
        public:
        
          /**
           * Constructor. Just set the theta expansion.
           */
          SimpleRBEvaluation(const Parallel::Communicator& comm)
            : RBEvaluation(comm)
          {
            set_rb_theta_expansion(cd_rb_theta_expansion);
          }
        
          /**
           * The coercivity constant is bounded below by 0.05.
           */
          virtual Real get_stability_lower_bound() { return 0.05; }
        
          /**
           * The object that stores the "theta" expansion of the parameter dependent PDE,
           * i.e. the set of parameter-dependent functions in the affine expansion of the PDE.
           */
          CDRBThetaExpansion cd_rb_theta_expansion;
        
        };
        
A simple subclass of Construction, which just needs to override build_rb_evaluation in order to build a SimpleRBEvaluation object, rather than an RBEvaluation object.
        class SimpleRBConstruction : public RBConstruction
        {
        public:
        
          SimpleRBConstruction (EquationSystems& es,
                                const std::string& name_in,
                                const unsigned int number_in)
          : Parent(es, name_in, number_in),
            dirichlet_bc(AutoPtr<DirichletBoundary>(NULL))
          {}
        
          /**
           * Destructor.
           */
          virtual ~SimpleRBConstruction () { }
        
          /**
           * The type of system.
           */
          typedef SimpleRBConstruction sys_type;
        
          /**
           * The type of the parent.
           */
          typedef RBConstruction Parent;
        
          /**
           * Initialize data structures.
           */
          virtual void init_data()
          {
            u_var = this->add_variable ("u", FIRST);
        
Generate a DirichletBoundary object
            dirichlet_bc = build_zero_dirichlet_boundary_object();
        
Set the Dirichet boundary IDs and the Dirichlet boundary variable numbers
            dirichlet_bc->b.insert(0);
            dirichlet_bc->b.insert(1);
            dirichlet_bc->b.insert(2);
            dirichlet_bc->b.insert(3);
            dirichlet_bc->variables.push_back(u_var);
        
Attach dirichlet_bc (must do this _before_ Parent::init_data)
            get_dof_map().add_dirichlet_boundary(*dirichlet_bc);
        
            Parent::init_data();
        
Set the rb_assembly_expansion for this Construction object. The theta expansion comes from the RBEvaluation object.
            set_rb_assembly_expansion(cd_rb_assembly_expansion);
        
We need to define an inner product matrix for this problem
            set_inner_product_assembly(cd_rb_assembly_expansion.A0_assembly);
          }
        
          /**
           * Pre-request all relevant element data.
           */
          virtual void init_context(FEMContext &c)
          {
For efficiency, we should prerequest all the data we will need to build the linear system before doing an element loop.
            c.element_fe_var[u_var]->get_JxW();
            c.element_fe_var[u_var]->get_phi();
            c.element_fe_var[u_var]->get_dphi();
          }
        
          /**
           * Variable number for u.
           */
          unsigned int u_var;
        
          /**
           * The object that stores the "assembly" expansion of the parameter dependent PDE,
           * i.e. the objects that define how to assemble the set of parameter-independent
           * operators in the affine expansion of the PDE.
           */
          CDRBAssemblyExpansion cd_rb_assembly_expansion;
        
          /**
           * The object that defines which degrees of freedom are on a Dirichlet boundary.
           */
          AutoPtr<DirichletBoundary> dirichlet_bc;
        
        };
        
        #endif



The source file reduced_basis_ex1.C with comments:

        #include <iostream>
        #include <algorithm>
        #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
        #include <cmath>
        #include <set>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/getpot.h"
        #include "libmesh/elem.h"
        
local includes
        #include "rb_classes.h"
        #include "assembly.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
The main program.
        int main (int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
        
        #if !defined(LIBMESH_HAVE_XDR)
We need XDR support to write out reduced bases
          libmesh_example_assert(false, "--enable-xdr");
        #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
XDR binary support requires double precision
          libmesh_example_assert(false, "--disable-singleprecision");
        #endif
FIXME: This example currently segfaults with Trilinos?
          libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
        
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Parse the input file (reduced_basis_ex1.in) using GetPot
          std::string parameters_filename = "reduced_basis_ex1.in";
          GetPot infile(parameters_filename);
        
          unsigned int n_elem = infile("n_elem", 1);       // Determines the number of elements in the "truth" mesh
          const unsigned int dim = 2;                      // The number of spatial dimensions
        
          bool store_basis_functions = infile("store_basis_functions", true); // Do we write the RB basis functions to disk?
        
Read the "online_mode" flag from the command line
          GetPot command_line (argc, argv);
          int online_mode = 0;
          if ( command_line.search(1, "-online_mode") )
            online_mode = command_line.next(online_mode);
        
Build a mesh on the default MPI communicator.
          Mesh mesh (init.comm(), dim);
          MeshTools::Generation::build_square (mesh,
                                               n_elem, n_elem,
                                               0., 1.,
                                               0., 1.,
                                               QUAD4);
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
We override RBConstruction with SimpleRBConstruction in order to specialize a few functions for this particular problem.
          SimpleRBConstruction & rb_con =
            equation_systems.add_system<SimpleRBConstruction> ("RBConvectionDiffusion");
        
Initialize the data structures for the equation system.
          equation_systems.init ();
        
Print out some information about the "truth" discretization
          equation_systems.print_info();
          mesh.print_info();
        
Build a new RBEvaluation object which will be used to perform Reduced Basis calculations. This is required in both the "Offline" and "Online" stages.
          SimpleRBEvaluation rb_eval(mesh.comm());
        
We need to give the RBConstruction object a pointer to our RBEvaluation object
          rb_con.set_rb_evaluation(rb_eval);
        
          if(!online_mode) // Perform the Offline stage of the RB method
          {
Read in the data that defines this problem from the specified text file
            rb_con.process_parameters_file(parameters_filename);
        
Print out info that describes the current setup of rb_con
            rb_con.print_info();
        
Prepare rb_con for the Construction stage of the RB method. This sets up the necessary data structures and performs initial assembly of the "truth" affine expansion of the PDE.
            rb_con.initialize_rb_construction();
        
Compute the reduced basis space by computing "snapshots", i.e. "truth" solves, at well-chosen parameter values and employing these snapshots as basis functions.
            rb_con.train_reduced_basis();
        
Write out the data that will subsequently be required for the Evaluation stage
            rb_con.get_rb_evaluation().write_offline_data_to_files();
        
If requested, write out the RB basis functions for visualization purposes
            if(store_basis_functions)
            {
Write out the basis functions
              rb_con.get_rb_evaluation().write_out_basis_functions(rb_con);
            }
          }
          else // Perform the Online stage of the RB method
          {
Read in the reduced basis data
            rb_eval.read_offline_data_from_files();
        
Read in online_N and initialize online parameters
            unsigned int online_N = infile("online_N",1);
            Real online_x_vel = infile("online_x_vel", 0.);
            Real online_y_vel = infile("online_y_vel", 0.);
            RBParameters online_mu;
            online_mu.set_value("x_vel", online_x_vel);
            online_mu.set_value("y_vel", online_y_vel);
            rb_eval.set_parameters(online_mu);
            rb_eval.print_parameters();
        
Now do the Online solve using the precomputed reduced basis
            rb_eval.rb_solve(online_N);
        
Print out outputs as well as the corresponding output error bounds.
            std::cout << "output 1, value = " << rb_eval.RB_outputs[0]
                      << ", bound = " << rb_eval.RB_output_error_bounds[0]
                      << std::endl;
            std::cout << "output 2, value = " << rb_eval.RB_outputs[1]
                      << ", bound = " << rb_eval.RB_output_error_bounds[1]
                      << std::endl;
            std::cout << "output 3, value = " << rb_eval.RB_outputs[2]
                      << ", bound = " << rb_eval.RB_output_error_bounds[2]
                      << std::endl;
            std::cout << "output 4, value = " << rb_eval.RB_outputs[3]
                      << ", bound = " << rb_eval.RB_output_error_bounds[3]
                      << std::endl << std::endl;
        
            if(store_basis_functions)
            {
Read in the basis functions
              rb_eval.read_in_basis_functions(rb_con);
        
Plot the solution
              rb_con.load_rb_solution();
        #ifdef LIBMESH_HAVE_EXODUS_API
              ExodusII_IO(mesh).write_equation_systems ("RB_sol.e",equation_systems);
        #endif
        
Plot the first basis function that was generated from the train_reduced_basis call in the Offline stage
              rb_con.load_basis_function(0);
        #ifdef LIBMESH_HAVE_EXODUS_API
              ExodusII_IO(mesh).write_equation_systems ("bf0.e",equation_systems);
        #endif
            }
          }
        
          return 0;
        }



The source file assembly.h without comments:

 
  #ifndef __assembly_h__
  #define __assembly_h__
  
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/fe.h"
  #include "libmesh/fe_interface.h"
  #include "libmesh/fe_base.h"
  #include "libmesh/elem_assembly.h"
  #include "libmesh/quadrature_gauss.h"
  
  #include "libmesh/rb_theta.h"
  #include "libmesh/rb_assembly_expansion.h"
  
  using libMesh::ElemAssembly;
  using libMesh::FEInterface;
  using libMesh::FEMContext;
  using libMesh::Number;
  using libMesh::Point;
  using libMesh::RBAssemblyExpansion;
  using libMesh::RBParameters;
  using libMesh::RBTheta;
  using libMesh::RBThetaExpansion;
  using libMesh::Real;
  using libMesh::RealGradient;
  
  struct ThetaA0 : RBTheta { virtual Number evaluate(const RBParameters& )   { return 0.05;  } };
  struct ThetaA1 : RBTheta { virtual Number evaluate(const RBParameters& mu) { return mu.get_value("x_vel"); } };
  struct ThetaA2 : RBTheta { virtual Number evaluate(const RBParameters& mu) { return mu.get_value("y_vel"); } };
  
  struct A0 : ElemAssembly
  {
    virtual void interior_assembly(FEMContext &c)
    {
      const unsigned int u_var = 0;
  
      const std::vector<Real> &JxW =
        c.element_fe_var[u_var]->get_JxW();
  
      const std::vector<std::vector<RealGradient> >& dphi =
        c.element_fe_var[u_var]->get_dphi();
  
      const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
      unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
      for (unsigned int qp=0; qp != n_qpoints; qp++)
        for (unsigned int i=0; i != n_u_dofs; i++)
          for (unsigned int j=0; j != n_u_dofs; j++)
            c.get_elem_jacobian()(i,j) += JxW[qp] * dphi[j][qp]*dphi[i][qp];
    }
  };
  
  
  struct A1 : ElemAssembly
  {
    virtual void interior_assembly(FEMContext &c)
    {
      const unsigned int u_var = 0;
  
      const std::vector<Real> &JxW =
        c.element_fe_var[u_var]->get_JxW();
  
      const std::vector<std::vector<Real> >& phi =
        c.element_fe_var[u_var]->get_phi();
  
      const std::vector<std::vector<RealGradient> >& dphi =
        c.element_fe_var[u_var]->get_dphi();
  
      const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
      unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
      for (unsigned int qp=0; qp != n_qpoints; qp++)
        for (unsigned int i=0; i != n_u_dofs; i++)
          for (unsigned int j=0; j != n_u_dofs; j++)
            c.get_elem_jacobian()(i,j) += JxW[qp] *dphi[j][qp](0)*phi[i][qp];
    }
  };
  
  struct A2 : ElemAssembly
  {
    virtual void interior_assembly(FEMContext &c)
    {
      const unsigned int u_var = 0;
  
      const std::vector<Real> &JxW =
        c.element_fe_var[u_var]->get_JxW();
  
      const std::vector<std::vector<Real> >& phi =
        c.element_fe_var[u_var]->get_phi();
  
      const std::vector<std::vector<RealGradient> >& dphi =
        c.element_fe_var[u_var]->get_dphi();
  
      const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
      unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
      for (unsigned int qp=0; qp != n_qpoints; qp++)
        for (unsigned int i=0; i != n_u_dofs; i++)
          for (unsigned int j=0; j != n_u_dofs; j++)
            c.get_elem_jacobian()(i,j) += JxW[qp] *dphi[j][qp](1)*phi[i][qp];
    }
  };
  
  struct F0 : ElemAssembly
  {
    virtual void interior_assembly(FEMContext &c)
    {
      const unsigned int u_var = 0;
  
      const std::vector<Real> &JxW =
        c.element_fe_var[u_var]->get_JxW();
  
      const std::vector<std::vector<Real> >& phi =
        c.element_fe_var[u_var]->get_phi();
  
      const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
      unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
      for (unsigned int qp=0; qp != n_qpoints; qp++)
        for (unsigned int i=0; i != n_u_dofs; i++)
          c.get_elem_residual()(i) += JxW[qp] * ( 1.*phi[i][qp] );
    }
  };
  
  struct OutputAssembly : ElemAssembly
  {
    OutputAssembly(Real min_x_in, Real max_x_in,
                   Real min_y_in, Real max_y_in)
                  :
                  min_x(min_x_in),
                  max_x(max_x_in),
                  min_y(min_y_in),
                  max_y(max_y_in)
    {}
  
    virtual void interior_assembly(FEMContext &c)
    {
      const unsigned int u_var = 0;
  
      const std::vector<Real> &JxW =
        c.element_fe_var[u_var]->get_JxW();
  
      const std::vector<std::vector<Real> >& phi =
        c.element_fe_var[u_var]->get_phi();
  
      const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
  
      unsigned int n_qpoints = (c.get_element_qrule())->n_points();
  
      Real output_area = (max_x-min_x) * (max_y-min_y);
  
      Point centroid = c.elem->centroid();
      if( (min_x <= centroid(0)) && (centroid(0) <= max_x) &&
          (min_y <= centroid(1)) && (centroid(1) <= max_y) )
        for (unsigned int qp=0; qp != n_qpoints; qp++)
          for (unsigned int i=0; i != n_u_dofs; i++)
            c.get_elem_residual()(i) += JxW[qp] * ( 1.*phi[i][qp] ) / output_area;
    }
  
    Real min_x, max_x, min_y, max_y;
  };
  
  struct CDRBThetaExpansion : RBThetaExpansion
  {
  
    /**
     * Constructor.
     */
    CDRBThetaExpansion()
    {
      attach_A_theta(&theta_a_0);   // Attach the lhs theta
      attach_A_theta(&theta_a_1);
      attach_A_theta(&theta_a_2);
  
      attach_F_theta(&rb_theta);    // Attach the rhs theta
  
      attach_output_theta(&rb_theta); // Attach output 0 theta
      attach_output_theta(&rb_theta); // Attach output 1 theta
      attach_output_theta(&rb_theta); // Attach output 2 theta
      attach_output_theta(&rb_theta); // Attach output 3 theta
    }
  
    ThetaA0 theta_a_0;
    ThetaA1 theta_a_1;
    ThetaA2 theta_a_2;
    RBTheta rb_theta; // Default RBTheta object, just returns 1.
  };
  
  struct CDRBAssemblyExpansion : RBAssemblyExpansion
  {
  
    /**
     * Constructor.
     */
    CDRBAssemblyExpansion()
      :
      L0(0.7,0.8,0.7,0.8),
      L1(0.2,0.3,0.7,0.8),
      L2(0.2,0.3,0.2,0.3),
      L3(0.7,0.8,0.2,0.3)
    {
      attach_A_assembly(&A0_assembly); // Attach the lhs assembly
      attach_A_assembly(&A1_assembly);
      attach_A_assembly(&A2_assembly);
  
      attach_F_assembly(&F0_assembly); // Attach the rhs assembly
  
      attach_output_assembly(&L0);       // Attach output 0 assembly
      attach_output_assembly(&L1);       // Attach output 1 assembly
      attach_output_assembly(&L2);       // Attach output 2 assembly
      attach_output_assembly(&L3);       // Attach output 3 assembly
    }
  
    A0 A0_assembly;
    A1 A1_assembly;
    A2 A2_assembly;
    F0 F0_assembly;
    OutputAssembly L0;
    OutputAssembly L1;
    OutputAssembly L2;
    OutputAssembly L3;
  };
  
  #endif



The source file rb_classes.h without comments:

 
  #ifndef __rb_classes_h__
  #define __rb_classes_h__
  
  #include "libmesh/rb_construction.h"
  #include "libmesh/fe_base.h"
  
  #include "assembly.h"
  
  using libMesh::AutoPtr;
  using libMesh::DirichletBoundary;
  using libMesh::EquationSystems;
  using libMesh::FEMContext;
  using libMesh::RBConstruction;
  using libMesh::RBEvaluation;
  using libMesh::Real;
  
  
  class SimpleRBEvaluation : public RBEvaluation
  {
  public:
  
    /**
     * Constructor. Just set the theta expansion.
     */
    SimpleRBEvaluation(const Parallel::Communicator& comm)
      : RBEvaluation(comm)
    {
      set_rb_theta_expansion(cd_rb_theta_expansion);
    }
  
    /**
     * The coercivity constant is bounded below by 0.05.
     */
    virtual Real get_stability_lower_bound() { return 0.05; }
  
    /**
     * The object that stores the "theta" expansion of the parameter dependent PDE,
     * i.e. the set of parameter-dependent functions in the affine expansion of the PDE.
     */
    CDRBThetaExpansion cd_rb_theta_expansion;
  
  };
  
  class SimpleRBConstruction : public RBConstruction
  {
  public:
  
    SimpleRBConstruction (EquationSystems& es,
                          const std::string& name_in,
                          const unsigned int number_in)
    : Parent(es, name_in, number_in),
      dirichlet_bc(AutoPtr<DirichletBoundary>(NULL))
    {}
  
    /**
     * Destructor.
     */
    virtual ~SimpleRBConstruction () { }
  
    /**
     * The type of system.
     */
    typedef SimpleRBConstruction sys_type;
  
    /**
     * The type of the parent.
     */
    typedef RBConstruction Parent;
  
    /**
     * Initialize data structures.
     */
    virtual void init_data()
    {
      u_var = this->add_variable ("u", FIRST);
  
      dirichlet_bc = build_zero_dirichlet_boundary_object();
  
      dirichlet_bc->b.insert(0);
      dirichlet_bc->b.insert(1);
      dirichlet_bc->b.insert(2);
      dirichlet_bc->b.insert(3);
      dirichlet_bc->variables.push_back(u_var);
  
      get_dof_map().add_dirichlet_boundary(*dirichlet_bc);
  
      Parent::init_data();
  
      set_rb_assembly_expansion(cd_rb_assembly_expansion);
  
      set_inner_product_assembly(cd_rb_assembly_expansion.A0_assembly);
    }
  
    /**
     * Pre-request all relevant element data.
     */
    virtual void init_context(FEMContext &c)
    {
      c.element_fe_var[u_var]->get_JxW();
      c.element_fe_var[u_var]->get_phi();
      c.element_fe_var[u_var]->get_dphi();
    }
  
    /**
     * Variable number for u.
     */
    unsigned int u_var;
  
    /**
     * The object that stores the "assembly" expansion of the parameter dependent PDE,
     * i.e. the objects that define how to assemble the set of parameter-independent
     * operators in the affine expansion of the PDE.
     */
    CDRBAssemblyExpansion cd_rb_assembly_expansion;
  
    /**
     * The object that defines which degrees of freedom are on a Dirichlet boundary.
     */
    AutoPtr<DirichletBoundary> dirichlet_bc;
  
  };
  
  #endif



The source file reduced_basis_ex1.C without comments:

 
  #include <iostream>
  #include <algorithm>
  #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
  #include <cmath>
  #include <set>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/getpot.h"
  #include "libmesh/elem.h"
  
  #include "rb_classes.h"
  #include "assembly.h"
  
  using namespace libMesh;
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
  #if !defined(LIBMESH_HAVE_XDR)
    libmesh_example_assert(false, "--enable-xdr");
  #elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
    libmesh_example_assert(false, "--disable-singleprecision");
  #endif
    libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    std::string parameters_filename = "reduced_basis_ex1.in";
    GetPot infile(parameters_filename);
  
    unsigned int n_elem = infile("n_elem", 1);       // Determines the number of elements in the "truth" mesh
    const unsigned int dim = 2;                      // The number of spatial dimensions
  
    bool store_basis_functions = infile("store_basis_functions", true); // Do we write the RB basis functions to disk?
  
    GetPot command_line (argc, argv);
    int online_mode = 0;
    if ( command_line.search(1, "-online_mode") )
      online_mode = command_line.next(online_mode);
  
    Mesh mesh (init.comm(), dim);
    MeshTools::Generation::build_square (mesh,
                                         n_elem, n_elem,
                                         0., 1.,
                                         0., 1.,
                                         QUAD4);
  
    EquationSystems equation_systems (mesh);
  
    SimpleRBConstruction & rb_con =
      equation_systems.add_system<SimpleRBConstruction> ("RBConvectionDiffusion");
  
    equation_systems.init ();
  
    equation_systems.print_info();
    mesh.print_info();
  
    SimpleRBEvaluation rb_eval(mesh.comm());
  
    rb_con.set_rb_evaluation(rb_eval);
  
    if(!online_mode) // Perform the Offline stage of the RB method
    {
      rb_con.process_parameters_file(parameters_filename);
  
      rb_con.print_info();
  
      rb_con.initialize_rb_construction();
  
      rb_con.train_reduced_basis();
  
      rb_con.get_rb_evaluation().write_offline_data_to_files();
  
      if(store_basis_functions)
      {
        rb_con.get_rb_evaluation().write_out_basis_functions(rb_con);
      }
    }
    else // Perform the Online stage of the RB method
    {
      rb_eval.read_offline_data_from_files();
  
      unsigned int online_N = infile("online_N",1);
      Real online_x_vel = infile("online_x_vel", 0.);
      Real online_y_vel = infile("online_y_vel", 0.);
      RBParameters online_mu;
      online_mu.set_value("x_vel", online_x_vel);
      online_mu.set_value("y_vel", online_y_vel);
      rb_eval.set_parameters(online_mu);
      rb_eval.print_parameters();
  
      rb_eval.rb_solve(online_N);
  
      std::cout << "output 1, value = " << rb_eval.RB_outputs[0]
                << ", bound = " << rb_eval.RB_output_error_bounds[0]
                << std::endl;
      std::cout << "output 2, value = " << rb_eval.RB_outputs[1]
                << ", bound = " << rb_eval.RB_output_error_bounds[1]
                << std::endl;
      std::cout << "output 3, value = " << rb_eval.RB_outputs[2]
                << ", bound = " << rb_eval.RB_output_error_bounds[2]
                << std::endl;
      std::cout << "output 4, value = " << rb_eval.RB_outputs[3]
                << ", bound = " << rb_eval.RB_output_error_bounds[3]
                << std::endl << std::endl;
  
      if(store_basis_functions)
      {
        rb_eval.read_in_basis_functions(rb_con);
  
        rb_con.load_rb_solution();
  #ifdef LIBMESH_HAVE_EXODUS_API
        ExodusII_IO(mesh).write_equation_systems ("RB_sol.e",equation_systems);
  #endif
  
        rb_con.load_basis_function(0);
  #ifdef LIBMESH_HAVE_EXODUS_API
        ExodusII_IO(mesh).write_equation_systems ("bf0.e",equation_systems);
  #endif
      }
    }
  
    return 0;
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/reduced_basis/reduced_basis_ex1'
***************************************************************
* Running Example reduced_basis_ex1:
*  mpirun -np 4 example-devel -online_mode 0 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
*** Warning, This code is untested, experimental, or likely to see future API changes: ../src/reduced_basis/rb_parametrized.C, line 41, compiled Apr 19 2013 at 11:42:51 ***
 EquationSystems
  n_systems()=1
   System #0, "RBConvectionDiffusion"
    Type "RBConstruction"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=2601
    n_local_dofs()=679
    n_constrained_dofs()=200
    n_local_constrained_dofs()=50
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 8.65013
      Average Off-Processor Bandwidth <= 0.256824
      Maximum  On-Processor Bandwidth <= 11
      Maximum Off-Processor Bandwidth <= 7
    DofMap Constraints
      Number of DoF Constraints = 200
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=2601
    n_local_nodes()=679
  n_elem()=2500
    n_local_elem()=625
    n_active_elem()=2500
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

Initializing training parameters with deterministic training set...
Parameter x_vel: log scaling = 0
Parameter y_vel: log scaling = 0


RBConstruction parameters:
system name: RBConvectionDiffusion
constrained_problem: 0
Nmax: 20
Aq operators attached: 3
Fq functions attached: 1
n_outputs: 4
output 0, Q_l = 1
output 1, Q_l = 1
output 2, Q_l = 1
output 3, Q_l = 1
Number of parameters: 2
Parameter x_vel: Min = -2, Max = 2, value = -2
Parameter y_vel: Min = -2, Max = 2, value = 0.79
n_training_samples: 100
single-matrix mode? 0
reuse preconditioner? 1
use a relative error bound in greedy? 0
write out data during basis training? 0
quiet mode? 1

Compute output dual inner products
output_dual_innerprods[0][0] = 0.323675
output_dual_innerprods[1][0] = 0.323675
output_dual_innerprods[2][0] = 0.323675
output_dual_innerprods[3][0] = 0.323675

---- Performing Greedy basis enrichment ----

---- Basis dimension: 0 ----
Performing RB solves on training set
Maximum (absolute) error bound is 3.74824

Performing truth solve at parameter:
x_vel: -2.000000e+00
y_vel: -2.000000e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 1 ----
Performing RB solves on training set
Maximum (absolute) error bound is 6.67625

Performing truth solve at parameter:
x_vel: 2.000000e+00
y_vel: 2.000000e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 2 ----
Performing RB solves on training set
Maximum (absolute) error bound is 5.18295

Performing truth solve at parameter:
x_vel: -2.000000e+00
y_vel: 2.000000e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 3 ----
Performing RB solves on training set
Maximum (absolute) error bound is 3.59028

Performing truth solve at parameter:
x_vel: 2.000000e+00
y_vel: -2.000000e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 4 ----
Performing RB solves on training set
Maximum (absolute) error bound is 2.71338

Performing truth solve at parameter:
x_vel: -2.222222e-01
y_vel: 2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 5 ----
Performing RB solves on training set
Maximum (absolute) error bound is 2.50784

Performing truth solve at parameter:
x_vel: 2.222222e-01
y_vel: -6.666667e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 6 ----
Performing RB solves on training set
Maximum (absolute) error bound is 2.15653

Performing truth solve at parameter:
x_vel: 6.666667e-01
y_vel: 2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 7 ----
Performing RB solves on training set
Maximum (absolute) error bound is 1.71352

Performing truth solve at parameter:
x_vel: -1.111111e+00
y_vel: -2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 8 ----
Performing RB solves on training set
Maximum (absolute) error bound is 1.61442

Performing truth solve at parameter:
x_vel: 2.222222e-01
y_vel: 1.555556e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 9 ----
Performing RB solves on training set
Maximum (absolute) error bound is 1.09712

Performing truth solve at parameter:
x_vel: -2.222222e-01
y_vel: -2.000000e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 10 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.9794

Performing truth solve at parameter:
x_vel: 2.000000e+00
y_vel: -2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 11 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.763958

Performing truth solve at parameter:
x_vel: -2.222222e-01
y_vel: -2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 12 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.616746

Performing truth solve at parameter:
x_vel: -1.111111e+00
y_vel: 1.111111e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 13 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.530181

Performing truth solve at parameter:
x_vel: -2.000000e+00
y_vel: 2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 14 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.475932

Performing truth solve at parameter:
x_vel: 2.222222e-01
y_vel: 6.666667e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 15 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.440399

Performing truth solve at parameter:
x_vel: -6.666667e-01
y_vel: -1.111111e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 16 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.423537

Performing truth solve at parameter:
x_vel: 1.111111e+00
y_vel: 1.111111e+00

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 17 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.408558

Performing truth solve at parameter:
x_vel: 1.111111e+00
y_vel: -6.666667e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 18 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.333764

Performing truth solve at parameter:
x_vel: 2.222222e-01
y_vel: -2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 19 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.196867

Performing truth solve at parameter:
x_vel: -6.666667e-01
y_vel: 2.222222e-01

Enriching the RB space
Updating RB matrices
Updating RB residual terms

---- Basis dimension: 20 ----
Performing RB solves on training set
Maximum (absolute) error bound is 0.191475

Maximum number of basis functions reached: Nmax = 20

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:54:40 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=2.24983, Active time=1.94564                                                      |
 -------------------------------------------------------------------------------------------------------------------
| 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()        1         0.0017      0.001680    0.0020      0.002013    0.09     0.10     |
|   build_constraint_matrix()           5625      0.0064      0.000001    0.0064      0.000001    0.33     0.33     |
|   build_sparsity()                    1         0.0013      0.001331    0.0060      0.006040    0.07     0.31     |
|   cnstrn_elem_mat_vec()               5625      0.0107      0.000002    0.0107      0.000002    0.55     0.55     |
|   create_dof_constraints()            1         0.0047      0.004745    0.0139      0.013935    0.24     0.72     |
|   distribute_dofs()                   1         0.0038      0.003832    0.0175      0.017514    0.20     0.90     |
|   dof_indices()                       14490     0.0532      0.000004    0.0532      0.000004    2.73     2.73     |
|   prepare_send_list()                 1         0.0000      0.000008    0.0000      0.000008    0.00     0.00     |
|   reinit()                            1         0.0072      0.007192    0.0072      0.007192    0.37     0.37     |
|                                                                                                                   |
| FE                                                                                                                |
|   compute_shape_functions()           12132     0.0413      0.000003    0.0413      0.000003    2.12     2.12     |
|   init_shape_functions()              900       0.0037      0.000004    0.0037      0.000004    0.19     0.19     |
|   inverse_map()                       1764      0.0049      0.000003    0.0049      0.000003    0.25     0.25     |
|                                                                                                                   |
| FEMap                                                                                                             |
|   compute_affine_map()                12132     0.0253      0.000002    0.0253      0.000002    1.30     1.30     |
|   compute_face_map()                  882       0.0039      0.000004    0.0090      0.000010    0.20     0.46     |
|   init_face_shape_functions()         18        0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|   init_reference_to_physical_map()    900       0.0030      0.000003    0.0030      0.000003    0.15     0.15     |
|                                                                                                                   |
| Mesh                                                                                                              |
|   find_neighbors()                    1         0.0037      0.003728    0.0074      0.007436    0.19     0.38     |
|   renumber_nodes_and_elem()           2         0.0004      0.000186    0.0004      0.000186    0.02     0.02     |
|                                                                                                                   |
| MeshCommunication                                                                                                 |
|   assign_global_indices()             1         0.0518      0.051763    0.0527      0.052665    2.66     2.71     |
|   compute_hilbert_indices()           2         0.0145      0.007253    0.0145      0.007253    0.75     0.75     |
|   find_global_indices()               2         0.0021      0.001049    0.0259      0.012948    0.11     1.33     |
|   parallel_sort()                     2         0.0007      0.000350    0.0086      0.004308    0.04     0.44     |
|                                                                                                                   |
| MeshTools::Generation                                                                                             |
|   build_cube()                        1         0.0011      0.001073    0.0011      0.001073    0.06     0.06     |
|                                                                                                                   |
| MetisPartitioner                                                                                                  |
|   partition()                         1         0.0212      0.021219    0.0350      0.034960    1.09     1.80     |
|                                                                                                                   |
| Parallel                                                                                                          |
|   allgather()                         14        0.0047      0.000333    0.0049      0.000349    0.24     0.25     |
|   barrier()                           1         0.0000      0.000019    0.0000      0.000019    0.00     0.00     |
|   broadcast()                         44        0.0003      0.000006    0.0003      0.000006    0.01     0.01     |
|   max(bool)                           5         0.0002      0.000037    0.0002      0.000037    0.01     0.01     |
|   max(scalar)                         119       0.0016      0.000013    0.0016      0.000013    0.08     0.08     |
|   max(vector)                         28        0.0004      0.000013    0.0012      0.000044    0.02     0.06     |
|   maxloc(scalar)                      21        0.0785      0.003737    0.0785      0.003737    4.03     4.03     |
|   min(bool)                           138       0.0015      0.000011    0.0015      0.000011    0.08     0.08     |
|   min(scalar)                         113       0.0362      0.000320    0.0362      0.000320    1.86     1.86     |
|   min(vector)                         28        0.0004      0.000015    0.0020      0.000072    0.02     0.10     |
|   probe()                             80        0.0024      0.000030    0.0024      0.000030    0.12     0.12     |
|   receive()                           64        0.0010      0.000016    0.0034      0.000053    0.05     0.18     |
|   send()                              52        0.0002      0.000003    0.0002      0.000003    0.01     0.01     |
|   send_receive()                      56        0.0002      0.000004    0.0030      0.000053    0.01     0.15     |
|   sum()                               44        0.0083      0.000188    0.0083      0.000190    0.43     0.43     |
|                                                                                                                   |
| Parallel::Request                                                                                                 |
|   wait()                              52        0.0001      0.000001    0.0001      0.000001    0.00     0.00     |
|                                                                                                                   |
| Partitioner                                                                                                       |
|   set_node_processor_ids()            1         0.0012      0.001181    0.0165      0.016520    0.06     0.85     |
|   set_parent_processor_ids()          1         0.0004      0.000403    0.0004      0.000403    0.02     0.02     |
|                                                                                                                   |
| PetscLinearSolver                                                                                                 |
|   solve()                             85        0.7486      0.008807    0.7486      0.008807    38.48    38.48    |
|                                                                                                                   |
| RBConstruction                                                                                                    |
|   add_scaled_matrix_and_vector()      9         0.1711      0.019015    0.3175      0.035280    8.80     16.32    |
|   clear()                             1         0.0005      0.000524    0.0005      0.000524    0.03     0.03     |
|   compute_Fq_representor_innerprods() 1         0.0013      0.001259    0.0101      0.010132    0.06     0.52     |
|   compute_max_error_bound()           21        0.0027      0.000129    0.1915      0.009121    0.14     9.84     |
|   compute_output_dual_innerprods()    1         0.0023      0.002252    0.0336      0.033627    0.12     1.73     |
|   enrich_RB_space()                   20        0.0286      0.001428    0.0286      0.001428    1.47     1.47     |
|   train_reduced_basis()               1         0.0027      0.002682    1.3991      1.399115    0.14     71.91    |
|   truth_assembly()                    20        0.0687      0.003433    0.0687      0.003433    3.53     3.53     |
|   truth_solve()                       20        0.0046      0.000229    0.2296      0.011482    0.24     11.80    |
|   update_RB_system_matrices()         20        0.1106      0.005533    0.1106      0.005533    5.69     5.69     |
|   update_residual_terms()             20        0.2402      0.012010    0.7922      0.039610    12.35    40.72    |
|                                                                                                                   |
| RBEvaluation                                                                                                      |
|   clear()                             1         0.0002      0.000218    0.0002      0.000218    0.01     0.01     |
|   compute_residual_dual_norm()        525       0.0952      0.000181    0.0952      0.000181    4.89     4.89     |
|   rb_solve()                          525       0.0145      0.000028    0.1098      0.000209    0.74     5.64     |
|   resize_data_structures()            1         0.0001      0.000055    0.0001      0.000055    0.00     0.00     |
|   write_offline_data_to_files()       1         0.0160      0.016031    0.0160      0.016031    0.82     0.82     |
|   write_out_basis_functions()         1         0.0003      0.000288    0.0874      0.087351    0.01     4.49     |
|   write_out_vectors()                 1         0.0332      0.033205    0.0871      0.087063    1.71     4.47     |
 -------------------------------------------------------------------------------------------------------------------
| Totals:                               56621     1.9456                                          100.00            |
 -------------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example reduced_basis_ex1:
*  mpirun -np 4 example-devel -online_mode 0 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
***************************************************************
* Running Example reduced_basis_ex1:
*  mpirun -np 4 example-devel -online_mode 1 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
*** Warning, This code is untested, experimental, or likely to see future API changes: ../src/reduced_basis/rb_parametrized.C, line 41, compiled Apr 19 2013 at 11:42:51 ***
 EquationSystems
  n_systems()=1
   System #0, "RBConvectionDiffusion"
    Type "RBConstruction"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=2601
    n_local_dofs()=679
    n_constrained_dofs()=200
    n_local_constrained_dofs()=50
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 8.65013
      Average Off-Processor Bandwidth <= 0.256824
      Maximum  On-Processor Bandwidth <= 11
      Maximum Off-Processor Bandwidth <= 7
    DofMap Constraints
      Number of DoF Constraints = 200
      Average DoF Constraint Length= 0
      Number of Node Constraints = 0

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=2601
    n_local_nodes()=679
  n_elem()=2500
    n_local_elem()=625
    n_active_elem()=2500
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

x_vel: -2.000000e+00
y_vel: 7.900000e-01

output 1, value = 0.119268, bound = 0.0646039
output 2, value = 0.371635, bound = 0.0646039
output 3, value = 0.255409, bound = 0.0646039
output 4, value = 0.115747, bound = 0.0646039


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:54:41 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=0.32433, Active time=0.310435                                                |
 --------------------------------------------------------------------------------------------------------------
| 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()   1         0.0017      0.001723    0.0021      0.002063    0.56     0.66     |
|   build_sparsity()               1         0.0014      0.001421    0.0062      0.006207    0.46     2.00     |
|   create_dof_constraints()       1         0.0039      0.003941    0.0112      0.011208    1.27     3.61     |
|   distribute_dofs()              1         0.0034      0.003444    0.0179      0.017870    1.11     5.76     |
|   dof_indices()                  4490      0.0155      0.000003    0.0155      0.000003    5.00     5.00     |
|   prepare_send_list()            1         0.0000      0.000009    0.0000      0.000009    0.00     0.00     |
|   reinit()                       1         0.0056      0.005615    0.0056      0.005615    1.81     1.81     |
|                                                                                                              |
| EquationSystems                                                                                              |
|   build_solution_vector()        2         0.0033      0.001657    0.0110      0.005502    1.07     3.54     |
|                                                                                                              |
| ExodusII_IO                                                                                                  |
|   write_nodal_data()             2         0.0792      0.039612    0.0792      0.039612    25.52    25.52    |
|                                                                                                              |
| Mesh                                                                                                         |
|   find_neighbors()               1         0.0039      0.003890    0.0072      0.007189    1.25     2.32     |
|   renumber_nodes_and_elem()      2         0.0004      0.000192    0.0004      0.000192    0.12     0.12     |
|                                                                                                              |
| MeshCommunication                                                                                            |
|   assign_global_indices()        1         0.0526      0.052629    0.0532      0.053176    16.95    17.13    |
|   compute_hilbert_indices()      2         0.0127      0.006371    0.0127      0.006371    4.10     4.10     |
|   find_global_indices()          2         0.0017      0.000833    0.0259      0.012961    0.54     8.35     |
|   parallel_sort()                2         0.0006      0.000289    0.0107      0.005331    0.19     3.43     |
|                                                                                                              |
| MeshOutput                                                                                                   |
|   write_equation_systems()       2         0.0001      0.000041    0.0905      0.045241    0.03     29.15    |
|                                                                                                              |
| MeshTools::Generation                                                                                        |
|   build_cube()                   1         0.0011      0.001100    0.0011      0.001100    0.35     0.35     |
|                                                                                                              |
| MetisPartitioner                                                                                             |
|   partition()                    1         0.0200      0.019981    0.0339      0.033854    6.44     10.91    |
|                                                                                                              |
| Parallel                                                                                                     |
|   allgather()                    14        0.0067      0.000478    0.0069      0.000493    2.15     2.22     |
|   barrier()                      1         0.0000      0.000020    0.0000      0.000020    0.01     0.01     |
|   broadcast()                    13        0.0001      0.000005    0.0001      0.000004    0.02     0.02     |
|   max(bool)                      1         0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|   max(scalar)                    135       0.0009      0.000007    0.0009      0.000007    0.29     0.29     |
|   max(vector)                    30        0.0002      0.000006    0.0006      0.000019    0.06     0.18     |
|   min(bool)                      156       0.0006      0.000004    0.0006      0.000004    0.20     0.20     |
|   min(scalar)                    129       0.0510      0.000395    0.0510      0.000395    16.43    16.43    |
|   min(vector)                    30        0.0003      0.000008    0.0016      0.000054    0.08     0.52     |
|   probe()                        64        0.0025      0.000038    0.0025      0.000038    0.79     0.79     |
|   receive()                      58        0.0002      0.000004    0.0027      0.000046    0.08     0.86     |
|   send()                         58        0.0002      0.000003    0.0002      0.000003    0.06     0.06     |
|   send_receive()                 56        0.0003      0.000005    0.0031      0.000055    0.09     0.99     |
|   sum()                          31        0.0105      0.000337    0.0221      0.000713    3.37     7.12     |
|                                                                                                              |
| Parallel::Request                                                                                            |
|   wait()                         58        0.0005      0.000009    0.0005      0.000009    0.16     0.16     |
|                                                                                                              |
| Partitioner                                                                                                  |
|   set_node_processor_ids()       1         0.0011      0.001108    0.0181      0.018054    0.36     5.82     |
|   set_parent_processor_ids()     1         0.0004      0.000361    0.0004      0.000361    0.12     0.12     |
|                                                                                                              |
| RBConstruction                                                                                               |
|   clear()                        1         0.0002      0.000211    0.0002      0.000211    0.07     0.07     |
|   load_basis_function()          1         0.0001      0.000122    0.0001      0.000122    0.04     0.04     |
|   load_rb_solution()             1         0.0002      0.000170    0.0002      0.000170    0.05     0.05     |
|                                                                                                              |
| RBEvaluation                                                                                                 |
|   clear()                        1         0.0001      0.000066    0.0001      0.000066    0.02     0.02     |
|   compute_residual_dual_norm()   1         0.0008      0.000787    0.0008      0.000787    0.25     0.25     |
|   rb_solve()                     1         0.0084      0.008378    0.0092      0.009165    2.70     2.95     |
|   read_in_basis_functions()      1         0.0000      0.000018    0.0822      0.082244    0.01     26.49    |
|   read_in_vectors()              1         0.0173      0.017349    0.0822      0.082225    5.59     26.49    |
|   read_offline_data_from_files() 1         0.0007      0.000662    0.0007      0.000735    0.21     0.24     |
|   resize_data_structures()       1         0.0001      0.000073    0.0001      0.000073    0.02     0.02     |
 --------------------------------------------------------------------------------------------------------------
| Totals:                          5360      0.3104                                          100.00            |
 --------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example reduced_basis_ex1:
*  mpirun -np 4 example-devel -online_mode 1 -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/reduced_basis/reduced_basis_ex1'

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

Hosted By:
SourceForge.net Logo