adjoints_ex6.C File Reference

Go to the source code of this file.

Functions
void	write_output (EquationSystems &es, unsigned int a_step, std::string solution_type)
void	set_system_parameters (PoissonSystem &system, FEMParameters &param)
std::unique_ptr< MeshRefinement >	build_mesh_refinement (MeshBase &mesh, FEMParameters &param)
std::unique_ptr< AdjointRefinementEstimator >	build_adjoint_refinement_error_estimator (QoISet &qois)
int	main (int argc, char **argv)

Function Documentation

build_adjoint_refinement_error_estimator()

std::unique_ptr<AdjointRefinementEstimator> build_adjoint_refinement_error_estimator

(

QoISet &

qois

)

Definition at line 190 of file adjoints_ex6.C.

{
  libMesh::out << "Computing the error estimate using the Adjoint Refinement Error Estimator\n" << std::endl;
 
  AdjointRefinementEstimator * adjoint_refinement_estimator = new AdjointRefinementEstimator;
 
  adjoint_refinement_estimator->qoi_set() = qois;
 
  // We enrich the FE space for the dual problem by doing 2 uniform h refinements
  adjoint_refinement_estimator->number_h_refinements = 2;
 
  return std::unique_ptr<AdjointRefinementEstimator>(adjoint_refinement_estimator);
}

References libMesh::AdjointRefinementEstimator::number_h_refinements, libMesh::out, and libMesh::AdjointRefinementEstimator::qoi_set().

Referenced by main().

build_mesh_refinement()

std::unique_ptr<MeshRefinement> build_mesh_refinement	(	MeshBase &	mesh,
		FEMParameters &	param
	)

Definition at line 170 of file adjoints_ex6.C.

{
  MeshRefinement * mesh_refinement = new MeshRefinement(mesh);
  mesh_refinement->coarsen_by_parents() = true;
  mesh_refinement->absolute_global_tolerance() = param.global_tolerance;
  mesh_refinement->nelem_target()      = param.nelem_target;
  mesh_refinement->refine_fraction()   = param.refine_fraction;
  mesh_refinement->coarsen_fraction()  = param.coarsen_fraction;
  mesh_refinement->coarsen_threshold() = param.coarsen_threshold;
 
  return std::unique_ptr<MeshRefinement>(mesh_refinement);
}

References libMesh::MeshRefinement::absolute_global_tolerance(), libMesh::MeshRefinement::coarsen_by_parents(), FEMParameters::coarsen_fraction, libMesh::MeshRefinement::coarsen_fraction(), FEMParameters::coarsen_threshold, libMesh::MeshRefinement::coarsen_threshold(), FEMParameters::global_tolerance, mesh, FEMParameters::nelem_target, libMesh::MeshRefinement::nelem_target(), FEMParameters::refine_fraction, and libMesh::MeshRefinement::refine_fraction().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 208 of file adjoints_ex6.C.

{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);
 
  // Skip adaptive examples on a non-adaptive libMesh build
#ifndef LIBMESH_ENABLE_AMR
  libmesh_example_requires(false, "--enable-amr");
#else
 
  // We use Dirichlet boundary conditions here
#ifndef LIBMESH_ENABLE_DIRICHLET
  libmesh_example_requires(false, "--enable-dirichlet");
#endif
 
  // This doesn't converge with Eigen BICGSTAB for some reason...
  libmesh_example_requires((libMesh::default_solver_package() != EIGEN_SOLVERS) &&
                           (libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE),
                           "--enable-petsc or --enable-trilinos");
 
  libMesh::out << "Started " << argv[0] << std::endl;
 
  // Make sure the general input file exists, and parse it
  {
    std::ifstream i("general.in");
    if (!i)
      libmesh_error_msg('[' << init.comm().rank() << "] Can't find general.in; exiting early.");
  }
  GetPot infile("general.in");
 
  // Read in parameters from the input file
  FEMParameters param(init.comm());
  param.read(infile);
 
  // Skip this default-2D example if libMesh was compiled as 1D-only.
  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
 
  // Create a mesh, with dimension to be overridden later, distributed
  // across the default MPI communicator.
  Mesh mesh(init.comm());
 
  // And an object to refine it
  std::unique_ptr<MeshRefinement> mesh_refinement =
    build_mesh_refinement(mesh, param);
 
  // And an EquationSystems to run on it
  EquationSystems equation_systems (mesh);
 
  libMesh::out << "Reading in and building the mesh" << std::endl;
 
  // Read in the mesh
  MeshTools::Generation::build_square
    (mesh, 2, 2,
     0.0, 1.0,
     0.0, 1.0,
     QUAD4);
 
  // Make all the elements of the mesh second order so we can compute
  // with a higher order basis
  mesh.all_second_order();
 
  // Create a mesh refinement object to do the initial uniform refinements
  // on the coarse grid read in from lshaped.xda
  MeshRefinement initial_uniform_refinements(mesh);
  initial_uniform_refinements.uniformly_refine(param.coarserefinements);
 
  libMesh::out << "Building system" << std::endl;
 
  // Build the FEMSystem
  PoissonSystem & system = equation_systems.add_system<PoissonSystem> ("PoissonSystem");
 
  // Set its parameters
  set_system_parameters(system, param);
 
  libMesh::out << "Initializing systems" << std::endl;
 
  equation_systems.init ();
 
  // Add an adjoint_solution0 vector to the system
  system.add_vector("adjoint_solution0", false, GHOSTED);
 
  // Print information about the mesh and system to the screen.
  mesh.print_info();
  equation_systems.print_info();
 
  // Get a pointer to the linear solver object to be able to reuse preconditioner
  LinearSolver<Number> * linear_solver = system.get_linear_solver();
 
  {
    // Adaptively solve the timestep
    unsigned int a_step = 0;
    for (; a_step != param.max_adaptivesteps; ++a_step)
      {
        // We can't adapt to both a tolerance and a
        // target mesh size
        if (param.global_tolerance != 0.)
          libmesh_assert_equal_to (param.nelem_target, 0);
        // If we aren't adapting to a tolerance we need a
        // target mesh size
        else
          libmesh_assert_greater (param.nelem_target, 0);
 
        // Dont reuse preconditioners before the primal solve
        linear_solver->reuse_preconditioner(false);
 
        // Solve the forward problem
        system.solve();
 
        // Write out the computed primal solution
        write_output(equation_systems, a_step, "primal");
 
        // Declare a QoISet object, we need this object to set weights for our QoI error contributions
        QoISet qois;
 
        // Declare a qoi_indices vector, each index will correspond to a QoI
        std::vector<unsigned int> qoi_indices;
        qoi_indices.push_back(0);
        qois.add_indices(qoi_indices);
 
        // Set weights for each index, these will weight the contribution of each QoI in the final error
        // estimate to be used for flagging elements for refinement
        qois.set_weight(0, 1.0);
 
        // Make sure we get the contributions to the adjoint RHS from the sides
        system.assemble_qoi_sides = true;
 
        // We are about to solve the adjoint system, but before we do this we see the same preconditioner
        // flag to reuse the preconditioner from the forward solver
        linear_solver->reuse_preconditioner(param.reuse_preconditioner);
 
        // Solve the adjoint system. This takes the transpose of the stiffness matrix and then
        // solves the resulting system
        system.adjoint_solve();
 
        //Now that we have solved the adjoint, set the adjoint_already_solved boolean to true, so we dont solve unnecessarily in the error estimator
        system.set_adjoint_already_solved(true);
 
        // 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, a_step, "adjoint_0");
 
        //Swap back
        primal_solution.swap(dual_solution_0);
 
        libMesh::out << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
                     << " active elements and "
                     << equation_systems.n_active_dofs()
                     << " active dofs." << std::endl ;
 
        // Postprocess, compute the approximate QoIs and write them out to the console
        libMesh::out << "Postprocessing: " << std::endl;
        system.postprocess_sides = true;
        system.postprocess();
 
        Number QoI_0_computed = system.get_QoI_value("computed", 0);
        Number QoI_0_exact = system.get_QoI_value("exact", 0);
 
        libMesh::out << "The computed QoI 0 is " << std::setprecision(17)
                     << QoI_0_computed << std::endl;
        libMesh::out << "The relative error in QoI 0 is " << std::setprecision(17)
                     << std::abs(QoI_0_computed - QoI_0_exact) << std::endl; // / std::abs(QoI_0_exact)
 
        // We will declare an error vector for passing to the adjoint refinement error estimator
        ErrorVector QoI_elementwise_error;
 
        // Build an adjoint refinement error estimator object
        std::unique_ptr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
          build_adjoint_refinement_error_estimator(qois);
 
        // Estimate the error in each element using the Adjoint Refinement estimator
        adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
 
        // Print out the computed error estimate, note that we access the global error estimates
        // using an accessor function, right now sum(QoI_elementwise_error) != global_QoI_error_estimate
        libMesh::out << "The computed relative error in QoI 0 is " << std::setprecision(17)
                     << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) << std::endl; // / std::abs(QoI_0_exact)
 
        // Also print out effectivity indices (estimated error/true error)
        libMesh::out << "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
                     << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
          std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
 
        // For refinement purposes we need to sort by error
        // *magnitudes*, but AdjointRefinement gives us signed errors.
        if (!param.refine_uniformly)
          for (std::size_t i=0; i<QoI_elementwise_error.size(); i++)
            if (QoI_elementwise_error[i] != 0.)
              QoI_elementwise_error[i] = std::abs(QoI_elementwise_error[i]);
 
        // We have to refine either based on reaching an error tolerance or
        // a number of elements target, which should be verified above
        // Otherwise we flag elements by error tolerance or nelem target
 
        // Uniform refinement
        if (param.refine_uniformly)
          {
            mesh_refinement->uniformly_refine(1);
          }
        // Adaptively refine based on reaching an error tolerance
        else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
          {
            mesh_refinement->flag_elements_by_error_tolerance (QoI_elementwise_error);
 
            mesh_refinement->refine_and_coarsen_elements();
          }
        // Adaptively refine based on reaching a target number of elements
        else
          {
            if (mesh.n_active_elem() >= param.nelem_target)
              {
                libMesh::out << "We reached the target number of elements.\n" << std::endl;
                break;
              }
 
            mesh_refinement->flag_elements_by_nelem_target (QoI_elementwise_error);
 
            mesh_refinement->refine_and_coarsen_elements();
          }
 
        // Dont forget to reinit the system after each adaptive refinement!
        equation_systems.reinit();
 
        libMesh::out << "Refined mesh to "
                     << mesh.n_active_elem()
                     << " active elements and "
                     << equation_systems.n_active_dofs()
                     << " active dofs." << std::endl;
      }
 
    // On the last adaptive step, dont refine elements and check regressions via asserts
    if (a_step == param.max_adaptivesteps)
      {
        linear_solver->reuse_preconditioner(false);
        system.solve();
 
        write_output(equation_systems, a_step, "primal");
 
        NumericVector<Number> & primal_solution = *system.solution;
 
        QoISet qois;
        std::vector<unsigned int> qoi_indices;
 
        qoi_indices.push_back(0);
        qois.add_indices(qoi_indices);
 
        qois.set_weight(0, 1.0);
 
        system.assemble_qoi_sides = true;
 
        linear_solver->reuse_preconditioner(param.reuse_preconditioner);
        system.adjoint_solve();
        system.set_adjoint_already_solved(true);
 
        NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
 
        primal_solution.swap(dual_solution_0);
        write_output(equation_systems, a_step, "adjoint_0");
 
        primal_solution.swap(dual_solution_0);
 
        libMesh::out << "Adaptive step " << a_step << ", we have " << mesh.n_active_elem()
                     << " active elements and "
                     << equation_systems.n_active_dofs()
                     << " active dofs." << std::endl ;
 
        libMesh::out << "Postprocessing: " << std::endl;
        system.postprocess_sides = true;
        system.postprocess();
 
        Number QoI_0_computed = system.get_QoI_value("computed", 0);
        Number QoI_0_exact = system.get_QoI_value("exact", 0);
 
        libMesh::out << "The computed QoI 0 is " << std::setprecision(17)
                     << QoI_0_computed << std::endl;
        libMesh::out << "The relative error in QoI 0 is " << std::setprecision(17)
                     << std::abs(QoI_0_computed - QoI_0_exact) << std::endl; // / std::abs(QoI_0_exact)
 
 
        // We will declare an error vector for passing to the adjoint
        // refinement error estimator Right now, only the first entry
        // of this vector will be filled (with the global QoI error
        // estimate) Later, each entry of the vector will contain
        // elementwise error that the user can sum to get the total
        // error
        ErrorVector QoI_elementwise_error;
 
        // Build an adjoint refinement error estimator object
        std::unique_ptr<AdjointRefinementEstimator> adjoint_refinement_error_estimator =
          build_adjoint_refinement_error_estimator(qois);
 
        // Estimate the error in each element using the Adjoint Refinement estimator
        adjoint_refinement_error_estimator->estimate_error(system, QoI_elementwise_error);
 
        // Print out the computed error estimate, note that we access
        // the global error estimates using an accessor function,
        // right now sum(QoI_elementwise_error) != global_QoI_error_estimate
        libMesh::out << "The computed relative error in QoI 0 is " << std::setprecision(17)
                     << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) << std::endl; // / std::abs(QoI_0_exact)
 
        // Also print out effectivity indices (estimated error/true error)
        libMesh::out << "The effectivity index for the computed error in QoI 0 is " << std::setprecision(17)
                     << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) /
          std::abs(QoI_0_computed - QoI_0_exact) << std::endl;
 
        // Hard coded assert to ensure that the actual numbers we are getting are what they should be
 
        // The effectivity index isn't exactly reproducible at single precision
        // libmesh_assert_less(std::abs(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact) - 0.84010976704434637), 1.e-5);
        // libmesh_assert_less(std::abs(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact) - 0.48294428289950514), 1.e-5);
 
        // But the effectivity indices should always be sane
        // libmesh_assert_less(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact), 2.5);
        // libmesh_assert_greater(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(0)) / std::abs(QoI_0_computed - QoI_0_exact), .4);
        // libmesh_assert_less(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact), 2.5);
        // libmesh_assert_greater(std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact), .4);
 
        // And the computed errors should still be low
        // libmesh_assert_less(std::abs(QoI_0_computed - QoI_0_exact), 2e-4);
        // libmesh_assert_less(std::abs(QoI_1_computed - QoI_1_exact), 2e-4);
      }
  }
 
  libMesh::err << '[' << mesh.processor_id()
               << "] Completing output." << std::endl;
 
#endif // #ifndef LIBMESH_ENABLE_AMR
 
  // All done.
  return 0;
}

References std::abs(), libMesh::QoISet::add_indices(), libMesh::EquationSystems::add_system(), libMesh::System::add_vector(), libMesh::DifferentiableSystem::adjoint_solve(), libMesh::MeshBase::all_second_order(), libMesh::DifferentiableQoI::assemble_qoi_sides, build_adjoint_refinement_error_estimator(), build_mesh_refinement(), libMesh::MeshTools::Generation::build_square(), libMesh::default_solver_package(), libMesh::EIGEN_SOLVERS, libMesh::err, libMesh::System::get_adjoint_solution(), libMesh::DifferentiableSystem::get_linear_solver(), PoissonSystem::get_QoI_value(), libMesh::GHOSTED, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::EquationSystems::n_active_dofs(), libMesh::MeshBase::n_active_elem(), libMesh::out, PoissonSystem::postprocess(), libMesh::DifferentiableSystem::postprocess_sides, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::ParallelObject::processor_id(), libMesh::QUAD4, FEMParameters::read(), libMesh::EquationSystems::reinit(), libMesh::LinearSolver< T >::reuse_preconditioner(), libMesh::System::set_adjoint_already_solved(), set_system_parameters(), libMesh::QoISet::set_weight(), libMesh::System::solution, libMesh::FEMSystem::solve(), libMesh::NumericVector< T >::swap(), libMesh::MeshRefinement::uniformly_refine(), and write_output().

set_system_parameters()

void set_system_parameters	(	PoissonSystem &	system,
		FEMParameters &	param
	)

Definition at line 118 of file adjoints_ex6.C.

{
  // Use analytical jacobians?
  system.analytic_jacobians() = param.analytic_jacobians;
 
  // Verify analytic jacobians against numerical ones?
  system.verify_analytic_jacobians = param.verify_analytic_jacobians;
 
  // Use the prescribed FE type
  system.fe_family() = param.fe_family[0];
  system.fe_order() = param.fe_order[0];
 
  // More desperate debugging options
  system.print_solution_norms = param.print_solution_norms;
  system.print_solutions      = param.print_solutions;
  system.print_residual_norms = param.print_residual_norms;
  system.print_residuals      = param.print_residuals;
  system.print_jacobian_norms = param.print_jacobian_norms;
  system.print_jacobians      = param.print_jacobians;
 
  // No transient time solver
  system.time_solver = libmesh_make_unique<SteadySolver>(system);
 
  // Nonlinear solver options
  {
    NewtonSolver * solver = new NewtonSolver(system);
    system.time_solver->diff_solver() = std::unique_ptr<DiffSolver>(solver);
 
    solver->quiet                       = param.solver_quiet;
    solver->max_nonlinear_iterations    = param.max_nonlinear_iterations;
    solver->minsteplength               = param.min_step_length;
    solver->relative_step_tolerance     = param.relative_step_tolerance;
    solver->relative_residual_tolerance = param.relative_residual_tolerance;
    solver->require_residual_reduction  = param.require_residual_reduction;
    solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
    if (system.time_solver->reduce_deltat_on_diffsolver_failure)
      {
        solver->continue_after_max_iterations = true;
        solver->continue_after_backtrack_failure = true;
      }
 
    // And the linear solver options
    solver->max_linear_iterations       = param.max_linear_iterations;
    solver->initial_linear_tolerance    = param.initial_linear_tolerance;
    solver->minimum_linear_tolerance    = param.minimum_linear_tolerance;
  }
}

References PoissonSystem::analytic_jacobians(), FEMParameters::analytic_jacobians, libMesh::DiffSolver::continue_after_backtrack_failure, libMesh::DiffSolver::continue_after_max_iterations, PoissonSystem::fe_family(), FEMParameters::fe_family, PoissonSystem::fe_order(), FEMParameters::fe_order, FEMParameters::initial_linear_tolerance, libMesh::DiffSolver::initial_linear_tolerance, FEMParameters::linear_tolerance_multiplier, libMesh::NewtonSolver::linear_tolerance_multiplier, FEMParameters::max_linear_iterations, libMesh::DiffSolver::max_linear_iterations, FEMParameters::max_nonlinear_iterations, libMesh::DiffSolver::max_nonlinear_iterations, FEMParameters::min_step_length, FEMParameters::minimum_linear_tolerance, libMesh::DiffSolver::minimum_linear_tolerance, libMesh::NewtonSolver::minsteplength, FEMParameters::print_jacobian_norms, libMesh::DifferentiableSystem::print_jacobian_norms, FEMParameters::print_jacobians, libMesh::DifferentiableSystem::print_jacobians, FEMParameters::print_residual_norms, libMesh::DifferentiableSystem::print_residual_norms, FEMParameters::print_residuals, libMesh::DifferentiableSystem::print_residuals, FEMParameters::print_solution_norms, libMesh::DifferentiableSystem::print_solution_norms, FEMParameters::print_solutions, libMesh::DifferentiableSystem::print_solutions, libMesh::DiffSolver::quiet, FEMParameters::relative_residual_tolerance, libMesh::DiffSolver::relative_residual_tolerance, FEMParameters::relative_step_tolerance, libMesh::DiffSolver::relative_step_tolerance, libMesh::NewtonSolver::require_residual_reduction, FEMParameters::require_residual_reduction, FEMParameters::solver_quiet, libMesh::DifferentiableSystem::time_solver, FEMParameters::verify_analytic_jacobians, and libMesh::FEMSystem::verify_analytic_jacobians.

Referenced by main().

write_output()

void write_output	(	EquationSystems &	es,
		unsigned int	a_step,
		std::string	solution_type
	)

Definition at line 97 of file adjoints_ex6.C.

{
#ifdef LIBMESH_HAVE_GMV
  MeshBase & mesh = es.get_mesh();
 
  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
}

References libMesh::EquationSystems::get_mesh(), mesh, and libMesh::MeshOutput< MT >::write_equation_systems().

Referenced by main().