24 #include <unordered_map> 25 #include <unordered_set> 28 #include "libmesh/dof_map.h" 29 #include "libmesh/elem.h" 30 #include "libmesh/equation_systems.h" 31 #include "libmesh/error_vector.h" 32 #include "libmesh/fe.h" 33 #include "libmesh/fe_interface.h" 34 #include "libmesh/libmesh_common.h" 35 #include "libmesh/libmesh_logging.h" 36 #include "libmesh/mesh_base.h" 37 #include "libmesh/mesh_refinement.h" 38 #include "libmesh/numeric_vector.h" 39 #include "libmesh/quadrature.h" 40 #include "libmesh/system.h" 41 #include "libmesh/diff_physics.h" 42 #include "libmesh/fem_system.h" 43 #include "libmesh/implicit_system.h" 44 #include "libmesh/partitioner.h" 45 #include "libmesh/adjoint_refinement_estimator.h" 46 #include "libmesh/enum_error_estimator_type.h" 47 #include "libmesh/enum_norm_type.h" 48 #include "libmesh/utility.h" 50 #ifdef LIBMESH_ENABLE_AMR 75 number_h_refinements(1),
76 number_p_refinements(0),
77 _residual_evaluation_physics(nullptr),
78 _adjoint_evaluation_physics(nullptr),
139 for (
unsigned int j=0,
140 n_qois = system.
n_qois(); j != n_qois; j++)
161 system.
assembly(
true,
false,
false,
true);
164 coarse_residual = &system.
get_vector(
"RHS Vector");
165 coarse_residual->
close();
168 std::ostringstream liftfunc_name;
169 liftfunc_name <<
"adjoint_lift_function" << j;
178 (system.
get_vector(liftfunc_name.str()), static_cast<unsigned int>(j));
181 std::cout<<
"The flux QoI "<<
static_cast<unsigned int>(j)<<
" is: "<<coarse_residual->
dot(system.
get_vector(liftfunc_name.str()))<<std::endl<<std::endl;
191 std::map<std::string, std::unique_ptr<NumericVector<Number>>> coarse_vectors;
193 coarse_vectors[var_name] = vec->
clone();
196 std::unique_ptr<NumericVector<Number>> coarse_solution = system.
solution->clone();
204 std::vector<bool> old_adjoints_projection_settings(system.
n_qois());
211 auto old_adjoint_vector_projection_setting = system.
vector_preservation(adjoint_vector_name);
214 old_adjoints_projection_settings[j] = old_adjoint_vector_projection_setting;
222 bool old_projection_setting;
229 bool old_project_with_constraints_setting;
235 std::unique_ptr<Partitioner> old_partitioner(
mesh.partitioner().release());
238 const bool old_renumbering_setting =
mesh.allow_renumbering();
239 mesh.allow_renumbering(
false);
242 if (solution_vector && solution_vector != system.
solution.get())
252 error_per_cell.clear();
253 error_per_cell.resize (
mesh.max_elem_id(), 0.);
261 const unsigned int sysnum = system.
number();
262 for (
const auto * node :
mesh.local_node_ptr_range())
263 for (
unsigned int v=0, nvars = node->n_vars(sysnum);
265 if (node->n_comp(sysnum, v))
267 local_dof_bearing_nodes++;
279 if(!swapping_adjoint_physics)
300 std::vector<std::unique_ptr<NumericVector<Number>>> coarse_adjoints;
313 coarse_adjoints.emplace_back(std::move(coarse_adjoint));
316 coarse_adjoints.emplace_back(
nullptr);
329 system.
assembly(
true,
false,
false,
true);
337 projected_residual.
close();
376 std::ostringstream liftfunc_name;
377 liftfunc_name <<
"adjoint_lift_function" << j;
410 std::ostringstream liftfunc_name;
411 liftfunc_name <<
"adjoint_lift_function" << j;
465 std::unordered_map<dof_id_type, unsigned int> shared_element_count;
472 std::unordered_set<dof_id_type> processed_node_ids;
476 for (
const auto & elem :
mesh.active_local_element_ptr_range())
477 for (
const Node & node : elem->node_ref_range())
483 if (processed_node_ids.find(node_id) == processed_node_ids.end())
486 std::set<const Elem *> fine_grid_neighbor_set;
489 elem->find_point_neighbors(node, fine_grid_neighbor_set);
492 std::vector<dof_id_type> coarse_grid_neighbors;
495 for (
const auto & fine_elem : fine_grid_neighbor_set)
498 const Elem * coarse_elem = fine_elem;
503 coarse_elem = coarse_elem->
parent();
512 for (std::size_t cgns = coarse_grid_neighbors.size(); j != cgns; j++)
513 if (coarse_grid_neighbors[j] == coarse_id)
518 if (j == coarse_grid_neighbors.size())
519 coarse_grid_neighbors.push_back(coarse_id);
524 shared_element_count[node_id] =
525 cast_int<unsigned int>(coarse_grid_neighbors.size());
528 processed_node_ids.insert(node_id);
536 std::vector<dof_id_type> dof_indices;
562 for (
const auto & elem :
mesh.active_local_element_ptr_range())
565 const Elem * coarse = elem;
571 coarse = coarse->
parent();
580 Number local_contribution = 0.;
583 for (
const auto & dof : dof_indices)
584 local_contribution += (*localized_projected_residual)(dof) * (*localized_adjoint_solution)(dof);
587 local_contribution *= error_weight;
592 (std::abs(local_contribution));
624 libmesh_assert_equal_to (n_coarse_elem,
mesh.n_active_elem());
630 for (
const auto * node :
mesh.local_node_ptr_range())
631 for (
unsigned int v=0, nvars = node->n_vars(sysnum);
633 if (node->n_comp(sysnum, v))
635 final_local_dof_bearing_nodes++;
638 libmesh_assert_equal_to (local_dof_bearing_nodes,
639 final_local_dof_bearing_nodes);
657 *system.
solution = *coarse_solution;
663 const std::string & var_name = pr.first;
667 if (
auto it = coarse_vectors.find(var_name);
668 it != coarse_vectors.end())
678 mesh.partitioner().reset(old_partitioner.release());
679 mesh.allow_renumbering(old_renumbering_setting);
690 #endif // #ifdef LIBMESH_ENABLE_AMR virtual ErrorEstimatorType type() const override
vectors_iterator vectors_end()
End of vectors container.
void uniformly_p_refine(unsigned int n=1)
Uniformly p refines the mesh n times.
This is the EquationSystems class.
const Elem * parent() const
A Node is like a Point, but with more information.
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Data structure for specifying which Quantities of Interest should be calculated in an adjoint or a pa...
DifferentiablePhysics * _adjoint_evaluation_physics
Pointer to object to use for adjoint assembly.
unsigned int n_qois() const
Number of currently active quantities of interest.
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
ErrorEstimatorType
Defines an enum for the different types of error estimators which are available.
const EquationSystems & get_equation_systems() const
This is the base class from which all geometric element types are derived.
void uniformly_p_coarsen(unsigned int n=1)
Attempts to uniformly p coarsen the mesh n times.
const Parallel::Communicator & comm() const
NumericVector< Number > & add_vector(std::string_view vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
The libMesh namespace provides an interface to certain functionality in the library.
vectors_iterator vectors_begin()
Beginning of vectors container.
QoISet _qoi_set
A QoISet to handle cases with multiple QoIs available.
virtual T dot(const NumericVector< T > &v) const =0
dof_id_type n_local_dofs() const
virtual std::pair< unsigned int, Real > adjoint_solve(const QoISet &qoi_indices=QoISet()) override
Assembles & solves the linear system (dR/du)^T*z = dq/du, for those quantities of interest q specifie...
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
dof_id_type n_dofs() const
This is the MeshBase class.
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogeneously constrains the numeric vector v, which represents an adjoint solution defined on the ...
This class provides a specific system class.
This class handles the numbering of degrees of freedom on a mesh.
bool has_index(std::size_t) const
Return whether or not this index is in the set to be calculated.
unsigned int number() const
virtual void assembly(bool, bool, bool=false, bool=false)
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
void remove_vector(std::string_view vec_name)
Removes the additional vector vec_name from this system.
Implements (adaptive) mesh refinement algorithms for a MeshBase.
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
Manages consistently variables, degrees of freedom, and coefficient vectors.
virtual void reinit()
Handle any mesh changes and reinitialize all the systems on the updated mesh.
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
void uniformly_coarsen(unsigned int n=1)
Attempts to uniformly coarsen the mesh n times.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
This class holds functions that will estimate the error in a finite element solution on a given mesh...
void reduce_error(std::vector< ErrorVectorReal > &error_per_cell, const Parallel::Communicator &comm) const
This method takes the local error contributions in error_per_cell from each processor and combines th...
void set_project_with_constraints(bool _project_with_constraints)
DifferentiablePhysics * _residual_evaluation_physics
Pointer to object to use for physics assembly evaluations.
bool is_adjoint_already_solved() const
Accessor for the adjoint_already_solved boolean.
virtual void close()=0
Calls the NumericVector's internal assembly routines, ensuring that the values are consistent across ...
virtual void estimate_error(const System &system, ErrorVector &error_per_cell, const NumericVector< Number > *solution_vector=nullptr, bool estimate_parent_error=false) override
This function does uniform refinements and an adjoint solve to get an adjoint solution on each cell...
Real weight(std::size_t) const
Get the weight for this index (default 1.0)
std::vector< Number > computed_global_QoI_errors
const std::string & vector_name(const unsigned int vec_num) const
virtual void update()
Update the local values to reflect the solution on neighboring processors.
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...
virtual void swap(NumericVector< T > &v)
Swaps the contents of this with v.
virtual std::unique_ptr< DifferentiableQoI > clone() override
We don't allow systems to be attached to each other.
bool & project_solution_on_reinit(void)
Tells the System whether or not to project the solution vector onto new grids when the system is rein...
const MeshBase & get_mesh() const
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
NumericVector< Number > & get_adjoint_solution(unsigned int i=0)
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
virtual void clear()
Restores the NumericVector<T> to a pristine state.
AdjointRefinementEstimator()
Constructor.
bool get_project_with_constraints()
Setter and getter functions for project_with_constraints boolean.
bool vector_preservation(std::string_view vec_name) const
const DofMap & get_dof_map() const
unsigned char number_h_refinements
How many h refinements to perform to get the fine grid.
const std::vector< dof_id_type > & get_send_list() const
Manages consistently variables, degrees of freedom, and coefficient vectors for explicit systems...
const NumericVector< Number > & get_vector(std::string_view vec_name) const
unsigned char number_p_refinements
How many p refinements to perform to get the fine grid.
void uniformly_refine(unsigned int n=1)
Uniformly refines the mesh n times.
Manages consistently variables, degrees of freedom, coefficient vectors, and matrices for implicit sy...
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.