reduced_basis_ex6.C File Reference

Go to the source code of this file.

Functions
void	transform_mesh_and_plot (EquationSystems &es, Real curvature, const std::string &filename)
int	main (int argc, char **argv)

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 111 of file reduced_basis_ex6.C.

References libMesh::EquationSystems::add_system(), libMesh::MeshTools::Generation::build_cube(), libMesh::ParallelObject::comm(), libMesh::command_line_next(), libMesh::default_solver_package(), libMesh::RBEIMConstruction::get_eim_assembly_objects(), libMesh::HEX8, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::RBEIMConstruction::initialize_eim_assembly_objects(), libMesh::RBEIMConstruction::initialize_eim_construction(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::RBConstruction::load_rb_solution(), mesh, libMesh::out, libMesh::RBEIMConstruction::print_info(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::RBEIMConstruction::process_parameters_file(), libMesh::RBDataDeserialization::RBEvaluationDeserialization::read_from_file(), libMesh::RBDataDeserialization::RBEIMEvaluationDeserialization::read_from_file(), libMesh::Real, libMesh::RBEIMConstruction::set_rb_eim_evaluation(), libMesh::RBConstruction::set_rb_evaluation(), libMesh::RBParameters::set_value(), libMesh::RBEIMConstruction::train_eim_approximation(), transform_mesh_and_plot(), libMesh::RBDataSerialization::RBEvaluationSerialization::write_to_file(), and libMesh::RBDataSerialization::RBEIMEvaluationSerialization::write_to_file().

{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");

#if !defined(LIBMESH_HAVE_XDR)
  // We need XDR support to write out reduced bases
  libmesh_example_requires(false, "--enable-xdr");
#elif defined(LIBMESH_DEFAULT_SINGLE_PRECISION)
  // XDR binary support requires double precision
  libmesh_example_requires(false, "--disable-singleprecision");
#elif defined(LIBMESH_DEFAULT_TRIPLE_PRECISION)
  // I have no idea why long double isn't working here... [RHS]
  libmesh_example_requires(false, "double precision");
#elif defined(LIBMESH_ENABLE_BLOCKED_STORAGE)
  // This example dies with "New nonzero at (0,2) caused a malloc"
  // when blocked storage is enabled.
  libmesh_example_requires(false, "--disable-blocked-storage");
#endif

  // This is a 3D example
  libmesh_example_requires(3 == LIBMESH_DIM, "3D support");

  // This example requires libmesh to be configured with both
  // DirichletBoundary and Cap'n Proto support.
#if !defined(LIBMESH_ENABLE_DIRICHLET) || !defined(LIBMESH_HAVE_CAPNPROTO)
  libmesh_example_requires(false, "--enable-dirichlet --enable-capnp");
#else

  // Parse the input file using GetPot
  std::string eim_parameters = "eim.in";
  std::string rb_parameters  = "rb.in";
  std::string main_parameters = "reduced_basis_ex6.in";
  GetPot infile(main_parameters);

  unsigned int n_elem_xy = infile("n_elem_xy", 1);
  unsigned int n_elem_z  = infile("n_elem_z", 1);

  // Read the "online_mode" flag from the command line
  const int online_mode =
    libMesh::command_line_next("-online_mode", 0);

  // Create a mesh, with dimension to be overridden by build_cube, on
  // the default MPI communicator.  We currently have to create a
  // ReplicatedMesh here due to a reduced_basis regression with
  // DistributedMesh
  ReplicatedMesh mesh(init.comm());

  MeshTools::Generation::build_cube (mesh,
                                     n_elem_xy, n_elem_xy, n_elem_z,
                                     -0.2, 0.2,
                                     -0.2, 0.2,
                                     0., 3.,
                                     HEX8);

  if (!online_mode)
    {
      // First run the Offline stage for the EIM approximation.
      {
        libMesh::out << "Perform EIM training" << std::endl << std::endl;

        // Initialize the EquationSystems object for this mesh and attach
        // the EIM and RB Construction objects
        EquationSystems equation_systems (mesh);

        SimpleEIMConstruction & eim_construction =
          equation_systems.add_system<SimpleEIMConstruction> ("EIM");

        // Initialize the data structures for the equation system.
        equation_systems.init ();

        // Print out some information about the "truth" discretization
        mesh.print_info();
        equation_systems.print_info();

        // Initialize the EIM RBEvaluation object
        SimpleEIMEvaluation eim_rb_eval(mesh.comm());

        // Set the rb_eval objects for the RBConstructions
        eim_construction.set_rb_eim_evaluation(eim_rb_eval);

        // Read data from input file and print state
        eim_construction.process_parameters_file(eim_parameters);
        eim_construction.print_info();

        // Perform the EIM Greedy and write out the data
        eim_construction.initialize_eim_construction();
        eim_construction.train_eim_approximation();

        RBDataSerialization::RBEIMEvaluationSerialization rb_eim_eval_writer(eim_rb_eval);
        rb_eim_eval_writer.write_to_file("rb_eim_eval.bin");

        // Write out the EIM basis functions
        eim_rb_eval.write_out_basis_functions("eim_data", /*binary=*/false);
      }

      {
        libMesh::out << std::endl << "Perform RB training" << std::endl << std::endl;

        // Create an equation systems object.
        EquationSystems equation_systems (mesh);

        SimpleEIMConstruction & eim_construction =
          equation_systems.add_system<SimpleEIMConstruction> ("EIM");
        SimpleRBConstruction & rb_construction =
          equation_systems.add_system<SimpleRBConstruction> ("RB");

        // 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();

        // Initialize the standard RBEvaluation object
        SimpleRBEvaluation rb_eval(mesh.comm());

        // Initialize the EIM RBEvaluation object
        SimpleEIMEvaluation eim_rb_eval(mesh.comm());

        // Set the rb_eval objects for the RBConstructions
        eim_construction.set_rb_eim_evaluation(eim_rb_eval);
        rb_construction.set_rb_evaluation(rb_eval);

        RBDataDeserialization::RBEIMEvaluationDeserialization rb_eim_eval_reader(eim_rb_eval);
        rb_eim_eval_reader.read_from_file("rb_eim_eval.bin");
        eim_rb_eval.read_in_basis_functions(rb_construction, "eim_data", /*binary=*/false);

        // Read data from input file and print state
        rb_construction.process_parameters_file(rb_parameters);

        // attach the EIM theta objects to the RBEvaluation
        eim_rb_eval.initialize_eim_theta_objects();
        rb_eval.get_rb_theta_expansion().attach_multiple_A_theta(eim_rb_eval.get_eim_theta_objects());

        // attach the EIM assembly objects to the RBConstruction
        eim_construction.initialize_eim_assembly_objects();
        rb_construction.get_rb_assembly_expansion().attach_multiple_A_assembly(eim_construction.get_eim_assembly_objects());

        // Print out the state of rb_construction now that the EIM objects have been attached
        rb_construction.print_info();

        // Need to initialize _after_ EIM greedy so that
        // the system knows how many affine terms there are
        rb_construction.initialize_rb_construction();
        rb_construction.train_reduced_basis();

        RBDataSerialization::RBEvaluationSerialization rb_eval_writer(rb_construction.get_rb_evaluation());
        rb_eval_writer.write_to_file("rb_eval.bin");

        // Write out the basis functions
        rb_construction.get_rb_evaluation().write_out_basis_functions(rb_construction,
                                                                      "rb_data");
      }
    }
  else
    {
      SimpleEIMEvaluation eim_rb_eval(mesh.comm());
      SimpleRBEvaluation rb_eval(mesh.comm());

      RBDataDeserialization::RBEIMEvaluationDeserialization rb_eim_eval_reader(eim_rb_eval);
      rb_eim_eval_reader.read_from_file("rb_eim_eval.bin");

      // attach the EIM theta objects to rb_eval objects
      eim_rb_eval.initialize_eim_theta_objects();
      rb_eval.get_rb_theta_expansion().attach_multiple_A_theta(eim_rb_eval.get_eim_theta_objects());

      RBDataDeserialization::RBEvaluationDeserialization rb_eval_reader(rb_eval);
      rb_eval_reader.read_from_file("rb_eval.bin", /*read_error_bound_data*/ true);

      // Get the parameters at which we will do a reduced basis solve
      Real online_curvature = infile("online_curvature", 0.);
      Real online_Bi        = infile("online_Bi", 0.);
      Real online_kappa     = infile("online_kappa", 0.);
      RBParameters online_mu;
      online_mu.set_value("curvature", online_curvature);
      online_mu.set_value("Bi", online_Bi);
      online_mu.set_value("kappa", online_kappa);
      rb_eval.set_parameters(online_mu);
      rb_eval.print_parameters();
      rb_eval.rb_solve(rb_eval.get_n_basis_functions());

      EquationSystems equation_systems (mesh);

      SimpleRBConstruction & rb_construction =
        equation_systems.add_system<SimpleRBConstruction> ("RB");

      equation_systems.init ();

      rb_construction.set_rb_evaluation(rb_eval);

      // read in the data from files
      rb_eval.read_in_basis_functions(rb_construction, "rb_data");
      rb_construction.load_rb_solution();

      transform_mesh_and_plot(equation_systems, online_curvature, "RB_sol.e");
    }

#endif // LIBMESH_ENABLE_DIRICHLET

  return 0;
}

transform_mesh_and_plot()

void transform_mesh_and_plot	(	EquationSystems &	es,
		Real	curvature,
		const std::string &	filename
	)

Definition at line 320 of file reduced_basis_ex6.C.

References libMesh::EquationSystems::get_mesh(), libMesh::libmesh_ignore(), mesh, libMesh::Real, and libMesh::ExodusII_IO::write_equation_systems().

Referenced by main().

{
  // Loop over the mesh nodes and move them!
  MeshBase & mesh = es.get_mesh();

  for (auto & node : mesh.node_ptr_range())
    {
      Real x = (*node)(0);
      Real z = (*node)(2);

      (*node)(0) = -1./curvature + (1./curvature + x)*cos(curvature*z);
      (*node)(2) = (1./curvature + x)*sin(curvature*z);
    }

#ifdef LIBMESH_HAVE_EXODUS_API
  ExodusII_IO(mesh).write_equation_systems(filename, es);
#endif

  // Avoid unused variable warnings in the --disable-exodus case
  libmesh_ignore(filename, es);
}