Commit e0ec98f5 authored by Ole Klein's avatar Ole Klein

Initial commit

parents
cmake_minimum_required(VERSION 2.8.12)
project(dune-randomfield CXX)
if(NOT (dune-common_DIR OR dune-common_ROOT OR
"${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
string(REPLACE ${CMAKE_PROJECT_NAME} dune-common dune-common_DIR
${PROJECT_BINARY_DIR})
endif()
#find dune-common and set the module path
find_package(dune-common REQUIRED)
list(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake/modules"
${dune-common_MODULE_PATH})
#include the dune macros
include(DuneMacros)
# start a dune project with information from dune.module
dune_project()
dune_enable_all_packages()
add_subdirectory(dune)
add_subdirectory(cmake/modules)
add_subdirectory(test)
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
Preparing the Sources
=========================
Additional to the software mentioned in README you'll need the
following programs installed on your system:
cmake >= 2.8.12
Getting started
---------------
If these preliminaries are met, you should run
dunecontrol all
which will find all installed dune modules as well as all dune modules
(not installed) which sources reside in a subdirectory of the current
directory. Note that if dune is not installed properly you will either
have to add the directory where the dunecontrol script resides (probably
./dune-common/bin) to your path or specify the relative path of the script.
Most probably you'll have to provide additional information to dunecontrol
(e. g. compilers, configure options) and/or make options.
The most convenient way is to use options files in this case. The files
define four variables:
CMAKE_FLAGS flags passed to cmake (during configure)
MAKE_FLAGS flags passed to make
An example options file might look like this:
#use this options to autogen, configure and make if no other options are given
CMAKE_FLAGS=" \
-DCMAKE_CXX_COMPILER=g++-4.9 \
-DCMAKE_CXX_FLAGS='-Wall -pedantic' \
-DCMAKE_INSTALL_PREFIX=/install/path" #Force g++-4.9 and set compiler flags
MAKE_FLAGS=install #Per default run make install instead of simply make
If you save this information into example.opts you can pass the opts file to
dunecontrol via the --opts option, e. g.
dunecontrol --opts=example.opts all
More info
---------
See
dunecontrol --help
for further options.
The full build-system is described in the dune-common/doc/buildsystem (Git version) or under share/doc/dune-common/buildsystem if you installed DUNE!
# module providing convenience methods for compiling binaries with FFTW3 support
#
# Provides the following functions:
#
# add_dune_fftw3_flags(target1 target2...)
#
# adds FFTW3 flags to the targets for compilation and linking
function(add_dune_fftw3_flags _targets)
if(FFTW3_FOUND)
foreach(_target ${_targets})
target_link_libraries(${_target} ${FFTW3_LIBRARY} -L${CMAKE_PREFIX_PATH}/lib -lfftw3_mpi)
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${FFTW3_INCLUDE_DIRECTORIES} )
endforeach(_target ${_targets})
endif(FFTW3_FOUND)
endfunction(add_dune_fftw3_flags)
# module providing convenience mehtods for compiling binries with HDF5 support
#
# Provides the following functions:
#
# add_dune_hdf5_flags(target1 target2...)
#
# adds HDF5 flags to the targets for compilation and linking
function(add_dune_hdf5_flags _targets)
if(HDF5_FOUND)
if(HDF5_IS_PARALLEL)
foreach(_target ${_targets})
target_link_libraries(${_target} ${HDF5_LIBRARIES} ${HDF5_C_LIBRARIES})
# set_property(TARGET ${_target} APPEND PROPERTY COMPILE_FLAGS ${HDF5_DEFINITIONS})
set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS H5_USE_16_API)
set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIRS})
endforeach(_target ${_targets})
else(HDF5_IS_PARALLEL)
message(FATAL_ERROR "HDF5 without parallel IO support found!")
endif(HDF5_IS_PARALLEL)
endif(HDF5_FOUND)
endfunction(add_dune_hdf5_flags)
set(modules "DuneRandomfieldMacros.cmake")
install(FILES ${modules} DESTINATION ${DUNE_INSTALL_MODULEDIR})
# File for module specific CMake tests.
find_package(FFTW3 REQUIRED)
include(AddFFTW3Flags)
find_package(HDF5 REQUIRED)
include(AddHDF5Flags)
# - Try to find FFTW3.
# Usage: find_package(FFTW3 [COMPONENTS [single double long-double threads]])
#
# Variables used by this module:
# FFTW3_ROOT_DIR - FFTW3 root directory
# Variables defined by this module:
# FFTW3_FOUND - system has FFTW3
# FFTW3_INCLUDE_DIR - the FFTW3 include directory (cached)
# FFTW3_INCLUDE_DIRS - the FFTW3 include directories
# (identical to FFTW3_INCLUDE_DIR)
# FFTW3[FL]?_LIBRARY - the FFTW3 library - double, single(F),
# long-double(L) precision (cached)
# FFTW3[FL]?_THREADS_LIBRARY - the threaded FFTW3 library - double, single(F),
# long-double(L) precision (cached)
# FFTW3_LIBRARIES - list of all FFTW3 libraries found
# Copyright (C) 2009-2010
# ASTRON (Netherlands Institute for Radio Astronomy)
# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
#
# This file is part of the LOFAR software suite.
# The LOFAR software suite is free software: you can redistribute it and/or
# modify it under the terms of the GNU General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# The LOFAR software suite is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
#
# $Id: FindFFTW3.cmake 15918 2010-06-25 11:12:42Z loose $
# Use double precision by default.
if(FFTW3_FIND_COMPONENTS MATCHES "^$")
set(_components double)
else()
set(_components ${FFTW3_FIND_COMPONENTS})
endif()
# Loop over each component.
set(_libraries)
foreach(_comp ${_components})
if(_comp STREQUAL "single")
list(APPEND _libraries fftw3f)
elseif(_comp STREQUAL "double")
list(APPEND _libraries fftw3)
elseif(_comp STREQUAL "long-double")
list(APPEND _libraries fftw3l)
elseif(_comp STREQUAL "threads")
set(_use_threads ON)
else(_comp STREQUAL "single")
message(FATAL_ERROR "FindFFTW3: unknown component `${_comp}' specified. "
"Valid components are `single', `double', `long-double', and `threads'.")
endif(_comp STREQUAL "single")
endforeach(_comp ${_components})
# If using threads, we need to link against threaded libraries as well.
if(_use_threads)
set(_thread_libs)
foreach(_lib ${_libraries})
list(APPEND _thread_libs ${_lib}_threads)
endforeach(_lib ${_libraries})
set(_libraries ${_thread_libs} ${_libraries})
endif(_use_threads)
# Keep a list of variable names that we need to pass on to
# find_package_handle_standard_args().
set(_check_list)
# Search for all requested libraries.
foreach(_lib ${_libraries})
string(TOUPPER ${_lib} _LIB)
find_library(${_LIB}_LIBRARY ${_lib}
HINTS ${FFTW3_ROOT_DIR} PATH_SUFFIXES lib)
mark_as_advanced(${_LIB}_LIBRARY)
list(APPEND FFTW3_LIBRARIES ${${_LIB}_LIBRARY})
list(APPEND _check_list ${_LIB}_LIBRARY)
endforeach(_lib ${_libraries})
# Search for the header file.
find_path(FFTW3_INCLUDE_DIR fftw3.h
HINTS ${FFTW3_ROOT_DIR} PATH_SUFFIXES include)
mark_as_advanced(FFTW3_INCLUDE_DIR)
list(APPEND _check_list FFTW3_INCLUDE_DIR)
# Handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if
# all listed variables are TRUE
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(FFTW3 DEFAULT_MSG ${_check_list})
/* begin dune-randomfield
put the definitions for config.h specific to
your project here. Everything above will be
overwritten
*/
/* begin private */
/* Name of package */
#define PACKAGE "@DUNE_MOD_NAME@"
/* Define to the address where bug reports for this package should be sent. */
#define PACKAGE_BUGREPORT "@DUNE_MAINTAINER@"
/* Define to the full name of this package. */
#define PACKAGE_NAME "@DUNE_MOD_NAME@"
/* Define to the full name and version of this package. */
#define PACKAGE_STRING "@DUNE_MOD_NAME@ @DUNE_MOD_VERSION@"
/* Define to the one symbol short name of this package. */
#define PACKAGE_TARNAME "@DUNE_MOD_NAME@"
/* Define to the home page for this package. */
#define PACKAGE_URL "@DUNE_MOD_URL@"
/* Define to the version of this package. */
#define PACKAGE_VERSION "@DUNE_MOD_VERSION@"
/* end private */
/* Define to the version of dune-randomfield */
#define DUNE_RANDOMFIELD_VERSION "@DUNE_RANDOMFIELD_VERSION@"
/* Define to the major version of dune-randomfield */
#define DUNE_RANDOMFIELD_VERSION_MAJOR @DUNE_RANDOMFIELD_VERSION_MAJOR@
/* Define to the minor version of dune-randomfield */
#define DUNE_RANDOMFIELD_VERSION_MINOR @DUNE_RANDOMFIELD_VERSION_MINOR@
/* Define to the revision of dune-randomfield */
#define DUNE_RANDOMFIELD_VERSION_REVISION @DUNE_RANDOMFIELD_VERSION_REVISION@
/* end dune-randomfield
Everything below here will be overwritten
*/
prefix=@prefix@
exec_prefix=@exec_prefix@
libdir=@libdir@
includedir=@includedir@
CXX=@CXX@
CC=@CC@
DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: dune-randomfield module
URL: http://dune-project.org/
Requires: dune-common
Libs: -L${libdir}
Cflags: -I${includedir}
################################
# Dune module information file #
################################
#Name of the module
Module: dune-randomfield
Version: 1.0
Maintainer: ole.klein@iwr.uni-heidelberg.de
#depending on
Depends: dune-common
add_subdirectory(randomfield)
#install headers
install(FILES fieldtraits.hh
io.hh
matrix.hh
randomfield.hh
stochastic.hh
trend.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/randomfield
)
This diff is collapsed.
// -*- tab-width: 2; indent-tabs-mode: nil -*-
#ifndef DUNE_RANDOMFIELD_IO_HH
#define DUNE_RANDOMFIELD_IO_HH
/// @todo get rid of this
#define HDF5_DATA_TYPE H5T_IEEE_F64LE //define for 64 bit machine
//#define HDF5_DATA_TYPE H5T_IEEE_F32LE //define for 32 bit machine
namespace Dune{
namespace RandomField{
/**
* @brief Check if file exists
*/
bool fileExists(std::string filename)
{
struct stat fileInfo;
int intStat;
// Attempt to get the file attributes
intStat = stat(filename.c_str(),&fileInfo);
// if attributes were found the file exists
return (intStat == 0);
}
/// @todo take care of missing files
/**
* @brief Read data from an HDF5 file (parallel)
*/
template<typename RF, unsigned int dim>
void readParallelFromHDF5(
std::vector<RF>& local_data,
const std::vector<unsigned int>& local_count,
const std::vector<unsigned int>& local_offset,
const MPI_Comm& communicator,
const std::string& data_name,
const std::string& data_filename)
{
// setup file access template with parallel IO access
hid_t access_pList = H5Pcreate(H5P_FILE_ACCESS);
assert(access_pList > -1);
herr_t status;
MPI_Info mpiInfo = MPI_INFO_NULL;
status = H5Pset_fapl_mpio(access_pList,communicator,mpiInfo);
assert(status > -1);
// open the file for reading
hid_t file_id = H5Fopen(data_filename.c_str(),H5F_ACC_RDONLY,access_pList);
assert(file_id > -1);
// Release file-access template
status = H5Pclose(access_pList);
assert(status > -1);
// open the dataset
hid_t dataset_id = H5Dopen(file_id,data_name.c_str());
assert(dataset_id > -1);
// get the dataspace
hid_t dataspace_id = H5Dget_space(dataset_id);
assert(dataspace_id > -1);
// some needed variables
hsize_t dimData;
hsize_t* dims;
// get the dimension (2D or 3D)
dimData = H5Sget_simple_extent_ndims(dataspace_id);
assert(dimData == dim);
// get the size of the data structure
dims = (hsize_t*)malloc(dim * sizeof (hsize_t));
status = H5Sget_simple_extent_dims(dataspace_id,dims,0);
assert(status > -1);
//set the local, offset, and count as hsize_t, which is needed by the HDF5 routines
hsize_t local_size = 1;
hsize_t offset[dim],count[dim];
for(unsigned int i=0; i < dim; i++ )
{
local_size *= local_count [i];
offset[dim-i-1] = local_offset[i];
count [dim-i-1] = local_count [i];
}
// create the memory space, if something needes to be read on this processor
hid_t memspace_id = 0;
if(local_size != 0)
memspace_id = H5Screate_simple(1,&local_size,NULL);
//select the hyperslab
status = H5Sselect_hyperslab(dataspace_id,H5S_SELECT_SET,offset,NULL,count,NULL);
assert(status > -1);
//resize the return data
local_data.resize(local_size);
// set up the collective transfer properties list
hid_t xferPropList = H5Pcreate(H5P_DATASET_XFER);
assert(xferPropList > -1);
// finally the reading from the file, only if something needes to be read
if(local_size != 0)
{
status = H5Dread(dataset_id,
H5T_NATIVE_DOUBLE,
memspace_id,
dataspace_id,
xferPropList,
&(local_data[0])
);
assert(status > -1);
}
// close the identifiers
H5Dclose(dataset_id);
H5Sclose(dataspace_id);
if(local_size != 0) //this identifier only exists if somethings needs to be read
H5Sclose(memspace_id);
free(dims);
status = H5Fclose( file_id );
assert( status > -1 );
}
/// @todo take care of missing files
/**
* @brief Write data to an HDF5 file (parallel)
*/
template<typename RF, unsigned int dim>
void writeParallelToHDF5(
const std::vector<unsigned int>& global_dim,
const std::vector<RF>& data,
const std::vector<unsigned int>& local_count,
const std::vector<unsigned int>& local_offset,
const MPI_Comm& communicator,
const std::string& data_name,
const std::string& data_filename)
{
//Info variable needed for HDF5
MPI_Info mpiInfo = MPI_INFO_NULL;
herr_t status;
// Set up file access property list with parallel I/O access
hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
H5Pset_fapl_mpio(plist_id,communicator,mpiInfo);
assert(plist_id > -1);
// Create a new file using default properties.
hid_t file_id = H5Fcreate(data_filename.c_str(),
H5F_ACC_TRUNC,
H5P_DEFAULT,
plist_id
);
assert(file_id > -1);
H5Pclose(plist_id);
// set the global size of the grid into a vector of type hsize_t (needed for HDF5 routines)
hsize_t global_dim_HDF5[dim];
for(unsigned int i = 0; i < dim; i++)
global_dim_HDF5[dim-i-1] = global_dim[i];
// set the count and offset in the different dimensions (determine the size of the hyperslab)
// (in hsize_t format, needed for HDF5 routines)
hsize_t count[dim], offset[dim];
for(unsigned int i = 0; i < dim; i++)
{
count[dim-i-1] = local_count [i];
offset[dim-i-1] = local_offset[i];
}
// define the total size of the local data
hsize_t nAllLocalCells = 1;
for (unsigned int i = 0; i < dim; i++)
nAllLocalCells *= count[i];
// Create the dataspace for the dataset.
hid_t filespace = H5Screate_simple(dim,global_dim_HDF5,NULL);
assert(filespace > -1);
// Create the dataset with default properties and close filespace.
hid_t dset_id = H5Dcreate(file_id,data_name.c_str(),HDF5_DATA_TYPE,filespace,H5P_DEFAULT);
H5Sclose(filespace);
assert(dset_id > -1);
//get the memoryspace (but only if something needs to be written on this processor!)
hid_t memspace_id;
if(nAllLocalCells != 0)
{ // -> otherwise HDF5 warning, because of writing nothing!
memspace_id = H5Screate_simple(dim,count,NULL);
assert(memspace_id > -1);
}
// Select hyperslab in the file
filespace = H5Dget_space(dset_id);
H5Sselect_hyperslab(filespace,H5S_SELECT_SET,offset,NULL,count,NULL);
// Create property list for collective dataset write.
plist_id = H5Pcreate(H5P_DATASET_XFER);
H5Pset_dxpl_mpio(plist_id,H5FD_MPIO_COLLECTIVE);
// finally write the data to the disk
// even if nothing should be written H5Dwrite needs to be called!
if(nAllLocalCells != 0)
{ // -> otherwise HDF5 warning, because of writing nothing!
status = H5Dwrite(dset_id,H5T_NATIVE_DOUBLE,memspace_id,filespace,plist_id,&(data[0]));
assert(status > -1);
}
else
{ // IMPORTANT. otherwise the H5Dwrite() blocks!
status = H5Dwrite(dset_id,H5T_NATIVE_DOUBLE,0,filespace,plist_id,&(data[0]));
assert(status > -1);
}
// Close the property list
status = H5Pclose(plist_id);
assert(status > -1);
// Close the filespace
status = H5Sclose(filespace);
assert(status > -1);
//if something written close the memspace
if(nAllLocalCells != 0)
{
status = H5Sclose(memspace_id);
assert(status > -1);
}
// Close the dataset
status = H5Dclose(dset_id);
assert(status > -1);
// Close the file
status = H5Fclose(file_id);
assert(status > -1);
//propably not needed. because the H5Dwrite blocks anyway
MPI_Barrier(communicator);
}
}
}
#endif // DUNE_RANDOMFIELD_IO_HH
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
dune_add_test(SOURCE fieldtest.cc)
target_link_dune_default_libraries(fieldtest)
add_dune_fftw3_flags(fieldtest)
add_dune_hdf5_flags(fieldtest)
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <iostream>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include<dune/common/fmatrix.hh>
#include<dune/common/fvector.hh>
#include<dune/randomfield/randomfield.hh>
template<typename DF, typename RF, unsigned int dimension>
class GridTraits
{
public:
enum {dim = dimension};
//typedef Dune::YaspGrid<dim> Grid;
//typedef typename Grid::LevelGridView GridView;
typedef RF RangeField;
typedef Dune::FieldVector<RF,1> Scalar;
typedef Dune::FieldVector<RF,dim> Vector;
typedef Dune::FieldMatrix<RF,dim,dim> Tensor;
typedef DF DomainField;
typedef Dune::FieldVector<DF,dim> Domain;
typedef Dune::FieldVector<DF,dim-1> IntersectionDomain;
//typedef typename GridView::Traits::template Codim<0>::Entity Element;
//typedef typename GridView::Intersection Intersection;
};
int main(int argc, char** argv)
{
try{
// Maybe initialize MPI
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
std::cout << "Hello World! This is dune-randomfield." << std::endl;
if(Dune::MPIHelper::isFake)
std::cout<< "This is a sequential program." << std::endl;
else
std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
<<" processes!"<<std::endl;
Dune::ParameterTree config;
Dune::ParameterTreeParser parser;
parser.readINITree("randomfield.ini",config);
typedef GridTraits<double,double,2> GridTraits;
typedef Dune::RandomField::ExponentialCovariance Covariance;
/*
GridTraits::Domain center;
center[0] = 0.5;
center[1] = 0.5;
randomFieldList.localize(center,0.15);
randomFieldList.writeToFile("sample");
Dune::RandomField::RandomFieldList<GridTraits,Covariance> copy(randomFieldList);
Dune::RandomField::RandomFieldList<GridTraits,Covariance> copy2(randomFieldList);
Dune::RandomField::RandomFieldList<GridTraits,Covariance> diff(randomFieldList);
copy.timesMatrixRoot();
copy.writeToFile("timesRoot");
copy.timesMatrixRoot();
copy.writeToFile("timesRootSquared");
diff = copy;
copy.timesInverseMatrix();
copy.writeToFile("rootRestored");
copy2.timesMatrix();
copy2.writeToFile("timesMatrix");
diff -= copy2;
diff.writeToFile("matrixDiff");
copy2.timesInverseMatrix();
copy2.writeToFile("matrixRestored");
randomFieldList -= copy;
randomFieldList.writeToFile("difference");
*/
std::cout << "------------------------------" << std::endl;
std::cout << "root inv root" << std::endl;
std::cout << "------------------------------" << std::endl;
for (unsigned int i = 0; i < 10; i++)
{
Dune::RandomField::RandomFieldList<GridTraits,Covariance> randomFieldList(config);
randomFieldList.generateUncorrelated();
const double fieldNorm = std::sqrt(randomFieldList * randomFieldList);
Dune::RandomField::RandomFieldList<GridTraits,Covariance> copy(randomFieldList);
copy.timesMatrixRoot();
copy.timesInverseMatrix();
copy.timesMatrixRoot();
randomFieldList -= copy;
const double diffNorm = std::sqrt(randomFieldList * randomFieldList);
std::cout << "norm: " << fieldNorm << " diffNorm: " << diffNorm << " opNorm: " << diffNorm/fieldNorm << std::endl;
}
std::cout << "------------------------------" << std::endl;
std::cout << "root root inv" << std::endl;
std::cout << "------------------------------" << std::endl;
for (unsigned int i = 0; i < 10; i++)
{
Dune::RandomField::RandomFieldList<GridTraits,Covariance> randomFieldList(config);
randomFieldList.generateUncorrelated();
const double fieldNorm = std::sqrt(randomFieldList * randomFieldList);
Dune::RandomField::RandomFieldList<GridTraits,Covariance> copy(randomFieldList);
copy.timesMatrixRoot();
copy.timesMatrixRoot();
copy.timesInverseMatrix();
randomFieldList -= copy;
const double diffNorm = std::sqrt(randomFieldList * randomFieldList);
std::cout << "norm: " << fieldNorm << " diffNorm: " << diffNorm << " opNorm: " << diffNorm/fieldNorm << std::endl;
}
std::cout << "------------------------------" << std::endl;
std::cout << "root inv root local" << std::endl;
std::cout << "------------------------------" << std::endl;
for (unsigned int i = 0; i < 10; i++)
{
Dune::RandomField::RandomFieldList<GridTraits,Covariance> randomFieldList(config);
randomFieldList.generateUncorrelated();
GridTraits::Domain center;
center[0] = 0.5;
center[1] = 0.5;
randomFieldList.localize(center,0.15);
const double fieldNorm = std::sqrt(randomFieldList * randomFieldList);
Dune::RandomField::RandomFieldList<GridTraits,Covariance> copy(randomFieldList);
copy.timesMatrixRoot();
copy.timesInverseMatrix();
copy.timesMatrixRoot();
randomFieldList -= copy;
const double diffNorm = std::sqrt(randomFieldList * randomFieldList);
std::cout << "norm: " << fieldNorm << " diffNorm: " << diffNorm << " opNorm: " << diffNorm/fieldNorm << std::endl;
}
std::cout << "------------------------------" << std::endl;
std::cout << "root root inv local" << std::endl;
std::cout << "------------------------------" << std::endl;
for (unsigned int i = 0; i < 10; i++)
{
Dune::RandomField::RandomFieldList<GridTraits,Covariance> randomFieldList(config);
randomFieldList.generateUncorrelated();
GridTraits::Domain center;
center[0] = 0.5;
center[1] = 0.5;
randomFieldList.localize(center,0.15);
const double fieldNorm = std::sqrt(randomFieldList * randomFieldList);
Dune::RandomField::RandomFieldList<GridTraits,Covariance> copy(randomFieldList);
copy.timesMatrixRoot();
copy.timesMatrixRoot();
copy.timesInverseMatrix();
randomFieldList -= copy;
const double diffNorm = std::sqrt(randomFieldList * randomFieldList);
std::cout << "norm: " << fieldNorm << " diffNorm: " << diffNorm << " opNorm: " << diffNorm/fieldNorm << std::endl;
}
return 0;
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;