Commit b210151d authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'feature/dune-randomfield' into 'master'

use dune-randomfield for computing the random parameter field

Closes #7

See merge request !6
parents f95e687f ef0170af
......@@ -93,6 +93,7 @@ Depending on your system configuration, there will be more packages necessary to
| [dune-typetree](https://gitlab.dune-project.org/staging/dune-typetree) | releases/2.5
| [dune-pdelab](https://gitlab.dune-project.org/pdelab/dune-pdelab) | releases/2.5
| [dune-testtools](https://gitlab.dune-project.org/quality/dune-testtools) | master
| [dune-randomfield](https://gitlab.dune-project.org/oklein/dune-randomfield) | master
### Optional Packages
......
......@@ -207,9 +207,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="N">
<definition> Suggested size (in cells) of the created random field. This value
may be disregarded by the random generator for small ``correlationLengths`` or
an extreme ``embeddingFactor``. </definition>
<definition> Suggested size (in cells) of the created random field. </definition>
<values> int &times; int (&times; int) </values>
<suggestion> 1000 1000 </suggestion>
</parameter>
......@@ -228,46 +226,12 @@ adding an empty line, make text **bold** or ``monospaced``.
<suggestion> 123456789 </suggestion>
</parameter>
<parameter name="beta">
<definition> Mean value of the resulting field, if ``millerSimilarity`` is used.
</definition>
<values> float </values>
<suggestion> .0 </suggestion>
</parameter>
<parameter name="embeddingFactor">
<definition> Defines how much the raw random field is 'zoomed'. </definition>
<values> float </values>
<suggestion> .01 </suggestion>
</parameter>
<parameter name="variogramModel">
<parameter name="covariance">
<definition> Defines which variogram function should be used to create the
random field. </definition>
<values> spherical, power, expomodulus, exponential, lorentz, gaussian </values>
<values> exponential, gaussian, spherical </values>
<suggestion> gaussian </suggestion>
<comment> spherical, power, expomodulus, exponential, lorentz, gaussian </comment>
</parameter>
<parameter name="newField">
<definition> If ``true``, force recalculation of the random field, even if a similar one
is already cached. </definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="newEV">
<definition> If ``true``, force recalculation of the eigenvalues of the circulant
matrix, even if they are already cached. </definition>
<values> true, false </values>
<suggestion> false </suggestion>
</parameter>
<parameter name="showEV">
<definition> If ``true``, display additional information on the calculated
eigenvalues. </definition>
<values> true, false </values>
<suggestion> false </suggestion>
<comment> exponential, gaussian, spherical </comment>
</parameter>
</category>
......
......@@ -7,5 +7,5 @@ Module: dorie
Version: 0.9
Maintainer: dorieteam@iup.uni-heidelberg.de
#depending on
Depends: dune-pdelab dune-uggrid
Depends: dune-pdelab dune-uggrid dune-randomfield
Suggests: dune-testtools
\ No newline at end of file
/*
* File: FieldData.hh
* Author: A. Ngo (09/2014)
*
* Data structure used to read the input data for the random field to be generated.
* Corresponds to "yfield.cfg".
*
*/
#ifndef FIELD_DATA_HH
#define FIELD_DATA_HH
namespace Dune {
namespace Dorie {
template<typename RF, const int dim>
class FieldData {
public:
Vector<RF> extensions;
Vector<unsigned int> nCells;
Vector<RF> gridsizes;
Vector<RF> correlations;
RF beta;
RF fieldVariance;
std::string variogramModel;
RF embeddingFactor;
bool newField;
bool newEV;
bool showEV;
int seed;
std::string location;
FieldData(){
extensions.resize(dim);
nCells.resize(dim);
gridsizes.resize(dim);
correlations.resize(dim);
};
void read(const Dune::ParameterTree& configuration){
Dune::FieldVector<RF,dim> extensions_ = configuration.get<Dune::FieldVector<RF,dim>>("generator.extensions");
extensions[0] = extensions_[0];
extensions[1] = extensions_[1];
if( dim == 3 )
extensions[2] = extensions_[2];
fieldVariance = configuration.get<RF>("generator.variance");
Dune::FieldVector<RF,dim> nCells_ = configuration.get<Dune::FieldVector<RF,dim>>("generator.fft.N");
nCells[0] = nCells_[0];
nCells[1] = nCells_[1];
if( dim == 3 )
nCells[2] = nCells_[2];
for(int i=0;i<dim;++i)
gridsizes[i] = extensions[i] / RF( nCells[i] );
Dune::FieldVector<RF,dim> correlations_ = configuration.get<Dune::FieldVector<RF,dim>>("generator.fft.correlationLengths");
correlations[0] = correlations_[0];
correlations[1] = correlations_[1];
if( dim == 3 )
correlations[2] = correlations_[2];
beta = configuration.get<RF>("generator.fft.beta");
embeddingFactor = configuration.get<RF>("generator.fft.embeddingFactor");
variogramModel = configuration.get<std::string>("generator.fft.variogramModel");
newField = ("true"==configuration.get<std::string>("generator.fft.newField"));
newEV = ("true"==configuration.get<std::string>("generator.fft.newEV"));
showEV = ("true"==configuration.get<std::string>("generator.fft.showEV"));
seed = configuration.get<int>("generator.fft.seed");
location = configuration.get<std::string>("generator.fft.outputPath");
userInputCheck();
};
void userInputCheck(){
RF nAllCells = nCells[0] * nCells[1];
if(dim == 3)
nAllCells *= nCells[2];
if( nCells[0] <= 0 || nCells[1] <= 0 || nAllCells <= 0 )
DUNE_THROW(IOError,"Random field extensions (N) need to be larger than zero!");
if(nAllCells > 2.5e8)
DUNE_THROW(IOError,"Random fields with more than 2.5e8 cells are not supported!");
//For more information look up https://www.hdfgroup.org/hdf5-quest.html#p2gb
if(correlations[0] <= .0 || correlations[1] <= .0 || (dim == 3 && correlations[2] <= .0))
DUNE_THROW(IOError,"Random field correlation lengths need to be larger than zero!");
};
void printInfos(){
std::cout << "Dimension : " << dim << std::endl;
std::cout << "Domain : " << Vector2CSV( extensions ) << std::endl;
std::cout << "Grid : " << Vector2CSV( nCells ) << std::endl;
std::cout << "cell-size : " << Vector2CSV( gridsizes ) << std::endl;
std::cout << "Variogram model : " << variogramModel << std::endl;
std::cout << "Correlation-lengths : " << Vector2CSV( correlations ) << std::endl;
std::cout << "Field trend : " << beta << std::endl;
std::cout << "Field variance : " << fieldVariance << std::endl;
std::cout << "embedding factor : " << embeddingFactor << std::endl;
std::cout << "seed : " << seed << std::endl;
};
};
}
}
#endif /* FIELD_DATA_HH */
This diff is collapsed.
This diff is collapsed.
/*
* File: Vector.hh
* Author: ngo
*
* Created on July 15, 2010, 2:32 PM
*/
#ifndef _VECTOR_HH
#define _VECTOR_HH
#include <vector>
#include <cmath>
#include <assert.h>
#include <iostream>
#include <sstream>
#include <iomanip>
template<typename ComponentType>
class Vector : public std::vector<ComponentType> {
typedef std::vector<ComponentType> basetype;
public:
Vector() : std::vector<ComponentType>() {
}
Vector( const std::vector<ComponentType>&x )
: std::vector<ComponentType>( x ) {
}
Vector( const Vector& orig ) : basetype(orig) {
}
Vector(const int size, const ComponentType defaultvalue_ = 0)
: std::vector<ComponentType>( size, defaultvalue_ ) {
}
/*
Vector(const int dim, const ComponentType* array_){
for(size_t i=0; i<dim; i++){
ComponentType a = array_[i];
this->push_back(a);
}
}
*/
void reset( const int dim ){
Vector &self = *this;
for( int i=0;i<dim;++i )
self[i] = 0;
}
// inner product with another vector:
ComponentType operator*(Vector & x) {
assert( x.size() == this->size() );
ComponentType sum( 0 );
Vector & self = *this;
for( size_t i = 0; i < this->size(); ++i )
sum += self[i] * x[i];
return sum;
}
// get the product of all the components
ComponentType componentsproduct() const {
ComponentType product( 1 );
const Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
product *= self[i];
return product;
}
// square of the euclidean norm:
ComponentType two_norm_sqr() const {
ComponentType sum( 0 );
const Vector & self = *this;
for (int i = 0; i < (int) this->size(); ++i)
sum += self[i] * self[i];
return sum;
}
// euclidean norm:
ComponentType two_norm() const {
return sqrt(two_norm_sqr());
}
// Manhattan norm:
ComponentType one_norm() const {
ComponentType sum( 0 );
const Vector & self = *this;
for (int i = 0; i < (int) this->size(); ++i)
sum += std::abs( self[i] );
return sum;
}
// Maximum norm:
ComponentType maximum_norm() const {
ComponentType maximum( 0 );
const Vector & self = *this;
for (int i = 0; i < (int) this->size(); ++i)
maximum = std::max( std::abs( self[i] ), maximum );
return maximum;
}
// Arithmetic mean:
ComponentType mean() const {
ComponentType sum( 0 );
const Vector & self = *this;
int N = (int) this->size();
for (int i = 0; i<N; ++i)
sum += self[i];
return REAL(sum) / REAL(N);
}
// Variance:
void mean_and_variance( REAL& mean_value, REAL& variance) const {
ComponentType sum( 0 );
const Vector & self = *this;
mean_value = mean();
int N = (int) this->size();
for (int i = 0; i<N; ++i){
ComponentType summand = self[i] - mean_value;
sum += summand * summand;
}
variance = REAL(sum) / REAL(N-1);
return;
}
// x[i] = value for all components
Vector &operator=(const ComponentType value) {
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] = value;
return *this;
}
// x[i] *= value
Vector &operator*=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] *= value;
return *this;
}
// x[i] /= value
Vector &operator/=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] /= value;
return *this;
}
// x[i] += value
Vector &operator+=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] += value;
return *this;
}
// x += y
Vector &operator+=(const Vector& b) {
Vector& self = *this;
assert( self.size() == b.size() );
for (size_t i = 0; i < this->size(); ++i)
self[i] += b[i];
return *this;
}
// x -= y
Vector &operator-=(const Vector& b) {
Vector& self = *this;
assert( self.size() == b.size() );
for (size_t i = 0; i < this->size(); ++i)
self[i] -= b[i];
return *this;
}
// x[i] -= value
Vector &operator-=(const ComponentType value) {
Vector& self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] -= value;
return *this;
}
Vector operator-(const Vector& b){
const Vector &a = *this;
assert( a.size() == b.size() );
Vector c(a);
assert( c.size() == b.size() );
c.axpy(-1,b);
return c;
};
// x = x + alpha * y :
Vector & axpy(const ComponentType alpha, const Vector & y) {
assert( this->size() == y.size());
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] += alpha * y[i];
return *this;
}
// x = x + alpha * y
Vector & scaled_add(const ComponentType alpha, const Vector & y) {
assert( this->size() == y.size());
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] += alpha * y[i];
return *this;
}
// x = alpha * x + y
Vector & a_times_x_plus_y(const ComponentType alpha, const Vector & y) {
assert( this->size() == y.size());
Vector & self = *this;
for (size_t i = 0; i < this->size(); ++i)
self[i] = alpha * self[i] + y[i];
return *this;
}
void print(void) const {
const Vector &self = *this;
std::cout << "( ";
for (size_t i = 0; i < self.size(); i++) {
std::cout << self[i] << " ";
}
std::cout << " )" << std::endl;
}
};
template<typename ComponentType>
std::ostream & operator <<(std::ostream & os, const Vector<ComponentType> & x) {
os << "( ";
for (int r = 0; r < (int) x.size(); ++r) {
os << x[r] << " ";
}
os << ")" << std::endl;
return os;
}
#endif /* _VECTOR_HH */
typedef double CTYPE;
typedef double REAL;
typedef unsigned int UINT;
typedef int INT;
template<int dim>
struct P0Layout
{
bool contains (Dune::GeometryType gt)
{
if (gt.dim()==dim) return true;
return false;
}
};
struct FEMType {
enum Type{ DG, CG };
};
// define the data type of the HDF5 files here. 32 or 64 bit!!!!
#define HDF5_DATA_TYPE H5T_IEEE_F64LE //define for 64 bit machine
//#define HDF5_DATA_TYPE H5T_IEEE_F32LE //define for 32 bit machine
......@@ -33,22 +33,24 @@
#pragma GCC diagnostic pop
// dorie-rfg includes
#include <dune/dorie-rfg/datatypes.hh>
#include <dune/dorie-rfg/Vector.hh>
#include <dune/dorie-rfg/tools.hh>
#if HDF5_PARALLEL
#include <dune/dorie-rfg/H5Tools_parallel.hh>
#else
#include <dune/dorie-rfg/H5Tools_sequential.hh>
#endif
#include <dune/dorie-rfg/FieldData.hh>
#include <dune/dorie-rfg/random_field_generator.hh>
#include <dune/randomfield/randomfield.hh>
/// Traits for the Random Field
template<unsigned int dimension>
struct GridTraits
{
enum {dim = dimension};
using RangeField = double;
using Scalar = Dune::FieldVector<RangeField,1>;
using DomainField = double;
using Domain = Dune::FieldVector<DomainField,dim>;
};
int main(int argc, char** argv)
{
try{
//Initialize Mpi
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
Dune::MPIHelper::instance(argc, argv);
if (argc!=2)
DUNE_THROW(Dune::IOError,"No parameter file specified!");
......@@ -72,28 +74,37 @@ int main(int argc, char** argv)
if (result != 0)
DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable");
const std::string outputFile = outputPath + "/yfield";
// standard values
inifile["stochastic.anisotropy"] = "axiparallel";
// comply to dune-randomfield config file scheme
inifile["grid.extensions"] = inifile["generator.extensions"];
inifile["grid.cells"] = inifile["generator.fft.N"];
inifile["stochastic.variance"] = inifile["generator.variance"];
inifile["stochastic.corrLength"] = inifile["generator.fft.correlationLengths"];
inifile["stochastic.covariance"] = inifile["generator.fft.covariance"];
// extract seed
const unsigned int seed = inifile.get<unsigned int>("generator.fft.seed");
// Create RFG objects
switch(dim){
case 2:
{
typedef Dune::Dorie::FieldData<REAL,2> FD;
FD fielddata;
fielddata.read(inifile);
fielddata.printInfos();
typedef Dune::Dorie::FFTFieldGenerator<FD,REAL,REAL*,2> YFG;
YFG yfg(fielddata, helper.getCommunicator());
yfg.init();
using Traits = GridTraits<2>;
Dune::RandomField::RandomField<Traits,false,false> field(inifile);
field.generate(seed);
field.writeToFile(outputFile);
}
break;
case 3:
{
typedef Dune::Dorie::FieldData<REAL,3> FD;
FD fielddata;
fielddata.read(inifile);
fielddata.printInfos();
typedef Dune::Dorie::FFTFieldGenerator<FD,REAL,REAL*,3> YFG;
YFG yfg(fielddata, helper.getCommunicator());
yfg.init();
using Traits = GridTraits<3>;
Dune::RandomField::RandomField<Traits,false,false> field(inifile);
field.generate(seed);
field.writeToFile(outputFile);
}
break;
default:
......
This diff is collapsed.
#ifndef DUNE_FUNCEP_FIELDGENTOOLS_HH
#define DUNE_FUNCEP_FIELDGENTOOLS_HH
namespace Dune {
namespace Dorie {
UINT indexconversion_to_1d( const Vector<UINT>& index,
const Vector<UINT>& nCells )
{
// "http://www.fftw.org/doc/Dynamic-Arrays-in-C.html#Dynamic-Arrays-in-C"!
if( nCells.size() > 2 )
return index[0] + nCells[0] * ( index[1] + nCells[1] * index[2] ); // 3d to 1d
else
return index[0] + nCells[0] * index[1]; // 2d to 1d
};
UINT indexconversion_3d_to_1d( const UINT iz
, const UINT iy
, const UINT ix
, const UINT Nz
, const UINT Ny
, const UINT Nx
) {
// "http://www.fftw.org/doc/Dynamic-Arrays-in-C.html#Dynamic-Arrays-in-C"
UINT l = ix + Nx * ( iy + Ny*iz );
return l;
};
// index conversion from 2d-index to 1d-index
UINT indexconversion_2d_to_1d(
const UINT iy
, const UINT ix
, const UINT Ny
, const UINT Nx
)
{
// UINT l = iy * Nx + ix;
UINT l = ix + Nx * iy;
return l;
};
template<typename BVEC, typename VEC>
void copy_to_std( const BVEC& backend_vector, VEC& standard_vector ){
for( auto it = backend_vector.begin(); it!=backend_vector.end(); ++it )
standard_vector.push_back(*it);
}
template<typename BVEC, typename VEC>
void copy_from_std( const VEC& standard_vector, BVEC& backend_vector ){
for( int i=0; i<standard_vector.size(); ++i )
backend_vector.base()[i] = standard_vector[i];
}
template<typename V>
std::string Vector2CSV( const V& v ){
std::stringstream output;
unsigned int j=0;
do{
output << v[j] << ", ";
++j;
} while(j<v.size()-1);
output << v[j];
return output.str();
}
void deleteFile( const std::string filename ){
if( 0 != remove( filename.c_str() ) ){
std::cout << "Error deleting file: " << filename << std::endl;
}
else{
std::cout << "File deleted: " << filename << std::endl;
}
};
/*
* This simple function checks if a file exists.
*/
bool bFileExists( const std::string& filename )
{
struct stat buf;
if( stat( filename.c_str(), &buf ) != -1 ){
UINT filesize = (UINT) buf.st_size;
std::cout << "File found: " << filename
<< ", file size = " << filesize
<< std::endl;
if( filesize > 0 ){
return true;
}
else{
deleteFile( filename );
return false;
}
}
else{
return false;
}
}