Dear Dumux developers,
I have a problem with random permeability simulations in parallel mode.
I’ve performed a test case similar to "tutorialproblem" with a random
permeability media*.* The intrinsic permeability field is provided by a
data file (here permeab.dat).
It works correctly for sequential simulations. My domain (2D,isotropic)
consists of 128000 cells (160 cells in X and 800 cells in Y direction)
corresponding to 128961 vertices.
Here is the code to distribute the permeability for cells:
###################
const Dune::FieldMatrix<Scalar, dim, dim> intrinsicPermeability(const
Element &element, /*@\label{tutorial-coupled:permeability}@*/
const FVElementGeometry
&fvGeometry,
const int scvIdx) const
{
Dune::FieldMatrix<Scalar, dim, dim> permLoc_;
Scalar permE = perm_[indexSet_.index(*(element.template
subEntity<dim> (scvIdx)))];
for (int i = 0; i < dim; i++){
permLoc_[i][i] = permE;
}
return permLoc_;
}
###################
where perm_ ( initialed by perm_(gridView.size(dim),0.0)) is read from
2x128000 array of the file permeab.dat.
Now, I try to test it in parallel mode shared by 4 processors.
The code will use local index and "local grid" (gridView.size returns the
size of the grid occupied by each processor) to read data from data file,
so only about ~1/4*128000=32000 first values of permeab.dat are read.
Therefore, 4 sub-domains have the same permeability field, it's not correct
because it leads to a bad description of permeability field of the whole
domain.
I guess that I must use ParallelIndexSet but since I am a new DUNE/Dumux
user, it's difficult for me to find the solution.
Do you have any documentations or examples concerning this problem?
Please find attached the tutorialspatialparams file.
Thank you in advance for your helps,
Regards,
Tri Dat
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program 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 2 of the License, or *
* (at your option) any later version. *
* *
* This program 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 this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief The spatial parameters for the fully coupled tutorial problem
* which uses the twophase box model.
*/
#ifndef DUMUX_TUTORIAL_SPATIAL_PARAMS_COUPLED_HH
#define DUMUX_TUTORIAL_SPATIAL_PARAMS_COUPLED_HH
// include parent spatial parameters
#include <dumux/material/spatialparams/implicitspatialparams.hh>
//#include <dumux/material/spatialparams/gstatrandompermeability.hh>
// include material laws
#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
/*@\label{tutorial-coupled:rawLawInclude}@*/
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
/*! \brief providing a random permeability field for 1-D, 2-D or 3-D using
gstat for gaussian simulation.
* \author Jochen Fritz
* gstat is an open source software tool which can (among other things)
generate
* geostatistical random fields. This functionality is used for this class.
See www.gstat.org.
* To use this class, unpack the zipped gstat tarball in the external
directory and build gstat
* according to the README file. You have to provide a control file for gstat
(examples can be found
* in the gstat tarball in the subdirectory DUMUX_stuff).
*/
#include<iostream>
#include <sstream>
#include<vector>
#include<set>
#include<map>
#include<string.h>
#include<stdio.h>
#include<stdlib.h>
namespace Dumux {
//forward declaration
template<class TypeTag>
class TutorialSpatialParamsCoupled;
namespace Properties
{
// The spatial parameters TypeTag
NEW_TYPE_TAG(TutorialSpatialParamsCoupled);/*@\label{tutorial-coupled:define-spatialparameters-typetag}@*/
// Set the spatial parameters
SET_TYPE_PROP(TutorialSpatialParamsCoupled, SpatialParams,
Dumux::TutorialSpatialParamsCoupled<TypeTag>);
/*@\label{tutorial-coupled:set-spatialparameters}@*/
// Set the material law
SET_PROP(TutorialSpatialParamsCoupled, MaterialLaw)
{
private:
// material law typedefs
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
// select material law to be used
typedef RegularizedBrooksCorey<Scalar> RawMaterialLaw;
/*@\label{tutorial-coupled:rawlaw}@*/
public:
// adapter for absolute law
typedef EffToAbsLaw<RawMaterialLaw> type;
/*@\label{tutorial-coupled:eff2abs}@*/
};
}
/*!
* \ingroup TwoPBoxModel
*
* \brief The spatial parameters for the fully coupled tutorial problem
* which uses the two phase box model.
*/
template<class TypeTag>
class TutorialSpatialParamsCoupled: public ImplicitSpatialParams<TypeTag>
/*@\label{tutorial-coupled:tutorialSpatialParameters}@*/
{
// Get informations for current implementation via property system
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GridView::ctype CoordScalar;
template<int dim>
struct ElementLayout
{
bool contains(Dune::GeometryType gt)
{
return gt.dim() == dim;
}
};
enum
{
dim = Grid::dimension,
dimWorld = Grid::dimensionworld
};
// Get object types for function arguments
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry)
FVElementGeometry;
typedef typename Grid::Traits::template Codim<0>::Entity Element;
// typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dim> LocalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
typedef Dune::FieldVector<CoordScalar,dimWorld> DimVector;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables)
ElementVolumeVariables;
typedef std::vector<Scalar> PermeabilityType;
typedef typename GET_PROP(TypeTag, ParameterTree) ParameterTree;
typedef Dune::FieldVector<Scalar, 1> PermType;
typedef Dune::BlockVector<Dune::FieldVector<Scalar,1> > RepresentationType;
typedef typename Grid::Traits::template Codim<0>::EntityPointer
ElementPointer;
typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, ElementLayout>
ElementMapper;
public:
// get material law from property system
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
// determine appropriate parameters depending on selected materialLaw
typedef typename MaterialLaw::Params MaterialLawParams;
/*@\label{tutorial-coupled:matLawObjectType}@*/
// NTD_randompermeability
typedef std::vector<MaterialLawParams> MaterialLawParamsVector;
/*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
* on the position in the domain
*
* \param element The finite volume element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* Alternatively, the function intrinsicPermeabilityAtPos(const
GlobalPosition& globalPos)
* could be defined, where globalPos is the vector including the global
coordinates
* of the finite volume.
*/
// Dumux V2.5 - with reference - NTD
/*
const Dune::FieldMatrix<Scalar, dim, dim> &intrinsicPermeability(const
Element &element,
const FVElementGeometry
&fvGeometry,
const int scvIdx) const
*/
// Dumux V2.6 - without_reference
const Dune::FieldMatrix<Scalar, dim, dim> intrinsicPermeability(const
Element &element, /*@\label{tutorial-coupled:permeability}@*/
const FVElementGeometry
&fvGeometry,
const int scvIdx) const
{
// return permeability_;
Dune::FieldMatrix<Scalar, dim, dim> permLoc_;
Scalar permE = perm_[indexSet_.index(*(element.template subEntity<dim>
(scvIdx)))];
for (int i = 0; i < dim; i++){
permLoc_[i][i] = permE;
}
return permLoc_;
}
/*! Defines the porosity \f$[-]\f$ of the porous medium depending
* on the position in the domain
*
* \param element The finite volume element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* Alternatively, the function porosityAtPos(const GlobalPosition&
globalPos)
* could be defined, where globalPos is the vector including the global
coordinates
* of the finite volume.
*/
Scalar porosity(const Element &element,
/*@\label{tutorial-coupled:porosity}@*/
const FVElementGeometry &fvGeometry,
const int scvIdx) const
{
return porosity_;
}
/*! Returns the parameter object for the material law (i.e. Brooks-Corey)
* depending on the position in the domain
*
* \param element The finite volume element
* \param fvGeometry The finite-volume geometry in the box scheme
* \param scvIdx The local vertex index
*
* Alternatively, the function materialLawParamsAtPos(const
GlobalPosition& globalPos)
* could be defined, where globalPos is the vector including the global
coordinates
* of the finite volume.
*/
const MaterialLawParams& materialLawParams(const Element &element,
/*@\label{tutorial-coupled:matLawParams}@*/
const FVElementGeometry
&fvGeometry,
const int scvIdx) const
{
return materialParams_;
}
// constructor
TutorialSpatialParamsCoupled(const GridView &gridView) :
ImplicitSpatialParams<TypeTag>(gridView),
indexSet_(gridView.indexSet()),
perm_(gridView.size(dim),0.0)
{
try
{
poro_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.Porosity);
}
catch (Dumux::ParameterException &e) {
std::cerr << e << ". Abort!\n";
exit(1) ;
}
catch (...) {
std::cerr << "Unknown exception thrown!\n";
exit(1);
}
//loadIntrinsicPermeability(const GridView& gridView);
loadperm(gridView);
porosity_ = poro_;
//set residual saturations
materialParams_.setSwr(0.0);
/*@\label{tutorial-coupled:setLawParams}@*/
materialParams_.setSnr(0.0);
//parameters of Brooks & Corey Law
materialParams_.setPe(0.0);
materialParams_.setLambda(2.0);
}
void loadperm(const GridView& gridView)
{
const unsigned size = gridView.size(dim);
std::cout << "Size=" << size << std::endl;
perm_.resize(size, 0.0);
std::string permeabilityFileName("permeab.dat");
if
(ParameterTree::tree().hasKey("SpatialParams.PermeabilityInputFileName"))
permeabilityFileName = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
std::string, SpatialParams,
PermeabilityInputFileName);
char concd[100];
strcpy(concd, "permeab.dat");
std::ifstream infile(permeabilityFileName);
std::cout << "Read permeability data from " << concd << std::endl;
std::string zeile;
ElementMapper elementMapper_=gridView;
ElementIterator eItEnd = gridView.template end<0> ();
Scalar totalVolume = 0;
Scalar meanPermeability = 0;
for (ElementIterator eIt = gridView.template begin<0>(); eIt !=
eItEnd; ++eIt)
{
int globalIdxI = elementMapper_.map(*eIt);
Scalar dummy1, dummy2, permi;
std::getline(infile, zeile);
std::istringstream ist(zeile);
if (dim == 1)
ist >> permi;
else if (dim == 2){
ist >> dummy1 >> permi;
}
else if (dim == 3)
ist >> dummy1 >> dummy2 >> permi;
else
DUNE_THROW(Dune::NotImplemented, "random permeability is not
implemented for dimension > 3");
perm_[indexSet_.index(*(eIt->template subEntity<dim> (3)))]
= pow(10.0, permi);
perm_[indexSet_.index(*(eIt->template subEntity<dim> (1)))]
= pow(10.0, permi);
perm_[indexSet_.index(*(eIt->template subEntity<dim> (2)))]
= pow(10.0, permi);
perm_[indexSet_.index(*(eIt->template subEntity<dim> (0)))]
= pow(10.0, permi);
Scalar volume = eIt->geometry().volume();
totalVolume += volume;
meanPermeability += volume / pow(10.0, permi);
}
for (int i=0; i != size; i++)
std::cout <<"perm_"<<i<<" = "<< perm_[i]<<"\n";
meanPermeability /= totalVolume;
meanPermeability = 1.0 / meanPermeability;
std::cout<<"meanPermeability = "<<meanPermeability<<"\n";
infile.close();
std::cout<<"Perm_[1] = "<<perm_[1]<<"\n";
std::cout<<"Size of perm_="<<perm_.size()<<"\n";
}
private:
Dune::FieldMatrix<Scalar, dim, dim> permeability_;
PermeabilityType perm_;
Scalar porosity_;
Scalar poro_;
const IndexSet& indexSet_;
// Object that holds the values/parameters of the selected material law.
MaterialLawParams materialParams_;
/*@\label{tutorial-coupled:matParamsObject}@*/
};
} // end namespace
#endif
_______________________________________________
Dumux mailing list
[email protected]
https://listserv.uni-stuttgart.de/mailman/listinfo/dumux
