Commit 7b29f2bd by Farid Mohammadi Committed by Katharina Heck

image
 ... ... @@ -5,7 +5,7 @@ This tutorial was copied from dumux/test/porousmediumflow/tracer/1ptracer. ## Problem set-up This example contains a contaminant transported by a base groundwater flow in a randomly distributed permeability field. The figure below shows the simulation set-up. The permeability values range between 6.12e-15 and 1.5 e-7 $m^2$. A pressure gradient between the top an the bottom boundary leads to a groundwater flux from the bottom to the top. Neumann no-flow boundaries are assigned to the left and right boundary. Initially, there is a contaminant concentration at the bottom of the domain. ![](./img/setup.png) ## Model description Two different models are applied to simulate the system: In a first step, the groundwater velocity is evaluated under stationary conditions. Therefore the single phase model is applied. ... ... @@ -14,36 +14,162 @@ In a second step, the contaminant gets transported based on the groundwater velo ### 1p Model The single phase model uses Darcy's law as the equation for the momentum conservation: $ \textbf v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho {\textbf g} \right) $ math \textbf v = - \frac{\textbf K}{\mu} \left(\textbf{grad}\, p - \varrho {\textbf g} \right)  With the darcy velocity $ \textbf v $, the permeability $ \textbf K$, the dynamic viscosity $ \mu$, the pressure $p$, the density $\rho$ and the gravity $\textbf g$. Darcy's law is inserted into the continuity equation: $ \phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0$ math \phi \frac{\partial \varrho}{\partial t} + \text{div} \textbf v = 0  with porosity $\phi$. The equation is discretized using a cell-centered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux handbook. The equation is discretized using a cell-centered finite volume scheme as spatial discretization for the pressure as primary variable. For details on the discretization scheme, have a look at the dumux [handbook](https://dumux.org/handbook). ### Tracer Model The transport of the contaminant component $\kappa$ is based on the previously evaluated velocity field $\textbf v$ with the help the following mass balance equation: $ \phi \frac{ \partial X^\kappa}{\partial t} - \text{div} \left\lbrace X^\kappa {\textbf v}+ D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa \right\rbrace = 0 $ math \phi \frac{ \partial X^\kappa}{\partial t} - \text{div} \left\lbrace X^\kappa {\textbf v}+ D^\kappa_\text{pm} \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa \right\rbrace = 0  With the porosity $\phi$, the mass fraction of the contaminant component $\kappa$: $X^\kappa$, the binary diffusion coefficient in the porous medium $ D^\kappa_\text{pm} $, the molar masses of the component $ M^\kappa $, the average molar mass of the phase $M_\alpha$ and the mole fraction $x$. With the porosity $\phi$, the mass fraction of the contaminant component $\kappa$: $X^\kappa$, the porous medium diffusivity $ D^\kappa_\text{pm} $, the molar masses of the component $ M^\kappa $, the average molar mass of the phase $M_\alpha$ and the mole fraction $x$. The porous medium diffusivity is yield out of the diffusion coefficient of the component, the porosity $\phi $ and the porous medium tortuosity $\tau$ by the following equation: The porous medium diffusivity is a function of the diffusion coefficient of the component $D^\kappa$, the porosity $\phi$ and the porous medium tortuosity $\tau$ by the following equation: $ math D^\kappa_\text{pm}= \phi \tau D^\kappa $  The primary variable of this model is the mass fraction $X^\kappa$. We apply the same spatial discretization as in the single pahse model and use the implicit Euler method for time discretization. For more information, have a look at the dumux handbook. In the following, we take a close look at the files containing the setup: At first, boundary conditions and spatially distributed parameters are set in problem_1p.hh and spatialparams_1p.hh, respectively, for the single phase model and subsequently in problem_tracer.hh and spatialparams_tracer.hh for the tracer model. Afterwards, we show the different steps for solving the model in the source file main.cc. At the end, we show some simulation results. ## The file spatialparams_1p.hh In this file, we generate the random permeability field in the constructor of the OnePTestSpatialParams class. Thereafter, spatial properties of the porous medium such as the permeability and the porosity are defined in various functions for the 1p problem. We want to generate a random permeability field. For this, we use a random number generation of the C++ standard library. cpp #include  In the file properties.hh all properties are declared. cpp #include  We include the spatial parameters for single-phase, finite volumes from which we will inherit. cpp #include namespace Dumux {  In the OnePTestSpatialParams class, we define all functions needed to describe the porous matrix, e.g. porosity and permeability for the 1p_problem. cpp template class OnePTestSpatialParams : public FVSpatialParamsOneP> {  We introduce using declarations that are derived from the property system, which we need in this class. cpp using GridView = typename FVGridGeometry::GridView; using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolume = typename FVElementGeometry::SubControlVolume; using Element = typename GridView::template Codim<0>::Entity; using ParentType = FVSpatialParamsOneP>; static constexpr int dimWorld = GridView::dimensionworld; using GlobalPosition = typename SubControlVolume::GlobalPosition; public: using PermeabilityType = Scalar; OnePTestSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry), K_(fvGridGeometry->gridView().size(0), 0.0) {  ### Generation of the random permeability field We get the permeability of the domain and the lens from the params.input file. cpp permeability_ = getParam("SpatialParams.Permeability"); permeabilityLens_ = getParam("SpatialParams.PermeabilityLens");  Further, we get the position of the lens, which is defined by the position of the lower left and the upper right corner. cpp lensLowerLeft_ = getParam("SpatialParams.LensLowerLeft"); lensUpperRight_ =getParam("SpatialParams.LensUpperRight");  We generate random fields for the permeability using a lognormal distribution, with the permeability_ as mean value and 10 % of it as standard deviation. A seperate distribution is used for the lens using permeabilityLens_. cpp std::mt19937 rand(0); std::lognormal_distribution K(std::log(permeability_), std::log(permeability_)*0.1); std::lognormal_distribution KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1); for (const auto& element : elements(fvGridGeometry->gridView())) { const auto eIdx = fvGridGeometry->elementMapper().index(element); const auto globalPos = element.geometry().center(); K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand); } }  ### Properties of the porous matrix We define the (intrinsic) permeability $[m^2]$ using the generated random permeability field. In this test, we use element-wise distributed permeabilities. cpp template const PermeabilityType& permeability(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const { return K_[scv.dofIndex()]; }  We set the porosity $[-]$ for the whole domain. cpp Scalar porosityAtPos(const GlobalPosition &globalPos) const { return 0.2; }  We reference to the permeability field. This is used in the main function to write an output for the permeability field. cpp const std::vector& getKField() const { return K_; } private:  We have a convenient definition of the position of the lens. cpp bool isInLens_(const GlobalPosition &globalPos) const { for (int i = 0; i < dimWorld; ++i) { if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_) return false; } return true; } GlobalPosition lensLowerLeft_; GlobalPosition lensUpperRight_; Scalar permeability_, permeabilityLens_; std::vector K_; static constexpr Scalar eps_ = 1.5e-7; }; } // end namespace Dumux  ## The file problem_1p.hh ... ... @@ -92,8 +218,7 @@ We enter the namespace Properties, which is a sub-namespace of the namespace Dum cpp namespace Properties {  A TypeTag for our simulation is created which inherits from the one-phase flow model and the cell centered, two-point-flux discretization scheme. A TypeTag for our simulation is created which inherits from the one-phase flow model and the cell centered, two-point-flux discretization scheme. cpp namespace TTag { struct IncompressibleTest { using InheritsFrom = std::tuple; }; ... ... @@ -115,12 +240,12 @@ template struct SpatialParams {  We define convenient shortcuts to the properties FVGridGeometry and Scalar: We define convenient shortcuts to the properties FVGridGeometry and Scalar: cpp using FVGridGeometry = GetPropType; using Scalar = GetPropType;  Finally we set the spatial parameters: Finally, we set the spatial parameters: cpp using type = OnePTestSpatialParams; }; ... ... @@ -146,12 +271,12 @@ description of water, which means we do not use tabulated values but more genera using type = FluidSystems::OnePLiquid >; };  We enable caching for the grid volume variables We enable caching for the grid volume variables. cpp template struct EnableGridVolumeVariablesCache { static constexpr bool value = true; };  We enable caching for the grid flux variables We enable caching for the grid flux variables. cpp template struct EnableGridFluxVariablesCache { static constexpr bool value = true; }; ... ... @@ -168,7 +293,7 @@ We leave the namespace Properties.  ### The problem class We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. As this is a porous medium problem, we inherit from the basic PorousMediumFlowProblem. As this is a porous medium problem, we inherit from the basic PorousMediumFlowProblem. cpp template class OnePTestProblem : public PorousMediumFlowProblem ... ... @@ -263,118 +388,96 @@ We leave the namespace Dumux. ## The file spatialparams_1p.hh ## The file spatialparams_tracer.hh In this file we generate the random permeability field in the constructor of the OnePTestSpatialParams class. Thereafter spatial properties of the porous medium like permeability and porosity are defined in various functions for the 1p problem. We want to generate a random permeability field. For this we use a random number generation of the C++ standard library. cpp #include  In the file properties.hh all properties are declared. In this file, we define spatial properties of the porous medium such as the permeability and the porosity in various functions for the tracer problem. Further, spatial dependent properties of the tracer fluid system are defined and in the end two functions handel the calculated volume fluxes from the solution of the 1p problem. In the file properties.hh, all properties are declared. cpp #include  We include the spatial parameters for single-phase, finite volumes from which we will inherit. As in the 1p spatialparams, we inherit from the spatial parameters for single-phase, finite volumes, which we include here. cpp #include  We enter the namespace Dumux cpp namespace Dumux {  In the OnePTestSpatialParams class, we define all functions needed to describe the porous matrix, e.g. define porosity and permeability for the 1p_problem. In the TracerTestSpatialParams class, we define all functions needed to describe spatially dependent parameters for the tracer_problem. cpp template class OnePTestSpatialParams class TracerTestSpatialParams : public FVSpatialParamsOneP> TracerTestSpatialParams> {  We introduce using declarations that are derived from the property system which we need in this class. cpp using GridView = typename FVGridGeometry::GridView; using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolume = typename FVElementGeometry::SubControlVolume; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using Element = typename GridView::template Codim<0>::Entity; using ParentType = FVSpatialParamsOneP>; TracerTestSpatialParams>; static constexpr int dimWorld = GridView::dimensionworld; using GlobalPosition = typename SubControlVolume::GlobalPosition; static const int dimWorld = GridView::dimensionworld; using GlobalPosition = typename Dune::FieldVector; public: using PermeabilityType = Scalar; OnePTestSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry), K_(fvGridGeometry->gridView().size(0), 0.0) {  ### Generation of the random permeability field We get the permeability of the domain and the lens from the params.input file. cpp permeability_ = getParam("SpatialParams.Permeability"); permeabilityLens_ = getParam("SpatialParams.PermeabilityLens");  Further, we get the position of the lens, which is defined by the position of the lower left and the upper right corner. cpp lensLowerLeft_ = getParam("SpatialParams.LensLowerLeft"); lensUpperRight_ =getParam("SpatialParams.LensUpperRight"); TracerTestSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) {}  We generate random fields for the permeability using a lognormal distribution, with the permeability_ as mean value and 10 % of it as standard deviation. A seperate distribution is used for the lens using permeabilityLens_. ### Properties of the porous matrix We define the same porosity for the whole domain as in the 1p spatialparams. cpp std::mt19937 rand(0); std::lognormal_distribution K(std::log(permeability_), std::log(permeability_)*0.1); std::lognormal_distribution KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1); for (const auto& element : elements(fvGridGeometry->gridView())) { const auto eIdx = fvGridGeometry->elementMapper().index(element); const auto globalPos = element.geometry().center(); K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand); } } Scalar porosityAtPos(const GlobalPosition& globalPos) const { return 0.2; }  ### Properties of the porous matrix We define the (intrinsic) permeability $[m^2]$ using the generated random permeability field. In this test, we use element-wise distributed permeabilities. We do not consider dispersivity for the tracer transport. So we set the dispersivity coefficient to zero. cpp template const PermeabilityType& permeability(const Element& element, const SubControlVolume& scv, const ElementSolution& elemSol) const { return K_[scv.dofIndex()]; } Scalar dispersivity(const Element &element, const SubControlVolume& scv, const ElementSolution& elemSol) const { return 0; }  We set the porosity $[-]$ for the whole domain. ### Properties of the fluid system In the following, we define fluid properties that are spatial parameters in the tracer model. They can possible vary with space but are usually constants. Further spatially constant values of the fluid system are defined in the TracerFluidSystem class in problem.hh. We define the fluid density to a constant value of 1000 $\frac{kg}{m^3}$. cpp Scalar porosityAtPos(const GlobalPosition &globalPos) const { return 0.2; } Scalar fluidDensity(const Element &element, const SubControlVolume& scv) const { return 1000; }  We reference to the permeability field. This is used in the main function to write an output for the permeability field. We define the fluid molar mass. cpp const std::vector& getKField() const { return K_; } Scalar fluidMolarMass(const Element &element, const SubControlVolume& scv) const { return 18.0; } private: Scalar fluidMolarMass(const GlobalPosition &globalPos) const { return 18.0; }  We have a convenience definition of the position of the lens. ### The volume fluxes We define a function which returns the field of volume fluxes. This is e.g. used to calculate the transport of the tracer. cpp bool isInLens_(const GlobalPosition &globalPos) const template Scalar volumeFlux(const Element &element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) const { for (int i = 0; i < dimWorld; ++i) { if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_) return false; } return true; return volumeFlux_[scvf.index()]; }  We define a function to set the volume flux. This is used in the main function to set the volume flux to the calculated value based on the solution of the 1p problem. cpp void setVolumeFlux(const std::vector& f) { volumeFlux_ = f; } GlobalPosition lensLowerLeft_; GlobalPosition lensUpperRight_; Scalar permeability_, permeabilityLens_; std::vector K_; static constexpr Scalar eps_ = 1.5e-7; private: std::vector volumeFlux_; }; } // end namespace Dumux ... ... @@ -388,7 +491,7 @@ We have a convenience definition of the position of the lens. Before we enter the problem class containing initial and boundary conditions, we include necessary files and introduce properties. ### Include files Again, we have to include the dune grid interphase: Again, we have to include the dune grid interface: cpp #include  ... ... @@ -426,7 +529,7 @@ We enter the namespace Properties, cpp namespace Properties {  A TypeTag for our simulation is created which inherits from the tracer model and the A TypeTag for our simulation is created which inherits from the tracer model and the cell centered, two-point-flux discretization scheme. cpp namespace TTag { ... ... @@ -468,20 +571,20 @@ We define that mass fractions are used to define the concentrations template struct UseMoles { static constexpr bool value = false; };  We don't use a solution dependent molecular diffusion coefficient: We do not use a solution dependent molecular diffusion coefficient: cpp template struct SolutionDependentMolecularDiffusion { static constexpr bool value = false; };  In the following we create a new tracer fluid system and derive it from the base fluid system. In the following, we create a new tracer fluid system and derive it from the base fluid system. cpp template class TracerFluidSystem : public FluidSystems::Base, TracerFluidSystem> {  We define convenient shortcuts to the properties Scalar, Problem, GridView, Element, FVElementGeometry and SubControlVolume: We define convenient shortcuts to the properties Scalar, Problem, GridView, Element, FVElementGeometry and SubControlVolume: cpp using Scalar = GetPropType; using Problem = GetPropType; ... ... @@ -506,22 +609,22 @@ and the number of components cpp static constexpr int numComponents = 1;  We set the component name for the component index (compIdx) for the vtk output: We set the component name for the component index (compIdx) for the vtk output: cpp static std::string componentName(int compIdx) { return "tracer_" + std::to_string(compIdx); }  We set the phase name for the phase index (phaseIdx) for velocity vtk output: We set the phase name for the phase index (phaseIdx) for velocity vtk output: cpp static std::string phaseName(int phaseIdx = 0) { return "Groundwater"; }  We set the molar mass of the tracer component with index compIdx We set the molar mass of the tracer component with index compIdx. cpp static Scalar molarMass(unsigned int compIdx) { return 0.300; }  And we set the value for the binary diffusion coefficient. This We set the value for the binary diffusion coefficient. This might depend on spatial parameters like pressure / temperature. But for our case it is 0.0: cpp static Scalar binaryDiffusionCoefficient(unsigned int compIdx, ... ... @@ -543,7 +646,7 @@ We leave the namespace Properties.  ### The problem class We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation. As this is a porous medium problem, we inherit from the basic PorousMediumFlowProblem. As this is a porous medium problem, we inherit from the basic PorousMediumFlowProblem. cpp template class TracerTestProblem : public PorousMediumFlowProblem ... ... @@ -634,104 +737,6 @@ We leave the namespace Dumux here. ## The file spatialparams_tracer.hh In this file we define spatial properties of the porous medium like permeability and porosity in various functions for the tracer problem. Further spatial dependent properties of the tracer fluid system are defined and in the end two functions handel the calculated volume fluxes from the solution of the 1p problem. In the file properties.hh all properties are declared. cpp #include  As in the 1p spatialparams we inherit from the spatial parameters for single-phase, finite volumes, which we include here. cpp #include  We enter the namespace Dumux cpp namespace Dumux {  In the TracerTestSpatialParams class, we define all functions needed to describe spatially dependent parameters for the tracer_problem. cpp template class TracerTestSpatialParams : public FVSpatialParamsOneP> { using GridView = typename FVGridGeometry::GridView; using FVElementGeometry = typename FVGridGeometry::LocalView; using SubControlVolume = typename FVElementGeometry::SubControlVolume; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; using Element = typename GridView::template Codim<0>::Entity; using ParentType = FVSpatialParamsOneP>; static const int dimWorld = GridView::dimensionworld; using GlobalPosition = typename Dune::FieldVector; public: TracerTestSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) {}  ### Properties of the porous matrix We define the same porosity for the whole domain as in the 1p spatialparams. cpp Scalar porosityAtPos(const GlobalPosition& globalPos) const { return 0.2; }  We don't consider dispersivity for the tracer transport. So we set the dispersivity coefficient to zero. cpp template Scalar dispersivity(const Element &element, const SubControlVolume& scv, const ElementSolution& elemSol) const { return 0; }  ### Properties of the fluid system In the following we define fluid properties that are spatial parameters in the tracer model. They can possible vary with space but are usually constants. Further spatially constant values of the fluid system are defined in the TracerFluidSystem class in problem.hh. We define the fluid density to a constant value of 1000 $\frac{kg}{m^3}$. cpp Scalar fluidDensity(const Element &element, const SubControlVolume& scv) const { return 1000; }  We define the fluid molar mass. cpp Scalar fluidMolarMass(const Element &element, const SubControlVolume& scv) const { return 18.0; } Scalar fluidMolarMass(const GlobalPosition &globalPos) const { return 18.0; }  ### The volume fluxes We define a function which returns the field of volume fluxes. This is e.g. used to calculate the transport of the tracer. cpp template Scalar volumeFlux(const Element &element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) const { return volumeFlux_[scvf.index()]; }  We define a function to set the volume flux. This is used in the main function to set the volume flux to the calculated value based on the solution of the 1p problem. cpp void setVolumeFlux(const std::vector& f) { volumeFlux_ = f; } private: std::vector volumeFlux_; }; } // end namespace Dumux  ## The file main.cc ... ... @@ -740,12 +745,12 @@ We look now at the main file for the tracer problem. We setup two problems in th cpp #include  We include both problems in the main file, the problem_1p.hh and the problem_tracer.hh. We include both problems in the main file, the problem_1p.hh and the problem_tracer.hh. cpp #include "problem_1p.hh" #include "problem_tracer.hh"  Further we include a standard header file for C++, to get time and date information Further, we include a standard header file for C++, to get time and date information cpp #include  ... ... @@ -760,7 +765,7 @@ Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which #include #include  In Dumux a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh. In Dumux, a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh. cpp #include  ... ... @@ -768,7 +773,7 @@ The following file contains the parameter class, which manages the definition of cpp #include `