44 #include "Teuchos_Array.hpp" 48 template <
typename value_type,
typename execution_space>
53 num_KL = solverParams.get<
int>(
"Number of KL Terms");
54 mean = solverParams.get<
double>(
"Mean");
55 std_dev = solverParams.get<
double>(
"Standard Deviation");
57 Teuchos::Array<double> domain_upper_bound_double;
58 Teuchos::Array<double> domain_lower_bound_double;
59 Teuchos::Array<double> correlation_length_double;
62 if (solverParams.isType<std::string>(
"Domain Upper Bounds"))
63 domain_upper_bound_double =
64 Teuchos::getArrayFromStringParameter<double>(
65 solverParams,
"Domain Upper Bounds");
67 domain_upper_bound_double =
68 solverParams.get< Teuchos::Array<double> >(
"Domain Upper Bounds");
71 Teuchos::Array<value_type> domain_upper_bound(domain_upper_bound_double.size());
72 for (
int i=0; i<domain_upper_bound.size(); i++)
73 domain_upper_bound[i]=domain_upper_bound_double[i];
76 if (solverParams.isType<std::string>(
"Domain Lower Bounds"))
77 domain_lower_bound_double =
78 Teuchos::getArrayFromStringParameter<double>(
79 solverParams,
"Domain Lower Bounds");
81 domain_lower_bound_double =
82 solverParams.get< Teuchos::Array<double> >(
"Domain Lower Bounds");
85 Teuchos::Array<value_type> domain_lower_bound(domain_lower_bound_double.size());
86 for (
int i=0; i<domain_lower_bound.size(); i++)
87 domain_lower_bound[i]=domain_lower_bound_double[i];
90 if (solverParams.isType<std::string>(
"Correlation Lengths"))
91 correlation_length_double =
92 Teuchos::getArrayFromStringParameter<double>(
93 solverParams,
"Correlation Lengths");
95 correlation_length_double =
96 solverParams.get< Teuchos::Array<double> >(
"Correlation Lengths");
99 Teuchos::Array<value_type> correlation_length(correlation_length_double.size());
100 for (
int i=0; i<correlation_length.size(); i++)
101 correlation_length[i]=correlation_length_double[i];
104 dim = domain_upper_bound.size();
105 Teuchos::Array< Teuchos::Array< one_d_eigen_pair_type > > eig_pairs(dim);
106 for (
int i=0; i<dim; i++) {
107 eig_pairs[i].resize(num_KL);
109 num_KL, domain_lower_bound[i], domain_upper_bound[i],
110 correlation_length[i], i, solverParams);
115 int num_prod =
static_cast<int>(
std::pow(static_cast<double>(num_KL),
116 static_cast<double>(dim)));
117 Teuchos::Array<product_eigen_pair_type> product_eig_pairs(num_prod);
118 Teuchos::Array<int>
index(dim, 0);
120 Teuchos::Array<value_type> vals(dim);
121 Teuchos::Array<one_d_eigen_pair_type> eigs(dim);
122 while (cnt < num_prod) {
123 for (
int i=0; i<dim; i++) {
124 eigs[i] = eig_pairs[i][
index[i]];
126 product_eig_pairs[cnt].set(eigs);
129 while (
j < dim-1 &&
index[
j] == num_KL) {
138 std::sort(product_eig_pairs.begin(), product_eig_pairs.end(),
142 product_eigen_funcs =
144 product_eigen_values =
146 typename eigen_func_array_type::HostMirror host_product_eigen_funcs =
148 typename eigen_value_array_type::HostMirror host_product_eigen_values =
150 for (
int i=0; i<num_prod; ++i) {
151 host_product_eigen_values(i) = 1.0;
152 for (
int j=0;
j<dim; ++
j) {
153 host_product_eigen_values(i) *= product_eig_pairs[i].eig_pairs[
j].eig_val;
154 host_product_eigen_funcs(i,
j) =
155 product_eig_pairs[i].eig_pairs[
j].eig_func;
162 template <
typename value_type,
typename execution_space>
163 template <
typename po
int_type,
typename rv_type>
164 KOKKOS_INLINE_FUNCTION
165 typename Teuchos::PromotionTraits<typename rv_type::value_type, value_type>::promote
167 evaluate(
const point_type& point,
const rv_type& random_variables)
const 169 typedef typename Teuchos::PromotionTraits<typename rv_type::value_type, value_type>::promote result_type;
170 result_type result = mean;
171 for (
int i=0; i<num_KL; ++i) {
172 result_type t = std_dev*
std::sqrt(product_eigen_values(i));
173 for (
int j=0;
j<dim; ++
j)
174 t *= product_eigen_funcs(i,
j).evaluate(point[
j]);
175 result += random_variables[i]*t;
180 template <
typename value_type,
typename execution_space>
181 template <
typename po
int_type>
182 KOKKOS_INLINE_FUNCTION
183 typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote
187 typedef typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote result_type;
188 result_type result = 0.0;
189 for (
int i=0; i<num_KL; i++) {
191 for (
int j=0;
j<dim; ++
j)
192 t *= product_eigen_funcs(i,
j).evaluate(point[
j]);
193 result += product_eigen_values(i).eig_val*t*t;
198 template <
typename value_type,
typename execution_space>
199 template <
typename po
int_type>
200 KOKKOS_INLINE_FUNCTION
201 typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote
205 typedef typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote result_type;
206 result_type t = std_dev*
std::sqrt(product_eigen_values(i));
207 for (
int j=0;
j<dim; ++
j)
208 t *= product_eigen_funcs(i,
j).evaluate(point[
j]);
212 template <
typename value_type,
typename execution_space>
217 os <<
"KL expansion using " << num_KL <<
" terms in " << dim
218 <<
" dimensions:" << std::endl;
219 os <<
"\tMean = " << mean <<
", standard deviation = " << std_dev
221 os <<
"\tEigenvalues, Eigenfunctions:" << std::endl;
222 for (
int i=0; i<num_KL; i++) {
223 os <<
"\t\t" << product_eigen_values(i) <<
", ";
224 for (
int j=0;
j<dim-1; ++
j) {
226 product_eigen_funcs(i,
j).print(os);
230 product_eigen_funcs(i,dim-1).print(os);
231 os <<
")" << std::endl;
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
SparseArrayIterator< index_iterator, value_iterator >::value_type index(const SparseArrayIterator< index_iterator, value_iterator > &it)
Kokkos::View< one_d_eigen_func_type **, Device > eigen_func_array_type
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
Kokkos::View< MeshScalar *, Device > eigen_value_array_type
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typename point_type::value_type, value_type >::promote evaluate_standard_deviation(const point_type &point) const
Evaluate standard deviation of random field at a point.
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typename point_type::value_type, value_type >::promote evaluate_eigenfunction(const point_type &point, int i) const
Evaluate given eigenfunction at a point.
const Teuchos::Array< eigen_pair_type > & getEigenPairs() const
Get eigenpairs.
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typename rv_type::value_type, value_type >::promote evaluate(const point_type &point, const rv_type &random_variables) const
Evaluate random field at a point.
Class representing an exponential covariance function and its KL eigevalues/eigenfunctions.
Predicate class for sorting product eigenfunctions based on eigenvalue.
ExponentialRandomField()
Default constructor.
void print(std::ostream &os) const
Print KL expansion.
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)