waves
Basic FE playground
wRandomField.h
Go to the documentation of this file.
1 #ifndef KATOPTRON_RANDOMFIELD_H
2 #define KATOPTRON_RANDOMFIELD_H
3 
4 #include "Stokhos_Sacado_Kokkos_MP_Vector.hpp"
5 #include "Stokhos_KL_ExponentialRandomField.hpp"
6 #include "EnsembleTraits.h"
7 #include "wSpectralApproach.h"
8 
9 template <typename Scalar, typename Device>
11 {
12 public:
13  typedef Kokkos::View<Scalar *, Kokkos::LayoutLeft, Device> RandomVariableView;
14  typedef typename RandomVariableView::size_type size_type;
15 
16  typedef typename Teuchos::ScalarTraits<Scalar>::coordinateType MeshScalar;
17  typedef Stokhos::KL::ExponentialRandomField<MeshScalar, Device> rf_type;
18  typedef Teuchos::ScalarTraits<MeshScalar> MST;
20 
21  bool isRandom = false;
22  bool isExpRandom = true;
23  rf_type m_rf; // Exponential random field
24  RandomVariableView m_rv; // KL random variables
25  mean_type g; // Mean
26  mean_type delta; // Dispersion level
28  int num_rv;
29  int ndim = 3; // Dimension of the physical domain
30  Teuchos::RCP<SpectralApproach<Scalar, Device>> sa;
31 
32  RandomField(Teuchos::RCP<Teuchos::ParameterList> randomParams, Kokkos::View<Scalar *, Kokkos::LayoutLeft, Device> _m_rv)
33  {
34  isRandom = randomParams->get("Is random", false);
35  if (isRandom)
36  {
37  g = randomParams->get("Mean", MST::one());
38  delta = randomParams->get("Dispersion level", MST::zero());
39  a = log(g / (sqrt(1 + delta * delta)));
40  b = sqrt(log(1 + delta * delta));
41 
42  isExpRandom = randomParams->get("Is exp", true);
43  if (isExpRandom)
44  {
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());
55 
56  if (correlation_lengthZ == 0)
57  ndim = 2;
58 
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);
62 
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);
66 
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);
78 
79  domain_upper[0] = EndX;
80  domain_upper[1] = EndY;
81  if (ndim == 3)
82  domain_upper[2] = EndZ;
83 
84  domain_lower[0] = BeginX;
85  domain_lower[1] = BeginY;
86  if (ndim == 3)
87  domain_lower[2] = BeginZ;
88 
89  correlation_lengths[0] = correlation_lengthX;
90  correlation_lengths[1] = correlation_lengthY;
91  if (ndim == 3)
92  correlation_lengths[2] = correlation_lengthZ;
93 
94  if (num_KL_X != 0 && num_KL_Y != 0 && num_KL_Z != 0)
95  {
96  num_KL_per_dim[0] = num_KL_X;
97  num_KL_per_dim[1] = num_KL_Y;
98  if (ndim == 3)
99  num_KL_per_dim[2] = num_KL_Z;
100  solverParams.set("Number of KL Terms per dimension", num_KL_per_dim);
101  }
102 
103  solverParams.set("Domain Upper Bounds", domain_upper);
104  solverParams.set("Domain Lower Bounds", domain_lower);
105  solverParams.set("Correlation Lengths", correlation_lengths);
106 
107  m_rf = rf_type(solverParams);
108  m_rv = _m_rv;
109 
110  std::ofstream outputfile_solverParams;
111  std::ofstream outputfile_m_rv;
112 
113  outputfile_solverParams.open("RandomFieldParams.txt");
114  outputfile_m_rv.open("RandomVariables.txt");
115 
116  outputfile_solverParams << solverParams << std::endl;
117 
118  for (size_t i = 0; i < m_rv.extent(0); ++i)
119  outputfile_m_rv << m_rv(i) << std::endl;
120 
121  outputfile_solverParams.close();
122  outputfile_m_rv.close();
123  }
124  else
125  {
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);
130 
131  unsigned long seed_Z = randomParams->get("Z seed", 42);
132  unsigned long seed_Phi = randomParams->get("Phi seed", 186);
133 
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);
136 
137  Kokkos::View<Scalar **, Kokkos::LayoutLeft, Device> Z;
138  Kokkos::View<Scalar **, Kokkos::LayoutLeft, Device> Phi;
139 
140  Z = construct_Z<Scalar, Device>(muw, seed_Z, wait_Z);
141  Phi = construct_Phi<Scalar, Device>(muw, seed_Phi, wait_Phi);
142 
143  sa = Teuchos::rcp(new SpectralApproach<Scalar, Device>(Z, Phi, a, r, muw, gauss));
144  }
145  }
146  }
147 
148  KOKKOS_INLINE_FUNCTION
149  Scalar operator()(const MeshScalar x, const MeshScalar y, const MeshScalar z) const
150  {
151  if (isRandom)
152  {
153  Scalar Xi = 0.;
154  if (isExpRandom)
155  {
156  if (ndim == 3)
157  {
158  const MeshScalar point[3] = {x, y, z};
159  Xi = m_rf.evaluate(point, m_rv);
160  }
161  else
162  {
163  const MeshScalar point[2] = {x, y};
164  Xi = m_rf.evaluate(point, m_rv);
165  }
166  }
167  else
168  Xi = sa->operator()(x, y);
169 
170  return exp(a + b * Xi);
171  }
172  else
173  return Teuchos::ScalarTraits<Scalar>::one();
174  }
175 };
176 
177 #endif //KATOPTRON_RANDOMFIELD_H
RandomField::MeshScalar
Teuchos::ScalarTraits< Scalar >::coordinateType MeshScalar
Definition: wRandomField.h:16
RandomField::size_type
RandomVariableView::size_type size_type
Definition: wRandomField.h:14
RandomField::RandomField
RandomField(Teuchos::RCP< Teuchos::ParameterList > randomParams, Kokkos::View< Scalar *, Kokkos::LayoutLeft, Device > _m_rv)
Definition: wRandomField.h:32
RandomField::rf_type
Stokhos::KL::ExponentialRandomField< MeshScalar, Device > rf_type
Definition: wRandomField.h:17
RandomField::sa
Teuchos::RCP< SpectralApproach< Scalar, Device > > sa
Definition: wRandomField.h:30
RandomField::ndim
int ndim
Definition: wRandomField.h:29
RandomField
Definition: wRandomField.h:10
wSpectralApproach.h
SpectralApproach
A class that constructs realization(s) of Gaussian or exponential scalar random field based on the sp...
Definition: wSpectralApproach.h:63
RandomField::g
mean_type g
Definition: wRandomField.h:25
EnsembleTraits
Definition: EnsembleTraits.h:8
RandomField::MST
Teuchos::ScalarTraits< MeshScalar > MST
Definition: wRandomField.h:18
RandomField::delta
mean_type delta
Definition: wRandomField.h:26
RandomField::RandomVariableView
Kokkos::View< Scalar *, Kokkos::LayoutLeft, Device > RandomVariableView
Definition: wRandomField.h:13
RandomField::operator()
KOKKOS_INLINE_FUNCTION Scalar operator()(const MeshScalar x, const MeshScalar y, const MeshScalar z) const
Definition: wRandomField.h:149
RandomField::m_rv
RandomVariableView m_rv
Definition: wRandomField.h:24
RandomField::isRandom
bool isRandom
Definition: wRandomField.h:21
EnsembleTraits.h
RandomField::mean_type
EnsembleTraits< Scalar >::value_type mean_type
Definition: wRandomField.h:19
RandomField::isExpRandom
bool isExpRandom
Definition: wRandomField.h:22
RandomField::num_rv
int num_rv
Definition: wRandomField.h:28
RandomField::m_rf
rf_type m_rf
Definition: wRandomField.h:23
RandomField::b
mean_type b
Definition: wRandomField.h:27
RandomField::a
mean_type a
Definition: wRandomField.h:27