Commit bcfa3405 authored by Philipp Nuske's avatar Philipp Nuske
Browse files

dumux/io/interfacemeshcreator

moved from devel: a meshcreator which gradually refines towards the (free flow - porous medium) interface

dumux/io/plotoverLine2d
moved from devel: a class for plotting variables between two (user-specified) points. Values are written to a *.dat file, that can be read in by an external plotting tool

dumux/dumx/common/dimensionlessnumbers
moved from devel: a class for calculating dimensionless numbers (Reynolds, Prandtl, ...) 

reviewed by Bernd


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11056 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 7cba3389
/*****************************************************************************
* Copyright (C) 2008 by Philipp Nuske *
* Institute for Modelling Hydraulic and Environmental Systems *
* University of Stuttgart, Germany *
* email: <givenname>.<name>@iws.uni-stuttgart.de *
* *
* 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 dimensionlessnumbers.hh
*
* \brief Collection of functions, calculating dimensionless numbers.
*
* All the input to the dimensionless numbers has to be provided as function arguments.
* Rendering this collection generic in the sense that it can be used by any model.
*/
#ifndef DIMENSIONLESS_NUMBERS_HH
#define DIMENSIONLESS_NUMBERS_HH
namespace Dumux
{
/*
* \brief Dimensionless Numbers
*
* Each number has it's own function.
* All the parameters for the calculation have to be handed over.
* Rendering this collection generic in the sense that it can be used by any model.
*/
template <class TypeTag>
class DimensionlessNumbers
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
/*!
* \brief Calculate the Reynolds Number [-] (Re).
*
* The Reynolds number is a measure for the relation of inertial to viscous forces.
* The bigger the value, the more important inertial (as compared to viscous) effects become.
* According to Bear [Dynamics of fluids in porous media (1972)] Darcy's law is valid for Re<1.
*
* source for Reynolds number definition: http://en.wikipedia.org/wiki/Reynolds_number
*
* \param darcyMagVelocity The absolute value of the darcy velocity. In the context of box models this
* leads to a problem: the velocities are defined on the faces while other things (storage, sources, output)
* are defined for the volume/vertex. Therefore, some sort of decision needs to be made which velocity to put
* into this function (e.g.: face-area weighted average). [m/s]
* \param charcteristicLength Typically, in the context of porous media flow, the mean grain size is taken as the characteristic length
* for calculation of Re. [m]
* \param kinematicViscosity Is defined as the dynamic (or absolute) viscosity divided by the density. http://en.wikipedia.org/wiki/Viscosity#Dynamic_viscosity. [m^2/s]
*
* \return The Reynolds Number as calculated from the input parameters
*/
static Scalar reynoldsNumber(const Scalar darcyMagVelocity,
const Scalar charcteristicLength,
const Scalar kinematicViscosity)
{
return darcyMagVelocity * charcteristicLength / kinematicViscosity ;
}
/*!
* \brief Calculate the Prandtl Number [-] (Pr).
*
* The Prandtl Number is a measure for the relation of viscosity and thermal diffusivity (temperaturleitfaehigkeit).
*
* It is defined as
* \f[
* \textnormal{Pr}= \frac{\nu}{\alpha} = \frac{c_p \mu}{\lambda}\, ,
* \f]
* with kinematic viscosity\f$\nu\f$, thermal diffusivity \f$\alpha\f$, heat capacity \f$c_p\f$, dynamic viscosity \f$\mu\f$ and thermal conductivity \f$\lambda\f$.
* Therefore, Pr is a material specific property (i.e.: not a function of flow directly but only of temperature, pressure and fluid).
*
* source for Prandtl number definition: http://en.wikipedia.org/wiki/Prandtl_number
*
* \param dynamicViscosity Dynamic (absolute) viscosity over density. http://en.wikipedia.org/wiki/Viscosity#Dynamic_viscosity [m^2/s]
* \param heatCapacity Heat capacity at constant pressure. Specifies the energy change for a given temperature change [J / (kg K)]
* \param thermalConductivity Conductivity to heat. Specifies how well matter transfers energy without moving. [W/(m K)]
* \return The Prandtl Number as calculated from the input parameters.
*/
static Scalar prandtlNumber(const Scalar dynamicViscosity,
const Scalar heatCapacity,
const Scalar thermalConductivity)
{
return dynamicViscosity * heatCapacity / thermalConductivity;
}
/*!
* \brief Calculate the Nusselt Number [-] (Nu).
*
* The Nusselt Number is a measure for the relation of convective- to conductive heat exchange.
*
* The Nusselt number is defined as Nu = h d / k,
* with h= heat transfer coefficient, d=characteristic length, k=heat conductivity(stagnant).
* However, the heat transfer coefficient from one phase to another is typically not known.
* Therefore, Nusselt numbers are usually given as *empirical* Nu(Reynolds, Prandtl) for a given flow
* field --forced convection-- and *empirical* Nu(Rayleigh, Prandtl) for flow caused by temperature
* differences --free convection--. The fluid characteristics enter via the Prandtl number.
*
* This function implements an *empirical* correlation for the case of porous media flow
* (packed bed flow as the chemical engineers call it).
*
* source for Nusselt number definition: http://en.wikipedia.org/wiki/Nusselt_number
* source for further empirical correlations for Nusselt Numbers: VDI-Gesellschaft, VDI-Waermeatlas, VDI-Verlag Duesseldorf, 2006
*
* \param prandtlNumber Dimensionless number relating viscosity and thermal diffusivity (temperaturleitfaehigkeit) [-].
* \param ReynoldsNumber Dimensionless number relating inertial and viscous forces [-].
* \return The Nusselt number as calculated from the input parameters [-].
*/
static Scalar nusseltNumberForced(const Scalar reynoldsNumber,
const Scalar prandtlNumber,
const Scalar porosity)
{
/* example: very common and simple case: flow straight circular pipe, only convection (no boiling), 10000<Re<120000, 0.7<Pr<120, far from pipe entrance, smooth surface of pipe ...
Dittus, F.W and Boelter, L.M.K, Heat Transfer in Automobile Radiators of the Tubular Type, Publications in Engineering, Vol. 2, pages 443-461, 1930
*/
// return 0.023 * pow(reynoldsNumber, 0.8) * pow(prandtlNumber,0.33);
/* example: flow through porous medium *single phase*, fit to many different data
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 293
*/
return 2. + 1.1 * pow(prandtlNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
/* example: VDI Waermeatlas 10. Auflage 2006, flow in packed beds, page Gj1, see also other sources and limitations therein.
* valid for 0.1<Re<10000, 0.6<Pr/Sc<10000, packed beds of perfect spheres.
*
*/
// Scalar numerator = 0.037 * pow(reynoldsNumber,0.8) * prandtlNumber ;
// Scalar reToMin01 = pow(reynoldsNumber,-0.1);
// Scalar prTo23 = pow(prandtlNumber, (2./3. ) ) ; // MIND THE pts! :-( otherwise the integer exponent version is chosen
// Scalar denominator = 1+ 2.443 * reToMin01 * (prTo23 -1.) ;
//
// Scalar nusseltTurbular = numerator / denominator;
// Scalar nusseltLaminar = 0.664 * sqrt(reynoldsNumber) * pow(prandtlNumber, (1./3.) );
// Scalar nusseltSingleSphere = 2 + sqrt( pow(nusseltLaminar,2.) + pow(nusseltTurbular,2.));
//
// Scalar funckyFactor = 1 + 1.5 * (1.-porosity); // for spheres of same size
// Scalar nusseltNumber = funckyFactor * nusseltSingleSphere ;
//
// return nusseltNumber;
}
/*!
* \brief Calculate the Schmidt Number [-] (Sc).
*
* The Schmidt Number is a measure for the relation of viscosity and mass diffusivity.
*
* It is defined as
* \f[
* \textnormal{Sc}= \frac{\nu}{D} = \frac{\mu}{\rho D}\, ,
* \f]
* with kinematic viscosity\f$\nu\f$, diffusion coefficient \f$D\f$, dynamic viscosity \f$\mu\f$ and mass density\f$\rho\f$.
* Therefore, Sc is a material specific property (i.e.: not a function of flow directly but only of temperature, pressure and fluid).
*
* source for Schmidt number definition: http://en.wikipedia.org/wiki/Schmidt_number
*
* \param dynamicViscosity Dynamic (absolute) viscosity over density. http://en.wikipedia.org/wiki/Viscosity#Dynamic_viscosity [m^2/s]
* \param massDensity Mass density of the considered phase. [kg / m^3]
* \param diffusionCoefficient Measure for how well a component can move through a phase due to a concentration gradient. [m^2/s]
* \return The Schmidt Number as calculated from the input parameters.
*/
static Scalar schmidtNumber(const Scalar dynamicViscosity,
const Scalar massDensity,
const Scalar diffusionCoefficient)
{
return dynamicViscosity / (massDensity * diffusionCoefficient);
}
/*!
* \brief Calculate the Sherwood Number [-] (Sh).
*
* The Sherwood Number is a measure for the relation of convective- to diffusive mass exchange.
*
* The Sherwood number is defined as Sh = K L/D,
* with K= mass transfer coefficient, L=characteristic length, D=mass diffusivity (stagnant).
*
* However, the mass transfer coefficient from one phase to another is typically not known.
* Therefore, Sherwood numbers are usually given as *empirical* Sh(Reynolds, Schmidt) for a given flow
* field (and fluid).
*
* Often, even the Sherwood number is not known. By means of the Chilton-Colburn analogy it can be deduced
* from the Nusselt number. According to the Chilton-Colburn analogy in a known Nusselt correltion one
* basically replaces Pr with Sc and Nu with Sh. For some very special cases this is actually accurate.
* (Source: Course Notes, Waerme- und Stoffuebertragung, Prof. Hans Hasse, Uni Stuttgart)
*
* This function implements an *empirical* correlation for the case of porous media flow
* (packed bed flow as the chemical engineers call it).
*
* source for Sherwood number definition: http://en.wikipedia.org/wiki/Sherwood_number
*
* \param schmidtNumber Dimensionless number relating viscosity and mass diffusivity [-].
* \param reynoldsNumber Dimensionless number relating inertial and viscous forces [-].
* \return The Nusselt number as calculated from the input parameters [-].
*/
static Scalar sherwoodNumber(const Scalar reynoldsNumber,
const Scalar schmidtNumber)
{
/* example: flow through porous medium *single phase*
* Wakao and Kaguei, Heat and mass Transfer in Packed Beds, Gordon and Breach Science Publishers, page 156
*/
return 2. + 1.1 * pow(schmidtNumber,(1./3.)) * pow(reynoldsNumber, 0.6);
}
/*!
* \brief Calculate the thermal diffusivity alpha [m^2/s].
*
* The thermal diffusivity is a measure for how fast "temperature (not heat!) spreads".
* It is defined as \alpha = k / (rho c_p)
* with \alpha: k: thermal conductivity [W/mK], rho: density [kg/m^3], c_p: cpecific heat capacity at constant pressure [J/kgK].
*
* Source for thermal diffusivity definition: http://en.wikipedia.org/wiki/Thermal_diffusivity
*
* \param thermalConductivity A material property defining how well heat is transported via conduction [W/(mK)].
* \param phaseDensity The density of the phase for which the thermal diffusivity is to be calculated [kg/m^3].
* \param heatCapacity A measure for how a much a material changes temperature for a given change of energy (at p=const.) [J/(kgm^3)].
* \return The thermal diffusivity as calculated from the input parameters [m^2/s].
*/
static Scalar thermalDiffusivity(const Scalar & thermalConductivity ,
const Scalar & phaseDensity ,
const Scalar & heatCapacity)
{
return thermalConductivity / (phaseDensity * heatCapacity);
}
}; // end class
} // end namespace Dumux
#endif // DIMENSIONLESS_NUMBERS_HH
#ifndef DUMUX_INTERFACEMESHCREATOR_HH
#define DUMUX_INTERFACEMESHCREATOR_HH
#include<dune/grid/common/gridfactory.hh>
#include <dune/grid/io/file/dgfparser/dgfs.hh>
namespace Dumux
{
template <class Grid>
class InterfaceMeshCreator
{
public:
enum {dim = Grid::dimension};
typedef typename Grid::ctype Scalar;
typedef typename Grid::LeafGridView GridView;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename Element::Geometry Geometry;
static Grid* create(const std::string& dgfName, const Dune::FieldVector<int, dim>& nElements,
const Scalar interfaceY, const Scalar gradingFactor, const bool refineTop = false)
{
typedef Dune::SGrid<dim,dim> HelperGrid;
Dune::GridPtr<HelperGrid> helperGridPtr(dgfName);
HelperGrid& helperGrid = *helperGridPtr;
typedef typename HelperGrid::LeafGridView HelperGridView;
typedef typename HelperGridView::template Codim<0>::Iterator HelperElementIterator;
typedef typename HelperGridView::Traits::template Codim<0>::Entity HelperElement;
typedef typename HelperElement::Geometry HelperGeometry;
HelperElementIterator helperElementIterator = helperGrid.leafView().template begin<0>();
const HelperElement& helperElement = *helperElementIterator;
const HelperGeometry& helperGeometry = helperElement.geometry();
Dune::FieldVector<Scalar,dim> refinePoint(0);
refinePoint[dim-1] = interfaceY;
Dune::FieldVector<Scalar,dim> gradingFactors(1);
gradingFactors[dim-1] = gradingFactor;
Dune::FieldVector<Scalar,dim> refinePointLocal(helperGeometry.local(refinePoint));
std::cout << "rglobal = " << refinePoint << ", rlocal = " << refinePointLocal << std::endl;
Dune::GridFactory<Grid> factory;
int nX = nElements[0];
int nY = nElements[1];
std::vector<std::vector<Scalar> > localPositions(dim);
for (int comp = 0; comp < dim; comp++)
{
Scalar lengthLeft = refinePointLocal[comp];
Scalar lengthRight = 1.0 - lengthLeft;
int nLeft, nRight;
Scalar hLeft, hRight;
if (lengthLeft < 1e-10)
{
nLeft = 0;
nRight = nElements[comp];
if (gradingFactors[comp] > 1.0)
hLeft = hRight = (1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nRight));
else
hLeft = hRight = 1.0/nElements[comp];
}
else if (lengthLeft > 1.0 - 1e-10)
{
nLeft = nElements[comp];
nRight = 0;
if (gradingFactors[comp] > 1.0)
hLeft = hRight = (1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nLeft));
else
hLeft = hRight = 1.0/nElements[comp];
}
else if (comp == dim - 1 && refineTop)
{
lengthLeft = refinePointLocal[comp];
lengthRight = (1 - refinePointLocal[comp])/2;
nLeft = nRight = nElements[comp]/3;
if (nElements[comp]%3 == 1)
nLeft += 1;
else if (nElements[comp]%3 == 2)
nRight += 1;
hLeft = lengthLeft*(1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nRight));
}
else if (lengthLeft > 0.5)
{
Scalar nLeftDouble = std::ceil(-log((1.0 + sqrt(1.0 + 4.0*pow(gradingFactors[comp], nElements[comp])*lengthRight/lengthLeft))
/(2.0*pow(gradingFactors[comp], nElements[comp])))/log(gradingFactors[comp]));
nLeft = std::min((int)std::ceil(nLeftDouble), nElements[comp]);
nRight = nElements[comp] - nLeft;
if (gradingFactors[comp] > 1.0)
{
hLeft = lengthLeft*(1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nRight));
}
else
hLeft = hRight = 1.0/nElements[comp];
}
else
{
Scalar nRightDouble = -log((1.0 + sqrt(1.0 + 4.0*pow(gradingFactors[comp], nElements[comp])*lengthLeft/lengthRight))
/(2.0*pow(gradingFactors[comp], nElements[comp])))/log(gradingFactors[comp]);
nRight = std::min((int)std::ceil(nRightDouble), nElements[comp]);
nLeft = nElements[comp] - nRight;
if (gradingFactors[comp] > 1.0)
{
hLeft = lengthLeft*(1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nLeft));
hRight = lengthRight*(1.0 - gradingFactors[comp])/(1.0 - pow(gradingFactors[comp], nRight));
}
else
hLeft = hRight = 1.0/nElements[comp];
}
std::cout << "lengthLeft = " << lengthLeft << ", lengthRight = " << lengthRight << ", hLeft = " << hLeft << ", hRight = " << hRight << ", nLeft = " << nLeft << ", nRight = " << nRight << std::endl;
int nVertices = nElements[comp] + 1;
localPositions[comp].resize(nVertices);
localPositions[comp][0] = 0.0;
for (int i = 0; i < nLeft; i++)
{
Scalar hI = hLeft*pow(gradingFactors[comp], nLeft-1-i);
localPositions[comp][i+1] = localPositions[comp][i] + hI;
}
for (int i = 0; i < nRight; i++)
{
Scalar hI = hRight*pow(gradingFactors[comp], i);
localPositions[comp][nLeft+i+1] = localPositions[comp][nLeft+i] + hI;
}
if (comp == dim - 1 && refineTop)
for (int i = 0; i < nRight; i++)
{
Scalar hI = hRight*pow(gradingFactors[comp], nRight-1-i);
localPositions[comp][nLeft+nRight+i+1] = localPositions[comp][nLeft+nRight+i] + hI;
}
if (localPositions[comp][nVertices-1] != 1.0)
{
for (int i = 0; i < nVertices; i++)
localPositions[comp][i] /= localPositions[comp][nVertices-1];
}
}
Dune::FieldVector<Scalar,dim> local;
for (int j = 0; j < nY + 1; j++)
{
local[1] = localPositions[1][j];
for (int i = 0; i < nX + 1; i++)
{
local[0] = localPositions[0][i];
Dune::FieldVector<Scalar,dim> position(helperGeometry.global(local));
factory.insertVertex(position);
}
}
for (int j = 0; j < nY; j++)
{
for (int i = 0; i < nX; i++)
{
std::vector<unsigned int> vertices(4);
vertices[0] = j*(nX+1) + i;
vertices[1] = j*(nX+1) + i+1;
vertices[2] = (j+1)*(nX+1) + i;
vertices[3] = (j+1)*(nX+1) + i+1;
factory.insertElement(Dune::GeometryType(Dune::GeometryType::cube,dim), vertices);
}
}
return factory.createGrid();
}
};
}
#endif
/*!
* \file
* \brief Plot variables over a line specified by two arguments.
* These output files are meant for visualization with another program (matlab, gnuplot...)
*
*/
#ifndef DUMUX_PLOTOVERLINE_2D_HH
#define DUMUX_PLOTOVERLINE_2D_HH
#include <dumux/common/valgrind.hh>
#include <dumux/common/propertysystem.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <iostream>
#include <fstream>
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(Scalar);
NEW_PROP_TAG(Problem);
NEW_PROP_TAG(GridView);
NEW_PROP_TAG(DofMapper);
NEW_PROP_TAG(ElementSolutionVector);
NEW_PROP_TAG(SolutionVector);
NEW_PROP_TAG(FVElementGeometry);
NEW_PROP_TAG(TwoPIAIndices);
NEW_PROP_TAG(NumEq);
NEW_PROP_TAG(MaterialLaw);
NEW_PROP_TAG(ElementVolumeVariables);
NEW_PROP_TAG(AwnSurface);
NEW_PROP_TAG(AwnSurfaceParams);
}
template<class TypeTag>
class PlotOverLine2D
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename MaterialLaw::Params aterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum {
wPhaseIdx = FluidSystem::wPhaseIdx,
nPhaseIdx = FluidSystem::nPhaseIdx,
sPhaseIdx = FluidSystem::sPhaseIdx,
wCompIdx = FluidSystem::wCompIdx,
nCompIdx = FluidSystem::nCompIdx,
// Grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
};
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
public:
/*
* \brief A function that writes results over a line (like paraview's plotOverline into a text file.)
*
* The writer needs to be called in postTimeStep().
*
* This function puts output variables (TemperaturePhase, Saturation, t, tIndex, ...) over space (1D, over a line) into a text file,
* so they can be read in by another program like matlab.
* The file can be found by the extension: dat
*/
void write(const Problem & problem,
const GlobalPosition & pointOne,
const GlobalPosition & pointTwo,
const std::string appendOutputName = "")
{
static_assert(dim==2, "this implements plot over Line: so far this works only for 2D");
FVElementGeometry fvGeometry;
ElementVolumeVariables elemVolVars;
// so many vertices in the domain
const unsigned int numGlobalVerts = problem.gridView().size(dim);
// check whether a vertex was already visited by true / false vector
bool isVisited[numGlobalVerts] ;
// loop all vertices of the main: initialize
for (int idx=0; idx < numGlobalVerts; idx++)
isVisited[idx] = false ;
// Looping over all elements of the domain
ElementIterator eEndIt = problem.gridView().template end<0>();
for (ElementIterator eIt = problem.gridView().template begin<0>() ; eIt not_eq eEndIt; ++eIt)
{
// updating the volume variables
fvGeometry.update(problem.gridView(), *eIt);
elemVolVars.update(problem, *eIt, fvGeometry, false);
// number of scv
const unsigned int numScv = fvGeometry.numScv;
for (int scvIdx=0; scvIdx < numScv; ++scvIdx) {
// find some global identification
const unsigned int globalIdx = problem.vertexMapper().map(*eIt, scvIdx, dim);
// only write out if the vertex was not already visited
if (isVisited[globalIdx])
continue;
isVisited[globalIdx] = true ;
// Getting the spatial coordinate
const GlobalPosition & globalPosCurrent
= fvGeometry.subContVol[scvIdx].global;
std::ofstream dataFile;
const unsigned int timeStepIndex = problem.timeManager().timeStepIndex() ;
// filename of the output file
std::string fileName = problem.name();
// filename consists of problem name + some function argument + .dat
// this way several plot overlines can be written from one simulation
fileName += appendOutputName;
fileName +=".dat";
// Writing a header into the output
if (timeStepIndex == 0)
{
dataFile.open(fileName.c_str());
dataFile << "# This is a DuMuX output file for further processing with the preferred graphics program of your choice. \n";
dataFile << "# This output file was written from "<< __FILE__ << ", line " <<__LINE__ << "\n";
dataFile << "# This output file was generated from code compiled at " << __TIME__ <<", "<< __DATE__<< "\n";
dataFile << "\n";
dataFile << "# Header\n";
dataFile << "#timestep\t time\t\t \t\t x \t\t y \t\tSw \t\t\t Tw\t\t Tn\t Ts \t xH2On \t xH2OnEquil \t xN2w \txN2wEquil " << std::endl;
dataFile.close();
}
// write output if the current location is between two specified points
if (isBetween(globalPosCurrent, pointOne, pointTwo))
{
const Scalar time = problem.timeManager().time();
const Scalar saturationW = elemVolVars[scvIdx].fluidState().saturation(wPhaseIdx);
const Scalar Tw = elemVolVars[scvIdx].fluidState().temperature(wPhaseIdx);
const Scalar Tn = elemVolVars[scvIdx].fluidState().temperature(nPhaseIdx);
const Scalar Ts = elemVolVars[scvIdx].fluidState().temperature(sPhaseIdx);
const Scalar xH2On = elemVolVars[scvIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
const Scalar xH2OnEquil = elemVolVars[scvIdx].xEquil(nPhaseIdx, wCompIdx);
const Scalar xN2w = elemVolVars[scvIdx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
const Scalar xN2wEquil = elemVolVars[scvIdx].xEquil(wPhaseIdx, nCompIdx);
// actual output into the text file
// This could be done much more efficiently
// if writing to the text file would be done all at once.
dataFile.open(fileName.c_str() , std::ios::app);