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

void	set_system_parameters (LaplaceSystem &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 236 of file adjoints_ex4.C.

References libMesh::out.

Referenced by main().

 {
   libMesh::out << "Computing the error estimate using the Adjoint Refinement Error Estimator" << std::endl << std::endl;
 
   auto adjoint_refinement_estimator = std::make_unique<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 adjoint_refinement_estimator;
 }

◆ build_mesh_refinement()

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

Definition at line 219 of file adjoints_ex4.C.

References FEMParameters::coarsen_fraction, FEMParameters::coarsen_threshold, FEMParameters::global_tolerance, mesh, FEMParameters::nelem_target, and FEMParameters::refine_fraction.

Referenced by main().

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

◆ main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 254 of file adjoints_ex4.C.

References std::abs(), libMesh::QoISet::add_indices(), libMesh::EquationSystems::add_system(), libMesh::DifferentiableSystem::adjoint_solve(), libMesh::MeshBase::all_second_order(), libMesh::DifferentiableQoI::assemble_qoi_sides, build_adjoint_refinement_error_estimator(), build_mesh_refinement(), libMesh::default_solver_package(), libMesh::err, libMesh::System::get_adjoint_solution(), libMesh::DifferentiableSystem::get_linear_solver(), LaplaceSystem::get_QoI_value(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, mesh, libMesh::EquationSystems::n_active_dofs(), libMesh::MeshBase::n_active_elem(), libMesh::out, LaplaceSystem::postprocess(), libMesh::DifferentiableSystem::postprocess_sides, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::ParallelObject::processor_id(), FEMParameters::read(), libMesh::MeshBase::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().

 {
   // 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");
 
   // Skip adaptive examples on a non-adaptive libMesh build
 #ifndef LIBMESH_ENABLE_AMR
   libmesh_example_requires(false, "--enable-amr");
 #else
 
   libMesh::out << "Started " << argv[0] << std::endl;
 
   // Make sure the general input file exists, and parse it
   {
     std::ifstream i("general.in");
     libmesh_error_msg_if(!i, '[' << init.comm().rank() << "] Can't find general.in; exiting early.");
   }
 
   // Read in parameters from the input file
   GetPot infile("general.in");
 
   // But allow the command line to override it.
   infile.parse_command_line(argc, argv);
 
   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
   mesh.read(param.domainfile.c_str());
   // Make all the elements of the mesh second order so we can compute
   // with a higher order basis
   mesh.all_second_order();
 
   // Create a mesh refinement object to do the initial uniform refinements
   // on the coarse grid read in from lshaped.xda
   MeshRefinement initial_uniform_refinements(mesh);
   initial_uniform_refinements.uniformly_refine(param.coarserefinements);
 
   libMesh::out << "Building system" << std::endl;
 
   // Build the FEMSystem
   LaplaceSystem & system = equation_systems.add_system<LaplaceSystem> ("LaplaceSystem");
 
   // Set its parameters
   set_system_parameters(system, param);
 
   libMesh::out << "Initializing systems" << std::endl;
 
   equation_systems.init ();
 
   // Print information about the mesh and system to the screen.
   mesh.print_info();
   equation_systems.print_info();
   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);
 
         linear_solver->reuse_preconditioner(false);
 
         // Solve the forward problem
         system.solve();
 
         // Write out the computed primal solution
         write_output(equation_systems, a_step, "primal", param);
 
         // Get a pointer to the primal solution vector
         NumericVector<Number> & primal_solution = *system.solution;
 
         // 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
         qois.add_indices({0,1});
 
         // 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, 0.5);
         qois.set_weight(1, 0.5);
 
         // 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 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", param);
 
         // Swap back
         primal_solution.swap(dual_solution_0);
 
         // Get a pointer to the solution vector of the adjoint problem for QoI 0
         NumericVector<Number> & dual_solution_1 = system.get_adjoint_solution(1);
 
         // Swap again
         primal_solution.swap(dual_solution_1);
         write_output(equation_systems, a_step, "adjoint_1", param);
 
         // Swap back again
         primal_solution.swap(dual_solution_1);
 
         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);
         Number QoI_1_computed = system.get_QoI_value("computed", 1);
         Number QoI_1_exact = system.get_QoI_value("exact", 1);
 
         libMesh::out << "The relative error in QoI 0 is "
                      << std::setprecision(17)
                      << std::abs(QoI_0_computed - QoI_0_exact) / std::abs(QoI_0_exact)
                      << std::endl;
 
         libMesh::out << "The relative error in QoI 1 is "
                      << std::setprecision(17)
                      << std::abs(QoI_1_computed - QoI_1_exact) / std::abs(QoI_1_exact)
                      << std::endl
                      << std::endl;
 
         // 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::abs(QoI_0_exact)
                      << std::endl;
 
         libMesh::out << "The computed relative error in QoI 1 is "
                      << std::setprecision(17)
                      << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_exact)
                      << std::endl
                      << std::endl;
 
         // 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;
 
         libMesh::out << "The effectivity index for the computed error in QoI 1 is "
                      << std::setprecision(17)
                      << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact)
                      << std::endl
                      << 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." << std::endl << 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;
       }
 
     // Do one last solve if necessary
     if (a_step == param.max_adaptivesteps)
       {
         linear_solver->reuse_preconditioner(false);
         system.solve();
 
         write_output(equation_systems, a_step, "primal", param);
 
         NumericVector<Number> & primal_solution = *system.solution;
 
         QoISet qois;
         qois.add_indices({0,1});
 
         qois.set_weight(0, 0.5);
         qois.set_weight(1, 0.5);
 
         system.assemble_qoi_sides = true;
         linear_solver->reuse_preconditioner(param.reuse_preconditioner);
         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);
 
         NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
 
         primal_solution.swap(dual_solution_0);
         write_output(equation_systems, a_step, "adjoint_0", param);
 
         primal_solution.swap(dual_solution_0);
 
         NumericVector<Number> & dual_solution_1 = system.get_adjoint_solution(1);
 
         primal_solution.swap(dual_solution_1);
         write_output(equation_systems, a_step, "adjoint_1", param);
 
         primal_solution.swap(dual_solution_1);
 
         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);
         Number QoI_1_computed = system.get_QoI_value("computed", 1);
         Number QoI_1_exact = system.get_QoI_value("exact", 1);
 
         libMesh::out << "The relative error in QoI 0 is "
                      << std::setprecision(17)
                      << std::abs(QoI_0_computed - QoI_0_exact) / std::abs(QoI_0_exact)
                      << std::endl;
 
         libMesh::out << "The relative error in QoI 1 is "
                      << std::setprecision(17)
                      << std::abs(QoI_1_computed - QoI_1_exact) / std::abs(QoI_1_exact)
                      << std::endl
                      << std::endl;
 
         // 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::abs(QoI_0_exact)
                      << std::endl;
 
         libMesh::out << "The computed relative error in QoI 1 is "
                      << std::setprecision(17)
                      << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_exact)
                      << std::endl
                      << std::endl;
 
         // 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;
 
         libMesh::out << "The effectivity index for the computed error in QoI 1 is "
                      << std::setprecision(17)
                      << std::abs(adjoint_refinement_error_estimator->get_global_QoI_error_estimate(1)) / std::abs(QoI_1_computed - QoI_1_exact)
                      << std::endl
                      << 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;
 }

◆ set_system_parameters()

void set_system_parameters	(	LaplaceSystem &	system,
		FEMParameters &	param
	)

Definition at line 156 of file adjoints_ex4.C.

Referenced by main().

 {
   // 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 = std::make_unique<SteadySolver>(system);
 
   // Nonlinear solver options
   if (param.use_petsc_snes)
   {
 #ifdef LIBMESH_HAVE_PETSC
     system.time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(system);
 #else
     libmesh_error_msg("This example requires libMesh to be compiled with PETSc support.");
 #endif
   }
   else
   {
     system.time_solver->diff_solver() = std::make_unique<NewtonSolver>(system);
     auto solver = cast_ptr<NewtonSolver*>(system.time_solver->diff_solver().get());
 
     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;
       }
     system.set_constrain_in_solver(param.constrain_in_solver);
 
     // 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;
   }
 }

◆ write_output()

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

Definition at line 95 of file adjoints_ex4.C.

References libMesh::EquationSystems::get_mesh(), libMesh::libmesh_ignore(), mesh, FEMParameters::output_exodus, FEMParameters::output_gmv, libMesh::MeshOutput< MT >::write_equation_systems(), and libMesh::ExodusII_IO::write_timestep().

Referenced by main().

 {
   // Ignore parameters when there are no output formats available.
   libmesh_ignore(es, a_step, solution_type, param);
 
 #ifdef LIBMESH_HAVE_GMV
   if (param.output_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
 
 #ifdef LIBMESH_HAVE_EXODUS_API
   if (param.output_exodus)
     {
       MeshBase & mesh = es.get_mesh();
 
       // We write out one file per adaptive step. The files are named in
       // the following way:
       // foo.e
       // foo.e-s002
       // foo.e-s003
       // ...
       // so that, if you open the first one with Paraview, it actually
       // opens the entire sequence of adapted files.
       std::ostringstream file_name_exodus;
 
       file_name_exodus << solution_type << ".e";
       if (a_step > 0)
         file_name_exodus << "-s"
                          << std::setw(3)
                          << std::setfill('0')
                          << std::right
                          << a_step + 1;
 
       // We write each adaptive step as a pseudo "time" step, where the
       // time simply matches the (1-based) adaptive step we are on.
       ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
                                        es,
                                        1,
                                        /*time=*/a_step + 1);
     }
 #endif
 }

Functions

Function Documentation

◆ build_adjoint_refinement_error_estimator()

◆ build_mesh_refinement()

◆ main()

◆ set_system_parameters()

◆ write_output()