dorie-rfg.cc 3.24 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// -*- tab-width: 4; indent-tabs-mode: nil -*-
/** \file

*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

// common includes
#include <random>
#include <fstream>
#include <vector>
#include <fftw3.h>
#include <fftw3-mpi.h>
#include <time.h>
#include <hdf5.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
Dion Haefner's avatar
Dion Haefner committed
20 21
#include <assert.h>
#include <sstream>
Dion Haefner's avatar
Dion Haefner committed
22 23

// DUNE includes
Dion Haefner's avatar
Dion Haefner committed
24 25 26
// Do not treat DUNE warnings as errors
#pragma GCC diagnostic push
#pragma GCC diagnostic warning "-Wall"
Dion Haefner's avatar
Dion Haefner committed
27 28 29 30 31
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
#include <dune/grid/common/gridenums.hh>
Dion Haefner's avatar
Dion Haefner committed
32 33
#include <dune/pdelab/common/geometrywrapper.hh>
#pragma GCC diagnostic pop
Dion Haefner's avatar
Dion Haefner committed
34 35

// dorie-rfg includes
36 37 38 39 40 41 42 43 44 45 46 47
#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>;
};
Dion Haefner's avatar
Dion Haefner committed
48 49 50 51 52

int main(int argc, char** argv)
{
	try{
		//Initialize Mpi
53
		Dune::MPIHelper::instance(argc, argv);
Dion Haefner's avatar
Dion Haefner committed
54 55 56 57 58

		if (argc!=2)
			DUNE_THROW(Dune::IOError,"No parameter file specified!");
		const std::string inifilename = argv[1];

59 60 61 62 63
		#if !(HDF5_PARALLEL)
			if (helper.size() > 1)
				DUNE_THROW(Dune::Exception,"Parallel HDF5 library is needed to run in parallel");
		#endif

Dion Haefner's avatar
Dion Haefner committed
64 65 66 67
		// Read ini file
		Dune::ParameterTree inifile;
		Dune::ParameterTreeParser ptreeparser;
		ptreeparser.readINITree(inifilename,inifile);
68
		const unsigned int dim = inifile.get<unsigned int>("generator.dimensions");
Dion Haefner's avatar
Dion Haefner committed
69 70 71 72 73 74 75 76

		// Attempt to create output directory
		const std::string outputPath = inifile.get<std::string>("generator.fft.outputPath");
		mkdir(outputPath.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
	  int result = access(outputPath.c_str(), W_OK);
	  if (result != 0)
	    DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable");

77 78
	  const std::string outputFile = outputPath + "YField";

79 80
		// standard values
		inifile["stochastic.anisotropy"] = "axiparallel";
81

82
	  // comply to dune-randomfield config file scheme
83 84 85
	  inifile["grid.extensions"] = inifile["generator.extensions"];
	  inifile["grid.cells"] = inifile["generator.fft.N"];
	  inifile["stochastic.variance"] = inifile["generator.variance"];
86
	  inifile["stochastic.corrLength"] = inifile["generator.fft.correlationLengths"];
87
	  inifile["stochastic.covariance"] = inifile["generator.fft.covariance"];
88

Dion Haefner's avatar
Dion Haefner committed
89 90 91 92
		// Create RFG objects
		switch(dim){
			case 2:
			{
93 94 95 96
				using Traits = GridTraits<2>;
				Dune::RandomField::RandomField<Traits,false,false> field(inifile);
				field.generate();
				field.writeToFile(outputFile);
Dion Haefner's avatar
Dion Haefner committed
97 98 99 100
			}
			break;
			case 3:
			{
101 102 103 104
				using Traits = GridTraits<3>;
				Dune::RandomField::RandomField<Traits,false,false> field(inifile);
  			field.generate();
  			field.writeToFile(outputFile);
Dion Haefner's avatar
Dion Haefner committed
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
			}
			break;
			default:
				DUNE_THROW(Dune::NotImplemented,"Only 2 and 3-dimensional fields are supported");
		}

		return 0;
	}

	catch (Dune::Exception &e)
	{
		std::cerr << "Dune reported error: " << e << std::endl;
		return 1;
  }
  catch (...)
	{
		std::cerr << "Unknown exception thrown!" << std::endl;
		throw;
		return 1;
  }
}