25 #include "libmesh/sparse_matrix.h" 26 #include "libmesh/numeric_vector.h" 27 #include "libmesh/dense_matrix.h" 28 #include "libmesh/dense_vector.h" 29 #include "libmesh/dof_map.h" 30 #include "libmesh/libmesh_logging.h" 31 #include "libmesh/equation_systems.h" 32 #include "libmesh/parallel.h" 33 #include "libmesh/parallel_algebra.h" 34 #include "libmesh/fe.h" 35 #include "libmesh/quadrature.h" 36 #include "libmesh/utility.h" 37 #include "libmesh/fe_interface.h" 38 #include "libmesh/fe_compute_data.h" 39 #include "libmesh/getpot.h" 40 #include "libmesh/exodusII_io.h" 41 #include "libmesh/fem_context.h" 42 #include "libmesh/elem.h" 43 #include "libmesh/int_range.h" 46 #include "libmesh/rb_construction_base.h" 47 #include "libmesh/rb_eim_construction.h" 48 #include "libmesh/rb_eim_evaluation.h" 49 #include "libmesh/rb_parameters.h" 50 #include "libmesh/rb_parametrized_function.h" 70 template <
typename DataMap>
71 void add(DataMap & u,
const Number k,
const DataMap & v)
73 for (
auto & [key, vec_vec_u] : u)
75 const std::vector<std::vector<Number>> & vec_vec_v = libmesh_map_find(v, key);
77 libmesh_error_msg_if (vec_vec_u.size() != vec_vec_v.size(),
"Size mismatch");
81 libmesh_error_msg_if (vec_vec_u[i].size() != vec_vec_v[i].size(),
"Size mismatch");
84 vec_vec_u[i][j] += k*vec_vec_v[i][j];
92 for (
auto & [key, vec_u] : u)
94 const std::vector<Number> & vec_v = libmesh_map_find(v, key);
96 libmesh_error_msg_if (vec_u.size() != vec_v.size(),
"Size mismatch");
100 vec_u[i] += k*vec_v[i];
108 const std::string & name_in,
109 const unsigned int number_in)
111 best_fit_type_flag(PROJECTION_BEST_FIT),
113 _set_Nmax_from_n_snapshots(false),
114 _Nmax_from_n_snapshots_increment(0),
115 _rel_training_tolerance(1.e-4),
116 _abs_training_tolerance(1.e-12),
117 _max_abs_value_in_training_set(0.),
118 _max_abs_value_in_training_set_index(0)
159 libmesh_error_msg_if(!
_rb_eim_eval,
"Error: RBEIMEvaluation object hasn't been initialized yet");
165 libmesh_error_msg_if(!
_rb_eim_eval,
"Error: RBEIMEvaluation object hasn't been initialized yet");
171 if (best_fit_type_string ==
"projection")
175 else if (best_fit_type_string ==
"eim")
179 else if (best_fit_type_string ==
"pod")
184 libmesh_error_msg(
"Error: invalid best_fit_type in input file");
190 libMesh::out << std::endl <<
"RBEIMConstruction parameters:" << std::endl;
195 libMesh::out <<
"Overruling Nmax based on number of snapshots, with increment set to " 216 libMesh::out <<
"EIM best fit type: projection" << std::endl;
234 GetPot infile(parameters_filename);
236 std::string best_fit_type_string = infile(
"best_fit_type",
"projection");
239 const unsigned int n_training_samples = infile(
"n_training_samples",0);
240 const bool deterministic_training = infile(
"deterministic_training",
false);
241 unsigned int training_parameters_random_seed_in =
242 static_cast<unsigned int>(-1);
243 training_parameters_random_seed_in = infile(
"training_parameters_random_seed",
244 training_parameters_random_seed_in);
245 const bool quiet_mode_in = infile(
"quiet_mode",
quiet_mode);
246 const unsigned int Nmax_in = infile(
"Nmax",
_Nmax);
247 const Real rel_training_tolerance_in = infile(
"rel_training_tolerance",
249 const Real abs_training_tolerance_in = infile(
"abs_training_tolerance",
253 unsigned int n_continuous_parameters = infile.vector_variable_size(
"parameter_names");
256 for (
unsigned int i=0; i<n_continuous_parameters; i++)
259 std::string param_name = infile(
"parameter_names",
"NONE", i);
262 Real min_val = infile(param_name, 0., 0);
263 mu_min_in.
set_value(param_name, min_val);
267 Real max_val = infile(param_name, 0., 1);
268 mu_max_in.
set_value(param_name, max_val);
272 std::map<std::string, std::vector<Real>> discrete_parameter_values_in;
274 unsigned int n_discrete_parameters = infile.vector_variable_size(
"discrete_parameter_names");
275 for (
unsigned int i=0; i<n_discrete_parameters; i++)
277 std::string param_name = infile(
"discrete_parameter_names",
"NONE", i);
279 unsigned int n_vals_for_param = infile.vector_variable_size(param_name);
280 std::vector<Real> vals_for_param(n_vals_for_param);
281 for (
auto j :
make_range(vals_for_param.size()))
282 vals_for_param[j] = infile(param_name, 0., j);
284 discrete_parameter_values_in[param_name] = vals_for_param;
287 std::map<std::string,bool> log_scaling_in;
291 for (
const auto & pr : mu_min_in)
292 log_scaling_in[pr.first] =
false;
296 deterministic_training,
297 training_parameters_random_seed_in,
300 rel_training_tolerance_in,
301 abs_training_tolerance_in,
304 discrete_parameter_values_in,
309 bool deterministic_training_in,
310 int training_parameters_random_seed_in,
312 unsigned int Nmax_in,
313 Real rel_training_tolerance_in,
314 Real abs_training_tolerance_in,
317 const std::map<std::string, std::vector<Real>> & discrete_parameter_values_in,
318 const std::map<std::string,bool> & log_scaling_in,
319 std::map<std::string, std::vector<RBParameter>> * training_sample_list)
337 const std::string & lookup_table_param_name =
340 libmesh_error_msg_if(!discrete_parameter_values_in.count(lookup_table_param_name),
341 "Lookup table parameter should be discrete");
344 std::map<std::string, std::vector<Real>> discrete_parameter_values_final(
345 discrete_parameter_values_in);
347 std::vector<Real> & lookup_table_param_values =
348 libmesh_map_find(discrete_parameter_values_final, lookup_table_param_name);
352 std::iota(lookup_table_param_values.begin(), lookup_table_param_values.end(), 0);
357 n_training_samples_in = lookup_table_param_values.size();
368 bool updated_deterministic_training = deterministic_training_in;
378 updated_deterministic_training =
false;
383 n_training_samples_in,
385 updated_deterministic_training);
387 if (training_sample_list)
399 const std::string & lookup_table_param_name =
405 std::vector<RBParameter> lookup_table_training_samples(n_training_samples_in, {val});
406 for (
auto & vec : lookup_table_training_samples)
434 LOG_SCOPE(
"train_eim_approximation_with_greedy()",
"RBEIMConstruction");
450 "Error: We currently only support EIM training starting from an empty basis");
452 libMesh::out << std::endl <<
"---- Performing Greedy EIM basis enrichment ----" << std::endl;
456 std::cout <<
"Maximum absolute value in the training set is " 464 Real greedy_error = -1.;
465 std::vector<RBParameters> greedy_param_list;
476 bool exit_on_next_iteration =
false;
481 bool bfs_equals_n_samples =
false;
488 <<
") equals number of training samples." << std::endl;
490 bfs_equals_n_samples =
true;
494 if (!exit_on_next_iteration)
498 libMesh::out <<
"Greedily selected parameter vector:" << std::endl;
502 libMesh::out <<
"Enriching the EIM approximation" << std::endl;
505 bool is_zero_bf = bfs_equals_n_samples || (greedy_error == 0.);
512 std::unique_ptr<EimPointData> eim_point_data;
520 !exit_on_next_iteration,
521 eim_point_data.get());
536 libMesh::out <<
"Computing EIM error on training set" << std::endl;
540 libMesh::out <<
"Maximum EIM error is " << greedy_error << std::endl << std::endl;
542 #ifdef LIBMESH_ENABLE_EXCEPTIONS 543 catch (
const std::exception & e)
547 if (exit_on_next_iteration)
549 std::cout <<
"Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
557 if (exit_on_next_iteration)
559 libMesh::out <<
"Extra EIM iteration for error indicator is complete, hence exiting EIM training now" << std::endl;
565 bool exit_condition_satisfied =
false;
569 libMesh::out <<
"Maximum number of basis functions reached: Nmax = " 571 exit_condition_satisfied =
true;
576 if (!exit_condition_satisfied)
579 libMesh::out <<
"Relative error tolerance reached." << std::endl;
580 exit_condition_satisfied =
true;
583 if (!exit_condition_satisfied)
586 libMesh::out <<
"Absolute error tolerance reached." << std::endl;
587 exit_condition_satisfied =
true;
591 if (!exit_condition_satisfied)
593 bool do_exit =
false;
597 for (
auto & param : greedy_param_list)
600 libMesh::out <<
"Exiting greedy because the same parameters were selected twice" 607 exit_condition_satisfied =
true;
610 if (exit_condition_satisfied)
618 exit_on_next_iteration =
true;
619 libMesh::out <<
"EIM error indicator is active, hence we will run one extra EIM iteration before exiting" 641 LOG_SCOPE(
"apply_normalization_to_solution_snapshots()",
"RBEIMConstruction");
643 libMesh::out <<
"Normalizing solution snapshots" << std::endl;
649 for (
unsigned int i=0; i<n_snapshots; i++)
657 apply_comp_scaling)));
667 apply_comp_scaling)));
677 apply_comp_scaling)));
699 libMesh::out <<
"Maximum absolute value in the training set after normalization: " 705 LOG_SCOPE(
"train_eim_approximation_with_POD()",
"RBEIMConstruction");
718 if (updated_Nmax > 0)
719 _Nmax =
static_cast<unsigned int>(updated_Nmax);
733 "Error: We currently only support EIM training starting from an empty basis");
735 libMesh::out << std::endl <<
"---- Performing POD EIM basis enrichment ----" << std::endl;
742 std::cout <<
"Start computing correlation matrix" << std::endl;
745 for (
unsigned int i=0; i<n_snapshots; i++)
747 for (
unsigned int j=0; j<=i; j++)
773 correlation_matrix(i,j) = inner_prod;
781 if ( (i+1) % 10 == 0)
782 std::cout <<
"Finished row " << (i+1) <<
" of " << n_snapshots << std::endl;
784 std::cout <<
"Finished computing correlation matrix" << std::endl;
802 correlation_matrix.
svd(sigma, U, VT );
807 bool exit_on_next_iteration =
false;
816 bool j_equals_n_snapshots =
false;
819 bool exit_condition_satisfied =
false;
821 if ((j == 0) && (sigma(0) == 0.))
823 libMesh::out <<
"Terminating EIM POD with empty basis because first singular value is zero" << std::endl;
824 exit_condition_satisfied =
true;
827 else if (j >= n_snapshots)
829 libMesh::out <<
"Number of basis functions equals number of training samples." << std::endl;
830 exit_condition_satisfied =
true;
831 j_equals_n_snapshots =
true;
843 rel_err = std::sqrt(sigma(j)) / std::sqrt(sigma(0));
846 if (exit_on_next_iteration)
848 libMesh::out <<
"Extra EIM iteration for error indicator is complete, POD error norm for extra iteration: " << rel_err << std::endl;
853 <<
", POD error norm: " << rel_err << std::endl;
855 if (!exit_condition_satisfied)
858 libMesh::out <<
"Maximum number of basis functions (" << j <<
") reached." << std::endl;
859 exit_condition_satisfied =
true;
862 if (!exit_condition_satisfied)
865 libMesh::out <<
"Training tolerance reached." << std::endl;
866 exit_condition_satisfied =
true;
869 if (exit_condition_satisfied)
878 exit_on_next_iteration =
true;
879 libMesh::out <<
"EIM error indicator is active, hence we will run one extra EIM iteration before exiting" 886 bool is_zero_bf = j_equals_n_snapshots || (rel_err == 0.);
896 for (
unsigned int i=0; i<n_snapshots; ++i )
899 Real norm_v = std::sqrt(sigma(j));
910 std::unique_ptr<EimPointData> eim_point_data;
918 !exit_on_next_iteration,
919 eim_point_data.get());
925 if (is_linearly_dependent)
930 exit_on_next_iteration =
true;
931 libMesh::out <<
"Linearly dependent data detected, finalizing iteration for the EIM error indicator before exiting" 938 if (is_linearly_dependent && !is_zero_bf)
949 !exit_on_next_iteration,
950 eim_point_data.get());
955 #ifdef LIBMESH_ENABLE_EXCEPTIONS 956 catch (
const std::exception & e)
960 if (exit_on_next_iteration)
962 std::cout <<
"Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
979 for (
unsigned int i=0; i<n_snapshots; ++i )
982 Real norm_v = std::sqrt(sigma(j));
993 std::unique_ptr<EimPointData> eim_point_data;
1001 !exit_on_next_iteration,
1002 eim_point_data.get());
1008 if (is_linearly_dependent)
1013 exit_on_next_iteration =
true;
1014 libMesh::out <<
"Linearly dependent data detected, finalizing iteration for the EIM error indicator before exiting" 1021 if (is_linearly_dependent && !is_zero_bf)
1032 !exit_on_next_iteration,
1033 eim_point_data.get());
1038 #ifdef LIBMESH_ENABLE_EXCEPTIONS 1039 catch (
const std::exception & e)
1043 if (exit_on_next_iteration)
1045 std::cout <<
"Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
1062 for (
unsigned int i=0; i<n_snapshots; ++i )
1065 Real norm_v = std::sqrt(sigma(j));
1076 std::unique_ptr<EimPointData> eim_point_data;
1084 !exit_on_next_iteration,
1085 eim_point_data.get());
1091 if (is_linearly_dependent)
1096 exit_on_next_iteration =
true;
1097 libMesh::out <<
"Linearly dependent data detected, finalizing iteration for the EIM error indicator before exiting" 1104 if (is_linearly_dependent && !is_zero_bf)
1115 !exit_on_next_iteration,
1116 eim_point_data.get());
1121 #ifdef LIBMESH_ENABLE_EXCEPTIONS 1122 catch (
const std::exception & e)
1126 if (exit_on_next_iteration)
1128 std::cout <<
"Exception occurred when enriching basis for error indicator hence we skip the error indicator in this case" << std::endl;
1142 std::cout <<
"Zero basis function encountered, hence exiting." << std::endl;
1175 if (
mesh.elem_dimensions().count(
dim))
1239 LOG_SCOPE(
"store_eim_solutions_for_training_set()",
"RBEIMConstruction");
1244 eim_solutions.clear();
1261 for (
unsigned int j=0; j<RB_size; j++)
1284 for (
unsigned int j=0; j<RB_size; j++)
1305 for (
unsigned int j=0; j<RB_size; j++)
1323 "Invalid index: " << training_index);
1330 "Invalid index: " << training_index);
1337 "Invalid index: " << training_index);
1373 LOG_SCOPE(
"compute_max_eim_error()",
"RBEIMConstruction");
1378 return std::make_pair(0.,0);
1382 unsigned int max_err_index = 0;
1386 "Error: Training samples should be the same on all procs");
1403 for (
unsigned int i=0; i<RB_size; i++)
1415 RB_inner_product_matrix_N.
lu_solve(best_fit_rhs, best_fit_coeffs);
1420 if (best_fit_error > max_err)
1422 max_err_index = training_index;
1423 max_err = best_fit_error;
1435 for (
unsigned int i=0; i<RB_size; i++)
1447 RB_inner_product_matrix_N.
lu_solve(best_fit_rhs, best_fit_coeffs);
1452 if (best_fit_error > max_err)
1454 max_err_index = training_index;
1455 max_err = best_fit_error;
1467 for (
unsigned int i=0; i<RB_size; i++)
1479 RB_inner_product_matrix_N.
lu_solve(best_fit_rhs, best_fit_coeffs);
1484 if (best_fit_error > max_err)
1486 max_err_index = training_index;
1487 max_err = best_fit_error;
1516 if (best_fit_error > max_err)
1518 max_err_index = training_index;
1519 max_err = best_fit_error;
1528 if (best_fit_error > max_err)
1530 max_err_index = training_index;
1531 max_err = best_fit_error;
1540 if (best_fit_error > max_err)
1542 max_err_index = training_index;
1543 max_err = best_fit_error;
1550 libmesh_error_msg(
"EIM best fit type not recognized");
1553 return std::make_pair(max_err,max_err_index);
1558 LOG_SCOPE(
"initialize_parametrized_functions_in_training_set()",
"RBEIMConstruction");
1561 "Error: We must have serial_training_set==true in " 1562 "RBEIMConstruction::initialize_parametrized_functions_in_training_set");
1564 libMesh::out <<
"Initializing parametrized functions in training set..." << std::endl;
1586 std::vector<Real> max_abs_value_per_component_in_training_set(n_comps);
1593 libMesh::out <<
"Initializing parametrized function for training sample " 1608 std::vector<std::vector<Number>> comps_and_qps(n_comps);
1609 for (
unsigned int comp=0; comp<n_comps; comp++)
1611 comps_and_qps[comp].resize(xyz_vector.size());
1616 elem_side_pair.first,
1617 elem_side_pair.second,
1619 comps_and_qps[comp][qp] =
value;
1628 if (abs_value > max_abs_value_per_component_in_training_set[comp])
1629 max_abs_value_per_component_in_training_set[comp] = abs_value;
1637 libMesh::out <<
"Parametrized functions in training set initialized" << std::endl;
1639 unsigned int max_id = 0;
1642 libMesh::out <<
"Maximum absolute value in the training set: " 1647 comm().
max(max_abs_value_per_component_in_training_set);
1656 max_abs_value_per_component_in_training_set[i] == 0.)
1667 libMesh::out <<
"Initializing parametrized function for training sample " 1679 const auto & node_id = pr.first;
1681 std::vector<Number> comps(n_comps);
1682 for (
unsigned int comp=0; comp<n_comps; comp++)
1687 comps[comp] =
value;
1696 if (abs_value > max_abs_value_per_component_in_training_set[comp])
1697 max_abs_value_per_component_in_training_set[comp] = abs_value;
1704 libMesh::out <<
"Parametrized functions in training set initialized" << std::endl;
1706 unsigned int max_id = 0;
1709 libMesh::out <<
"Maximum absolute value in the training set: " 1714 comm().
max(max_abs_value_per_component_in_training_set);
1723 max_abs_value_per_component_in_training_set[i] == 0.)
1734 libMesh::out <<
"Initializing parametrized function for training sample " 1747 std::vector<std::vector<Number>> comps_and_qps(n_comps);
1748 for (
unsigned int comp=0; comp<n_comps; comp++)
1750 comps_and_qps[comp].resize(xyz_vector.size());
1755 comps_and_qps[comp][qp] =
value;
1764 if (abs_value > max_abs_value_per_component_in_training_set[comp])
1765 max_abs_value_per_component_in_training_set[comp] = abs_value;
1773 libMesh::out <<
"Parametrized functions in training set initialized" << std::endl;
1775 unsigned int max_id = 0;
1778 libMesh::out <<
"Maximum absolute value in the training set: " 1783 comm().
max(max_abs_value_per_component_in_training_set);
1792 max_abs_value_per_component_in_training_set[i] == 0.)
1806 LOG_SCOPE(
"initialize_qp_data()",
"RBEIMConstruction");
1810 libMesh::out <<
"Initializing quadrature point locations" << std::endl;
1814 libMesh::out <<
"Initializing quadrature point and perturbation locations" << std::endl;
1825 const std::set<boundary_id_type> & parametrized_function_boundary_ids =
1827 libmesh_error_msg_if (parametrized_function_boundary_ids.empty(),
1828 "Need to have non-empty boundary IDs to initialize side data");
1839 const auto & binfo =
mesh.get_boundary_info();
1840 std::vector<boundary_id_type> side_boundary_ids;
1842 for (
const auto & elem :
mesh.active_local_element_ptr_range())
1850 const std::vector<Real> & JxW = elem_fe->get_JxW();
1851 const std::vector<Point> & xyz = elem_fe->get_xyz();
1854 auto side_fe = context.
get_side_fe(0, elem->dim());
1855 const std::vector<Real> & JxW_side = side_fe->get_JxW();
1856 const std::vector< Point > & xyz_side = side_fe->get_xyz();
1858 for (context.
side = 0;
1865 binfo.boundary_ids(elem, context.
side, side_boundary_ids);
1867 bool has_side_boundary_id =
false;
1870 if(parametrized_function_boundary_ids.count(side_boundary_id))
1872 has_side_boundary_id =
true;
1873 matching_boundary_id = side_boundary_id;
1877 if(has_side_boundary_id)
1881 auto elem_side_pair = std::make_pair(elem_id, context.
side);
1895 std::vector<std::vector<Point>> xyz_perturb_vec_at_qps;
1897 for (
const Point & xyz_qp : xyz_side)
1899 std::vector<Point> xyz_perturb_vec;
1900 if (elem->dim() == 3)
1911 std::unique_ptr<const Elem> elem_side;
1912 elem->build_side_ptr(elem_side, context.
side);
1924 Point xi_eta_perturb = xi_eta;
1926 xi_eta_perturb(0) += fd_delta;
1927 Point xyz_perturb_0 =
1931 xi_eta_perturb(0) -= fd_delta;
1933 xi_eta_perturb(1) += fd_delta;
1934 Point xyz_perturb_1 =
1938 xi_eta_perturb(1) -= fd_delta;
1943 Point unit_0 = (xyz_perturb_0-xyz_qp).unit();
1944 Point unit_1 = (xyz_perturb_1-xyz_qp).unit();
1946 xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_0);
1947 xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_1);
1956 xyz_perturb_vec_at_qps.emplace_back(xyz_perturb_vec);
1966 if (elem->dim() == 2)
1967 for (
unsigned int shellface_index=0; shellface_index<2; shellface_index++)
1969 binfo.shellface_boundary_ids(elem, shellface_index, side_boundary_ids);
1971 bool has_side_boundary_id =
false;
1974 if(parametrized_function_boundary_ids.count(side_boundary_id))
1976 has_side_boundary_id =
true;
1977 matching_boundary_id = side_boundary_id;
1981 if(has_side_boundary_id)
1987 auto elem_side_pair = std::make_pair(elem_id, shellface_index);
2001 std::vector<std::vector<Point>> xyz_perturb_vec_at_qps;
2003 for (
const Point & xyz_qp : xyz)
2005 std::vector<Point> xyz_perturb_vec;
2019 Point xi_eta_perturb = xi_eta;
2021 xi_eta_perturb(0) += fd_delta;
2022 Point xyz_perturb_0 =
2026 xi_eta_perturb(0) -= fd_delta;
2028 xi_eta_perturb(1) += fd_delta;
2029 Point xyz_perturb_1 =
2033 xi_eta_perturb(1) -= fd_delta;
2038 Point unit_0 = (xyz_perturb_0-xyz_qp).unit();
2039 Point unit_1 = (xyz_perturb_1-xyz_qp).unit();
2041 xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_0);
2042 xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_1);
2045 xyz_perturb_vec_at_qps.emplace_back(xyz_perturb_vec);
2056 const std::set<boundary_id_type> & parametrized_function_boundary_ids =
2058 libmesh_error_msg_if (parametrized_function_boundary_ids.empty(),
2059 "Need to have non-empty boundary IDs to initialize node data");
2064 const auto & binfo =
mesh.get_boundary_info();
2070 std::set<dof_id_type> nodes_with_nodesets;
2071 for (
const auto & t : binfo.build_node_list())
2072 nodes_with_nodesets.insert(std::get<0>(t));
2075 std::vector<boundary_id_type> node_boundary_ids;
2079 const Node * node =
mesh.node_ptr(node_id);
2084 binfo.boundary_ids(node, node_boundary_ids);
2086 bool has_node_boundary_id =
false;
2091 has_node_boundary_id =
true;
2096 if(has_node_boundary_id)
2111 for (
const auto & elem :
mesh.active_local_element_ptr_range())
2114 const std::vector<Real> & JxW = elem_fe->get_JxW();
2115 const std::vector<Point> & xyz = elem_fe->get_xyz();
2130 std::vector<std::vector<Point>> xyz_perturb_vec_at_qps;
2132 for (
const Point & xyz_qp : xyz)
2134 std::vector<Point> xyz_perturb_vec;
2135 if (elem->dim() == 3)
2137 Point xyz_perturb = xyz_qp;
2139 xyz_perturb(0) += fd_delta;
2140 xyz_perturb_vec.emplace_back(xyz_perturb);
2141 xyz_perturb(0) -= fd_delta;
2143 xyz_perturb(1) += fd_delta;
2144 xyz_perturb_vec.emplace_back(xyz_perturb);
2145 xyz_perturb(1) -= fd_delta;
2147 xyz_perturb(2) += fd_delta;
2148 xyz_perturb_vec.emplace_back(xyz_perturb);
2149 xyz_perturb(2) -= fd_delta;
2151 else if (elem->dim() == 2)
2172 Point xi_eta_perturb = xi_eta;
2174 xi_eta_perturb(0) += fd_delta;
2175 Point xyz_perturb_0 =
2179 xi_eta_perturb(0) -= fd_delta;
2181 xi_eta_perturb(1) += fd_delta;
2182 Point xyz_perturb_1 =
2186 xi_eta_perturb(1) -= fd_delta;
2191 Point unit_0 = (xyz_perturb_0-xyz_qp).unit();
2192 Point unit_1 = (xyz_perturb_1-xyz_qp).unit();
2194 xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_0);
2195 xyz_perturb_vec.emplace_back(xyz_qp + fd_delta*unit_1);
2205 xyz_perturb_vec_at_qps.emplace_back(xyz_perturb_vec);
2217 LOG_SCOPE(
"inner_product()",
"RBEIMConstruction");
2221 for (
const auto & [elem_id, v_comp_and_qp] : v)
2223 const auto & w_comp_and_qp = libmesh_map_find(w, elem_id);
2226 for (
const auto & comp :
index_range(v_comp_and_qp))
2228 const std::vector<Number> & v_qp = v_comp_and_qp[comp];
2229 const std::vector<Number> & w_qp = w_comp_and_qp[comp];
2231 Real comp_scaling = 1.;
2232 if (apply_comp_scaling)
2240 val += JxW[qp] * comp_scaling * v_qp[qp] *
libmesh_conj(w_qp[qp]);
2251 LOG_SCOPE(
"side_inner_product()",
"RBEIMConstruction");
2255 for (
const auto & [elem_and_side, v_comp_and_qp] : v)
2257 const auto & w_comp_and_qp = libmesh_map_find(w, elem_and_side);
2260 for (
const auto & comp :
index_range(v_comp_and_qp))
2262 const std::vector<Number> & v_qp = v_comp_and_qp[comp];
2263 const std::vector<Number> & w_qp = w_comp_and_qp[comp];
2265 Real comp_scaling = 1.;
2266 if (apply_comp_scaling)
2274 val += JxW[qp] * comp_scaling * v_qp[qp] *
libmesh_conj(w_qp[qp]);
2285 LOG_SCOPE(
"node_inner_product()",
"RBEIMConstruction");
2289 for (
const auto & [node_id, v_comps] : v)
2291 const auto & w_comps = libmesh_map_find(w, node_id);
2298 Real comp_scaling = 1.;
2299 if (apply_comp_scaling)
2306 val += comp_scaling * v_comps[comp] *
libmesh_conj(w_comps[comp]);
2316 Real max_value = 0.;
2318 for (
const auto & pr : v)
2320 const auto & values = pr.second;
2323 const auto &
value = values[comp];
2325 Real comp_scaling = 1.;
2330 "Invalid vector index");
2334 max_value = std::max(max_value, std::abs(
value * comp_scaling));
2343 bool add_basis_function,
2346 LOG_SCOPE(
"enrich_eim_approximation()",
"RBEIMConstruction");
2369 bool add_basis_function,
2388 for (
unsigned int i=0; i<RB_size; i++)
2409 Number optimal_value = 0.;
2410 Point optimal_point;
2411 unsigned int optimal_comp = 0;
2413 unsigned int optimal_side_index = 0;
2416 unsigned int optimal_qp = 0;
2417 std::vector<Point> optimal_point_perturbs;
2418 std::vector<Real> optimal_point_phi_i_qp;
2421 Real largest_abs_value = -1.;
2427 for (
const auto & [elem_and_side, comp_and_qp] : local_pf)
2430 unsigned int side_index = elem_and_side.second;
2437 std::vector<std::vector<Real>> phi;
2446 side_fe->reinit(&elem_ref, side_index);
2448 phi = side_fe->get_phi();
2450 else if (side_type == 1)
2455 phi = elem_fe->get_phi();
2458 libmesh_error_msg (
"Unrecognized side_type: " << side_type);
2460 for (
const auto & comp :
index_range(comp_and_qp))
2462 const std::vector<Number> & qp_values = comp_and_qp[comp];
2469 bool update_optimal_point =
false;
2470 if (!eim_point_data)
2471 update_optimal_point = (abs_value > largest_abs_value);
2473 update_optimal_point = (elem_id == eim_point_data->
elem_id) &&
2474 (side_index == eim_point_data->
side_index) &&
2478 if (update_optimal_point)
2480 largest_abs_value = abs_value;
2481 optimal_value =
value;
2482 optimal_comp = comp;
2483 optimal_elem_id = elem_id;
2484 optimal_side_index = side_index;
2487 optimal_point_phi_i_qp.resize(phi.size());
2489 optimal_point_phi_i_qp[i] = phi[i][qp];
2491 const auto & point_list =
2494 libmesh_error_msg_if(qp >= point_list.size(),
"Error: Invalid qp");
2496 optimal_point = point_list[qp];
2503 const auto & perturb_list =
2506 libmesh_error_msg_if(qp >= perturb_list.size(),
"Error: Invalid qp");
2508 optimal_point_perturbs = perturb_list[qp];
2517 unsigned int proc_ID_index;
2518 this->
comm().
maxloc(largest_abs_value, proc_ID_index);
2525 this->
comm().
broadcast(optimal_subdomain_id, proc_ID_index);
2528 this->
comm().
broadcast(optimal_point_perturbs, proc_ID_index);
2529 this->
comm().
broadcast(optimal_point_phi_i_qp, proc_ID_index);
2533 if (add_basis_function)
2535 if (optimal_value == 0.)
2537 libMesh::out <<
"Encountered linearly dependent data in EIM enrichment, hence skip adding new basis function" << std::endl;
2553 optimal_subdomain_id,
2554 optimal_boundary_id,
2556 optimal_point_perturbs,
2557 optimal_point_phi_i_qp);
2564 bool add_basis_function,
2583 for (
unsigned int i=0; i<RB_size; i++)
2602 Number optimal_value = 0.;
2603 Point optimal_point;
2604 unsigned int optimal_comp = 0;
2609 Real largest_abs_value = -1.;
2611 for (
const auto & [node_id, values] : local_pf)
2618 bool update_optimal_point =
false;
2619 if (!eim_point_data)
2620 update_optimal_point = (abs_value > largest_abs_value);
2622 update_optimal_point = (node_id == eim_point_data->
node_id) &&
2625 if (update_optimal_point)
2627 largest_abs_value = abs_value;
2628 optimal_value =
value;
2629 optimal_comp = comp;
2630 optimal_node_id = node_id;
2641 unsigned int proc_ID_index;
2642 this->
comm().
maxloc(largest_abs_value, proc_ID_index);
2652 if (add_basis_function)
2654 if (optimal_value == 0.)
2656 libMesh::out <<
"Encountered linearly dependent data in EIM enrichment, hence skip adding new basis function" << std::endl;
2671 optimal_boundary_id);
2678 bool add_basis_function,
2697 for (
unsigned int i=0; i<RB_size; i++)
2717 Number optimal_value = 0.;
2718 Point optimal_point;
2719 unsigned int optimal_comp = 0;
2722 unsigned int optimal_qp = 0;
2723 std::vector<Point> optimal_point_perturbs;
2724 std::vector<Real> optimal_point_phi_i_qp;
2726 std::vector<Real> optimal_JxW_all_qp;
2727 std::vector<std::vector<Real>> optimal_phi_i_all_qp;
2729 Point optimal_dxyzdxi_elem_center;
2730 Point optimal_dxyzdeta_elem_center;
2733 Real largest_abs_value = -1.;
2746 for (
const auto & [elem_id, comp_and_qp] : local_pf)
2753 const std::vector<std::vector<Real>> & phi = elem_fe->get_phi();
2754 const auto & JxW = elem_fe->get_JxW();
2755 const auto & dxyzdxi = elem_fe->get_dxyzdxi();
2756 const auto & dxyzdeta = elem_fe->get_dxyzdeta();
2758 elem_fe->reinit(&elem_ref);
2760 for (
const auto & comp :
index_range(comp_and_qp))
2762 const std::vector<Number> & qp_values = comp_and_qp[comp];
2772 bool update_optimal_point =
false;
2773 if (!eim_point_data)
2774 update_optimal_point = (abs_value > largest_abs_value);
2776 update_optimal_point = (elem_id == eim_point_data->
elem_id) &&
2780 if (update_optimal_point)
2782 largest_abs_value = abs_value;
2783 optimal_value =
value;
2784 optimal_comp = comp;
2785 optimal_elem_id = elem_id;
2787 optimal_elem_type = elem_ref.
type();
2789 optimal_point_phi_i_qp.resize(phi.size());
2791 optimal_point_phi_i_qp[i] = phi[i][qp];
2793 const auto & point_list =
2796 libmesh_error_msg_if(qp >= point_list.size(),
"Error: Invalid qp");
2798 optimal_point = point_list[qp];
2804 const auto & perturb_list =
2807 libmesh_error_msg_if(qp >= perturb_list.size(),
"Error: Invalid qp");
2809 optimal_point_perturbs = perturb_list[qp];
2814 optimal_JxW_all_qp = JxW;
2815 optimal_phi_i_all_qp = phi;
2823 elem_fe->reinit (&elem_ref, &nodes);
2825 Point dxyzdxi_pt, dxyzdeta_pt;
2827 dxyzdxi_pt = dxyzdxi[0];
2829 dxyzdeta_pt = dxyzdeta[0];
2831 optimal_dxyzdxi_elem_center = dxyzdxi_pt;
2832 optimal_dxyzdeta_elem_center = dxyzdeta_pt;
2834 elem_fe->reinit(&elem_ref);
2843 unsigned int proc_ID_index;
2844 this->
comm().
maxloc(largest_abs_value, proc_ID_index);
2850 this->
comm().
broadcast(optimal_subdomain_id, proc_ID_index);
2852 this->
comm().
broadcast(optimal_point_perturbs, proc_ID_index);
2853 this->
comm().
broadcast(optimal_point_phi_i_qp, proc_ID_index);
2855 this->
comm().
broadcast(optimal_phi_i_all_qp, proc_ID_index);
2856 this->
comm().
broadcast(optimal_dxyzdxi_elem_center, proc_ID_index);
2857 this->
comm().
broadcast(optimal_dxyzdeta_elem_center, proc_ID_index);
2861 int optimal_elem_type_int =
static_cast<int>(optimal_elem_type);
2862 this->
comm().
broadcast(optimal_elem_type_int, proc_ID_index);
2863 optimal_elem_type =
static_cast<ElemType>(optimal_elem_type_int);
2868 int optimal_qrule_order_int =
static_cast<int>(optimal_qrule_order);
2869 this->
comm().
broadcast(optimal_qrule_order_int, proc_ID_index);
2870 optimal_qrule_order =
static_cast<Order>(optimal_qrule_order_int);
2875 if (add_basis_function)
2877 if (optimal_value == 0.)
2879 libMesh::out <<
"Encountered linearly dependent data in EIM enrichment, hence skip adding new basis function" << std::endl;
2894 optimal_subdomain_id,
2896 optimal_point_perturbs,
2897 optimal_point_phi_i_qp,
2900 optimal_phi_i_all_qp,
2901 optimal_qrule_order,
2902 optimal_dxyzdxi_elem_center,
2903 optimal_dxyzdeta_elem_center);
2911 LOG_SCOPE(
"update_eim_matrices()",
"RBEIMConstruction");
2916 libmesh_assert_msg(RB_size >= 1,
"Must have at least 1 basis function.");
2918 if (set_eim_error_indicator)
2929 for (
unsigned int j=0; j<RB_size; j++)
2939 extra_point_row(j) =
value;
2945 for (
unsigned int j=0; j<RB_size; j++)
2953 extra_point_row(j) =
value;
2959 for (
unsigned int j=0; j<RB_size; j++)
2968 extra_point_row(j) =
value;
2980 for (
unsigned int i=(RB_size-1); i<RB_size; i++)
2982 for (
unsigned int j=0; j<RB_size; j++)
2998 for (
unsigned int j=0; j<RB_size; j++)
3015 for (
unsigned int i=(RB_size-1); i<RB_size; i++)
3017 for (
unsigned int j=0; j<RB_size; j++)
3033 for (
unsigned int j=0; j<RB_size; j++)
3048 for (
unsigned int i=(RB_size-1); i<RB_size; i++)
3050 for (
unsigned int j=0; j<RB_size; j++)
3066 for (
unsigned int j=0; j<RB_size; j++)
3083 for (
auto & pr : local_pf)
3085 auto & values = pr.second;
3086 for (
auto &
value : values)
3087 value *= scaling_factor;
3098 std::default_random_engine gen;
3099 std::uniform_int_distribution<> dist{0,
static_cast<int>(n)};
3114 if (
comm().size() > 1)
3129 bool error_finding_new_element =
false;
3130 if (
comm().rank() == 0)
3135 std::set<dof_id_type> previous_elem_ids(vec_eval_input.
elem_ids.begin(), vec_eval_input.
elem_ids.end());
3151 std::set<dof_id_type> new_elem_ids;
3152 for (
const auto & v_pair : *v_ptr)
3153 if (previous_elem_ids.count(v_pair.first) == 0)
3154 new_elem_ids.insert(v_pair.first);
3160 error_finding_new_element = (new_elem_ids.empty());
3162 if (!error_finding_new_element)
3166 auto item = new_elem_ids.begin();
3167 std::advance(item, random_elem_idx);
3168 eim_point_data.
elem_id = *item;
3172 if (!error_finding_new_element)
3175 const auto & vars_and_qps = libmesh_map_find(*v_ptr,eim_point_data.
elem_id);
3180 const auto & qps = libmesh_map_find(*v_ptr,eim_point_data.
elem_id)[eim_point_data.
comp_index];
3187 libmesh_error_msg_if(error_finding_new_element,
"Could not find new element in get_random_point()");
3194 return eim_point_data;
3208 if (
comm().size() > 1)
3223 bool error_finding_new_element_and_side =
false;
3224 if (
comm().rank() == 0)
3228 std::pair<dof_id_type,unsigned int> elem_and_side;
3230 std::set<std::pair<dof_id_type,unsigned int>> previous_elem_and_side_ids;
3233 previous_elem_and_side_ids.insert(
3240 std::set<std::pair<dof_id_type,unsigned int>> new_elem_and_side_ids;
3241 for (
const auto & v_pair : *v_ptr)
3242 if (previous_elem_and_side_ids.count(v_pair.first) == 0)
3243 new_elem_and_side_ids.insert(v_pair.first);
3249 error_finding_new_element_and_side = (new_elem_and_side_ids.empty());
3251 if (!error_finding_new_element_and_side)
3255 auto item = new_elem_and_side_ids.begin();
3256 std::advance(item, random_elem_and_side_idx);
3257 elem_and_side = *item;
3258 eim_point_data.
elem_id = elem_and_side.first;
3259 eim_point_data.
side_index = elem_and_side.second;
3263 if (!error_finding_new_element_and_side)
3266 const auto & vars_and_qps = libmesh_map_find(*v_ptr,elem_and_side);
3271 const auto & qps = libmesh_map_find(*v_ptr,elem_and_side)[eim_point_data.
comp_index];
3278 libmesh_error_msg_if(error_finding_new_element_and_side,
"Could not find new (element,side) in get_random_point()");
3286 return eim_point_data;
3300 if (
comm().size() > 1)
3315 bool error_finding_new_node =
false;
3316 if (
comm().rank() == 0)
3321 std::set<dof_id_type> previous_node_ids(vec_eval_input.
node_ids.begin(), vec_eval_input.
node_ids.end());
3325 std::set<dof_id_type> new_node_ids;
3326 for (
const auto & v_pair : *v_ptr)
3327 if (previous_node_ids.count(v_pair.first) == 0)
3328 new_node_ids.insert(v_pair.first);
3334 error_finding_new_node = (new_node_ids.empty());
3336 if (!error_finding_new_node)
3340 auto item = new_node_ids.begin();
3341 std::advance(item, random_node_idx);
3342 eim_point_data.
node_id = *item;
3346 if (!error_finding_new_node)
3348 const auto & vars = libmesh_map_find(*v_ptr,eim_point_data.
node_id);
3354 libmesh_error_msg_if(error_finding_new_node,
"Could not find new node in get_random_point()");
3360 return eim_point_data;
virtual unsigned int get_n_components() const =0
Specify the number of components in this parametrized function.
ElemType
Defines an enum for geometric element types.
Real get_max_abs_value(const DataMap &v) const
Get the maximum absolute value from a vector stored in the format that we use for basis functions...
virtual void preevaluate_parametrized_function_on_mesh_sides(const RBParameters &mu, const std::map< std::pair< dof_id_type, unsigned int >, std::vector< Point >> &side_all_xyz, const std::map< std::pair< dof_id_type, unsigned int >, subdomain_id_type > &sbd_ids, const std::map< std::pair< dof_id_type, unsigned int >, boundary_id_type > &side_boundary_ids, const std::map< std::pair< dof_id_type, unsigned int >, unsigned int > &side_types, const std::map< std::pair< dof_id_type, unsigned int >, std::vector< std::vector< Point >> > &side_all_xyz_perturb, const System &sys)
Same as preevaluate_parametrized_function_on_mesh() except for mesh sides.
This is the EquationSystems class.
std::map< std::pair< dof_id_type, unsigned int >, std::vector< std::vector< Number > > > SideQpDataMap
Type of the data structure used to map from (elem id, side index) -> [n_vars][n_qp] data...
Order
defines an enum for polynomial orders.
bool quiet_mode
Flag to indicate whether we print out extra information during the Offline stage. ...
A Node is like a Point, but with more information.
Real get_max_abs_value_in_training_set() const
Get the maximum value (across all processors) from the parametrized functions in the training set...
std::unordered_map< dof_id_type, boundary_id_type > _local_node_boundary_ids
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
boost::multiprecision::float128 real(const boost::multiprecision::float128 in)
const SideQpDataMap & get_side_basis_function(unsigned int i) const
Get a reference to the i^th side basis function.
void add_basis_function(const QpDataMap &bf)
Add bf to our EIM basis.
bool enrich_eim_approximation_on_nodes(const NodeDataMap &node_pf, bool add_basis_function, EimPointData *eim_point_data)
Implementation of enrich_eim_approximation() for the case of element nodes.
void set_interpolation_matrix_entry(unsigned int i, unsigned int j, Number value)
Set entry of the EIM interpolation matrix.
std::map< dof_id_type, std::vector< Number > > NodeDataMap
Type of the data structure used to map from (node id) -> [n_vars] data.
void initialize_eim_construction()
Perform initialization of this object to prepare for running train_eim_approximation().
bool _normalize_solution_snapshots
Set this boolean to true if we want to normalize solution snapshots used in training to have norm of ...
RBEIMEvaluation & get_rb_eim_evaluation()
Get a reference to the RBEvaluation object.
const std::set< unsigned int > & scale_components_in_enrichment() const
Get _scale_components_in_enrichment.
const std::set< boundary_id_type > & get_parametrized_function_boundary_ids() const
For RBParametrizedFunctions defined on element sides or nodes, we get/set the boundary IDs that this ...
void enable_set_Nmax_from_n_snapshots(int increment)
Call this method to set _set_Nmax_from_n_snapshots=true and _Nmax_from_n_snapshots_increment=incremen...
void rb_eim_solves(const std::vector< RBParameters > &mus, unsigned int N)
Perform rb_eim_solves at each mu in mus and store the results in _rb_eim_solutions.
Number get_eim_basis_function_value(unsigned int basis_function_index, dof_id_type elem_id, unsigned int comp, unsigned int qp) const
Same as above, except that we just return the value at the qp^th quadrature point.
static constexpr Real TOLERANCE
unsigned char get_elem_dim() const
const Elem & get_elem() const
Accessor for current Elem object.
virtual void reinit()
Reinitializes degrees of freedom and other required data on the current mesh.
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
virtual Real train_eim_approximation()
Generate the EIM approximation for the specified parametrized function using either POD or the Greedy...
std::vector< SideQpDataMap > _local_side_parametrized_functions_for_training
Same as _local_parametrized_functions_for_training except for side data.
Real get_parameter_min(const std::string ¶m_name) const
Get minimum allowable value of parameter param_name.
unsigned int get_n_basis_functions() const
Return the current number of EIM basis functions.
void add_node_basis_function(const NodeDataMap &node_bf)
Add node_bf to our EIM basis.
Number side_inner_product(const SideQpDataMap &v, const SideQpDataMap &w, bool apply_comp_scaling)
Same as inner_product() except for side data.
Real get_parameter_max(const std::string ¶m_name) const
Get maximum allowable value of parameter param_name.
std::map< std::pair< dof_id_type, unsigned int >, unsigned int > _local_side_quad_point_side_types
For side data, we also store "side type" info.
const RBParameters & get_parameters_max() const
Get an RBParameters object that specifies the maximum allowable value for each parameter.
std::unordered_map< dof_id_type, std::vector< Real > > _local_quad_point_JxW
This struct is used to encapsulate the arguments required to specify an EIM point that we may add to ...
This is the base class from which all geometric element types are derived.
virtual ~RBEIMConstruction()
void add_interpolation_data(Point p, unsigned int comp, dof_id_type elem_id, subdomain_id_type subdomain_id, unsigned int qp, const std::vector< Point > &perturbs, const std::vector< Real > &phi_i_qp, ElemType elem_type, const std::vector< Real > &JxW_all_qp, const std::vector< std::vector< Real >> &phi_i_all_qp, Order qrule_order, const Point &dxyz_dxi_elem_center, const Point &dxyz_deta_elem_center)
Add interpolation data associated with a new basis function.
void set_quiet_mode(bool quiet_mode_in)
Set the quiet_mode flag.
const std::vector< DenseVector< Number > > & get_rb_eim_solutions() const
Return the EIM solution coefficients from the most recent call to rb_eim_solves().
virtual void set_best_fit_type_flag(const std::string &best_fit_type_string)
Specify which type of "best fit" we use to guide the EIM greedy algorithm.
void reinit_eim_projection_matrix()
Zero the _eim_projection_matrix and resize it to be get_Nmax() x get_Nmax().
virtual void process_parameters_file(const std::string ¶meters_filename)
Read parameters in from file and set up this system accordingly.
Real _max_abs_value_in_training_set
Maximum value in _local_parametrized_functions_for_training across all processors.
const Parallel::Communicator & comm() const
bool is_quiet() const
Is the system in quiet mode?
RBEIMEvaluation * _rb_eim_eval
The RBEIMEvaluation object that we use to perform the EIM training.
const QpDataMap & get_parametrized_function_from_training_set(unsigned int training_index) const
Get a const reference to the specified parametrized function from the training set.
RBEIMConstruction(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
unsigned int get_interpolation_points_side_index(unsigned int index) const
virtual void clear() override
Clear this object.
void decrement_vector(QpDataMap &v, const DenseVector< Number > &coeffs)
Subtract coeffs[i]*basis_function[i] from v.
EimPointData get_random_point(const QpDataMap &v)
Helper function that identifies a random EIM point from v.
void set_training_parameter_values(const std::string ¶m_name, const std::vector< RBParameter > &values)
Overwrite the local training samples for param_name using values.
const NodeDataMap & get_node_basis_function(unsigned int i) const
Get a reference to the i^th node basis function.
The libMesh namespace provides an interface to certain functionality in the library.
unsigned int get_interpolation_points_comp(unsigned int index) const
const std::vector< DenseVector< Number > > & get_eim_solutions_for_training_set() const
Return a const reference to the EIM solutions for the parameters in the training set.
void set_abs_training_tolerance(Real new_training_tolerance)
Get/set the absolute tolerance for the basis training.
void update_eim_matrices(bool set_eim_error_indicator)
Update the matrices used in training the EIM approximation.
virtual Real train_eim_approximation_with_greedy()
Generate the EIM approximation for the specified parametrized function using the Greedy Algorithm...
const MeshBase & get_mesh() const
Number node_inner_product(const NodeDataMap &v, const NodeDataMap &w, bool apply_comp_scaling)
Same as inner_product() except for node data.
virtual void set_Nmax(unsigned int Nmax)
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
virtual void preevaluate_parametrized_function_on_mesh(const RBParameters &mu, const std::unordered_map< dof_id_type, std::vector< Point >> &all_xyz, const std::unordered_map< dof_id_type, subdomain_id_type > &sbd_ids, const std::unordered_map< dof_id_type, std::vector< std::vector< Point >> > &all_xyz_perturb, const System &sys)
Store the result of vectorized_evaluate.
void add_side_basis_function(const SideQpDataMap &side_bf)
Add side_bf to our EIM basis.
virtual Number lookup_preevaluated_value_on_mesh(unsigned int comp, dof_id_type elem_id, unsigned int qp) const
Look up the preevaluate values of the parametrized function for component comp, element elem_id...
static Number get_parametrized_function_value(const Parallel::Communicator &comm, const QpDataMap &pf, dof_id_type elem_id, unsigned int comp, unsigned int qp)
Same as above, except that we just return the value at the qp^th quadrature point.
This is the MeshBase class.
virtual void print_info()
Print out info that describes the current setup of this RBConstruction.
void side_decrement_vector(SideQpDataMap &v, const DenseVector< Number > &coeffs)
Same as decrement_vector() except for Side data.
std::map< std::pair< dof_id_type, unsigned int >, subdomain_id_type > _local_side_quad_point_subdomain_ids
virtual void elem_fe_reinit(const std::vector< Point > *const pts=nullptr)
Reinitializes interior FE objects on the current geometric element.
virtual void load_training_set(const std::map< std::string, std::vector< RBParameter >> &new_training_set)
Overwrite the training parameters with new_training_set.
Real get_node_max_abs_value(const NodeDataMap &v) const
Get the maximum absolute value from a vector stored in the format that we use for basis functions...
const std::unordered_map< dof_id_type, std::vector< Real > > & get_local_quad_point_JxW()
Get the interior and side quadrature weights.
virtual void initialize_training_parameters(const RBParameters &mu_min, const RBParameters &mu_max, const unsigned int n_global_training_samples, const std::map< std::string, bool > &log_param_scale, const bool deterministic=true)
Initialize the parameter ranges and indicate whether deterministic or random training parameters shou...
Real get_rel_training_tolerance()
void set_training_random_seed(int seed)
Set the seed that is used to randomly generate training parameters.
Real get_abs_training_tolerance()
unsigned int n_parameters() const
Get the number of parameters that have been added.
bool is_lookup_table
Boolean to indicate if this parametrized function is defined based on a lookup table or not...
std::string lookup_table_param_name
If this is a lookup table, then lookup_table_param_name specifies the parameter that is used to index...
const QpDataMap & get_basis_function(unsigned int i) const
Get a reference to the i^th basis function.
std::vector< QpDataMap > _local_parametrized_functions_for_training
The parametrized functions that are used for training.
static const boundary_id_type invalid_id
Number used for internal use.
Real _abs_training_tolerance
const RBParameters & get_parameters_min() const
Get an RBParameters object that specifies the minimum allowable value for each parameter.
virtual T el(const unsigned int i, const unsigned int j) const override final
const System & get_system() const
Accessor for associated system.
std::vector< NodeDataMap > _local_node_parametrized_functions_for_training
Same as _local_parametrized_functions_for_training except for node data.
void resize_data_structures(const unsigned int Nmax)
Resize the data structures for storing data associated with this object.
void svd(DenseVector< Real > &sigma)
Compute the singular value decomposition of the matrix.
const boundary_id_type node_boundary_id
unsigned int _Nmax
Maximum number of EIM basis functions we are willing to use.
Manages consistently variables, degrees of freedom, and coefficient vectors.
std::map< dof_id_type, std::vector< std::vector< Number > > > QpDataMap
Type of the data structure used to map from (elem id) -> [n_vars][n_qp] data.
Number get_eim_basis_function_side_value(unsigned int basis_function_index, dof_id_type elem_id, unsigned int side_index, unsigned int comp, unsigned int qp) const
Same as get_eim_basis_function_value() except for side data.
std::map< std::pair< dof_id_type, unsigned int >, std::vector< Real > > _local_side_quad_point_JxW
virtual Number lookup_preevaluated_side_value_on_mesh(unsigned int comp, dof_id_type elem_id, unsigned int side_index, unsigned int qp) const
Look up the preevaluated values of the parametrized function for component comp, element elem_id...
virtual void init_context(FEMContext &)
Pre-request FE data needed for calculations.
numeric_index_type get_local_n_training_samples() const
Get the total number of training samples local to this processor.
void store_eim_solutions_for_training_set()
Get the EIM solution vector at all parametrized functions in the training set.
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
void add_side_interpolation_data(Point p, unsigned int comp, dof_id_type elem_id, unsigned int side_index, subdomain_id_type subdomain_id, boundary_id_type boundary_id, unsigned int qp, const std::vector< Point > &perturbs, const std::vector< Real > &phi_i_qp)
Add interpolation data associated with a new basis function.
This class provides all data required for a physics package (e.g.
bool on_mesh_nodes() const
const Elem * reference_elem() const
void add_node_interpolation_data(Point p, unsigned int comp, dof_id_type node_id, boundary_id_type boundary_id)
Add interpolation data associated with a new basis function.
void node_decrement_vector(NodeDataMap &v, const DenseVector< Number > &coeffs)
Same as decrement_vector() except for node data.
const RBParameters & get_parameters() const
Get the current parameters.
void print_parameters() const
Print the current parameters.
unsigned int _max_abs_value_in_training_set_index
The training sample index at which we found _max_abs_value_in_training_set.
bool requires_all_elem_qp_data
Boolean to indicate whether this parametrized function requires data from all qps on the current elem...
void get_principal_submatrix(unsigned int sub_m, unsigned int sub_n, DenseMatrix< T > &dest) const
Put the sub_m x sub_n principal submatrix into dest.
void maxloc(T &r, unsigned int &max_id) const
This class is part of the rbOOmit framework.
bool requires_all_elem_center_data
Boolean to indicate whether this parametrized function requires data from the center on the current e...
void apply_normalization_to_solution_snapshots()
Rescale solution snapshots so that they all have unity norm.
numeric_index_type get_n_training_samples() const
Get the number of global training samples.
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
virtual std::unique_ptr< ElemAssembly > build_eim_assembly(unsigned int bf_index)=0
Build an element assembly object that will access basis function bf_index.
BEST_FIT_TYPE best_fit_type_flag
Enum that indicates which type of "best fit" algorithm we should use.
DenseMatrix< Number > _eim_projection_matrix
The matrix we use in order to perform L2 projections of parametrized functions as part of EIM trainin...
EimPointData get_random_point_from_training_sample()
Get a random point using the 0^th training sample as input to get_random_point(). ...
static void scale_parametrized_function(DataMap &local_pf, Number scaling_factor)
Scale all values in pf by scaling_factor.
void set_params_from_training_set(unsigned int global_index)
Set parameters to the RBParameters stored in index global_index of the global training set...
virtual void preevaluate_parametrized_function_on_mesh_nodes(const RBParameters &mu, const std::unordered_map< dof_id_type, Point > &all_xyz, const std::unordered_map< dof_id_type, boundary_id_type > &node_boundary_ids, const System &sys)
Same as preevaluate_parametrized_function_on_mesh() except for mesh nodes.
int _Nmax_from_n_snapshots_increment
virtual unsigned int n_sides() const =0
bool _set_Nmax_from_n_snapshots
If _set_Nmax_from_n_snapshots=true, then we overrule Nmax to be Nmax += _Nmax_from_n_snapshots_increm...
bool set_parameters(const RBParameters ¶ms)
Set the current parameters to params The parameters are checked for validity; an error is thrown if t...
dof_id_type get_interpolation_points_node_id(unsigned int index) const
virtual bool use_eim_error_indicator() const
Virtual function to indicate if we use the EIM error indicator in this case.
void set_rb_eim_evaluation(RBEIMEvaluation &rb_eim_eval_in)
Set the RBEIMEvaluation object.
const Elem * neighbor_ptr(unsigned int i) const
DenseVector< Number > rb_eim_solve(DenseVector< Number > &EIM_rhs)
Calculate the EIM approximation for the given right-hand side vector EIM_rhs.
void lu_solve(const DenseVector< T > &b, DenseVector< T > &x)
Solve the system Ax=b given the input vector b.
Number get_eim_basis_function_node_value(unsigned int basis_function_index, dof_id_type node_id, unsigned int var) const
Same as get_eim_basis_function_value() except for node data.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< std::pair< dof_id_type, unsigned int >, boundary_id_type > _local_side_quad_point_boundary_ids
void initialize_rb_property_map()
Initialize the rb_property_map of RBEIMEvaluation (this) from the rb_property_map stored in RBParamet...
std::unordered_map< dof_id_type, Point > _local_node_locations
Same as above except for node data.
const VectorizedEvalInput & get_vec_eval_input() const
Get the VectorizedEvalInput data.
const SideQpDataMap & get_side_parametrized_function_from_training_set(unsigned int training_index) const
virtual void initialize_eim_assembly_objects()
Build a vector of ElemAssembly objects that accesses the basis functions stored in this RBEIMConstruc...
unsigned char side
Current side for side_* to examine.
void max(const T &r, T &o, Request &req) const
static Number get_parametrized_function_side_value(const Parallel::Communicator &comm, const SideQpDataMap &pf, dof_id_type elem_id, unsigned int side_index, unsigned int comp, unsigned int qp)
Same as get_parametrized_function_value() except for side data.
virtual unsigned short dim() const =0
void set_rel_training_tolerance(Real new_training_tolerance)
Get/set the relative tolerance for the basis training.
std::vector< std::unique_ptr< ElemAssembly > > & get_eim_assembly_objects()
void set_value(const std::string ¶m_name, Real value)
Set the value of the specified parameter.
static Point map(const unsigned int dim, const Elem *elem, const Point &reference_point)
Number inner_product(const QpDataMap &v, const QpDataMap &w, bool apply_comp_scaling)
Evaluate the inner product of vec1 and vec2 which specify values at quadrature points.
void set_error_indicator_interpolation_row(const DenseVector< Number > &error_indicator_row)
bool enrich_eim_approximation_on_sides(const SideQpDataMap &side_pf, bool add_basis_function, EimPointData *eim_point_data)
Implementation of enrich_eim_approximation() for the case of element sides.
virtual const Elem & elem_ref(const dof_id_type i) const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
const std::map< std::pair< dof_id_type, unsigned int >, std::vector< Real > > & get_local_side_quad_point_JxW()
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...
void initialize_parametrized_functions_in_training_set()
Compute and store the parametrized function for each parameter in the training set at all the stored ...
dof_id_type get_interpolation_points_elem_id(unsigned int index) const
unsigned int get_Nmax() const
Get/set Nmax, the maximum number of RB functions we are willing to compute.
void disable_set_Nmax_from_n_snapshots()
Call this method to set _set_Nmax_from_n_snapshots=false and reset _Nmax_from_n_snapshots_increment t...
std::pair< Real, unsigned int > compute_max_eim_error()
Find the training sample that has the largest EIM approximation error based on the current EIM approx...
RBEIMEvaluation::NodeDataMap NodeDataMap
Type of the data structure used to map from node id -> [n_vars] data.
bool serial_training_set
This boolean flag indicates whether or not the training set should be the same on all processors...
virtual Number lookup_preevaluated_node_value_on_mesh(unsigned int comp, dof_id_type node_id) const
Look up the preevaluate values of the parametrized function for component comp, node node_id...
RBParameters get_params_from_training_set(unsigned int global_index)
Return the RBParameters in index global_index of the global training set.
This class is part of the rbOOmit framework.
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
std::vector< Real > _component_scaling_in_training_set
Keep track of a scaling factor for each component of the parametrized functions in the training set w...
bool enrich_eim_approximation_on_interiors(const QpDataMap &interior_pf, bool add_basis_function, EimPointData *eim_point_data)
Implementation of enrich_eim_approximation() for the case of element interiors.
void enrich_eim_approximation(unsigned int training_index, bool add_basis_function, EimPointData *eim_point_data)
Add a new basis function to the EIM approximation.
Real _rel_training_tolerance
Relative and absolute tolerances for training the EIM approximation.
RBEIMEvaluation::QpDataMap QpDataMap
Type of the data structure used to map from (elem id) -> [n_vars][n_qp] data.
unsigned int get_n_parametrized_functions_for_training() const
Get the number of parametrized functions used for training.
bool on_mesh_sides() const
Defines a dense vector for use in Finite Element-type computations.
virtual Real train_eim_approximation_with_POD()
Generate the EIM approximation for the specified parametrized function using Proper Orthogonal Decomp...
const std::string & name() const
std::unordered_map< dof_id_type, std::vector< Point > > _local_quad_point_locations
The quadrature point locations, quadrature point weights (JxW), and subdomain IDs on every element lo...
unsigned int n_vars() const
const NodeDataMap & get_node_parametrized_function_from_training_set(unsigned int training_index) const
std::map< std::pair< dof_id_type, unsigned int >, std::vector< Point > > _local_side_quad_point_locations
Same as above except for side data.
RBParametrizedFunction & get_parametrized_function()
Get a reference to the parametrized function.
virtual void initialize_lookup_table()
If this parametrized function is defined based on a lookup table then we can call this function to in...
void initialize_parameters(const RBParameters &mu_min_in, const RBParameters &mu_max_in, const std::map< std::string, std::vector< Real >> &discrete_parameter_values)
Initialize the parameter ranges and set current_parameters.
std::vector< std::unique_ptr< ElemAssembly > > _rb_eim_assembly_objects
The vector of assembly objects that are created to point to this RBEIMConstruction.
RBEIMEvaluation::SideQpDataMap SideQpDataMap
Type of the data structure used to map from (elem id,side_index) -> [n_vars][n_qp] data...
std::map< std::pair< dof_id_type, unsigned int >, std::vector< std::vector< Point > > > _local_side_quad_point_locations_perturbations
This class enables evaluation of an Empirical Interpolation Method (EIM) approximation.
processor_id_type processor_id() const
void initialize_qp_data()
Initialize the data associated with each quad point (location, JxW, etc.) so that we can use this in ...
std::unordered_map< dof_id_type, std::vector< std::vector< Point > > > _local_quad_point_locations_perturbations
EIM approximations often arise when applying a geometric mapping to a Reduced Basis formulation...
virtual ElemType type() const =0
unsigned int get_n_params() const
Get the number of parameters.
A Point defines a location in LIBMESH_DIM dimensional Real space.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
unsigned int get_interpolation_points_qp(unsigned int index) const
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
static void scale_node_parametrized_function(NodeDataMap &local_pf, Number scaling_factor)
Scale all values in pf by scaling_factor The templated function above handles the elem and side cases...
virtual void clear()
Clear all the data structures associated with the system.
bool is_discrete_parameter(const std::string &mu_name) const
Is parameter mu_name discrete?
Point vertex_average() const
void set_rb_construction_parameters(unsigned int n_training_samples_in, bool deterministic_training_in, int training_parameters_random_seed_in, bool quiet_mode_in, unsigned int Nmax_in, Real rel_training_tolerance_in, Real abs_training_tolerance_in, const RBParameters &mu_min_in, const RBParameters &mu_max_in, const std::map< std::string, std::vector< Real >> &discrete_parameter_values_in, const std::map< std::string, bool > &log_scaling, std::map< std::string, std::vector< RBParameter >> *training_sample_list=nullptr)
Set the state of this RBConstruction object based on the arguments to this function.
const std::set< unsigned char > & elem_dimensions() const
std::unordered_map< dof_id_type, subdomain_id_type > _local_quad_point_subdomain_ids
void print_discrete_parameter_values() const
Print out all the discrete parameter values.
static unsigned int get_random_int_0_to_n(unsigned int n)
Static helper function that is used by get_random_point().
void set_union(T &data, const unsigned int root_id) const
Real fd_delta
The finite difference step size in the case that this function in the case that this function uses fi...
static Number get_parametrized_function_node_value(const Parallel::Communicator &comm, const NodeDataMap &pf, dof_id_type node_id, unsigned int comp)
Same as get_parametrized_function_value() except for node data.