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

void	set_system_parameters (HeatSystem &system, FEMParameters &param)

int	main (int argc, char **argv)

Function Documentation

◆ main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 337 of file adjoints_ex5.C.

 {
   // Skip adaptive examples on a non-adaptive libMesh build
 #ifndef LIBMESH_ENABLE_AMR
   libmesh_ignore(argc, argv);
   libmesh_example_requires(false, "--enable-amr");
 #else
   // Skip this 2D example if libMesh was compiled as 1D-only.
   libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
 
   // We use Dirichlet boundary conditions here
 #ifndef LIBMESH_ENABLE_DIRICHLET
   libmesh_example_requires(false, "--enable-dirichlet");
 #endif
 
   // Initialize libMesh.
   LibMeshInit init (argc, argv);
 
   // This doesn't converge with Trilinos for some reason...
   libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
 
   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.");
   }
   GetPot infile("general.in");
 
   // But allow the command line to override it.
   infile.parse_command_line(argc, argv);
 
   // Read in parameters from the input file
   FEMParameters param(init.comm());
   param.read(infile);
 
   // Create a mesh with the given dimension, distributed
   // across the default MPI communicator.
   Mesh mesh(init.comm(), cast_int<unsigned char>(param.dimension));
 
   // And an object to refine it
   auto mesh_refinement = std::make_unique<MeshRefinement>(mesh);
 
   // And an EquationSystems to run on it
   EquationSystems equation_systems (mesh);
 
   libMesh::out << "Building mesh" << std::endl;
 
   // Build a unit square
   ElemType elemtype;
 
   if (param.elementtype == "tri" ||
       param.elementtype == "unstructured")
     elemtype = TRI3;
   else
     elemtype = QUAD4;
 
   MeshTools::Generation::build_square (mesh, param.coarsegridx, param.coarsegridy,
                                        param.domain_xmin, param.domain_xmin + param.domain_edge_width,
                                        param.domain_ymin, param.domain_ymin + param.domain_edge_length,
                                        elemtype);
 
   libMesh::out << "Building system" << std::endl;
 
   HeatSystem & system = equation_systems.add_system<HeatSystem> ("HeatSystem");
 
   set_system_parameters(system, param);
 
   libMesh::out << "Initializing systems" << std::endl;
 
   // Initialize the system
   equation_systems.init ();
 
   // Refine the grid again if requested
   for (unsigned int i=0; i != param.extrarefinements; ++i)
     {
       mesh_refinement->uniformly_refine(1);
       equation_systems.reinit();
     }
 
   libMesh::out << "Setting primal initial conditions" << std::endl;
 
   read_initial_parameters();
 
   system.project_solution(initial_value, initial_grad,
                           equation_systems.parameters);
 
   // Output the H1 norm of the initial conditions
   libMesh::out << "|U("
                << system.time
                << ")|= "
                << system.calculate_norm(*system.solution, 0, H1)
                << std::endl
                << std::endl;
 
   // Tell system that we have two qois
   system.init_qois(1);
 
   // Add the adjoint vectors (and the old adjoint vectors) to system
   system.time_solver->init_adjoints();
 
   // Close up any resources initial.C needed
   finish_initialization();
 
   // Plot the initial conditions
   write_output(equation_systems, 0, "primal", param);
 
   // Print information about the mesh and system to the screen.
   mesh.print_info();
   equation_systems.print_info();
 
   // In optimized mode we catch any solver errors, so that we can
   // write the proper footers before closing.  In debug mode we just
   // let the exception throw so that gdb can grab it.
 #if defined(NDEBUG) && defined(LIBMESH_ENABLE_EXCEPTIONS)
   try
 #endif
     {
       // Now we begin the timestep loop to compute the time-accurate
       // solution of the equations.
       for (unsigned int t_step=param.initial_timestep;
            t_step != param.initial_timestep + param.n_timesteps; ++t_step)
         {
           // A pretty update message
           libMesh::out << " Solving time step "
                        << t_step
                        << ", time = "
                        << system.time
                        << std::endl;
 
           // Solve the forward problem at time t, to obtain the solution at time t + dt
           system.solve();
 
 
           // Output the H1 norm of the computed solution
           libMesh::out << "|U("
                       << system.time + system.time_solver->last_completed_timestep_size()
                       << ")|= "
                       << system.calculate_norm(*system.solution, 0, H1)
                       << std::endl;
 
           // Advance to the next timestep in a transient problem
           libMesh::out << "Advancing timestep" << std::endl << std::endl;
           system.time_solver->advance_timestep();
 
           // Write out this timestep
           write_output(equation_systems, t_step+1, "primal", param);
         }
       // End timestep loop
 
 
       // Now we will solve the backwards in time adjoint problem
       libMesh::out << std::endl << "Solving the adjoint problem" << std::endl;
 
       // We need to tell the library that it needs to project the adjoint, so
       // MemorySolutionHistory knows it has to save it
 
       const std::string & adjoint_solution_name0 = "adjoint_solution0";
       const std::string & old_adjoint_solution_name0 = "_old_adjoint_solution0";
       // Tell the library to project the adjoint vector, and hence, memory solution history to
       // save it
       system.set_vector_preservation(adjoint_solution_name0, true);
       system.set_vector_preservation(old_adjoint_solution_name0, true);
 
       libMesh::out << "Setting adjoint initial conditions Z("
                    << system.time
                    << ")"
                    <<std::endl;
 
       // The first thing we have to do is to apply the adjoint initial
       // condition. The user should supply these. Here they are specified
       // in the functions adjoint_initial_value and adjoint_initial_gradient
       system.project_vector(adjoint_initial_value,
                             adjoint_initial_grad,
                             equation_systems.parameters,
                             system.get_adjoint_solution(0));
 
       // Since we have specified an adjoint solution for the current
       // time (T), set the adjoint_already_solved boolean to true, so
       // we dont solve unnecessarily in the adjoint sensitivity method
       system.set_adjoint_already_solved(true);
 
       Real Z_norm = system.calculate_norm(system.get_adjoint_solution(), 0, H1);
 
       libMesh::out << "|Z("
                    << system.time
                    << ")|= "
                    << Z_norm
                    << std::endl
                    << std::endl;
 
       // Call the adjoint advance timestep here to save the adjoint
       // initial conditions to disk
       system.time_solver->adjoint_advance_timestep();
 
       Real Z_old_norm = system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1);
 
       // Make sure that we have copied the old adjoint solution properly
       libmesh_assert(Z_norm == Z_old_norm);
 
       libMesh::out << "|Z_old("
                        << system.time + (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
                        << ")|= "
                        << Z_old_norm
                        << std::endl
                        << std::endl;
 
       write_output(equation_systems, param.n_timesteps, "dual", param);
 
       // Now that the adjoint initial condition is set, we will start the
       // backwards in time adjoint integration
 
       // For loop stepping backwards in time
       for (unsigned int t_step=param.initial_timestep;
            t_step != param.initial_timestep + param.n_timesteps; ++t_step)
         {
           //A pretty update message
           libMesh::out << " Solving adjoint time step "
                        << t_step
                        << ", time = "
                        << system.time
                        << std::endl;
 
           // Output the H1 norm of the retrieved primal solution from the last call
           // to adjoint_advance_timestep
           libMesh::out << "|U("
                        << system.time
                        << ")|= "
                        << system.calculate_norm(*system.solution, 0, H1)
                        << std::endl;
 
           libMesh::out << "|U("
                          << system.time - (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
                          << ")|= "
                          << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
                          << std::endl;
 
           system.set_adjoint_already_solved(false);
 
           system.adjoint_solve();
 
           Z_norm = system.calculate_norm(system.get_adjoint_solution(), 0, H1);
 
           libMesh::out << "|Z("
                        << system.time
                        << ")|= "
                        << Z_norm
                        << std::endl
                        << std::endl;
 
           // 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);
 
           libMesh::out << "Saving adjoint and retrieving primal solutions at time t=" << system.time - system.deltat << std::endl;
 
           // The adjoint_advance_timestep function calls the retrieve and store
           // function of the memory_solution_history class via the
           // memory_solution_history object we declared earlier.  The
           // retrieve function sets the system primal vectors to their
           // values at the current time instance.
           system.time_solver->adjoint_advance_timestep();
 
           Z_old_norm = system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1);
 
           // Make sure that we have copied the old adjoint solution properly
           libmesh_assert(Z_norm == Z_old_norm);
 
           libMesh::out << "|Z_old("
                        << system.time + (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
                        << ")|= "
                        << Z_old_norm
                        << std::endl
                        << std::endl;
 
           // Get a pointer to the primal solution vector
           NumericVector<Number> & primal_solution = *system.solution;
 
           // Get a pointer to the solution vector of the adjoint problem for QoI 0
           NumericVector<Number> & dual_solution_0 = system.get_adjoint_solution(0);
 
           // Swap the primal and dual solutions so we can write out the adjoint solution
           primal_solution.swap(dual_solution_0);
 
           write_output(equation_systems, param.n_timesteps - (t_step + 1), "dual", param);
 
           // Swap back
           primal_solution.swap(dual_solution_0);
         }
       // End adjoint timestep loop
 
       // Now that we have computed both the primal and adjoint solutions, we compute the sensitivities to the parameter p
       // dQ/dp = int_{0}^{T} partialQ/partialp - partialR/partialp(u,z;p) dt
       // The quantity partialQ/partialp - partialR/partialp(u,z;p) is evaluated internally by the ImplicitSystem::adjoint_qoi_parameter_sensitivity function.
       // This sensitivity evaluation is called internally by an overloaded TimeSolver::integrate_adjoint_sensitivity method which we call below.
 
       // Reset the time
       system.time = 0.0;
 
       // Prepare the quantities we need to pass to TimeSolver::integrate_adjoint_sensitivity
       QoISet qois;
 
       qois.add_indices({0});
       qois.set_weight(0, 1.0);
 
       // A SensitivityData object
       SensitivityData sensitivities(qois, system, system.get_parameter_vector());
 
       // Accumulator for the time integrated total sensitivity
       Number total_sensitivity = 0.0;
 
       // Retrieve the primal and adjoint solutions at the current timestep
       system.time_solver->retrieve_timestep();
 
       Z_old_norm = system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1);
 
       // Assert that the old adjoint values match hard coded numbers for 1st timestep or 1st half time step from the adjoint solve
       if(param.timesolver_tolerance)
       {
         libmesh_error_msg_if(std::abs(Z_old_norm - (2.23366)) >= 2.e-4,
                              "Mismatch in expected Z0_old norm for the 1st half timestep!");
        }
       else
       {
         libmesh_error_msg_if(std::abs(Z_old_norm - (2.23627)) >= 2.e-4,
                              "Mismatch in expected Z0_old norm for the 1st timestep!");
       }
 
       // A pretty update message
       libMesh::out << "Retrieved, "
                << "time = "
                << system.time
                << std::endl;
 
       libMesh::out << "|U("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(*system.solution, 0, H1)
                    << std::endl;
 
       libMesh::out << "|U_old("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
                    << std::endl;
 
       libMesh::out << "|Z("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
                    << std::endl
                    << std::endl;
 
       libMesh::out << "|Z_old("
                        << system.time + (system.time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
                        << ")|= "
                        << system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1)
                        << std::endl
                        << std::endl;
 
       // Now we begin the timestep loop to compute the time-accurate
       // adjoint sensitivities
       for (unsigned int t_step=param.initial_timestep;
            t_step != param.initial_timestep + param.n_timesteps; ++t_step)
         {
           // Call the postprocess function which we have overloaded to compute
           // accumulate the perturbed residuals
           system.time_solver->integrate_adjoint_sensitivity(qois, system.get_parameter_vector(), sensitivities);
 
           // A pretty update message
           libMesh::out << "Retrieved, "
                << "time = "
                << system.time
                << std::endl;
 
           libMesh::out << "|U("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(*system.solution, 0, H1)
                    << std::endl;
 
           libMesh::out << "|U_old("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(system.get_vector("_old_nonlinear_solution"), 0, H1)
                    << std::endl;
 
           libMesh::out << "|Z("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(system.get_adjoint_solution(0), 0, H1)
                    << std::endl
                    << std::endl;
 
           libMesh::out << "|Z_old("
                    << system.time
                    << ")|= "
                    << system.calculate_norm(system.get_vector(old_adjoint_solution_name0), 0, H1)
                    << std::endl
                    << std::endl;
 
           // Add the contribution of the current timestep to the total sensitivity
           total_sensitivity += sensitivities[0][0];
         }
 
       // Print it out
       libMesh::out << "Sensitivity of QoI 0 w.r.t parameter 0 is: "
                    << total_sensitivity
                    << std::endl;
 
       // Hard coded test to ensure that the actual numbers we are
       // getting are what they should be
       // The 2e-4 tolerance is chosen to ensure success even with
       // 32-bit floats
       if(param.timesolver_tolerance)
       {
         libmesh_error_msg_if(std::abs(system.time - (1.0089)) >= 2.e-4,
                              "Mismatch in end time reached by adaptive timestepper!");
 
         libmesh_error_msg_if(std::abs(total_sensitivity - 4.87767) >= 3.e-3,
                              "Mismatch in sensitivity gold value!");
       }
       else
       {
         libmesh_error_msg_if(std::abs(total_sensitivity - 4.83551) >= 2.e-4,
                              "Mismatch in sensitivity gold value!");
       }
     }
 #if defined(NDEBUG) && defined(LIBMESH_ENABLE_EXCEPTIONS)
   catch (...)
     {
       libMesh::err << '[' << mesh.processor_id()
                    << "] Caught exception; exiting early." << std::endl;
     }
 #endif
 
   libMesh::err << '[' << mesh.processor_id()
                << "] Completing output."
                << std::endl;
 
   // All done.
   return 0;
 
 #endif // LIBMESH_ENABLE_AMR
 }

◆ set_system_parameters()

void set_system_parameters	(	HeatSystem &	system,
		FEMParameters &	param
	)

Definition at line 179 of file adjoints_ex5.C.

Referenced by main().

 {
   // Use the prescribed FE type
   system.fe_family() = param.fe_family[0];
   system.fe_order() = param.fe_order[0];
 
   // Use analytical jacobians?
   system.analytic_jacobians() = param.analytic_jacobians;
 
   // Verify analytic jacobians against numerical ones?
   system.verify_analytic_jacobians = param.verify_analytic_jacobians;
   system.numerical_jacobian_h = param.numerical_jacobian_h;
 
   // More desperate debugging options
   system.print_solution_norms    = param.print_solution_norms;
   system.print_solutions         = param.print_solutions;
   system.print_residual_norms    = param.print_residual_norms;
   system.print_residuals         = param.print_residuals;
   system.print_jacobian_norms    = param.print_jacobian_norms;
   system.print_jacobians         = param.print_jacobians;
   system.print_element_jacobians = param.print_element_jacobians;
   system.print_element_residuals = param.print_element_residuals;
 
   // Solve this as a time-dependent or steady system
   if (param.transient)
     {
       std::unique_ptr<UnsteadySolver> innersolver;
       if (param.timesolver_core == "euler2")
         {
           auto euler2solver = std::make_unique<Euler2Solver>(system);
           euler2solver->theta = param.timesolver_theta;
           innersolver = std::move(euler2solver);
         }
       else if (param.timesolver_core == "euler")
         {
           auto eulersolver = std::make_unique<EulerSolver>(system);
           eulersolver->theta = param.timesolver_theta;
           innersolver = std::move(eulersolver);
         }
       else
         libmesh_error_msg("This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods.");
 
       if (param.timesolver_tolerance)
         {
           system.time_solver = std::make_unique<TwostepTimeSolver>(system);
           auto timesolver = cast_ptr<TwostepTimeSolver *>(system.time_solver.get());
 
           timesolver->max_growth       = param.timesolver_maxgrowth;
           timesolver->target_tolerance = param.timesolver_tolerance;
           timesolver->upper_tolerance  = param.timesolver_upper_tolerance;
           timesolver->component_norm   = SystemNorm(param.timesolver_norm);
           timesolver->core_time_solver = std::move(innersolver);
 
           // The Memory/File Solution History object we will set the system SolutionHistory object to
           if(param.solution_history_type == "memory")
           {
             MemorySolutionHistory heatsystem_solution_history(system);
             system.time_solver->set_solution_history(heatsystem_solution_history);
           }
           else if (param.solution_history_type == "file")
           {
             FileSolutionHistory heatsystem_solution_history(system);
             system.time_solver->set_solution_history(heatsystem_solution_history);
           }
           else
             libmesh_error_msg("Unrecognized solution history type: " << param.solution_history_type);
         }
       else
       {
         system.time_solver = std::move(innersolver);
 
         // The Memory/File Solution History object we will set the system SolutionHistory object to
         if(param.solution_history_type == "memory")
         {
           MemorySolutionHistory heatsystem_solution_history(system);
           system.time_solver->set_solution_history(heatsystem_solution_history);
         }
         else if (param.solution_history_type == "file")
         {
           FileSolutionHistory heatsystem_solution_history(system);
           system.time_solver->set_solution_history(heatsystem_solution_history);
         }
         else
           libmesh_error_msg("Unrecognized solution history type: " << param.solution_history_type);
       }
     }
   else // !param.transient
     system.time_solver = std::make_unique<SteadySolver>(system);
 
   system.time_solver->reduce_deltat_on_diffsolver_failure = param.deltat_reductions;
   system.time_solver->quiet = param.time_solver_quiet;
 
  #ifdef LIBMESH_ENABLE_DIRICHLET
   // Create any Dirichlet boundary conditions
   typedef
     std::map<boundary_id_type, FunctionBase<Number> *>::
     const_iterator Iter;
 
   for (Iter i = param.dirichlet_conditions.begin();
        i != param.dirichlet_conditions.end(); ++i)
     {
       boundary_id_type b = i->first;
       FunctionBase<Number> *f = i->second;
 
       system.get_dof_map().add_dirichlet_boundary(DirichletBoundary({b},
                                                                     param.dirichlet_condition_variables[b],
                                                                     f));
 
       libMesh::out << "Added Dirichlet boundary " << b << " for variables ";
       for (std::size_t vi=0; vi != param.dirichlet_condition_variables[b].size(); ++vi)
         libMesh::out << param.dirichlet_condition_variables[b][vi];
       libMesh::out << std::endl;
     }
  #endif // LIBMESH_ENABLE_DIRICHLET
 
   // Set the time stepping options
   system.deltat = param.deltat;
 
   // And the integration options
   system.extra_quadrature_order = param.extra_quadrature_order;
 
   // And the nonlinear solver options
   if (param.use_petsc_snes)
     {
 #ifdef LIBMESH_HAVE_PETSC
       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->verbose                     = param.solver_verbose;
       solver->max_nonlinear_iterations    = param.max_nonlinear_iterations;
       solver->minsteplength               = param.min_step_length;
       solver->relative_step_tolerance     = param.relative_step_tolerance;
       solver->relative_residual_tolerance = param.relative_residual_tolerance;
       solver->require_residual_reduction  = param.require_residual_reduction;
       solver->linear_tolerance_multiplier = param.linear_tolerance_multiplier;
       if (system.time_solver->reduce_deltat_on_diffsolver_failure)
         {
           solver->continue_after_max_iterations = true;
           solver->continue_after_backtrack_failure = true;
         }
       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	t_step,
		std::string	solution_type,
		FEMParameters &	param
	)

Definition at line 122 of file adjoints_ex5.C.

References libMesh::libmesh_ignore(), mesh, FEMParameters::output_exodus, and FEMParameters::output_gmv.

Referenced by main().

 {
   // Ignore parameters when there are no output formats available.
   libmesh_ignore(es, t_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
                     << t_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 timestep. 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 (t_step > 0)
         file_name_exodus << "-s"
                          << std::setw(3)
                          << std::setfill('0')
                          << std::right
                          << t_step + 1;
 
       // TODO: Get the current time from the System...
       ExodusII_IO(mesh).write_timestep(file_name_exodus.str(),
                                        es,
                                        1,
                                        /*time=*/t_step + 1);
     }
 #endif
 }

Functions

Function Documentation

◆ main()

◆ set_system_parameters()

◆ write_output()