Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_KL_ExponentialRandomFieldImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include <cmath>
43#include <algorithm>
44#include "Teuchos_Array.hpp"
45
47
48template <typename value_type, typename execution_space>
50ExponentialRandomField(Teuchos::ParameterList& solverParams)
51{
52 // Get required parameters
53 num_KL = solverParams.get<int>("Number of KL Terms");
54 mean = solverParams.get<double>("Mean");
55 std_dev = solverParams.get<double>("Standard Deviation");
56
57 Teuchos::Array<double> domain_upper_bound_double;
58 Teuchos::Array<double> domain_lower_bound_double;
59 Teuchos::Array<double> correlation_length_double;
60
61 // Get Domain Upper Bounds
62 if (solverParams.isType<std::string>("Domain Upper Bounds"))
63 domain_upper_bound_double =
64 Teuchos::getArrayFromStringParameter<double>(
65 solverParams, "Domain Upper Bounds");
66 else
67 domain_upper_bound_double =
68 solverParams.get< Teuchos::Array<double> >("Domain Upper Bounds");
69
70 // (Convert each element from double to value_type)
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];
74
75 // Get Domain Lower Bounds
76 if (solverParams.isType<std::string>("Domain Lower Bounds"))
77 domain_lower_bound_double =
78 Teuchos::getArrayFromStringParameter<double>(
79 solverParams, "Domain Lower Bounds");
80 else
81 domain_lower_bound_double =
82 solverParams.get< Teuchos::Array<double> >("Domain Lower Bounds");
83
84 // (Convert each element from double to value_type)
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];
88
89 // Get Correlation Lengths
90 if (solverParams.isType<std::string>("Correlation Lengths"))
91 correlation_length_double =
92 Teuchos::getArrayFromStringParameter<double>(
93 solverParams, "Correlation Lengths");
94 else
95 correlation_length_double =
96 solverParams.get< Teuchos::Array<double> >("Correlation Lengths");
97
98 // (Convert each element from double to value_type)
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];
102
103 // Compute 1-D eigenfunctions for each dimension
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);
108 OneDExponentialCovarianceFunction<value_type> cov_func(
109 num_KL, domain_lower_bound[i], domain_upper_bound[i],
110 correlation_length[i], i, solverParams);
111 eig_pairs[i] = cov_func.getEigenPairs();
112 }
113
114 // Compute all possible tensor product combinations of 1-D eigenfunctions
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);
119 int cnt = 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]];
125 }
126 product_eig_pairs[cnt].set(eigs);
127 ++index[0];
128 int j = 0;
129 while (j < dim-1 && index[j] == num_KL) {
130 index[j] = 0;
131 ++j;
132 ++index[j];
133 }
134 ++cnt;
135 }
136
137 // Sort product eigenfunctions based on product eigenvalue
138 std::sort(product_eig_pairs.begin(), product_eig_pairs.end(),
139 ProductEigenPairGreater<one_d_eigen_func_type, execution_space>());
141 // Copy eigenpairs into view
142 product_eigen_funcs =
143 eigen_func_array_type("product eigen functions", num_prod, dim);
144 product_eigen_values =
145 eigen_value_array_type("product eigen vvalues", num_prod);
146 typename eigen_func_array_type::HostMirror host_product_eigen_funcs =
147 Kokkos::create_mirror_view(product_eigen_funcs);
148 typename eigen_value_array_type::HostMirror host_product_eigen_values =
149 Kokkos::create_mirror_view(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;
156 }
157 }
158 Kokkos::deep_copy(product_eigen_funcs, host_product_eigen_funcs);
159 Kokkos::deep_copy(product_eigen_values, host_product_eigen_values);
160}
161
162template <typename value_type, typename execution_space>
163template <typename point_type, typename rv_type>
164KOKKOS_INLINE_FUNCTION
165typename Teuchos::PromotionTraits<typename rv_type::value_type, value_type>::promote
167evaluate(const point_type& point, const rv_type& random_variables) const
168{
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;
176 }
177 return result;
178}
179
180template <typename value_type, typename execution_space>
181template <typename point_type>
182KOKKOS_INLINE_FUNCTION
183typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote
185evaluate_standard_deviation(const point_type& point) const
186{
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++) {
190 result_type t = 1.0;
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;
194 }
195 return std::sqrt(result);
196}
197
198template <typename value_type, typename execution_space>
199template <typename point_type>
200KOKKOS_INLINE_FUNCTION
201typename Teuchos::PromotionTraits<typename point_type::value_type, value_type>::promote
203evaluate_eigenfunction(const point_type& point, int i) const
204{
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]);
209 return t;
210}
211
212template <typename value_type, typename execution_space>
213void
215print(std::ostream& os) const
216{
217 os << "KL expansion using " << num_KL << " terms in " << dim
218 << " dimensions:" << std::endl;
219 os << "\tMean = " << mean << ", standard deviation = " << std_dev
220 << std::endl;
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) {
225 os << "(";
226 product_eigen_funcs(i,j).print(os);
227 os << ") * ";
228 }
229 os << "(";
230 product_eigen_funcs(i,dim-1).print(os);
231 os << ")" << std::endl;
232 }
233}
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typenamepoint_type::value_type, value_type >::promote evaluate_eigenfunction(const point_type &point, int i) const
Evaluate given eigenfunction at a point.
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typenamerv_type::value_type, value_type >::promote evaluate(const point_type &point, const rv_type &random_variables) const
Evaluate random field at a point.
Kokkos::View< one_d_eigen_func_type **, Device > eigen_func_array_type
void print(std::ostream &os) const
Print KL expansion.
KOKKOS_INLINE_FUNCTION Teuchos::PromotionTraits< typenamepoint_type::value_type, value_type >::promote evaluate_standard_deviation(const point_type &point) const
Evaluate standard deviation of random field at a point.
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)