31 #include "libmesh/rb_construction_base.h" 32 #include "libmesh/rb_parameters.h" 35 #include "libmesh/id_types.h" 36 #include "libmesh/libmesh_common.h" 37 #include "libmesh/numeric_vector.h" 38 #include "libmesh/equation_systems.h" 39 #include "libmesh/parallel.h" 40 #include "libmesh/condensed_eigen_system.h" 41 #include "libmesh/linear_implicit_system.h" 42 #include "libmesh/int_range.h" 43 #include "libmesh/utility.h" 46 #include "libmesh/dof_map.h" 47 #include "libmesh/shell_matrix.h" 48 #include "libmesh/sparse_matrix.h" 59 std::pair<unsigned int, unsigned int> calculate_n_local_samples_and_index(
61 const unsigned int n_global_samples,
64 unsigned int n_local_samples = n_global_samples;
65 unsigned int first_local_index = 0;
68 return {n_local_samples, first_local_index};
71 unsigned int quotient = n_global_samples/
communicator.size();
72 unsigned int remainder = n_global_samples%
communicator.size();
75 n_local_samples = (quotient + 1);
80 n_local_samples = quotient;
81 first_local_index =
communicator.rank()*quotient + remainder;
83 return {n_local_samples, first_local_index};
96 const std::string & name_in,
97 const unsigned int number_in)
98 : Base(es, name_in, number_in),
100 serial_training_set(false),
101 _training_parameters_initialized(false),
102 _first_local_index(0),
103 _n_local_training_samples(0),
104 _n_global_training_samples(0),
105 _training_parameters_random_seed(-1)
109 template <
class Base>
112 template <
class Base>
117 RBParametrized::clear();
118 _training_parameters.clear();
121 template <
class Base>
129 inner_product_storage_vector->init (this->n_dofs(), this->n_local_dofs(),
false,
PARALLEL);
132 template <
class Base>
134 std::pair<numeric_index_type, Real> & error_pair)
138 unsigned int proc_ID_index;
142 communicator.broadcast(error_pair.first, proc_ID_index);
145 template <
class Base>
148 libmesh_error_msg_if(!_training_parameters_initialized,
149 "Error: training parameters must first be initialized.");
156 if (_training_parameters.empty())
158 if (serial_training_set)
161 return this->comm().size();
164 return _n_global_training_samples;
167 template <
class Base>
170 libmesh_error_msg_if(!_training_parameters_initialized,
171 "Error: training parameters must first be initialized.");
177 if (_training_parameters.empty())
180 return _n_local_training_samples;
183 template <
class Base>
186 libmesh_error_msg_if(!_training_parameters_initialized,
187 "Error: training parameters must first be initialized.");
194 if (_training_parameters.empty())
196 if (serial_training_set)
199 return this->comm().rank();
202 return _first_local_index;
205 template <
class Base>
208 libmesh_error_msg_if(!_training_parameters_initialized,
209 "Error: training parameters must first be initialized.");
211 if (_training_parameters.empty())
214 return _first_local_index + _n_local_training_samples;
217 template <
class Base>
220 set_parameters(get_params_from_training_set(global_index));
223 template <
class Base>
226 libmesh_error_msg_if(!_training_parameters_initialized,
227 "Error: training parameters must first be initialized.");
232 if (!_training_parameters.empty())
234 libmesh_error_msg_if((global_index < this->get_first_local_training_index()) ||
235 (global_index >= this->get_last_local_training_index()),
238 <<
" must be within range: " 239 << this->get_first_local_training_index()
241 << this->get_last_local_training_index());
243 const numeric_index_type local_index = global_index - get_first_local_training_index();
244 for (
const auto & [param_name, sample_vector] : _training_parameters)
245 params.
set_value(param_name, sample_vector[local_index]);
250 const auto & mine = get_parameters();
251 for (
const auto & [key, extra_sample_vector] :
252 as_range(mine.extra_begin(), mine.extra_end()))
262 template <
class Base>
265 libmesh_error_msg_if(!_training_parameters_initialized,
266 "Error: training parameters must first be initialized.");
269 if ((this->get_first_local_training_index() <= global_index) &&
270 (global_index < this->get_last_local_training_index()))
273 set_params_from_training_set(global_index);
276 root_id = this->processor_id();
280 this->comm().max(root_id);
281 broadcast_parameters(root_id);
284 template <
class Base>
287 const unsigned int n_global_training_samples,
288 const std::map<std::string,bool> & log_param_scale,
289 const bool deterministic)
294 libMesh::out <<
"Initializing training parameters with " 295 << (deterministic ?
"deterministic " :
"random " )
296 <<
"training set..." << std::endl;
298 for (
const auto & pr : log_param_scale)
301 <<
": log scaling = " 310 const auto [first_local_index, last_local_index] =
311 generate_training_parameters_deterministic(this->comm(),
313 _training_parameters,
314 n_global_training_samples,
317 serial_training_set);
318 _first_local_index = first_local_index;
319 _n_local_training_samples = last_local_index-first_local_index;
324 const auto [first_local_index, last_local_index] =
325 generate_training_parameters_random(this->comm(),
327 _training_parameters,
328 n_global_training_samples,
331 this->_training_parameters_random_seed,
332 serial_training_set);
333 _first_local_index = first_local_index;
334 _n_local_training_samples = last_local_index-first_local_index;
336 _n_global_training_samples = _n_local_training_samples;
338 if (!serial_training_set)
339 this->comm().sum(_n_global_training_samples);
343 if (get_n_discrete_params() > 0)
345 for (
auto & [param_name, sample_vector] : _training_parameters)
347 if (is_discrete_parameter(param_name))
349 std::vector<Real> discrete_values =
350 get_discrete_parameter_values().find(param_name)->second;
351 for (
const auto sample_idx :
index_range(sample_vector))
354 std::vector<Real> discretized_vector(sample_vector[sample_idx].size());
355 std::transform(sample_vector[sample_idx].cbegin(),
356 sample_vector[sample_idx].cend(),
357 discretized_vector.begin(),
358 [&discrete_values](
const Real & val) {
359 return get_closest_value(val, discrete_values);
361 sample_vector[sample_idx] = discretized_vector;
367 _training_parameters_initialized =
true;
370 template <
class Base>
374 libmesh_parallel_only(this->comm());
377 libmesh_error_msg_if(!_training_parameters_initialized,
378 "Error: load_training_set cannot be used to initialize parameters");
381 const unsigned int n_params = get_n_params();
382 libmesh_error_msg_if(new_training_set.size() > n_params,
383 "Error: new_training_set should not have more than get_n_params() parameters.");
387 const bool size_matches = (new_training_set.size() == n_params);
388 this->comm().verify(size_matches);
396 _first_local_index = 0;
397 _n_local_training_samples =
398 cast_int<numeric_index_type>(new_training_set.begin()->second.size());
399 _n_global_training_samples = _n_local_training_samples;
401 if (!serial_training_set)
403 this->comm().sum(_n_global_training_samples);
406 std::vector<numeric_index_type> local_sizes (this->n_processors(), 0);
407 local_sizes[this->processor_id()] = _n_local_training_samples;
408 this->comm().sum(local_sizes);
412 for (
auto p :
make_range(this->processor_id()))
413 _first_local_index += local_sizes[p];
417 for (
const auto & pr : _training_parameters)
418 libmesh_error_msg_if(!new_training_set.count(pr.first),
419 "Parameters must be identical in order to overwrite dataset.");
422 _training_parameters = new_training_set;
430 for (
auto & [param_name, sample_vector]: _training_parameters)
432 if (new_training_set.count(param_name))
434 for (
const auto i :
make_range(get_local_n_training_samples()))
436 const unsigned int num_new_samples = libmesh_map_find(new_training_set,param_name).size();
437 libmesh_error_msg_if (num_new_samples==0,
"new_training_set set should not be empty");
439 const unsigned int new_training_set_index = i % num_new_samples;
440 sample_vector[i] = libmesh_map_find(new_training_set,param_name)[new_training_set_index];
447 template <
class Base>
449 const std::string & param_name,
const std::vector<RBParameter> & values)
451 libmesh_error_msg_if(!_training_parameters_initialized,
452 "Training parameters must be initialized before calling set_training_parameter_values");
453 libmesh_error_msg_if(values.size() != get_local_n_training_samples(),
454 "Inconsistent sizes");
457 auto & training_vector = libmesh_map_find(_training_parameters, param_name);
458 training_vector = values;
462 template <
class Base>
463 std::pair<std::size_t, std::size_t>
465 const std::map<std::string, bool> & log_param_scale,
466 std::map<std::string, std::vector<RBParameter>> & local_training_parameters_in,
467 const unsigned int n_global_training_samples_in,
470 const int training_parameters_random_seed,
471 const bool serial_training_set)
473 const unsigned int num_params = min_parameters.
n_parameters();
474 libmesh_error_msg_if(num_params!=max_parameters.
n_parameters(),
475 "Number of parameters must be identical for min/max.");
478 local_training_parameters_in.clear();
483 if (training_parameters_random_seed < 0)
485 if (!serial_training_set)
490 std::srand( static_cast<unsigned>( std::time(0)*(1+
communicator.rank()) ));
499 unsigned int current_time =
static_cast<unsigned>( std::time(0) );
501 std::srand(current_time);
506 if (!serial_training_set)
511 std::srand( static_cast<unsigned>( training_parameters_random_seed*(1+
communicator.rank()) ));
517 std::srand( static_cast<unsigned>( training_parameters_random_seed ));
528 const auto & [n_local_training_samples, first_local_index] =
529 calculate_n_local_samples_and_index(
communicator, n_global_training_samples_in,
530 serial_training_set);
531 for (
const auto & pr : min_parameters)
532 local_training_parameters_in[pr.first] = std::vector<RBParameter>(n_local_training_samples);
535 for (
auto & [param_name, sample_vector] : local_training_parameters_in)
537 for (
auto i :
make_range(n_local_training_samples))
539 Real random_number =
static_cast<Real>(std::rand()) / RAND_MAX;
542 if (libmesh_map_find(log_param_scale, param_name))
544 Real log_min = std::log10(min_parameters.get_value(param_name));
545 Real log_range = std::log10(max_parameters.
get_value(param_name) / min_parameters.get_value(param_name));
547 sample_vector[i] = {
std::pow(10., log_min + random_number*log_range )};
553 random_number * (max_parameters.
get_value(param_name) -
554 min_parameters.get_value(param_name)) +
555 min_parameters.get_value(param_name)};
559 return {first_local_index, first_local_index+n_local_training_samples};
562 template <
class Base>
563 std::pair<std::size_t, std::size_t>
565 const std::map<std::string, bool> & log_param_scale,
566 std::map<std::string, std::vector<RBParameter>> & local_training_parameters_in,
567 const unsigned int n_global_training_samples_in,
570 const bool serial_training_set)
573 const unsigned int num_params = min_parameters.
n_parameters();
579 libmesh_not_implemented_msg(
"ERROR: Deterministic training sample generation " 580 "not implemented for more than three parameters.");
589 const auto &[n_local_training_samples, first_local_index] =
590 calculate_n_local_samples_and_index(
communicator, n_global_training_samples_in,
591 serial_training_set);
592 const auto last_local_index = first_local_index + n_local_training_samples;
593 for (
const auto & pr : min_parameters)
594 local_training_parameters_in[pr.first] = std::vector<RBParameter>(n_local_training_samples);
600 std::vector<unsigned int> n_training_samples_per_param(3);
601 for (
unsigned int param=0; param<3; param++)
603 if (param < num_params)
605 n_training_samples_per_param[param] =
606 static_cast<unsigned int>( std::round(
std::pow(static_cast<Real>(n_global_training_samples_in), 1./num_params)) );
610 n_training_samples_per_param[param] = 1;
618 unsigned int total_samples_check = 1;
619 for (
unsigned int n_samples : n_training_samples_per_param)
621 total_samples_check *= n_samples;
624 libmesh_error_msg_if(total_samples_check != n_global_training_samples_in,
625 "Error: Number of training samples = " 626 << n_global_training_samples_in
627 <<
" does not enable a uniform grid of samples with " 628 << num_params <<
" parameters. Try " 629 << total_samples_check <<
" samples instead?");
634 std::vector<std::vector<Real>> training_samples_per_param(num_params);
637 for (
const auto & pr : min_parameters)
639 const std::string & param_name = pr.first;
640 const bool use_log_scaling = libmesh_map_find(log_param_scale, param_name);
641 Real min_param = min_parameters.get_value(param_name);
644 training_samples_per_param[i].resize(n_training_samples_per_param[i]);
646 for (
unsigned int j=0; j<n_training_samples_per_param[i]; j++)
651 Real epsilon = 1.e-6;
652 Real log_min = std::log10(min_param + epsilon);
653 Real log_range = std::log10( (max_param-epsilon) / (min_param+epsilon) );
654 Real step_size = log_range /
655 std::max((
unsigned int)1,(n_training_samples_per_param[i]-1));
657 if (j<(n_training_samples_per_param[i]-1))
659 training_samples_per_param[i][j] =
std::pow(10., log_min + j*step_size );
665 training_samples_per_param[i][j] = max_param;
671 Real step_size = (max_param - min_param) /
672 std::max((
unsigned int)1,(n_training_samples_per_param[i]-1));
673 training_samples_per_param[i][j] = j*step_size + min_param;
683 std::vector<unsigned int> indices(3);
684 unsigned int index_count = 0;
685 for (indices[0]=0; indices[0]<n_training_samples_per_param[0]; indices[0]++)
687 for (indices[1]=0; indices[1]<n_training_samples_per_param[1]; indices[1]++)
689 for (indices[2]=0; indices[2]<n_training_samples_per_param[2]; indices[2]++)
691 unsigned int param_count = 0;
692 for (
const auto & pr : min_parameters)
694 std::vector<RBParameter> & training_vector =
695 libmesh_map_find(local_training_parameters_in, pr.first);
696 if (first_local_index <= index_count && index_count < last_local_index)
697 training_vector[index_count - first_local_index] =
698 {training_samples_per_param[param_count][indices[param_count]]};
707 return {first_local_index, first_local_index+n_local_training_samples};
711 template <
class Base>
714 libmesh_assert_less (proc_id, this->n_processors());
718 libmesh_error_msg_if(current_parameters.
n_samples()!=1,
719 "Only single-sample RBParameter objects can be broadcast.");
724 const std::size_t nparams = current_parameters.
n_parameters();
725 const std::size_t nsamples = current_parameters.
n_samples();
728 std::vector<std::size_t> param_value_sizes;
729 param_value_sizes.reserve(nparams);
730 for (
const auto & pr : current_parameters)
731 param_value_sizes.push_back(pr.second[0].size());
734 this->comm().broadcast(param_value_sizes, proc_id);
735 std::size_t buffsize = std::accumulate(param_value_sizes.cbegin(), param_value_sizes.cend(), 0ul);
736 std::vector<Real> serialized_parameters;
737 serialized_parameters.reserve(buffsize);
740 for (
const auto & pr : current_parameters)
742 for (
const auto sample_idx :
make_range(nsamples))
743 serialized_parameters.insert(serialized_parameters.end(),
744 pr.second[sample_idx].cbegin(),
745 pr.second[sample_idx].cend());
749 this->comm().broadcast(serialized_parameters, proc_id);
752 std::size_t param_idx = 0;
753 auto val_idx = serialized_parameters.cbegin();
754 for (
const auto & pr : current_parameters)
756 const std::size_t param_value_size = param_value_sizes[param_idx];
757 for (
const auto sample_idx:
make_range(nsamples))
759 auto end_val_idx =
std::next(val_idx,param_value_size);
761 current_parameters.set_value(pr.first, sample_idx, sample_val);
762 val_idx = end_val_idx;
768 set_parameters(current_parameters);
771 template <
class Base>
774 this->_training_parameters_random_seed = seed;
780 #if defined(LIBMESH_HAVE_SLEPC) Real get_value(const std::string ¶m_name) const
Get the value of the specified parameter, throw an error if it does not exist.
This is the EquationSystems class.
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
The libMesh namespace provides an interface to certain functionality in the library.
uint8_t processor_id_type
unsigned int n_parameters() const
Get the number of parameters that have been added.
dof_id_type numeric_index_type
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
RBConstructionBase(EquationSystems &es, const std::string &name, const unsigned int number)
Constructor.
std::vector< Real > RBParameter
Typedef for an individual RB parameter.
This class is part of the rbOOmit framework.
unsigned int n_samples() const
Returns the number of samples stored for all parameters.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void set_extra_value(const std::string ¶m_name, Real value)
Set the value of the specified extra parameter.
DIE A HORRIBLE DEATH HERE typedef MPI_Comm communicator
void set_value(const std::string ¶m_name, Real value)
Set the value of the specified parameter.
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...
This class is part of the rbOOmit framework.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...