Commit 12977e86 authored by Ole Klein's avatar Ole Klein

extend core functionality to arbitrary dimensions (no overlap above 3D)

parent 7188cd0d
......@@ -2,8 +2,9 @@
dune-randomfield provides Gaussian random fields based on
circulant embedding, with the following features:
- support for 1D (processes), 2D and full 3D realizations
of random fields
- support for random fields of arbitrary dimensionality
- data redistribution and parallel overlap for 1D (processes),
2D and 3D realizations of random fields
- exponential, Gaussian, Matérn, spherical and cubic
covariance functions, among others
- axiparallel and full geometric anisotropy as options
......
......@@ -306,20 +306,13 @@ namespace Dune {
template<typename T>
void getFFTData(T& allocLocal, T& localN0, T& local0Start) const
{
if (dim == 3)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0],
(ptrdiff_t)extendedCells[1],(ptrdiff_t)extendedCells[2]};
allocLocal = fftw_mpi_local_size_3d(n[2] , n[1], n[0], comm, &localN0, &local0Start);
}
else if (dim == 2)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0],(ptrdiff_t)extendedCells[1]};
allocLocal = fftw_mpi_local_size_2d(n[1], n[0], comm, &localN0, &local0Start);
}
else if (dim == 1)
ptrdiff_t n[dim];
for (unsigned int i = 0; i < dim; i++)
n[i] = extendedCells[dim-1-i];
if (dim == 1)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0]};
ptrdiff_t localN02, local0Start2;
allocLocal = fftw_mpi_local_size_1d(n[0], comm, FFTW_FORWARD, FFTW_ESTIMATE,
&localN0, &local0Start, &localN02, &local0Start2);
......@@ -327,59 +320,38 @@ namespace Dune {
DUNE_THROW(Dune::Exception,"1d size / offset results don't match");
}
else
DUNE_THROW(Dune::Exception,"dimension of field has to be 1, 2 or 3");
allocLocal = fftw_mpi_local_size(dim, n, comm, &localN0, &local0Start);
}
/**
* @brief Convert an index tuple into a one dimensional encoding
*/
unsigned int indicesToIndex(
template<unsigned int currentDim = 0>
static unsigned int indicesToIndex(
const std::array<unsigned int,dim>& indices,
const std::array<unsigned int,dim>& bound
) const
)
{
if (dim == 3)
{
return indices[0] + bound[0] * (indices[1] + bound[1]*indices[2]);
}
else if (dim == 2)
{
return indices[1] * bound[0] + indices[0];
}
else if (dim == 1)
{
return indices[0];
}
if constexpr(currentDim == dim-1)
return indices[currentDim];
else
DUNE_THROW(Dune::Exception,"dimension of field has to be 1, 2 or 3");
return indices[currentDim]
+ bound[currentDim] * indicesToIndex<currentDim+1>(indices,bound);
}
/**
* @brief Convert a one dimensional encoding into the original index tuple
*/
void indexToIndices(
template<unsigned int currentDim = 0>
static void indexToIndices(
const unsigned int index,
std::array<unsigned int,dim>& indices,
const std::array<unsigned int,dim>& bound
) const
)
{
if (dim == 3)
{
indices[0] = index % bound[0];
indices[1] = (index / bound[0]) % bound[1];
indices[2] = (index / bound[0]) / bound[1];
}
else if (dim == 2)
{
indices[0] = index % bound[0];
indices[1] = index / bound[0];
}
else if (dim == 1)
{
indices[0] = index;
}
else
DUNE_THROW(Dune::Exception,"dimension of field has to be 1, 2 or 3");
indices[currentDim] = index % bound[currentDim];
if constexpr(currentDim != dim-1)
indexToIndices<currentDim+1>(index/bound[currentDim],indices,bound);
}
/**
......
......@@ -68,7 +68,7 @@ namespace Dune{
hsize_t dimData;
hsize_t* dims;
// get the dimension (2D or 3D)
// get the dimension
dimData = H5Sget_simple_extent_ndims(dataspace_id);
assert(dimData == dim);
......@@ -87,7 +87,7 @@ namespace Dune{
count [dim-i-1] = local_count [i];
}
// create the memory space, if something needes to be read on this processor
// create the memory space, if something needs to be read on this processor
hid_t memspace_id = 0;
if(local_size != 0)
memspace_id = H5Screate_simple(1,&local_size,NULL);
......@@ -245,6 +245,127 @@ namespace Dune{
}
#endif //HAVE_HDF5
template<typename RF, unsigned int dim>
bool writeToXDMF(
const std::array<unsigned int,dim>& cells,
const std::array<RF,dim>& extensions,
const std::string& fileName
)
{
if (dim > 4)
return false;
std::ofstream file(fileName+".xdmf",std::ofstream::trunc);
file
<< "<?xml version=\"1.0\" ?>\n"
<< "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n"
<< "<Xdmf Version=\"2.0\">\n"
<< " <Domain>\n";
if (dim < 4)
{
file
<< " <Grid Name=\"StructuredGrid\">\n"
<< " <Topology TopologyType=\"3DRectMesh\" NumberOfElements=\"";
for (unsigned int i = 0; i < dim; i++)
file << cells[dim-(i+1)] << " ";
file
<< "\"/>\n"
<< " <Geometry GeometryType=\"origin_dxdydz\">\n"
<< " <DataItem Dimensions=\"3\">\n"
<< " 0. 0. 0.\n"
<< " </DataItem>\n"
<< " <DataItem Dimensions=\"3\">\n"
<< " ";
for (unsigned int i = 0; i < 3-dim; i++)
file << extensions[0]/cells[0] << " "; // additional entries to visualize 1D and 2D files
for (unsigned int i = 0; i < dim; i++)
file << extensions[i]/cells[i] << " ";
file
<< "\n"
<< " </DataItem>\n"
<< " </Geometry>\n"
<< " <Attribute Name=\""
<< fileName
<< "\" Center=\"Cell\">\n"
<< " <DataItem Dimensions=\"";
for (unsigned int i = 0; i < dim; i++)
file << cells[dim-(i+1)] << " ";
file
<< "\" Format=\"HDF\">\n"
<< " " << fileName+".stoch.h5" << ":/stochastic\n"
<< " </DataItem>\n"
<< " </Attribute>\n"
<< " </Grid>\n";
}
else if (dim == 4)
{
file
<< " <Grid Name=\"GridTime\" GridType=\"Collection\""
<< " CollectionType=\"Temporal\">\n";
for (unsigned int j = 0; j < cells[dim-1]; j++)
{
file
<< " <Grid Name=\"StructuredGrid\">\n"
<< " <Time Value=\"" << j << "\"/>\n"
<< " <Topology TopologyType=\"3DRectMesh\" NumberOfElements=\"";
for (unsigned int i = 0; i < dim-1; i++)
file << cells[dim-(i+2)] << " ";
file
<< "\"/>\n"
<< " <Geometry GeometryType=\"origin_dxdydz\">\n"
<< " <DataItem Dimensions=\"3\">\n"
<< " 0. 0. 0.\n"
<< " </DataItem>\n"
<< " <DataItem Dimensions=\"3\">\n"
<< " ";
for (unsigned int i = 0; i < dim-1; i++)
file << extensions[i]/cells[i] << " ";
file
<< "\n"
<< " </DataItem>\n"
<< " </Geometry>\n"
<< " <Attribute Name=\""
<< fileName
<< "\" Center=\"Cell\">\n"
<< " <DataItem ItemType=\"HyperSlab\" Dimensions=\"";
for (unsigned int i = 0; i < dim-1; i++)
file << cells[dim-(i+2)] << " ";
file
<< "1\">\n"
<< " <DataItem Dimensions=\"3 4\">\n"
<< " " << j << " 0 0 0\n"
<< " 1 1 1 1\n"
<< " 1 ";
for (unsigned int i = 0; i < dim-1; i++)
file << cells[dim-(i+2)] << " ";
file
<< "\n"
<< " </DataItem>\n"
<< " <DataItem Dimensions=\"";
for (unsigned int i = 0; i < dim; i++)
file << cells[dim-(i+1)] << " ";
file
<< "\" Format=\"HDF\">\n"
<< " " << fileName+".stoch.h5" << ":/stochastic\n"
<< " </DataItem>\n"
<< " </DataItem>\n"
<< " </Attribute>\n"
<< " </Grid>\n";
}
file
<< " </Grid>\n";
}
file
<< " </Domain>\n"
<< "</Xdmf>"
<< std::endl;
return true;
}
}
}
......
......@@ -325,7 +325,7 @@ namespace Dune {
for (unsigned int index = 0; index < localExtendedDomainSize; index++)
{
(*traits).indexToIndices(index,indices,localExtendedCells);
Traits::indexToIndices(index,indices,localExtendedCells);
for (unsigned int i = 0; i < dim; i++)
{
......@@ -349,27 +349,12 @@ namespace Dune {
{
fftw_plan plan_forward;
if (dim == 3)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0],
(ptrdiff_t)extendedCells[1],(ptrdiff_t)extendedCells[2]};
plan_forward = fftw_mpi_plan_dft_3d(n[2], n[1], n[0], vector, vector,
(*traits).comm, FFTW_FORWARD, FFTW_ESTIMATE);
}
else if (dim == 2)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0],(ptrdiff_t)extendedCells[1]};
plan_forward = fftw_mpi_plan_dft_2d(n[1], n[0], vector, vector,
(*traits).comm, FFTW_FORWARD, FFTW_ESTIMATE);
}
else if (dim == 1)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0]};
plan_forward = fftw_mpi_plan_dft_1d(n[0], vector, vector,
(*traits).comm, FFTW_FORWARD, FFTW_ESTIMATE);
}
else
DUNE_THROW(Dune::Exception,"dimension of field has to be 1, 2 or 3");
ptrdiff_t n[dim];
for (unsigned int i = 0; i < dim; i++)
n[i] = extendedCells[dim-1-i];
plan_forward = fftw_mpi_plan_dft(dim,n,vector,vector,
(*traits).comm, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(plan_forward);
fftw_destroy_plan(plan_forward);
......@@ -383,27 +368,12 @@ namespace Dune {
{
fftw_plan plan_backward;
if (dim == 3)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0],
(ptrdiff_t)extendedCells[1],(ptrdiff_t)extendedCells[2]};
plan_backward = fftw_mpi_plan_dft_3d(n[2], n[1], n[0], vector, vector,
(*traits).comm, FFTW_BACKWARD, FFTW_ESTIMATE);
}
else if (dim == 2)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0],(ptrdiff_t)extendedCells[1]};
plan_backward = fftw_mpi_plan_dft_2d(n[1], n[0], vector, vector,
(*traits).comm, FFTW_BACKWARD, FFTW_ESTIMATE);
}
else if (dim == 1)
{
ptrdiff_t n[] = {(ptrdiff_t)extendedCells[0]};
plan_backward = fftw_mpi_plan_dft_1d(n[0], vector, vector,
(*traits).comm, FFTW_BACKWARD, FFTW_ESTIMATE);
}
else
DUNE_THROW(Dune::Exception,"dimension of field has to be 1, 2 or 3");
ptrdiff_t n[dim];
for (unsigned int i = 0; i < dim; i++)
n[i] = extendedCells[dim-1-i];
plan_backward = fftw_mpi_plan_dft(dim,n,vector,vector,
(*traits).comm, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(plan_backward);
fftw_destroy_plan(plan_backward);
......@@ -543,8 +513,8 @@ namespace Dune {
std::array<unsigned int,dim> indices;
for (unsigned int index = 0; index < localDomainSize; index++)
{
(*traits).indexToIndices(index,indices,localCells);
const unsigned int extIndex = (*traits).indicesToIndex(indices,localExtendedCells);
Traits::indexToIndices(index,indices,localCells);
const unsigned int extIndex = Traits::indicesToIndex(indices,localExtendedCells);
extendedField[extIndex][0] = field[index];
}
......@@ -571,10 +541,10 @@ namespace Dune {
for (unsigned int index = 0; index < localDomainSize; index++)
{
(*traits).indexToIndices(index,indices,localCells);
Traits::indexToIndices(index,indices,localCells);
const unsigned int offset = i * localExtendedDomainSize/embeddingFactor;
const unsigned int extIndex
= (*traits).indicesToIndex(indices,localExtendedCells) + offset;
= Traits::indicesToIndex(indices,localExtendedCells) + offset;
extendedField[extIndex][0] = localCopy[index];
}
......@@ -600,8 +570,8 @@ namespace Dune {
std::array<unsigned int,dim> indices;
for (unsigned int index = 0; index < localDomainSize; index++)
{
(*traits).indexToIndices(index,indices,localCells);
const unsigned int extIndex = (*traits).indicesToIndex(indices,localExtendedCells);
Traits::indexToIndices(index,indices,localCells);
const unsigned int extIndex = Traits::indicesToIndex(indices,localExtendedCells);
field[index] = extendedField[extIndex][0];
}
......@@ -625,9 +595,9 @@ namespace Dune {
localCopy[i].resize(localDomainSize);
for (unsigned int index = 0; index < localDomainSize; index++)
{
(*traits).indexToIndices(index,indices,localCells);
Traits::indexToIndices(index,indices,localCells);
const unsigned int offset = i * localExtendedDomainSize/embeddingFactor;
const unsigned int extIndex = (*traits).indicesToIndex(indices,localExtendedCells);
const unsigned int extIndex = Traits::indicesToIndex(indices,localExtendedCells);
localCopy[i][index] = extendedField[extIndex + offset][0];
}
......
......@@ -109,7 +109,7 @@ namespace Dune {
const typename Traits::DomainType& coords = elem.geometry().center();
(*traits).coordsToIndices(coords,evalIndices,localEvalOffset);
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
evalVector[index] = vLocal[0];
}
......@@ -129,7 +129,7 @@ namespace Dune {
const typename Traits::DomainType& coords = referenceElement(elem.geometry()).position(0,0);
const typename Traits::DomainType& global = elem.geometry().global(coords);
(*traits).coordsToIndices(global,evalIndices,localEvalOffset);
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
Dune::FieldVector<RF,1> value;
dgf.evaluate(elem,coords,value);
......@@ -177,8 +177,8 @@ namespace Dune {
{
localEvalOffset[0] = rank * localEvalCells[0];
}
else
DUNE_THROW(Dune::Exception,"dimension of field has to be 1, 2 or 3");
else if ((*traits).verbose)
std::cout << "Note: dimension of field has to be 1, 2 or 3 for data redistribution and overlap" << std::endl;
dataVector.resize(localDomainSize);
evalVector.resize(localDomainSize);
......@@ -200,54 +200,12 @@ namespace Dune {
#if HAVE_HDF5
if ((*traits).verbose && rank == 0)
std::cout << "writing random field to file " << fileName << std::endl;
writeParallelToHDF5<RF,dim>((*traits).cells, dataVector, localCells, localOffset,
(*traits).comm, "/stochastic", fileName+".stoch.h5");
if (rank == 0)
{
std::ofstream file(fileName+".xdmf",std::ofstream::trunc);
file
<< "<?xml version=\"1.0\" ?>\n"
<< "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []>\n"
<< "<Xdmf Version=\"2.0\">\n"
<< " <Domain>\n"
<< " <Grid Name=\"StructuredGrid\" GridType=\"Uniform\">\n"
<< " <Topology TopologyType=\"3DRectMesh\" NumberOfElements=\"";
for (unsigned int i = 0; i < dim; i++)
file << cells[dim-(i+1)] << " ";
file
<< "\"/>\n"
<< " <Geometry GeometryType=\"origin_dxdydz\">\n"
<< " <DataItem Dimensions=\"3\" NumberType=\"Float\" Precision=\"4\" Format=\"XML\">\n"
<< " 0. 0. 0.\n"
<< " </DataItem>\n"
<< " <DataItem Dimensions=\"3\" NumberType=\"Float\" Precision=\"4\" Format=\"XML\">\n"
<< " ";
for (unsigned int i = 0; i < 3 - dim; i++)
file << extensions[0]/cells[0] << " "; // additional entries to visualize 1D and 2D files
for (unsigned int i = 0; i < dim; i++)
file << extensions[i]/cells[i] << " ";
file
<< "\n"
<< " </DataItem>\n"
<< " </Geometry>\n"
<< " <Attribute Name=\""
<< fileName
<< "\" AttributeType=\"Scalar\" Center=\"Cell\">\n"
<< " <DataItem Dimensions=\"";
for (unsigned int i = 0; i < dim; i++)
file << cells[dim-(i+1)] << " ";
file
<< "\" NumberType=\"Float\" Precision=\"4\" Format=\"HDF\">\n"
<< " " << fileName+".stoch.h5" << ":/stochastic\n"
<< " </DataItem>\n"
<< " </Attribute>\n"
<< " </Grid>\n"
<< " </Domain>\n"
<< "</Xdmf>"
<< std::endl;
}
writeToXDMF<RF,dim>((*traits).cells,(*traits).extensions,fileName);
#else //HAVE_HDF5
DUNE_THROW(Dune::NotImplemented,"Writing and reading field files requires parallel HDF5 support");
#endif //HAVE_HDF5
......@@ -402,7 +360,7 @@ namespace Dune {
evalIndices[i]--;
}
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
output[0] = evalVector[index];
}
}
......@@ -426,7 +384,7 @@ namespace Dune {
evalIndices[i]--;
}
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
output[0] = evalVector[index];
}
}
......@@ -446,7 +404,7 @@ namespace Dune {
evalIndices[i]--;
}
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
output[0] = evalVector[index];
}
}
......@@ -497,8 +455,8 @@ namespace Dune {
newIndices[1] = 2*oldIndices[1];
newIndices[2] = 2*oldIndices[2];
const unsigned int oldIndex = (*traits).indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = (*traits).indicesToIndex(newIndices,localCells);
const unsigned int oldIndex = Traits::indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = Traits::indicesToIndex(newIndices,localCells);
const RF oldValue = oldData[oldIndex];
dataVector[newIndex ] = oldValue;
......@@ -519,8 +477,8 @@ namespace Dune {
newIndices[0] = 2*oldIndices[0];
newIndices[1] = 2*oldIndices[1];
const unsigned int oldIndex = (*traits).indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = (*traits).indicesToIndex(newIndices,localCells);
const unsigned int oldIndex = Traits::indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = Traits::indicesToIndex(newIndices,localCells);
const RF oldValue = oldData[oldIndex];
dataVector[newIndex ] = oldValue;
......@@ -570,8 +528,8 @@ namespace Dune {
oldIndices[1] = 2*newIndices[1];
oldIndices[2] = 2*newIndices[2];
const unsigned int oldIndex = (*traits).indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = (*traits).indicesToIndex(newIndices,localCells);
const unsigned int oldIndex = Traits::indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = Traits::indicesToIndex(newIndices,localCells);
RF newValue = 0.;
newValue += oldData[oldIndex ];
......@@ -593,8 +551,8 @@ namespace Dune {
oldIndices[0] = 2*newIndices[0];
oldIndices[1] = 2*newIndices[1];
const unsigned int oldIndex = (*traits).indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = (*traits).indicesToIndex(newIndices,localCells);
const unsigned int oldIndex = Traits::indicesToIndex(oldIndices,oldLocalCells);
const unsigned int newIndex = Traits::indicesToIndex(newIndices,localCells);
RF newValue = 0.;
newValue += oldData[oldIndex ];
......@@ -651,8 +609,8 @@ namespace Dune {
for (unsigned int i = 0; i < localDomainSize; i++)
{
(*traits).indexToIndices(i,cellIndices,localCells);
(*traits).indicesToCoords(cellIndices,localOffset,location);
Traits::indexToIndices(i,cellIndices,localCells);
Traits::indicesToCoords(cellIndices,localOffset,location);
distSquared = 0.;
for (unsigned int j = 0; j < dim; j++)
......@@ -851,12 +809,12 @@ namespace Dune {
evalIndices[iNextNext] < localEvalCells[iNextNext]; evalIndices[iNextNext]++)
{
evalIndices[i] = 0;
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
extract[2*i ][evalIndices[iNext] + evalIndices[iNextNext] * localEvalCells[iNext]]
= evalVector[index];
evalIndices[i] = localEvalCells[i] - 1;
const unsigned int index2 = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index2 = Traits::indicesToIndex(evalIndices,localEvalCells);
extract[2*i+1][evalIndices[iNext] + evalIndices[iNextNext] * localEvalCells[iNext]]
= evalVector[index2];
}
......@@ -881,11 +839,11 @@ namespace Dune {
for (evalIndices[iNext] = 0; evalIndices[iNext] < localEvalCells[iNext]; evalIndices[iNext]++)
{
evalIndices[i] = 0;
const unsigned int index = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index = Traits::indicesToIndex(evalIndices,localEvalCells);
extract[2*i ][evalIndices[iNext]] = evalVector[index];
evalIndices[i] = localEvalCells[i] - 1;
const unsigned int index2 = (*traits).indicesToIndex(evalIndices,localEvalCells);
const unsigned int index2 = Traits::indicesToIndex(evalIndices,localEvalCells);
extract[2*i+1][evalIndices[iNext]] = evalVector[index2];
}
}
......
......@@ -259,8 +259,10 @@ void generateFields(const Dune::MPIHelper& helper, const std::string& configFile
generate<2>(config);
else if (extensions.size() == 3)
generate<3>(config);
else if (extensions.size() == 4)
generate<4>(config);
else
DUNE_THROW(Dune::NotImplemented,"dimension (size of grid.extensions) has to be 1, 2 or 3");
DUNE_THROW(Dune::NotImplemented,"dimension (size of grid.extensions) has to be 1, 2, 3 or 4");
}
else
{
......@@ -272,8 +274,10 @@ void generateFields(const Dune::MPIHelper& helper, const std::string& configFile
generateList<2>(config);
else if (extensions.size() == 3)
generateList<3>(config);
else if (extensions.size() == 4)
generateList<4>(config);
else
DUNE_THROW(Dune::NotImplemented,"dimension (size of grid.extensions) has to be 1, 2 or 3");
DUNE_THROW(Dune::NotImplemented,"dimension (size of grid.extensions) has to be 1, 2, 3 or 4");
}
}
......
dune_symlink_to_source_files(FILES randomfield1d.ini randomfield2d.ini randomfield3d.ini)
dune_symlink_to_source_files(FILES randomfield1d.ini randomfield2d.ini
randomfield3d.ini randomfield4d.ini)
foreach(dim 1 2 3)
foreach(dim 1 2 3 4)
set(name fieldtest_${dim})
dune_add_test(SOURCES fieldtest.cc
NAME ${name}
......
......@@ -143,56 +143,6 @@ void runTests(Dune::ParameterTree config, std::string covariance)
diffMultInv(randomField);
}
/**
* @brief 1D version of tests
*/
void test1d()
{
Dune::ParameterTree config;
Dune::ParameterTreeParser parser;
parser.readINITree("randomfield1d.ini",config);
using GridTraits = GridTraits<double,double,1>;
std::cout << "--------------" << std::endl;
std::cout << "1D Exponential" << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"exponential");
std::cout << "--------------" << std::endl;
std::cout << "1D Gaussian " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"gaussian");
std::cout << "--------------" << std::endl;
std::cout << "1D Spherical " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"spherical");
}
/**
* @brief 2D version of tests
*/
void test2d()
{
Dune::ParameterTree config;
Dune::ParameterTreeParser parser;
parser.readINITree("randomfield2d.ini",config);
using GridTraits = GridTraits<double,double,2>;
std::cout << "--------------" << std::endl;
std::cout << "2D Exponential" << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"exponential");
std::cout << "--------------" << std::endl;
std::cout << "2D Gaussian " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"gaussian");
std::cout << "--------------" << std::endl;
std::cout << "2D Spherical " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"spherical");
}
template<unsigned int dim>
void test()
{
......@@ -210,10 +160,13 @@ void test()
std::cout << dim << "D Gaussian " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"gaussian");
std::cout << "--------------" << std::endl;
std::cout << dim << "D Spherical " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"spherical");
if (dim < 4)
{
std::cout << "--------------" << std::endl;
std::cout << dim << "D Spherical " << std::endl;
std::cout << "--------------" << std::endl;
runTests<GridTraits>(config,"spherical");
}
}
int main(int argc, char** argv)
......@@ -227,7 +180,7 @@ int main(int argc, char** argv)
std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
<<" processes!"<<std::endl;
static_assert(DIMENSION >= 1 && DIMENSION <= 3, "only dimension 1, 2 and 3 supported!");
static_assert(DIMENSION >= 1 && DIMENSION <= 4, "only dimension 1, 2, 3 and 4 supported!");
test<DIMENSION>();
......
[grid]
extensions = 1 1 1 1
cells = 8 8 8 8
[stochastic]
variance = 1