1 #ifndef KATOPTRON_RANDOMFIELD_H
2 #define KATOPTRON_RANDOMFIELD_H
4 #include "Stokhos_Sacado_Kokkos_MP_Vector.hpp"
5 #include "Stokhos_KL_ExponentialRandomField.hpp"
9 template <
typename Scalar,
typename Device>
14 typedef typename RandomVariableView::size_type
size_type;
16 typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType
MeshScalar;
17 typedef Stokhos::KL::ExponentialRandomField<MeshScalar, Device>
rf_type;
18 typedef Teuchos::ScalarTraits<MeshScalar>
MST;
30 Teuchos::RCP<SpectralApproach<Scalar, Device>>
sa;
32 RandomField(Teuchos::RCP<Teuchos::ParameterList> randomParams, Kokkos::View<Scalar *, Kokkos::LayoutLeft, Device> _m_rv)
34 isRandom = randomParams->get(
"Is random",
false);
37 g = randomParams->get(
"Mean", MST::one());
38 delta = randomParams->get(
"Dispersion level", MST::zero());
45 const MeshScalar BeginX = randomParams->get(
"Begin X", MST::zero());
46 const MeshScalar BeginY = randomParams->get(
"Begin Y", MST::zero());
47 const MeshScalar BeginZ = randomParams->get(
"Begin Z", MST::zero());
48 const MeshScalar EndX = randomParams->get(
"End X", MST::one());
49 const MeshScalar EndY = randomParams->get(
"End Y", MST::one());
50 const MeshScalar EndZ = randomParams->get(
"End Z", MST::one());
51 num_rv = randomParams->get(
"Number random variables", 0);
52 const MeshScalar correlation_lengthX = randomParams->get(
"Correlation length X", MST::one());
53 const MeshScalar correlation_lengthY = randomParams->get(
"Correlation length Y", MST::one());
54 const MeshScalar correlation_lengthZ = randomParams->get(
"Correlation length Z", MST::one());
56 if (correlation_lengthZ == 0)
59 const int num_KL_X = randomParams->get(
"Number of KL Terms X", 0);
60 const int num_KL_Y = randomParams->get(
"Number of KL Terms Y", 0);
61 const int num_KL_Z = randomParams->get(
"Number of KL Terms Z", 0);
63 mean_type eps = randomParams->get(
"Bound Perturbation Size", 1e-6);
64 mean_type tol = randomParams->get(
"Nonlinear Solver Tolerance", 1e-10);
65 int max_it = randomParams->get(
"Maximum Nonlinear Solver Iterations", 100);
67 Teuchos::ParameterList solverParams;
68 solverParams.set(
"Number of KL Terms",
num_rv);
69 solverParams.set(
"Mean", MST::zero());
70 solverParams.set(
"Standard Deviation", MST::one());
71 solverParams.set(
"Bound Perturbation Size", eps);
72 solverParams.set(
"Nonlinear Solver Tolerance", tol);
73 solverParams.set(
"Maximum Nonlinear Solver Iterations", max_it);
74 Teuchos::Array<MeshScalar> domain_upper(
ndim, 1.0),
75 domain_lower(
ndim, 0.0),
76 correlation_lengths(
ndim, correlation_lengthX);
77 Teuchos::Array<int> num_KL_per_dim(
ndim, 1.0);
79 domain_upper[0] = EndX;
80 domain_upper[1] = EndY;
82 domain_upper[2] = EndZ;
84 domain_lower[0] = BeginX;
85 domain_lower[1] = BeginY;
87 domain_lower[2] = BeginZ;
89 correlation_lengths[0] = correlation_lengthX;
90 correlation_lengths[1] = correlation_lengthY;
92 correlation_lengths[2] = correlation_lengthZ;
94 if (num_KL_X != 0 && num_KL_Y != 0 && num_KL_Z != 0)
96 num_KL_per_dim[0] = num_KL_X;
97 num_KL_per_dim[1] = num_KL_Y;
99 num_KL_per_dim[2] = num_KL_Z;
100 solverParams.set(
"Number of KL Terms per dimension", num_KL_per_dim);
103 solverParams.set(
"Domain Upper Bounds", domain_upper);
104 solverParams.set(
"Domain Lower Bounds", domain_lower);
105 solverParams.set(
"Correlation Lengths", correlation_lengths);
110 std::ofstream outputfile_solverParams;
111 std::ofstream outputfile_m_rv;
113 outputfile_solverParams.open(
"RandomFieldParams.txt");
114 outputfile_m_rv.open(
"RandomVariables.txt");
116 outputfile_solverParams << solverParams << std::endl;
118 for (
size_t i = 0; i <
m_rv.extent(0); ++i)
119 outputfile_m_rv <<
m_rv(i) << std::endl;
121 outputfile_solverParams.close();
122 outputfile_m_rv.close();
126 MeshScalar a = randomParams->get(
"Correlation length", MST::one());
127 MeshScalar r = randomParams->get(
"Wavenumber cutoff", MST::one());
128 size_t muw = randomParams->get(
"Wavenumber discretization points", 40);
129 bool gauss = randomParams->get(
"Gaussian field",
true);
131 unsigned long seed_Z = randomParams->get(
"Z seed", 42);
132 unsigned long seed_Phi = randomParams->get(
"Phi seed", 186);
134 size_t wait_Z = randomParams->get(
"Number Z of previously drawn samples", 0);
135 size_t wait_Phi = randomParams->get(
"Number Phi of previously drawn samples", 0);
137 Kokkos::View<Scalar **, Kokkos::LayoutLeft, Device> Z;
138 Kokkos::View<Scalar **, Kokkos::LayoutLeft, Device> Phi;
140 Z = construct_Z<Scalar, Device>(muw, seed_Z, wait_Z);
141 Phi = construct_Phi<Scalar, Device>(muw, seed_Phi, wait_Phi);
148 KOKKOS_INLINE_FUNCTION
168 Xi =
sa->operator()(x, y);
170 return exp(
a +
b * Xi);
173 return Teuchos::ScalarTraits<Scalar>::one();
177 #endif //KATOPTRON_RANDOMFIELD_H