ug4
Loading...
Searching...
No Matches
lognormal_random_field_impl.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 * Authors: Ivo Muha, Martin Rupp
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33/*
34 * Method described in Comput Visual Sci (2006) 9: 1�10, DOI 10.1007/s00791-006-0012-2
35 * Simulation of lognormal random fields with varying resolution scale and local average for Darcy flow
36 * For the generation of the d-dimensional Gaussian random field we choose
37 * a simple spectral method which can sim- ulate anisotropicly correlated fields.
38 * The Gaussian random field f (x) is realized as a superposition of a large number
39 * of randomly chosen harmonic modes following the method introduced by Kraichnan
40 *
41 */
42#ifndef __H__UG__LIB_DISC__SPATIAL_DISC__LOGNORMAL_RANDOM_FIELD_IMPL__
43#define __H__UG__LIB_DISC__SPATIAL_DISC__LOGNORMAL_RANDOM_FIELD_IMPL__
44
47#include "common/math/misc/math_util.h" // urand
48
49#ifndef M_PI
50#define M_PI 3.14159265358979323846264338327950288 /* pi */
51#endif
52
53namespace ug{
54
55template <typename TData, int dim, typename TRet>
56TRet LognormalRandomField<TData,dim,TRet>::evaluate(TData& D, const MathVector<dim>& x, number time, int si) const
57{
58 double k = eval_K(x);
59 for(size_t i = 0; i < dim; ++i)
60 {
61 for(size_t j = 0; j < dim; ++j)
62 {
63 D[i][j] = 0.0;
64 if (i==j)
65 D[i][j] = k;
66 }
67 }
68 return;
69}
70
71template <typename TData, int dim, typename TRet>
73{
74 // from Numerical Recipes
75 int iset;
76 static double gset;
77 double fac, rsq, v1, v2, x1, x2;
78
79 iset = 0;
80 if (iset == 0)
81 {
82 do
83 {
84 x1 = urand(0.0, 1.0);
85 x2 = urand(0.0, 1.0);
86 v1 = 2.0 * x1 - 1.0;
87 v2 = 2.0 * x2 - 1.0;
88 rsq = v1 * v1 + v2 * v2;
89 } while (rsq >= 1.0 || rsq == 0.0);
90 fac = sqrt(-2.0 * log(rsq) / rsq);
91 gset = v1 * fac;
92 iset = 1;
93 return v2 * fac;
94 }
95 else
96 {
97 iset = 0;
98 return gset;
99 }
100
101}
102
103template<typename TData, int dim, typename TRet>
105{
106 return urand(0.0, 1.0);
107}
108
109template <typename TData, int dim, typename TRet>
111{
112
113 double result = 0.0;
114 for(int i = 0; i < m_N; i++)
115 result += cos(VecDot(m_vRandomQvec[i], x) + m_vRandomAlpha[i]);
116
117
118 double f = m_dMean_f + sqrt(2*m_dSigma_f*m_dSigma_f/m_N)*result;
119
120 if(m_bNoExp)
121 return f;
122 else
123 return exp(f);
124}
125
126template <typename TData, int dim, typename TRet>
127void LognormalRandomField<TData,dim,TRet>::set_config(size_t N, double mean_f, double sigma_f, double sigma)
128{
129 m_N = N;
130 for(int j=0; j<dim; j++)
131 m_sigma[j] = sqrt(1.0/sigma);
132 m_dMean_f = mean_f;
133 m_dSigma_f = sigma_f;
134 m_dSigma = sigma;
135
136 m_vRandomQvec.clear();
137 m_vRandomAlpha.clear();
138
140
141 m_vRandomQvec.resize(N);
142 for(int j=0; j<dim; j++)
143 for(int i = 0; i < m_N; i++)
144 m_vRandomQvec[i][j] = m_sigma[j]*gasdev();
145
146
147 for(int i = 0; i < m_N; i++)
148 m_vRandomAlpha.push_back(undev()*2*M_PI);
149
150// UG_LOG("corrx = " << m_sigma[0] << " corry = " << m_sigma[1] << "\n");
151// PRINT_VECTOR(m_vRandomQvec, "m_vRandomQvec");
152// PRINT_VECTOR(m_vRandomAlpha, "m_vRandomAlpha");
153
154}
155
156template <typename TData, int dim, typename TRet>
158{
159 std::stringstream ss;
160
161 //ss << "LognormalRandomField < TData = " << TypeName<TData>() << ", dim = " << dim << ", TRet = " << TypeName<TRet>() << " >\n";
162 ss << " LognormalRandomField<" << dim << "d> ( m_N = " << m_N << ", m_dMean_f = " << m_dMean_f <<
163 ", m_dSigma_f = " << m_dSigma_f << ", sigma = " << m_dSigma << ", m_bNoExp = " << TrueFalseString(m_bNoExp) << ")";
164 return ss.str();
165
166}
167
168
169} // ug
170
171#endif // __H__UG__LIB_DISC__SPATIAL_DISC__LOGNORMAL_RANDOM_FIELD_IMPL__
double gasdev()
Definition lognormal_random_field_impl.h:72
std::string config_string() const
Definition lognormal_random_field_impl.h:157
TRet evaluate(TData &D, const MathVector< dim > &x, number time, int si) const
Definition lognormal_random_field_impl.h:56
double undev()
Definition lognormal_random_field_impl.h:104
void set_config(size_t N, double mean_f, double sigma_f, double sigma)
Definition lognormal_random_field_impl.h:127
double eval_K(const MathVector< dim > &x) const
Definition lognormal_random_field_impl.h:110
a mathematical Vector with N entries.
Definition math_vector.h:97
double number
Definition types.h:124
TNumber urand(TNumber lowerBound, TNumber upperBound)
uniform distributed random numbers in [lowerBound, upperBound[. Use srand to set a seed.
Definition math_util_impl.hpp:67
vector_t::value_type VecDot(const vector_t &v1, const vector_t &v2)
returns the dot-product of two vector_ts
Definition math_vector_functions_common_impl.hpp:385
the ug namespace
const char * TrueFalseString(bool b)
Definition string_util.h:437