81 #include "adjoint_initial.h" 82 #include "femparameters.h" 83 #include "heatsystem.h" 86 #include "libmesh/equation_systems.h" 87 #include "libmesh/dirichlet_boundaries.h" 88 #include "libmesh/system_norm.h" 89 #include "libmesh/numeric_vector.h" 90 #include "libmesh/enum_solver_package.h" 92 #include "libmesh/mesh.h" 93 #include "libmesh/mesh_generation.h" 94 #include "libmesh/mesh_refinement.h" 96 #include "libmesh/petsc_diff_solver.h" 97 #include "libmesh/steady_solver.h" 98 #include "libmesh/euler_solver.h" 99 #include "libmesh/euler2_solver.h" 100 #include "libmesh/newton_solver.h" 101 #include "libmesh/petsc_diff_solver.h" 102 #include "libmesh/twostep_time_solver.h" 104 #include "libmesh/getpot.h" 105 #include "libmesh/tecplot_io.h" 106 #include "libmesh/gmv_io.h" 107 #include "libmesh/exodusII_io.h" 109 #include "libmesh/sensitivity_data.h" 112 #include "libmesh/solution_history.h" 113 #include "libmesh/memory_solution_history.h" 114 #include "libmesh/file_solution_history.h" 118 #include <sys/time.h> 124 std::string solution_type,
130 #ifdef LIBMESH_HAVE_GMV 133 MeshBase &
mesh = es.get_mesh();
135 std::ostringstream file_name_gmv;
136 file_name_gmv << solution_type
143 GMVIO(
mesh).write_equation_systems(file_name_gmv.str(), es);
147 #ifdef LIBMESH_HAVE_EXODUS_API 150 MeshBase &
mesh = es.get_mesh();
160 std::ostringstream file_name_exodus;
162 file_name_exodus << solution_type <<
".e";
164 file_name_exodus <<
"-s" 171 ExodusII_IO(
mesh).write_timestep(file_name_exodus.str(),
205 std::unique_ptr<UnsteadySolver> innersolver;
208 auto euler2solver = std::make_unique<Euler2Solver>(system);
210 innersolver = std::move(euler2solver);
214 auto eulersolver = std::make_unique<EulerSolver>(system);
216 innersolver = std::move(eulersolver);
219 libmesh_error_msg(
"This example (and unsteady adjoints in libMesh) only support Backward Euler and explicit methods.");
223 system.
time_solver = std::make_unique<TwostepTimeSolver>(system);
224 auto timesolver = cast_ptr<TwostepTimeSolver *>(system.
time_solver.get());
230 timesolver->core_time_solver = std::move(innersolver);
235 MemorySolutionHistory heatsystem_solution_history(system);
236 system.
time_solver->set_solution_history(heatsystem_solution_history);
240 FileSolutionHistory heatsystem_solution_history(system);
241 system.
time_solver->set_solution_history(heatsystem_solution_history);
253 MemorySolutionHistory heatsystem_solution_history(system);
254 system.
time_solver->set_solution_history(heatsystem_solution_history);
258 FileSolutionHistory heatsystem_solution_history(system);
259 system.
time_solver->set_solution_history(heatsystem_solution_history);
266 system.
time_solver = std::make_unique<SteadySolver>(system);
271 #ifdef LIBMESH_ENABLE_DIRICHLET 274 std::map<boundary_id_type, FunctionBase<Number> *>::
281 FunctionBase<Number> *f = i->second;
287 libMesh::out <<
"Added Dirichlet boundary " << b <<
" for variables ";
292 #endif // LIBMESH_ENABLE_DIRICHLET 303 #ifdef LIBMESH_HAVE_PETSC 304 system.
time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(system);
306 libmesh_error_msg(
"This example requires libMesh to be compiled with PETSc support.");
311 system.
time_solver->diff_solver() = std::make_unique<NewtonSolver>(system);
312 auto solver = cast_ptr<NewtonSolver *>(system.
time_solver->diff_solver().get());
322 if (system.
time_solver->reduce_deltat_on_diffsolver_failure)
324 solver->continue_after_max_iterations =
true;
325 solver->continue_after_backtrack_failure =
true;
337 int main (
int argc,
char ** argv)
340 #ifndef LIBMESH_ENABLE_AMR 342 libmesh_example_requires(
false,
"--enable-amr");
345 libmesh_example_requires(2 <= LIBMESH_DIM,
"2D support");
348 #ifndef LIBMESH_ENABLE_DIRICHLET 349 libmesh_example_requires(
false,
"--enable-dirichlet");
353 LibMeshInit
init (argc, argv);
362 std::ifstream i(
"general.in");
363 libmesh_error_msg_if(!i,
'[' <<
init.comm().rank() <<
"] Can't find general.in; exiting early.");
365 GetPot infile(
"general.in");
368 infile.parse_command_line(argc, argv);
376 Mesh
mesh(
init.comm(), cast_int<unsigned char>(param.dimension));
379 auto mesh_refinement = std::make_unique<MeshRefinement>(
mesh);
382 EquationSystems equation_systems (
mesh);
389 if (param.elementtype ==
"tri" ||
390 param.elementtype ==
"unstructured")
396 param.domain_xmin, param.domain_xmin + param.domain_edge_width,
397 param.domain_ymin, param.domain_ymin + param.domain_edge_length,
409 equation_systems.init ();
412 for (
unsigned int i=0; i != param.extrarefinements; ++i)
414 mesh_refinement->uniformly_refine(1);
415 equation_systems.reinit();
418 libMesh::out <<
"Setting primal initial conditions" << std::endl;
423 equation_systems.parameters);
447 equation_systems.print_info();
452 #if defined(NDEBUG) && defined(LIBMESH_ENABLE_EXCEPTIONS) 458 for (
unsigned int t_step=param.initial_timestep;
459 t_step != param.initial_timestep + param.n_timesteps; ++t_step)
480 libMesh::out <<
"Advancing timestep" << std::endl << std::endl;
484 write_output(equation_systems, t_step+1,
"primal", param);
491 libMesh::out << std::endl <<
"Solving the adjoint problem" << std::endl;
496 const std::string & adjoint_solution_name0 =
"adjoint_solution0";
497 const std::string & old_adjoint_solution_name0 =
"_old_adjoint_solution0";
513 equation_systems.parameters,
540 << system.
time + (system.
time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
546 write_output(equation_systems, param.n_timesteps,
"dual", param);
552 for (
unsigned int t_step=param.initial_timestep;
553 t_step != param.initial_timestep + param.n_timesteps; ++t_step)
571 << system.
time - (system.
time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
594 libMesh::out <<
"Saving adjoint and retrieving primal solutions at time t=" << system.
time - system.
deltat << std::endl;
609 << system.
time + (system.
time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
622 primal_solution.
swap(dual_solution_0);
624 write_output(equation_systems, param.n_timesteps - (t_step + 1),
"dual", param);
627 primal_solution.swap(dual_solution_0);
642 qois.add_indices({0});
643 qois.set_weight(0, 1.0);
649 Number total_sensitivity = 0.0;
657 if(param.timesolver_tolerance)
659 libmesh_error_msg_if(
std::abs(Z_old_norm - (2.23366)) >= 2.e-4,
660 "Mismatch in expected Z0_old norm for the 1st half timestep!");
664 libmesh_error_msg_if(
std::abs(Z_old_norm - (2.23627)) >= 2.e-4,
665 "Mismatch in expected Z0_old norm for the 1st timestep!");
694 << system.
time + (system.
time_solver->last_completed_timestep_size())/((param.timesolver_tolerance) ? 2.0 : 1.0)
702 for (
unsigned int t_step=param.initial_timestep;
703 t_step != param.initial_timestep + param.n_timesteps; ++t_step)
742 total_sensitivity += sensitivities[0][0];
746 libMesh::out <<
"Sensitivity of QoI 0 w.r.t parameter 0 is: " 754 if(param.timesolver_tolerance)
756 libmesh_error_msg_if(
std::abs(system.
time - (1.0089)) >= 2.e-4,
757 "Mismatch in end time reached by adaptive timestepper!");
759 libmesh_error_msg_if(
std::abs(total_sensitivity - 4.87767) >= 3.e-3,
760 "Mismatch in sensitivity gold value!");
764 libmesh_error_msg_if(
std::abs(total_sensitivity - 4.83551) >= 2.e-4,
765 "Mismatch in sensitivity gold value!");
768 #if defined(NDEBUG) && defined(LIBMESH_ENABLE_EXCEPTIONS) 772 <<
"] Caught exception; exiting early." << std::endl;
777 <<
"] Completing output." 783 #endif // LIBMESH_ENABLE_AMR
libMesh::Real timesolver_upper_tolerance
Real time
For time-dependent problems, this is the time t at the beginning of the current timestep.
bool print_solution_norms
Number adjoint_initial_value(const Point &p, const Parameters &, const std::string &, const std::string &)
ElemType
Defines an enum for geometric element types.
void write_output(EquationSystems &es, unsigned int t_step, std::string solution_type, FEMParameters ¶m)
bool print_element_residuals
std::string solution_history_type
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
libMesh::Real timesolver_tolerance
void finish_initialization()
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
virtual void set_constrain_in_solver(bool enable)
set_constrain_in_solver to false to apply constraints only via residual terms in the systems to be so...
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
Gradient adjoint_initial_grad(const Point &p, const Parameters &, const std::string &, const std::string &)
int extra_quadrature_order
A member int that can be employed to indicate increased or reduced quadrature order.
Real numerical_jacobian_h
If calculating numeric jacobians is required, the FEMSystem will perturb each solution vector entry b...
std::string timesolver_core
bool require_residual_reduction
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we're going to use.
bool print_residual_norms
void init_qois(unsigned int n_qois)
Accessors for qoi and qoi_error_estimates vectors.
libMesh::Real timesolver_maxgrowth
int main(int argc, char **argv)
libMesh::Real relative_step_tolerance
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
std::vector< unsigned int > fe_order
bool print_element_residuals
Set print_element_residuals to true to print each R_elem contribution.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
SolverPackage default_solver_package()
void project_solution(FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr) const
Projects arbitrary functions onto the current solution.
bool print_solution_norms
Set print_residual_norms to true to print |U| whenever it is used in an assembly() call...
void libmesh_ignore(const Args &...)
unsigned int max_linear_iterations
bool print_residual_norms
Set print_residual_norms to true to print |F| whenever it is assembled.
double minimum_linear_tolerance
ParameterVector & get_parameter_vector()
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
std::vector< std::string > fe_family
libMesh::Real verify_analytic_jacobians
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Real calculate_norm(const NumericVector< Number > &v, unsigned int var, FEMNormType norm_type, std::set< unsigned int > *skip_dimensions=nullptr) const
bool print_residuals
Set print_residuals to true to print F whenever it is assembled.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
double initial_linear_tolerance
int extra_quadrature_order
void read(GetPot &input, const std::vector< std::string > *other_variable_names=nullptr)
bool print_solutions
Set print_solutions to true to print U whenever it is used in an assembly() call. ...
libMesh::Real relative_residual_tolerance
void set_system_parameters(HeatSystem &system, FEMParameters ¶m)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void set_vector_preservation(const std::string &vec_name, bool preserve)
Allows one to set the boolean controlling whether the vector identified by vec_name should be "preser...
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
libMesh::Real min_step_length
std::string & fe_family()
libMesh::Real timesolver_theta
Gradient initial_grad(const Point &, const Parameters &, const std::string &, const std::string &)
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
bool print_element_jacobians
unsigned int max_nonlinear_iterations
double linear_tolerance_multiplier
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
std::map< libMesh::boundary_id_type, std::vector< unsigned int > > dirichlet_condition_variables
void project_vector(NumericVector< Number > &new_vector, FunctionBase< Number > *f, FunctionBase< Gradient > *g=nullptr, int is_adjoint=-1) const
Projects arbitrary functions onto a vector of degree of freedom values for the current system...
virtual void solve() override
Invokes the solver associated with the system.
bool & analytic_jacobians()
void read_initial_parameters()
libMesh::Real numerical_jacobian_h
std::map< libMesh::boundary_id_type, libMesh::FunctionBase< libMesh::Number > * > dirichlet_conditions
const DofMap & get_dof_map() const
std::vector< libMesh::FEMNormType > timesolver_norm
template class LIBMESH_EXPORT NumericVector< Number >
Number initial_value(const Point &, const Parameters &, const std::string &, const std::string &)
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) override
This function sets the _is_adjoint boolean member of TimeSolver to true and then calls the adjoint_so...
unsigned int & fe_order()
const NumericVector< Number > & get_vector(std::string_view vec_name) const
unsigned int deltat_reductions