66 #include "libmesh/eigen_sparse_linear_solver.h" 67 #include "libmesh/enum_solver_package.h" 68 #include "libmesh/enum_solver_type.h" 69 #include "libmesh/equation_systems.h" 70 #include "libmesh/error_vector.h" 71 #include "libmesh/mesh.h" 72 #include "libmesh/mesh_refinement.h" 73 #include "libmesh/newton_solver.h" 74 #include "libmesh/numeric_vector.h" 75 #include "libmesh/petsc_diff_solver.h" 76 #include "libmesh/steady_solver.h" 77 #include "libmesh/system_norm.h" 80 #include "libmesh/parameter_vector.h" 81 #include "libmesh/sensitivity_data.h" 84 #include "libmesh/kelly_error_estimator.h" 85 #include "libmesh/patch_recovery_error_estimator.h" 88 #include "libmesh/adjoint_residual_error_estimator.h" 89 #include "libmesh/qoi_set.h" 92 #include "libmesh/getpot.h" 93 #include "libmesh/gmv_io.h" 94 #include "libmesh/exodusII_io.h" 97 #include "femparameters.h" 112 std::string solution_type,
118 #ifdef LIBMESH_HAVE_GMV 123 std::ostringstream file_name_gmv;
124 file_name_gmv << solution_type
135 #ifdef LIBMESH_HAVE_EXODUS_API 148 std::ostringstream file_name_exodus;
150 file_name_exodus << solution_type <<
".e";
152 file_name_exodus <<
"-s" 173 #ifdef LIBMESH_HAVE_EIGEN_SPARSE 177 if (eigen_linear_solver)
189 auto solver = cast_ptr<NewtonSolver*>(diff_solver);
222 system.
time_solver = std::make_unique<SteadySolver>(system);
227 #ifdef LIBMESH_HAVE_PETSC 228 system.
time_solver->diff_solver() = std::make_unique<PetscDiffSolver>(system);
230 libmesh_error_msg(
"This example requires libMesh to be compiled with PETSc support.");
235 system.
time_solver->diff_solver() = std::make_unique<NewtonSolver>(system);
236 auto solver = cast_ptr<NewtonSolver*>(system.
time_solver->diff_solver().get());
245 if (system.
time_solver->reduce_deltat_on_diffsolver_failure)
247 solver->continue_after_max_iterations =
true;
248 solver->continue_after_backtrack_failure =
true;
262 #ifdef LIBMESH_ENABLE_AMR 267 auto mesh_refinement = std::make_unique<MeshRefinement>(
mesh);
268 mesh_refinement->coarsen_by_parents() =
true;
275 return mesh_refinement;
278 #endif // LIBMESH_ENABLE_AMR 291 libMesh::out <<
"Using Kelly Error Estimator" << std::endl;
293 return std::make_unique<KellyErrorEstimator>();
297 libMesh::out <<
"Using Adjoint Residual Error Estimator with Patch Recovery Weights" << std::endl << std::endl;
299 auto adjoint_residual_estimator = std::make_unique<AdjointResidualErrorEstimator>();
301 adjoint_residual_estimator->error_plot_suffix =
"error.gmv";
303 adjoint_residual_estimator->primal_error_estimator() = std::make_unique<PatchRecoveryErrorEstimator>();
304 adjoint_residual_estimator->primal_error_estimator()->error_norm.set_type(0,
H1_SEMINORM);
306 adjoint_residual_estimator->dual_error_estimator() = std::make_unique<PatchRecoveryErrorEstimator>();
307 adjoint_residual_estimator->dual_error_estimator()->error_norm.set_type(0,
H1_SEMINORM);
309 return adjoint_residual_estimator;
312 libmesh_error_msg(
"Unknown indicator_type = " << param.
indicator_type);
316 int main (
int argc,
char ** argv)
323 "--enable-petsc, --enable-trilinos, or --enable-eigen");
326 #ifndef LIBMESH_ENABLE_AMR 327 libmesh_example_requires(
false,
"--enable-amr");
334 std::ifstream i(
"general.in");
335 libmesh_error_msg_if(!i,
'[' <<
init.comm().rank() <<
"] Can't find general.in; exiting early.");
339 GetPot infile(
"general.in");
341 infile.parse_command_line(argc, argv);
347 libmesh_example_requires(2 <= LIBMESH_DIM,
"2D support");
354 std::unique_ptr<MeshRefinement> mesh_refinement =
360 libMesh::out <<
"Reading in and building the mesh" << std::endl;
363 mesh.
read(param.domainfile.c_str());
382 qois.set_weight(0, 0.5);
395 equation_systems.
init ();
403 unsigned int a_step = 0;
404 for (; a_step != param.max_adaptivesteps; ++a_step)
408 if (param.global_tolerance != 0.)
409 libmesh_assert_equal_to (param.nelem_target, 0);
413 libmesh_assert_greater (param.nelem_target, 0);
419 write_output(equation_systems, a_step,
"primal", param);
440 GetPot infile_l_shaped(
"l-shaped.in");
442 Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
443 Number sensitivity_QoI_0_0_exact = infile_l_shaped(
"sensitivity_0_0", 0.0);
444 Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
445 Number sensitivity_QoI_0_1_exact = infile_l_shaped(
"sensitivity_0_1", 0.0);
451 <<
" active elements and " 456 libMesh::out <<
"Sensitivity of QoI one to Parameter one is " 457 << sensitivity_QoI_0_0_computed
459 libMesh::out <<
"Sensitivity of QoI one to Parameter two is " 460 << sensitivity_QoI_0_1_computed
463 libMesh::out <<
"The relative error in sensitivity QoI_0_0 is " 464 << std::setprecision(17)
465 <<
std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact) /
std::abs(sensitivity_QoI_0_0_exact)
468 libMesh::out <<
"The relative error in sensitivity QoI_0_1 is " 469 << std::setprecision(17)
470 <<
std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact) /
std::abs(sensitivity_QoI_0_1_exact)
478 primal_solution.
swap(dual_solution_0);
479 write_output(equation_systems, a_step,
"adjoint_0", param);
482 primal_solution.
swap(dual_solution_0);
489 if (param.refine_uniformly)
491 libMesh::out <<
"Refining Uniformly" << std::endl << std::endl;
493 mesh_refinement->uniformly_refine(1);
496 else if (param.global_tolerance >= 0. && param.nelem_target == 0.)
502 std::unique_ptr<ErrorEstimator> error_estimator =
506 error_estimator->estimate_error(system, error);
508 mesh_refinement->flag_elements_by_error_tolerance (error);
510 mesh_refinement->refine_and_coarsen_elements();
519 std::unique_ptr<ErrorEstimator> error_estimator =
523 error_estimator->estimate_error(system, error);
527 libMesh::out<<
"We reached the target number of elements."<<std::endl <<std::endl;
531 mesh_refinement->flag_elements_by_nelem_target (error);
533 mesh_refinement->refine_and_coarsen_elements();
537 equation_systems.
reinit();
544 <<
" active elements and " 551 if (a_step == param.max_adaptivesteps)
555 write_output(equation_systems, a_step,
"primal", param);
570 GetPot infile_l_shaped(
"l-shaped.in");
572 Number sensitivity_QoI_0_0_computed = sensitivities[0][0];
573 Number sensitivity_QoI_0_0_exact = infile_l_shaped(
"sensitivity_0_0", 0.0);
574 Number sensitivity_QoI_0_1_computed = sensitivities[0][1];
575 Number sensitivity_QoI_0_1_exact = infile_l_shaped(
"sensitivity_0_1", 0.0);
581 <<
" active elements and " 586 libMesh::out <<
"Sensitivity of QoI one to Parameter one is " 587 << sensitivity_QoI_0_0_computed
590 libMesh::out <<
"Sensitivity of QoI one to Parameter two is " 591 << sensitivity_QoI_0_1_computed
595 << std::setprecision(17)
596 <<
std::abs(sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact
600 << std::setprecision(17)
601 <<
std::abs(sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact
606 libmesh_assert_less(
std::abs((sensitivity_QoI_0_0_computed - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
607 libmesh_assert_less(
std::abs((sensitivity_QoI_0_1_computed - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
612 const bool forward_sensitivity = infile(
"--forward_sensitivity",
true);
617 if (forward_sensitivity)
625 libmesh_assert_less(
std::abs((forward_sensitivities[0][0] - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact), 2.e-4);
626 libmesh_assert_less(
std::abs((forward_sensitivities[0][1] - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact), 2.e-4);
631 (
std::abs((forward_sensitivities[0][0] - sensitivity_QoI_0_0_computed)/sensitivity_QoI_0_0_computed),
TOLERANCE);
633 (
std::abs((forward_sensitivities[0][1] - sensitivity_QoI_0_1_computed)/sensitivity_QoI_0_1_computed),
TOLERANCE);
635 libMesh::out <<
"The error in forward calculation of sensitivity QoI_0_0 is " 636 << std::setprecision(17)
637 <<
std::abs(forward_sensitivities[0][0] - sensitivity_QoI_0_0_exact)/sensitivity_QoI_0_0_exact
640 libMesh::out <<
"The error in forward calculation of sensitivity QoI_0_1 is " 641 << std::setprecision(17)
642 <<
std::abs(forward_sensitivities[0][1] - sensitivity_QoI_0_1_exact)/sensitivity_QoI_0_1_exact
651 primal_solution.
swap(dual_solution_0);
652 write_output(equation_systems, a_step,
"adjoint_0", param);
654 primal_solution.
swap(dual_solution_0);
659 <<
"] Completing output." 662 #endif // #ifndef LIBMESH_ENABLE_AMR unsigned int nelem_target
bool print_solution_norms
This is the EquationSystems class.
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
virtual dof_id_type n_active_elem() const =0
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
int main(int argc, char **argv)
void set_adjoint_already_solved(bool setting)
Setter for the adjoint_already_solved boolean.
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
std::unique_ptr< ErrorEstimator > build_error_estimator(FEMParameters ¶m)
static constexpr Real TOLERANCE
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...
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
bool print_jacobian_norms
Set print_jacobian_norms to true to print |J| whenever it is assembled.
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
void adjust_linear_solver(LinearSolver< Number > &linear_solver)
std::string & fe_family()
bool require_residual_reduction
std::unique_ptr< TimeSolver > time_solver
A pointer to the solver object we're going to use.
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
This class provides an interface to Eigen iterative solvers that is compatible with the libMesh Linea...
bool print_residual_norms
libMesh::Real relative_step_tolerance
std::string indicator_type
bool & analytic_jacobians()
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
The libMesh namespace provides an interface to certain functionality in the library.
This class implements writing meshes in the GMV format.
libMesh::Real refine_fraction
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
virtual std::unique_ptr< DiffSolver > & diff_solver()
An implicit linear or nonlinear solver to use at each timestep.
This is the MeshBase class.
std::vector< unsigned int > fe_order
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
virtual void forward_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector ¶meters, SensitivityData &sensitivities) override
Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with r...
SolverPackage default_solver_package()
libMesh::Real coarsen_threshold
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 &...)
Implements (adaptive) mesh refinement algorithms for a MeshBase.
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
Data structure for holding completed parameter sensitivity calculations.
virtual void adjoint_qoi_parameter_sensitivity(const QoISet &qoi_indices, const ParameterVector ¶meters, SensitivityData &sensitivities) override
Solves for the derivative of each of the system's quantities of interest q in qoi[qoi_indices] with r...
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
std::vector< std::string > fe_family
libMesh::Real verify_analytic_jacobians
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
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
libMesh::Real global_tolerance
void attach_qoi(DifferentiableQoI *qoi_in)
Attach external QoI object.
void set_solver_type(const SolverType st)
Sets the type of solver to use.
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
unsigned int & fe_order()
void add_command_line_names(const GetPot &getpot)
Merge a GetPot object's requested names into the set of queried command-line names.
virtual LinearSolver< Number > * get_linear_solver() const override
bool print_jacobian_norms
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
void adjust_linear_solvers(LaplaceSystem &system)
bool assemble_qoi_sides
If assemble_qoi_sides is true (it is false by default), the assembly loop for a quantity of interest ...
void write_timestep(const std::string &fname, const EquationSystems &es, const int timestep, const Real time, const std::set< std::string > *system_names=nullptr)
Writes out the solution at a specific timestep.
void write_output(EquationSystems &es, unsigned int a_step, std::string solution_type, FEMParameters ¶m)
void add_indices(const std::vector< unsigned int > &indices)
Add this indices to the set to be calculated.
const MeshBase & get_mesh() const
libMesh::Real min_step_length
std::size_t n_active_dofs() const
void all_second_order(const bool full_ordered=true)
Calls the range-based version of this function with a range consisting of all elements in the mesh...
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
unsigned int max_nonlinear_iterations
double linear_tolerance_multiplier
virtual void solve() override
Invokes the solver associated with the system.
virtual void init()
Initialize all the systems.
std::unique_ptr< MeshRefinement > build_mesh_refinement(MeshBase &mesh, FEMParameters ¶m)
virtual System & add_system(std::string_view system_type, std::string_view name)
Add the system of type system_type named name to the systems array.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
processor_id_type processor_id() const
ParameterVector & get_parameter_vector()
void set_system_parameters(LaplaceSystem &system, FEMParameters ¶m)
libMesh::Real coarsen_fraction
TimeSolver & get_time_solver()
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.